Actual source code: precond.c
slepc-3.13.3 2020-06-14
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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: PetscBool usermat; /* the user has set a matrix */
19: Mat mat; /* user-provided matrix */
20: } ST_PRECOND;
22: static PetscErrorCode STSetDefaultKSP_Precond(ST st)
23: {
25: PC pc;
26: PCType pctype;
27: PetscBool t0,t1;
30: KSPGetPC(st->ksp,&pc);
31: PCGetType(pc,&pctype);
32: if (!pctype && st->A && st->A[0]) {
33: if (st->matmode == ST_MATMODE_SHELL) {
34: PCSetType(pc,PCJACOBI);
35: } else {
36: MatHasOperation(st->A[0],MATOP_DUPLICATE,&t0);
37: if (st->nmat>1) {
38: MatHasOperation(st->A[0],MATOP_AXPY,&t1);
39: } else t1 = PETSC_TRUE;
40: PCSetType(pc,(t0 && t1)?PCBJACOBI:PCNONE);
41: }
42: }
43: KSPSetErrorIfNotConverged(st->ksp,PETSC_FALSE);
44: return(0);
45: }
47: PetscErrorCode STPostSolve_Precond(ST st)
48: {
50: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
53: if (st->matmode == ST_MATMODE_INPLACE && !(ctx->mat || (PetscAbsScalar(st->sigma)>=PETSC_MAX_REAL && st->nmat>1))) {
54: if (st->nmat>1) {
55: MatAXPY(st->A[0],st->sigma,st->A[1],st->str);
56: } else {
57: MatShift(st->A[0],st->sigma);
58: }
59: st->Astate[0] = ((PetscObject)st->A[0])->state;
60: st->state = ST_STATE_INITIAL;
61: st->opready = PETSC_FALSE;
62: }
63: return(0);
64: }
66: /*
67: Operator (precond):
68: Op P M
69: if nmat=1: --- A-sI NULL
70: if nmat=2: --- A-sB NULL
71: */
72: PetscErrorCode STComputeOperator_Precond(ST st)
73: {
75: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
78: /* if the user did not set the shift, use the target value */
79: if (!st->sigma_set) st->sigma = st->defsigma;
80: st->M = NULL;
82: /* P = A-sigma*B */
83: if (ctx->mat) {
84: PetscObjectReference((PetscObject)ctx->mat);
85: MatDestroy(&st->P);
86: st->P = ctx->mat;
87: } else {
88: PetscObjectReference((PetscObject)st->A[1]);
89: MatDestroy(&st->T[0]);
90: st->T[0] = st->A[1];
91: if (!(PetscAbsScalar(st->sigma) < PETSC_MAX_REAL) && st->nmat>1) {
92: PetscObjectReference((PetscObject)st->T[0]);
93: MatDestroy(&st->P);
94: st->P = st->T[0];
95: } else if (PetscAbsScalar(st->sigma)<PETSC_MAX_REAL) {
96: STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),&st->T[1]);
97: PetscObjectReference((PetscObject)st->T[1]);
98: MatDestroy(&st->P);
99: st->P = st->T[1];
100: }
101: }
102: return(0);
103: }
105: PetscErrorCode STSetUp_Precond(ST st)
106: {
108: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
111: if (st->P) {
112: KSPSetOperators(st->ksp,ctx->ksphasmat?st->P:NULL,st->P);
113: /* NOTE: we do not call KSPSetUp() here because some eigensolvers such as JD require a lazy setup */
114: }
115: return(0);
116: }
118: PetscErrorCode STSetShift_Precond(ST st,PetscScalar newshift)
119: {
121: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
124: if (st->transform && !ctx->mat) {
125: STMatMAXPY_Private(st,-newshift,-st->sigma,0,NULL,PETSC_FALSE,&st->T[1]);
126: if (st->P!=st->T[1]) {
127: PetscObjectReference((PetscObject)st->T[1]);
128: MatDestroy(&st->P);
129: st->P = st->T[1];
130: }
131: }
132: if (st->P) {
133: if (!st->ksp) { STGetKSP(st,&st->ksp); }
134: KSPSetOperators(st->ksp,ctx->ksphasmat?st->P:NULL,st->P);
135: }
136: return(0);
137: }
139: static PetscErrorCode STPrecondGetMatForPC_Precond(ST st,Mat *mat)
140: {
141: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
144: *mat = ctx->mat;
145: return(0);
146: }
148: /*@
149: STPrecondGetMatForPC - Returns the matrix previously set by STPrecondSetMatForPC().
151: Not Collective, but the Mat is shared by all processors that share the ST
153: Input Parameter:
154: . st - the spectral transformation context
156: Output Parameter:
157: . mat - the matrix that will be used in constructing the preconditioner or
158: NULL if no matrix was set by STPrecondSetMatForPC().
160: Level: advanced
162: .seealso: STPrecondSetMatForPC()
163: @*/
164: PetscErrorCode STPrecondGetMatForPC(ST st,Mat *mat)
165: {
171: PetscUseMethod(st,"STPrecondGetMatForPC_C",(ST,Mat*),(st,mat));
172: return(0);
173: }
175: static PetscErrorCode STPrecondSetMatForPC_Precond(ST st,Mat mat)
176: {
177: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
181: STCheckNotSeized(st,1);
182: PetscObjectReference((PetscObject)mat);
183: MatDestroy(&ctx->mat);
184: ctx->mat = mat;
185: ctx->usermat = mat? PETSC_TRUE: PETSC_FALSE;
186: st->state = ST_STATE_INITIAL;
187: st->opready = PETSC_FALSE;
188: return(0);
189: }
191: /*@
192: STPrecondSetMatForPC - Sets the matrix that must be used to build the preconditioner.
194: Logically Collective on st
196: Input Parameter:
197: + st - the spectral transformation context
198: - mat - the matrix that will be used in constructing the preconditioner
200: Level: advanced
202: Notes:
203: This matrix will be passed to the KSP object (via KSPSetOperators) as
204: the matrix to be used when constructing the preconditioner.
205: If no matrix is set or mat is set to NULL, A - sigma*B will
206: be used to build the preconditioner, being sigma the value set by STSetShift().
208: .seealso: STPrecondSetMatForPC(), STSetShift()
209: @*/
210: PetscErrorCode STPrecondSetMatForPC(ST st,Mat mat)
211: {
218: PetscTryMethod(st,"STPrecondSetMatForPC_C",(ST,Mat),(st,mat));
219: return(0);
220: }
222: static PetscErrorCode STPrecondSetKSPHasMat_Precond(ST st,PetscBool ksphasmat)
223: {
224: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
227: if (ctx->ksphasmat != ksphasmat) {
228: ctx->ksphasmat = ksphasmat;
229: st->state = ST_STATE_INITIAL;
230: }
231: return(0);
232: }
234: /*@
235: STPrecondSetKSPHasMat - Sets a flag indicating that during STSetUp the coefficient
236: matrix of the KSP linear solver (A) must be set to be the same matrix as the
237: preconditioner (P).
239: Collective on st
241: Input Parameter:
242: + st - the spectral transformation context
243: - ksphasmat - the flag
245: Notes:
246: Often, the preconditioner matrix is used only in the PC object, but
247: in some solvers this matrix must be provided also as the A-matrix in
248: the KSP object.
250: Level: developer
252: .seealso: STPrecondGetKSPHasMat(), STSetShift()
253: @*/
254: PetscErrorCode STPrecondSetKSPHasMat(ST st,PetscBool ksphasmat)
255: {
261: PetscTryMethod(st,"STPrecondSetKSPHasMat_C",(ST,PetscBool),(st,ksphasmat));
262: return(0);
263: }
265: static PetscErrorCode STPrecondGetKSPHasMat_Precond(ST st,PetscBool *ksphasmat)
266: {
267: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
270: *ksphasmat = ctx->ksphasmat;
271: return(0);
272: }
274: /*@
275: STPrecondGetKSPHasMat - Returns the flag indicating if the coefficient
276: matrix of the KSP linear system (A) is set to be the same matrix as the
277: preconditioner (P).
279: Not Collective
281: Input Parameter:
282: . st - the spectral transformation context
284: Output Parameter:
285: . ksphasmat - the flag
287: Level: developer
289: .seealso: STPrecondSetKSPHasMat(), STSetShift()
290: @*/
291: PetscErrorCode STPrecondGetKSPHasMat(ST st,PetscBool *ksphasmat)
292: {
298: PetscUseMethod(st,"STPrecondGetKSPHasMat_C",(ST,PetscBool*),(st,ksphasmat));
299: return(0);
300: }
302: PetscErrorCode STDestroy_Precond(ST st)
303: {
304: ST_PRECOND *ctx = (ST_PRECOND*)st->data;
308: MatDestroy(&ctx->mat);
309: PetscFree(st->data);
310: PetscObjectComposeFunction((PetscObject)st,"STPrecondGetMatForPC_C",NULL);
311: PetscObjectComposeFunction((PetscObject)st,"STPrecondSetMatForPC_C",NULL);
312: PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",NULL);
313: PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",NULL);
314: return(0);
315: }
317: SLEPC_EXTERN PetscErrorCode STCreate_Precond(ST st)
318: {
320: ST_PRECOND *ctx;
323: PetscNewLog(st,&ctx);
324: st->data = (void*)ctx;
326: st->usesksp = PETSC_TRUE;
328: st->ops->apply = STApply_Generic;
329: st->ops->applytrans = STApplyTranspose_Generic;
330: st->ops->setshift = STSetShift_Precond;
331: st->ops->getbilinearform = STGetBilinearForm_Default;
332: st->ops->setup = STSetUp_Precond;
333: st->ops->computeoperator = STComputeOperator_Precond;
334: st->ops->postsolve = STPostSolve_Precond;
335: st->ops->destroy = STDestroy_Precond;
336: st->ops->setdefaultksp = STSetDefaultKSP_Precond;
338: PetscObjectComposeFunction((PetscObject)st,"STPrecondGetMatForPC_C",STPrecondGetMatForPC_Precond);
339: PetscObjectComposeFunction((PetscObject)st,"STPrecondSetMatForPC_C",STPrecondSetMatForPC_Precond);
340: PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",STPrecondGetKSPHasMat_Precond);
341: PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",STPrecondSetKSPHasMat_Precond);
342: return(0);
343: }