Actual source code: precond.c

slepc-3.18.3 2023-03-24
Report Typos and Errors
  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:    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: {
 22:   PC             pc;
 23:   PCType         pctype;
 24:   PetscBool      t0,t1;

 26:   KSPGetPC(st->ksp,&pc);
 27:   PCGetType(pc,&pctype);
 28:   if (!pctype && st->A && st->A[0]) {
 29:     if (st->matmode == ST_MATMODE_SHELL) PCSetType(pc,PCJACOBI);
 30:     else {
 31:       MatHasOperation(st->A[0],MATOP_DUPLICATE,&t0);
 32:       if (st->nmat>1) MatHasOperation(st->A[0],MATOP_AXPY,&t1);
 33:       else t1 = PETSC_TRUE;
 34:       PCSetType(pc,(t0 && t1)?PCBJACOBI:PCNONE);
 35:     }
 36:   }
 37:   KSPSetErrorIfNotConverged(st->ksp,PETSC_FALSE);
 38:   return 0;
 39: }

 41: PetscErrorCode STPostSolve_Precond(ST st)
 42: {
 43:   if (st->matmode == ST_MATMODE_INPLACE && !(st->Pmat || (PetscAbsScalar(st->sigma)>=PETSC_MAX_REAL && st->nmat>1))) {
 44:     if (st->nmat>1) MatAXPY(st->A[0],st->sigma,st->A[1],st->str);
 45:     else MatShift(st->A[0],st->sigma);
 46:     st->Astate[0] = ((PetscObject)st->A[0])->state;
 47:     st->state   = ST_STATE_INITIAL;
 48:     st->opready = PETSC_FALSE;
 49:   }
 50:   return 0;
 51: }

 53: /*
 54:    Operator (precond):
 55:                Op        P         M
 56:    if nmat=1:  ---       A-sI      NULL
 57:    if nmat=2:  ---       A-sB      NULL
 58: */
 59: PetscErrorCode STComputeOperator_Precond(ST st)
 60: {
 61:   /* if the user did not set the shift, use the target value */
 62:   if (!st->sigma_set) st->sigma = st->defsigma;
 63:   st->M = NULL;

 65:   /* build custom preconditioner from the split matrices */
 66:   if (st->Psplit) {
 67:     if (!(PetscAbsScalar(st->sigma) < PETSC_MAX_REAL) && st->nmat>1) {
 68:       PetscObjectReference((PetscObject)st->Psplit[0]);
 69:       MatDestroy(&st->Pmat);
 70:       st->Pmat = st->Psplit[0];
 71:     } else if (PetscAbsScalar(st->sigma)<PETSC_MAX_REAL) STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PETSC_TRUE,PETSC_TRUE,&st->Pmat);
 72:   }

 74:   /* P = A-sigma*B */
 75:   if (st->Pmat) {
 76:     PetscObjectReference((PetscObject)st->Pmat);
 77:     MatDestroy(&st->P);
 78:     st->P = st->Pmat;
 79:   } else {
 80:     PetscObjectReference((PetscObject)st->A[1]);
 81:     MatDestroy(&st->T[0]);
 82:     st->T[0] = st->A[1];
 83:     if (!(PetscAbsScalar(st->sigma) < PETSC_MAX_REAL) && st->nmat>1) {
 84:       PetscObjectReference((PetscObject)st->T[0]);
 85:       MatDestroy(&st->P);
 86:       st->P = st->T[0];
 87:     } else if (PetscAbsScalar(st->sigma)<PETSC_MAX_REAL) {
 88:       STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[1]);
 89:       PetscObjectReference((PetscObject)st->T[1]);
 90:       MatDestroy(&st->P);
 91:       st->P = st->T[1];
 92:     }
 93:   }
 94:   return 0;
 95: }

 97: PetscErrorCode STSetUp_Precond(ST st)
 98: {
 99:   ST_PRECOND     *ctx = (ST_PRECOND*)st->data;

101:   if (st->P) {
102:     ST_KSPSetOperators(st,ctx->ksphasmat?st->P:NULL,st->P);
103:     /* NOTE: we do not call KSPSetUp() here because some eigensolvers such as JD require a lazy setup */
104:   }
105:   return 0;
106: }

108: PetscErrorCode STSetShift_Precond(ST st,PetscScalar newshift)
109: {
110:   ST_PRECOND     *ctx = (ST_PRECOND*)st->data;

112:   if (st->Psplit) { /* update custom preconditioner from the split matrices */
113:     if (PetscAbsScalar(st->sigma)<PETSC_MAX_REAL || st->nmat==1) STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PETSC_FALSE,PETSC_TRUE,&st->Pmat);
114:   }
115:   if (st->transform && !st->Pmat) {
116:     STMatMAXPY_Private(st,-newshift,-st->sigma,0,NULL,PETSC_FALSE,PETSC_FALSE,&st->T[1]);
117:     if (st->P!=st->T[1]) {
118:       PetscObjectReference((PetscObject)st->T[1]);
119:       MatDestroy(&st->P);
120:       st->P = st->T[1];
121:     }
122:   }
123:   if (st->P) ST_KSPSetOperators(st,ctx->ksphasmat?st->P:NULL,st->P);
124:   return 0;
125: }

127: static PetscErrorCode STPrecondSetKSPHasMat_Precond(ST st,PetscBool ksphasmat)
128: {
129:   ST_PRECOND *ctx = (ST_PRECOND*)st->data;

131:   if (ctx->ksphasmat != ksphasmat) {
132:     ctx->ksphasmat = ksphasmat;
133:     st->state      = ST_STATE_INITIAL;
134:   }
135:   return 0;
136: }

138: /*@
139:    STPrecondSetKSPHasMat - Sets a flag indicating that during STSetUp the coefficient
140:    matrix of the KSP linear solver (A) must be set to be the same matrix as the
141:    preconditioner (P).

143:    Collective on st

145:    Input Parameters:
146: +  st - the spectral transformation context
147: -  ksphasmat - the flag

149:    Notes:
150:    Often, the preconditioner matrix is used only in the PC object, but
151:    in some solvers this matrix must be provided also as the A-matrix in
152:    the KSP object.

154:    Level: developer

156: .seealso: STPrecondGetKSPHasMat(), STSetShift()
157: @*/
158: PetscErrorCode STPrecondSetKSPHasMat(ST st,PetscBool ksphasmat)
159: {
162:   PetscTryMethod(st,"STPrecondSetKSPHasMat_C",(ST,PetscBool),(st,ksphasmat));
163:   return 0;
164: }

166: static PetscErrorCode STPrecondGetKSPHasMat_Precond(ST st,PetscBool *ksphasmat)
167: {
168:   ST_PRECOND *ctx = (ST_PRECOND*)st->data;

170:   *ksphasmat = ctx->ksphasmat;
171:   return 0;
172: }

174: /*@
175:    STPrecondGetKSPHasMat - Returns the flag indicating if the coefficient
176:    matrix of the KSP linear system (A) is set to be the same matrix as the
177:    preconditioner (P).

179:    Not Collective

181:    Input Parameter:
182: .  st - the spectral transformation context

184:    Output Parameter:
185: .  ksphasmat - the flag

187:    Level: developer

189: .seealso: STPrecondSetKSPHasMat(), STSetShift()
190: @*/
191: PetscErrorCode STPrecondGetKSPHasMat(ST st,PetscBool *ksphasmat)
192: {
195:   PetscUseMethod(st,"STPrecondGetKSPHasMat_C",(ST,PetscBool*),(st,ksphasmat));
196:   return 0;
197: }

199: PetscErrorCode STDestroy_Precond(ST st)
200: {
201:   PetscFree(st->data);
202:   PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",NULL);
203:   PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",NULL);
204:   return 0;
205: }

207: SLEPC_EXTERN PetscErrorCode STCreate_Precond(ST st)
208: {
209:   ST_PRECOND     *ctx;

211:   PetscNew(&ctx);
212:   st->data = (void*)ctx;

214:   st->usesksp = PETSC_TRUE;

216:   st->ops->apply           = STApply_Generic;
217:   st->ops->applymat        = STApplyMat_Generic;
218:   st->ops->applytrans      = STApplyTranspose_Generic;
219:   st->ops->setshift        = STSetShift_Precond;
220:   st->ops->getbilinearform = STGetBilinearForm_Default;
221:   st->ops->setup           = STSetUp_Precond;
222:   st->ops->computeoperator = STComputeOperator_Precond;
223:   st->ops->postsolve       = STPostSolve_Precond;
224:   st->ops->destroy         = STDestroy_Precond;
225:   st->ops->setdefaultksp   = STSetDefaultKSP_Precond;

227:   PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",STPrecondGetKSPHasMat_Precond);
228:   PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",STPrecondSetKSPHasMat_Precond);
229:   return 0;
230: }