Actual source code: precond.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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>          /*I "slepcst.h" I*/

 16: typedef struct {
 17:   PetscBool setmat;
 18: } ST_PRECOND;

 20: static PetscErrorCode STSetDefaultKSP_Precond(ST st)
 21: {
 23:   PC             pc;
 24:   PCType         pctype;
 25:   Mat            P;
 26:   PetscBool      t0,t1;

 29:   KSPGetPC(st->ksp,&pc);
 30:   PetscObjectGetType((PetscObject)pc,&pctype);
 31:   STPrecondGetMatForPC(st,&P);
 32:   if (!pctype && st->A && st->A[0]) {
 33:     if (P || st->shift_matrix == 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)?PCJACOBI:PCNONE);
 41:     }
 42:   }
 43:   KSPSetErrorIfNotConverged(st->ksp,PETSC_FALSE);
 44:   return(0);
 45: }

 47: PetscErrorCode STSetUp_Precond(ST st)
 48: {
 49:   Mat            P;
 50:   PC             pc;
 51:   PetscBool      t0,setmat,destroyP=PETSC_FALSE,builtP;

 55:   /* if the user did not set the shift, use the target value */
 56:   if (!st->sigma_set) st->sigma = st->defsigma;

 58:   /* If either pc is none and no matrix has to be set, or pc is shell , exit */
 59:   if (!st->ksp) { STGetKSP(st,&st->ksp); }
 60:   KSPGetPC(st->ksp,&pc);
 61:   PetscObjectTypeCompare((PetscObject)pc,PCSHELL,&t0);
 62:   if (t0) return(0);
 63:   PetscObjectTypeCompare((PetscObject)pc,PCNONE,&t0);
 64:   STPrecondGetKSPHasMat(st,&setmat);
 65:   if (t0 && !setmat) return(0);

 67:   /* Check if a user matrix is set */
 68:   STPrecondGetMatForPC(st,&P);

 70:   /* If not, create A - shift*B */
 71:   if (P) {
 72:     builtP = PETSC_FALSE;
 73:     destroyP = PETSC_TRUE;
 74:     PetscObjectReference((PetscObject)P);
 75:   } else {
 76:     builtP = PETSC_TRUE;

 78:     if (!(PetscAbsScalar(st->sigma) < PETSC_MAX_REAL) && st->nmat>1) {
 79:       P = st->A[1];
 80:       destroyP = PETSC_FALSE;
 81:     } else if (st->sigma == 0.0) {
 82:       P = st->A[0];
 83:       destroyP = PETSC_FALSE;
 84:     } else if (PetscAbsScalar(st->sigma) < PETSC_MAX_REAL && st->shift_matrix != ST_MATMODE_SHELL) {
 85:       if (st->shift_matrix == ST_MATMODE_INPLACE) {
 86:         P = st->A[0];
 87:         destroyP = PETSC_FALSE;
 88:       } else {
 89:         MatDuplicate(st->A[0],MAT_COPY_VALUES,&P);
 90:         destroyP = PETSC_TRUE;
 91:       }
 92:       if (st->nmat>1) {
 93:         MatAXPY(P,-st->sigma,st->A[1],st->str);
 94:       } else {
 95:         MatShift(P,-st->sigma);
 96:       }
 97:       /* TODO: in case of ST_MATMODE_INPLACE should keep the Hermitian flag of st->A and restore at the end */
 98:       STMatSetHermitian(st,P);
 99:     } else builtP = PETSC_FALSE;
100:   }

102:   /* If P was not possible to obtain, set pc to PCNONE */
103:   if (!P) {
104:     PCSetType(pc,PCNONE);

106:     /* If some matrix has to be set to ksp, a shell matrix is created */
107:     if (setmat) {
108:       STMatShellCreate(st,-st->sigma,0,NULL,NULL,&P);
109:       STMatSetHermitian(st,P);
110:       destroyP = PETSC_TRUE;
111:     }
112:   }

114:   KSPSetOperators(st->ksp,setmat?P:NULL,P);

116:   if (destroyP) {
117:     MatDestroy(&P);
118:   } else if (st->shift_matrix == ST_MATMODE_INPLACE && builtP) {
119:     if (st->sigma != 0.0 && PetscAbsScalar(st->sigma) < PETSC_MAX_REAL) {
120:       if (st->nmat>1) {
121:         MatAXPY(st->A[0],st->sigma,st->A[1],st->str);
122:       } else {
123:         MatShift(st->A[0],st->sigma);
124:       }
125:     }
126:   }
127:   return(0);
128: }

130: PetscErrorCode STSetShift_Precond(ST st,PetscScalar newshift)
131: {

135:   /* Nothing to be done if STSetUp has not been called yet */
136:   if (!st->state) return(0);
137:   st->sigma = newshift;
138:   if (st->shift_matrix != ST_MATMODE_SHELL) {
139:     STSetUp_Precond(st);
140:   }
141:   return(0);
142: }

144: static PetscErrorCode STPrecondGetMatForPC_Precond(ST st,Mat *mat)
145: {
147:   PC             pc;
148:   PetscBool      flag;

151:   if (!st->ksp) { STGetKSP(st,&st->ksp); }
152:   KSPGetPC(st->ksp,&pc);
153:   PCGetOperatorsSet(pc,NULL,&flag);
154:   if (flag) {
155:     PCGetOperators(pc,NULL,mat);
156:   } else *mat = NULL;
157:   return(0);
158: }

160: /*@
161:    STPrecondGetMatForPC - Returns the matrix previously set by STPrecondSetMatForPC().

163:    Not Collective, but the Mat is shared by all processors that share the ST

165:    Input Parameter:
166: .  st - the spectral transformation context

168:    Output Parameter:
169: .  mat - the matrix that will be used in constructing the preconditioner or
170:    NULL if no matrix was set by STPrecondSetMatForPC().

172:    Level: advanced

174: .seealso: STPrecondSetMatForPC()
175: @*/
176: PetscErrorCode STPrecondGetMatForPC(ST st,Mat *mat)
177: {

183:   PetscUseMethod(st,"STPrecondGetMatForPC_C",(ST,Mat*),(st,mat));
184:   return(0);
185: }

187: static PetscErrorCode STPrecondSetMatForPC_Precond(ST st,Mat mat)
188: {
189:   PC             pc;
190:   Mat            A;
191:   PetscBool      flag;

195:   if (!st->ksp) { STGetKSP(st,&st->ksp); }
196:   KSPGetPC(st->ksp,&pc);
197:   /* Yes, all these lines are needed to safely set mat as the preconditioner
198:      matrix in pc */
199:   PCGetOperatorsSet(pc,&flag,NULL);
200:   if (flag) {
201:     PCGetOperators(pc,&A,NULL);
202:     PetscObjectReference((PetscObject)A);
203:   } else A = NULL;
204:   PetscObjectReference((PetscObject)mat);
205:   PCSetOperators(pc,A,mat);
206:   MatDestroy(&A);
207:   MatDestroy(&mat);
208:   STPrecondSetKSPHasMat(st,PETSC_TRUE);
209:   return(0);
210: }

212: /*@
213:    STPrecondSetMatForPC - Sets the matrix that must be used to build the preconditioner.

215:    Logically Collective on ST and Mat

217:    Input Parameter:
218: +  st - the spectral transformation context
219: -  mat - the matrix that will be used in constructing the preconditioner

221:    Level: advanced

223:    Notes:
224:    This matrix will be passed to the KSP object (via KSPSetOperators) as
225:    the matrix to be used when constructing the preconditioner.
226:    If no matrix is set or mat is set to NULL, A - sigma*B will
227:    be used to build the preconditioner, being sigma the value set by STSetShift().

229: .seealso: STPrecondSetMatForPC(), STSetShift()
230: @*/
231: PetscErrorCode STPrecondSetMatForPC(ST st,Mat mat)
232: {

239:   PetscTryMethod(st,"STPrecondSetMatForPC_C",(ST,Mat),(st,mat));
240:   return(0);
241: }

243: static PetscErrorCode STPrecondSetKSPHasMat_Precond(ST st,PetscBool setmat)
244: {
245:   ST_PRECOND *data = (ST_PRECOND*)st->data;

248:   data->setmat = setmat;
249:   return(0);
250: }

252: /*@
253:    STPrecondSetKSPHasMat - Sets a flag indicating that during STSetUp the coefficient
254:    matrix of the KSP linear system (A) must be set to be the same matrix as the
255:    preconditioner (P).

257:    Collective on ST

259:    Input Parameter:
260: +  st - the spectral transformation context
261: -  setmat - the flag

263:    Notes:
264:    In most cases, the preconditioner matrix is used only in the PC object, but
265:    in external solvers this matrix must be provided also as the A-matrix in
266:    the KSP object.

268:    Level: developer

270: .seealso: STPrecondGetKSPHasMat(), STSetShift()
271: @*/
272: PetscErrorCode STPrecondSetKSPHasMat(ST st,PetscBool setmat)
273: {

279:   PetscTryMethod(st,"STPrecondSetKSPHasMat_C",(ST,PetscBool),(st,setmat));
280:   return(0);
281: }

283: static PetscErrorCode STPrecondGetKSPHasMat_Precond(ST st,PetscBool *setmat)
284: {
285:   ST_PRECOND *data = (ST_PRECOND*)st->data;

288:   *setmat = data->setmat;
289:   return(0);
290: }

292: /*@
293:    STPrecondGetKSPHasMat - Returns the flag indicating if the coefficient
294:    matrix of the KSP linear system (A) is set to be the same matrix as the
295:    preconditioner (P).

297:    Not Collective

299:    Input Parameter:
300: .  st - the spectral transformation context

302:    Output Parameter:
303: .  setmat - the flag

305:    Level: developer

307: .seealso: STPrecondSetKSPHasMat(), STSetShift()
308: @*/
309: PetscErrorCode STPrecondGetKSPHasMat(ST st,PetscBool *setmat)
310: {

316:   PetscUseMethod(st,"STPrecondGetKSPHasMat_C",(ST,PetscBool*),(st,setmat));
317:   return(0);
318: }

320: PetscErrorCode STDestroy_Precond(ST st)
321: {

325:   PetscFree(st->data);
326:   PetscObjectComposeFunction((PetscObject)st,"STPrecondGetMatForPC_C",NULL);
327:   PetscObjectComposeFunction((PetscObject)st,"STPrecondSetMatForPC_C",NULL);
328:   PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",NULL);
329:   PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",NULL);
330:   return(0);
331: }

333: SLEPC_EXTERN PetscErrorCode STCreate_Precond(ST st)
334: {
336:   ST_PRECOND     *ctx;

339:   PetscNewLog(st,&ctx);
340:   st->data = (void*)ctx;

342:   st->ops->getbilinearform = STGetBilinearForm_Default;
343:   st->ops->setup           = STSetUp_Precond;
344:   st->ops->setshift        = STSetShift_Precond;
345:   st->ops->destroy         = STDestroy_Precond;
346:   st->ops->setdefaultksp   = STSetDefaultKSP_Precond;

348:   PetscObjectComposeFunction((PetscObject)st,"STPrecondGetMatForPC_C",STPrecondGetMatForPC_Precond);
349:   PetscObjectComposeFunction((PetscObject)st,"STPrecondSetMatForPC_C",STPrecondSetMatForPC_Precond);
350:   PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",STPrecondGetKSPHasMat_Precond);
351:   PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",STPrecondSetKSPHasMat_Precond);

353:   STPrecondSetKSPHasMat_Precond(st,PETSC_TRUE);
354:   return(0);
355: }