Actual source code: precond.c

slepc-3.14.0 2020-09-30
Report Typos and Errors
  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->applymat        = STApplyMat_Generic;
330:   st->ops->applytrans      = STApplyTranspose_Generic;
331:   st->ops->setshift        = STSetShift_Precond;
332:   st->ops->getbilinearform = STGetBilinearForm_Default;
333:   st->ops->setup           = STSetUp_Precond;
334:   st->ops->computeoperator = STComputeOperator_Precond;
335:   st->ops->postsolve       = STPostSolve_Precond;
336:   st->ops->destroy         = STDestroy_Precond;
337:   st->ops->setdefaultksp   = STSetDefaultKSP_Precond;

339:   PetscObjectComposeFunction((PetscObject)st,"STPrecondGetMatForPC_C",STPrecondGetMatForPC_Precond);
340:   PetscObjectComposeFunction((PetscObject)st,"STPrecondSetMatForPC_C",STPrecondSetMatForPC_Precond);
341:   PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",STPrecondGetKSPHasMat_Precond);
342:   PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",STPrecondSetKSPHasMat_Precond);
343:   return(0);
344: }