Actual source code: shift.c
slepc-3.18.3 2023-03-24
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: Shift spectral transformation, applies (A + sigma I) as operator, or
12: inv(B)(A + sigma B) for generalized problems
13: */
15: #include <slepc/private/stimpl.h>
17: PetscErrorCode STBackTransform_Shift(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
18: {
19: PetscInt j;
21: for (j=0;j<n;j++) {
22: eigr[j] += st->sigma;
23: }
24: return 0;
25: }
27: PetscErrorCode STPostSolve_Shift(ST st)
28: {
29: if (st->matmode == ST_MATMODE_INPLACE) {
30: if (st->nmat>1) MatAXPY(st->A[0],st->sigma,st->A[1],st->str);
31: else MatShift(st->A[0],st->sigma);
32: st->Astate[0] = ((PetscObject)st->A[0])->state;
33: st->state = ST_STATE_INITIAL;
34: st->opready = PETSC_FALSE;
35: }
36: return 0;
37: }
39: /*
40: Operator (shift):
41: Op P M
42: if nmat=1: A-sI NULL A-sI
43: if nmat=2: B^-1 (A-sB) B A-sB
44: */
45: PetscErrorCode STComputeOperator_Shift(ST st)
46: {
47: st->usesksp = (st->nmat>1)? PETSC_TRUE: PETSC_FALSE;
48: PetscObjectReference((PetscObject)st->A[1]);
49: MatDestroy(&st->T[1]);
50: st->T[1] = st->A[1];
51: STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[0]);
52: if (st->nmat>1) PetscObjectReference((PetscObject)st->T[1]);
53: MatDestroy(&st->P);
54: st->P = (st->nmat>1)? st->T[1]: NULL;
55: st->M = st->T[0];
56: if (st->nmat>1 && st->Psplit) { /* build custom preconditioner from the split matrices */
57: MatDestroy(&st->Pmat);
58: PetscObjectReference((PetscObject)st->Psplit[1]);
59: st->Pmat = st->Psplit[1];
60: }
61: return 0;
62: }
64: PetscErrorCode STSetUp_Shift(ST st)
65: {
66: PetscInt k,nc,nmat=st->nmat;
67: PetscScalar *coeffs=NULL;
69: if (nmat>1) STSetWorkVecs(st,1);
70: if (nmat>2) { /* set-up matrices for polynomial eigenproblems */
71: if (st->transform) {
72: nc = (nmat*(nmat+1))/2;
73: PetscMalloc1(nc,&coeffs);
74: /* Compute coeffs */
75: STCoeffs_Monomial(st,coeffs);
76: /* T[n] = A_n */
77: k = nmat-1;
78: PetscObjectReference((PetscObject)st->A[k]);
79: MatDestroy(&st->T[k]);
80: st->T[k] = st->A[k];
81: for (k=0;k<nmat-1;k++) STMatMAXPY_Private(st,nmat>2?st->sigma:-st->sigma,0.0,k,coeffs?coeffs+((nmat-k)*(nmat-k-1))/2:NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[k]);
82: PetscFree(coeffs);
83: PetscObjectReference((PetscObject)st->T[nmat-1]);
84: MatDestroy(&st->P);
85: st->P = st->T[nmat-1];
86: if (st->Psplit) { /* build custom preconditioner from the split matrices */
87: STMatMAXPY_Private(st,st->sigma,0.0,nmat-1,coeffs?coeffs:NULL,PETSC_TRUE,PETSC_TRUE,&st->Pmat);
88: }
89: ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P);
90: } else {
91: for (k=0;k<nmat;k++) {
92: PetscObjectReference((PetscObject)st->A[k]);
93: MatDestroy(&st->T[k]);
94: st->T[k] = st->A[k];
95: }
96: }
97: }
98: if (st->P) KSPSetUp(st->ksp);
99: return 0;
100: }
102: PetscErrorCode STSetShift_Shift(ST st,PetscScalar newshift)
103: {
104: PetscInt k,nc,nmat=PetscMax(st->nmat,2);
105: PetscScalar *coeffs=NULL;
107: if (st->transform) {
108: if (st->matmode == ST_MATMODE_COPY && nmat>2) {
109: nc = (nmat*(nmat+1))/2;
110: PetscMalloc1(nc,&coeffs);
111: /* Compute coeffs */
112: STCoeffs_Monomial(st,coeffs);
113: }
114: for (k=0;k<nmat-1;k++) STMatMAXPY_Private(st,nmat>2?newshift:-newshift,nmat>2?st->sigma:-st->sigma,k,coeffs?coeffs+((nmat-k)*(nmat-k-1))/2:NULL,PETSC_FALSE,PETSC_FALSE,&st->T[k]);
115: if (st->matmode == ST_MATMODE_COPY && nmat>2) PetscFree(coeffs);
116: if (st->nmat<=2) st->M = st->T[0];
117: }
118: return 0;
119: }
121: SLEPC_EXTERN PetscErrorCode STCreate_Shift(ST st)
122: {
123: st->usesksp = PETSC_TRUE;
125: st->ops->apply = STApply_Generic;
126: st->ops->applytrans = STApplyTranspose_Generic;
127: st->ops->backtransform = STBackTransform_Shift;
128: st->ops->setshift = STSetShift_Shift;
129: st->ops->getbilinearform = STGetBilinearForm_Default;
130: st->ops->setup = STSetUp_Shift;
131: st->ops->computeoperator = STComputeOperator_Shift;
132: st->ops->postsolve = STPostSolve_Shift;
133: st->ops->setdefaultksp = STSetDefaultKSP_Default;
134: return 0;
135: }