Actual source code: precond.c
slepc-3.16.0 2021-09-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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: Implements the ST class for preconditioned eigenvalue methods
12: */
14: #include <slepc/private/stimpl.h>
16: typedef struct {
17: PetscBool ksphasmat; /* the KSP must have the same matrix as PC */
18: } ST_PRECOND;
20: static PetscErrorCode STSetDefaultKSP_Precond(ST st)
21: {
23: PC pc;
24: PCType pctype;
25: PetscBool t0,t1;
28: KSPGetPC(st->ksp,&pc);
29: PCGetType(pc,&pctype);
30: if (!pctype && st->A && st->A[0]) {
31: if (st->matmode == ST_MATMODE_SHELL) {
32: PCSetType(pc,PCJACOBI);
33: } else {
34: MatHasOperation(st->A[0],MATOP_DUPLICATE,&t0);
35: if (st->nmat>1) {
36: MatHasOperation(st->A[0],MATOP_AXPY,&t1);
37: } else t1 = PETSC_TRUE;
38: PCSetType(pc,(t0 && t1)?PCBJACOBI:PCNONE);
39: }
40: }
41: KSPSetErrorIfNotConverged(st->ksp,PETSC_FALSE);
42: return(0);
43: }
45: PetscErrorCode STPostSolve_Precond(ST st)
46: {
50: if (st->matmode == ST_MATMODE_INPLACE && !(st->Pmat || (PetscAbsScalar(st->sigma)>=PETSC_MAX_REAL && st->nmat>1))) {
51: if (st->nmat>1) {
52: MatAXPY(st->A[0],st->sigma,st->A[1],st->str);
53: } else {
54: MatShift(st->A[0],st->sigma);
55: }
56: st->Astate[0] = ((PetscObject)st->A[0])->state;
57: st->state = ST_STATE_INITIAL;
58: st->opready = PETSC_FALSE;
59: }
60: return(0);
61: }
63: /*
64: Operator (precond):
65: Op P M
66: if nmat=1: --- A-sI NULL
67: if nmat=2: --- A-sB NULL
68: */
69: PetscErrorCode STComputeOperator_Precond(ST st)
70: {
74: /* if the user did not set the shift, use the target value */
75: if (!st->sigma_set) st->sigma = st->defsigma;
76: st->M = NULL;
78: /* P = A-sigma*B */
79: if (st->Pmat) {
80: PetscObjectReference((PetscObject)st->Pmat);
81: MatDestroy(&st->P);
82: st->P = st->Pmat;
83: } else {
84: PetscObjectReference((PetscObject)st->A[1]);
85: MatDestroy(&st->T[0]);
86: st->T[0] = st->A[1];
87: if (!(PetscAbsScalar(st->sigma) < PETSC_MAX_REAL) && st->nmat>1) {
88: PetscObjectReference((PetscObject)st->T[0]);
89: MatDestroy(&st->P);
90: st->P = st->T[0];
91: } else if (PetscAbsScalar(st->sigma)<PETSC_MAX_REAL) {
92: STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),&st->T[1]);
93: PetscObjectReference((PetscObject)st->T[1]);
94: MatDestroy(&st->P);
95: st->P = st->T[1];
96: }
97: }
98: return(0);
99: }
101: PetscErrorCode STSetUp_Precond(ST st)
102: {
104: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
107: if (st->P) {
108: STKSPSetOperators(st,ctx->ksphasmat?st->P:NULL,st->P);
109: /* NOTE: we do not call KSPSetUp() here because some eigensolvers such as JD require a lazy setup */
110: }
111: return(0);
112: }
114: PetscErrorCode STSetShift_Precond(ST st,PetscScalar newshift)
115: {
117: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
120: if (st->transform && !st->Pmat) {
121: STMatMAXPY_Private(st,-newshift,-st->sigma,0,NULL,PETSC_FALSE,&st->T[1]);
122: if (st->P!=st->T[1]) {
123: PetscObjectReference((PetscObject)st->T[1]);
124: MatDestroy(&st->P);
125: st->P = st->T[1];
126: }
127: }
128: if (st->P) {
129: STKSPSetOperators(st,ctx->ksphasmat?st->P:NULL,st->P);
130: }
131: return(0);
132: }
134: static PetscErrorCode STPrecondSetKSPHasMat_Precond(ST st,PetscBool ksphasmat)
135: {
136: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
139: if (ctx->ksphasmat != ksphasmat) {
140: ctx->ksphasmat = ksphasmat;
141: st->state = ST_STATE_INITIAL;
142: }
143: return(0);
144: }
146: /*@
147: STPrecondSetKSPHasMat - Sets a flag indicating that during STSetUp the coefficient
148: matrix of the KSP linear solver (A) must be set to be the same matrix as the
149: preconditioner (P).
151: Collective on st
153: Input Parameters:
154: + st - the spectral transformation context
155: - ksphasmat - the flag
157: Notes:
158: Often, the preconditioner matrix is used only in the PC object, but
159: in some solvers this matrix must be provided also as the A-matrix in
160: the KSP object.
162: Level: developer
164: .seealso: STPrecondGetKSPHasMat(), STSetShift()
165: @*/
166: PetscErrorCode STPrecondSetKSPHasMat(ST st,PetscBool ksphasmat)
167: {
173: PetscTryMethod(st,"STPrecondSetKSPHasMat_C",(ST,PetscBool),(st,ksphasmat));
174: return(0);
175: }
177: static PetscErrorCode STPrecondGetKSPHasMat_Precond(ST st,PetscBool *ksphasmat)
178: {
179: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
182: *ksphasmat = ctx->ksphasmat;
183: return(0);
184: }
186: /*@
187: STPrecondGetKSPHasMat - Returns the flag indicating if the coefficient
188: matrix of the KSP linear system (A) is set to be the same matrix as the
189: preconditioner (P).
191: Not Collective
193: Input Parameter:
194: . st - the spectral transformation context
196: Output Parameter:
197: . ksphasmat - the flag
199: Level: developer
201: .seealso: STPrecondSetKSPHasMat(), STSetShift()
202: @*/
203: PetscErrorCode STPrecondGetKSPHasMat(ST st,PetscBool *ksphasmat)
204: {
210: PetscUseMethod(st,"STPrecondGetKSPHasMat_C",(ST,PetscBool*),(st,ksphasmat));
211: return(0);
212: }
214: PetscErrorCode STDestroy_Precond(ST st)
215: {
219: PetscFree(st->data);
220: PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",NULL);
221: PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",NULL);
222: return(0);
223: }
225: SLEPC_EXTERN PetscErrorCode STCreate_Precond(ST st)
226: {
228: ST_PRECOND *ctx;
231: PetscNewLog(st,&ctx);
232: st->data = (void*)ctx;
234: st->usesksp = PETSC_TRUE;
236: st->ops->apply = STApply_Generic;
237: st->ops->applymat = STApplyMat_Generic;
238: st->ops->applytrans = STApplyTranspose_Generic;
239: st->ops->setshift = STSetShift_Precond;
240: st->ops->getbilinearform = STGetBilinearForm_Default;
241: st->ops->setup = STSetUp_Precond;
242: st->ops->computeoperator = STComputeOperator_Precond;
243: st->ops->postsolve = STPostSolve_Precond;
244: st->ops->destroy = STDestroy_Precond;
245: st->ops->setdefaultksp = STSetDefaultKSP_Precond;
247: PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",STPrecondGetKSPHasMat_Precond);
248: PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",STPrecondSetKSPHasMat_Precond);
249: return(0);
250: }