Actual source code: kspimpl.h

petsc-3.11.2 2019-05-18
Report Typos and Errors

  2: #ifndef _KSPIMPL_H
  3: #define _KSPIMPL_H

  5:  #include <petscksp.h>
  6:  #include <petsc/private/petscimpl.h>

  8: PETSC_EXTERN PetscBool KSPRegisterAllCalled;
  9: PETSC_EXTERN PetscErrorCode KSPRegisterAll(void);
 10: PETSC_EXTERN PetscErrorCode KSPGuessRegisterAll(void);
 11: PETSC_EXTERN PetscErrorCode KSPMatRegisterAll(void);

 13: typedef struct _KSPOps *KSPOps;

 15: struct _KSPOps {
 16:   PetscErrorCode (*buildsolution)(KSP,Vec,Vec*);       /* Returns a pointer to the solution, or
 17:                                                           calculates the solution in a
 18:                                                           user-provided area. */
 19:   PetscErrorCode (*buildresidual)(KSP,Vec,Vec,Vec*);   /* Returns a pointer to the residual, or
 20:                                                           calculates the residual in a
 21:                                                           user-provided area.  */
 22:   PetscErrorCode (*solve)(KSP);                        /* actual solver */
 23:   PetscErrorCode (*setup)(KSP);
 24:   PetscErrorCode (*setfromoptions)(PetscOptionItems*,KSP);
 25:   PetscErrorCode (*publishoptions)(KSP);
 26:   PetscErrorCode (*computeextremesingularvalues)(KSP,PetscReal*,PetscReal*);
 27:   PetscErrorCode (*computeeigenvalues)(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt *);
 28:   PetscErrorCode (*computeritz)(KSP,PetscBool,PetscBool,PetscInt*,Vec[],PetscReal*,PetscReal*);
 29:   PetscErrorCode (*destroy)(KSP);
 30:   PetscErrorCode (*view)(KSP,PetscViewer);
 31:   PetscErrorCode (*reset)(KSP);
 32:   PetscErrorCode (*load)(KSP,PetscViewer);
 33: };

 35: typedef struct _KSPGuessOps *KSPGuessOps;

 37: struct _KSPGuessOps {
 38:   PetscErrorCode (*formguess)(KSPGuess,Vec,Vec); /* Form initial guess */
 39:   PetscErrorCode (*update)(KSPGuess,Vec,Vec);    /* Update database */
 40:   PetscErrorCode (*setfromoptions)(KSPGuess);
 41:   PetscErrorCode (*setup)(KSPGuess);
 42:   PetscErrorCode (*destroy)(KSPGuess);
 43:   PetscErrorCode (*view)(KSPGuess,PetscViewer);
 44:   PetscErrorCode (*reset)(KSPGuess);
 45: };

 47: /*
 48:    Defines the KSPGuess data structure.
 49: */
 50: struct _p_KSPGuess {
 51:   PETSCHEADER(struct _KSPGuessOps);
 52:   KSP              ksp;       /* the parent KSP */
 53:   Mat              A;         /* the current linear operator */
 54:   PetscObjectState omatstate; /* previous linear operator state */
 55:   void             *data;     /* pointer to the specific implementation */
 56: };

 58: PETSC_EXTERN PetscErrorCode KSPGuessCreate_Fischer(KSPGuess);
 59: PETSC_EXTERN PetscErrorCode KSPGuessCreate_POD(KSPGuess);

 61: /*
 62:      Maximum number of monitors you can run with a single KSP
 63: */
 64: #define MAXKSPMONITORS 5
 65: typedef enum {KSP_SETUP_NEW, KSP_SETUP_NEWMATRIX, KSP_SETUP_NEWRHS} KSPSetUpStage;

 67: /*
 68:    Defines the KSP data structure.
 69: */
 70: struct _p_KSP {
 71:   PETSCHEADER(struct _KSPOps);
 72:   DM              dm;
 73:   PetscBool       dmAuto;       /* DM was created automatically by KSP */
 74:   PetscBool       dmActive;     /* KSP should use DM for computing operators */
 75:   /*------------------------- User parameters--------------------------*/
 76:   PetscInt        max_it;                     /* maximum number of iterations */
 77:   KSPGuess        guess;
 78:   PetscBool       guess_zero,                  /* flag for whether initial guess is 0 */
 79:                   calc_sings,                  /* calculate extreme Singular Values */
 80:                   calc_ritz,                   /* calculate (harmonic) Ritz pairs */
 81:                   guess_knoll;                /* use initial guess of PCApply(ksp->B,b */
 82:   PCSide          pc_side;                  /* flag for left, right, or symmetric preconditioning */
 83:   PetscInt        normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */
 84:   PetscReal       rtol,                     /* relative tolerance */
 85:                   abstol,                     /* absolute tolerance */
 86:                   ttol,                     /* (not set by user)  */
 87:                   divtol;                   /* divergence tolerance */
 88:   PetscReal       rnorm0;                   /* initial residual norm (used for divergence testing) */
 89:   PetscReal       rnorm;                    /* current residual norm */
 90:   KSPConvergedReason    reason;
 91:   PetscBool             errorifnotconverged; /* create an error if the KSPSolve() does not converge */

 93:   Vec vec_sol,vec_rhs;            /* pointer to where user has stashed
 94:                                       the solution and rhs, these are
 95:                                       never touched by the code, only
 96:                                       passed back to the user */
 97:   PetscReal     *res_hist;            /* If !0 stores residual at iterations*/
 98:   PetscReal     *res_hist_alloc;      /* If !0 means user did not provide buffer, needs deallocation */
 99:   PetscInt      res_hist_len;         /* current size of residual history array */
100:   PetscInt      res_hist_max;         /* actual amount of data in residual_history */
101:   PetscBool     res_hist_reset;       /* reset history to size zero for each new solve */

103:   PetscInt      chknorm;             /* only compute/check norm if iterations is great than this */
104:   PetscBool     lagnorm;             /* Lag the residual norm calculation so that it is computed as part of the
105:                                         MPI_Allreduce() for computing the inner products for the next iteration. */
106:   /* --------User (or default) routines (most return -1 on error) --------*/
107:   PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */
108:   PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void**);         /* */
109:   void *monitorcontext[MAXKSPMONITORS];                  /* residual calculation, allows user */
110:   PetscInt  numbermonitors;                                   /* to, for instance, print residual norm, etc. */

112:   PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
113:   PetscErrorCode (*convergeddestroy)(void*);
114:   void       *cnvP;

116:   void       *user;             /* optional user-defined context */

118:   PC         pc;

120:   void       *data;                      /* holder for misc stuff associated
121:                                    with a particular iterative solver */

123:   PetscBool         view,   viewPre,   viewReason,   viewMat,   viewPMat,   viewRhs,   viewSol,   viewMatExp,   viewEV,   viewSV,   viewEVExp,   viewFinalRes,   viewPOpExp,   viewDScale;
124:   PetscViewer       viewer, viewerPre, viewerReason, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp, viewerDScale;
125:   PetscViewerFormat format, formatPre, formatReason, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp, formatDScale;

127:   /* ----------------Default work-area management -------------------- */
128:   PetscInt       nwork;
129:   Vec            *work;

131:   KSPSetUpStage  setupstage;
132:   PetscBool      setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */

134:   PetscInt       its;       /* number of iterations so far computed in THIS linear solve*/
135:   PetscInt       totalits;   /* number of iterations used by this KSP object since it was created */

137:   PetscBool      transpose_solve;    /* solve transpose system instead */

139:   KSPNormType    normtype;          /* type of norm used for convergence tests */

141:   PCSide         pc_side_set;   /* PC type set explicitly by user */
142:   KSPNormType    normtype_set;  /* Norm type set explicitly by user */

144:   /*   Allow diagonally scaling the matrix before computing the preconditioner or using
145:        the Krylov method. Note this is NOT just Jacobi preconditioning */

147:   PetscBool    dscale;       /* diagonal scale system; used with KSPSetDiagonalScale() */
148:   PetscBool    dscalefix;    /* unscale system after solve */
149:   PetscBool    dscalefix2;   /* system has been unscaled */
150:   Vec          diagonal;     /* 1/sqrt(diag of matrix) */
151:   Vec          truediagonal;

153:   PetscInt     setfromoptionscalled;
154:   PetscBool    skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */

156:   PetscViewer  eigviewer;   /* Viewer where computed eigenvalues are displayed */

158:   PetscErrorCode (*presolve)(KSP,Vec,Vec,void*);
159:   PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*);
160:   void           *prectx,*postctx;
161: };

163: typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
164:   PetscReal coef;
165:   PetscReal bnrm;
166: } KSPDynTolCtx;

168: typedef struct {
169:   PetscBool  initialrtol;    /* default relative residual decrease is computing from initial residual, not rhs */
170:   PetscBool  mininitialrtol; /* default relative residual decrease is computing from min of initial residual and rhs */
171:   Vec        work;
172: } KSPConvergedDefaultCtx;

174: PETSC_STATIC_INLINE PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm)
175: {

179:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
180:   if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) {
181:     ksp->res_hist[ksp->res_hist_len++] = norm;
182:   }
183:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
184:   return(0);
185: }

187: PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP,PetscBool,KSPNormType*,PCSide*);

189: PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP,PetscInt,const PetscReal*,const PetscReal*);

191: typedef struct _p_DMKSP *DMKSP;
192: typedef struct _DMKSPOps *DMKSPOps;
193: struct _DMKSPOps {
194:   PetscErrorCode (*computeoperators)(KSP,Mat,Mat,void*);
195:   PetscErrorCode (*computerhs)(KSP,Vec,void*);
196:   PetscErrorCode (*computeinitialguess)(KSP,Vec,void*);
197:   PetscErrorCode (*destroy)(DMKSP*);
198:   PetscErrorCode (*duplicate)(DMKSP,DMKSP);
199: };

201: struct _p_DMKSP {
202:   PETSCHEADER(struct _DMKSPOps);
203:   void *operatorsctx;
204:   void *rhsctx;
205:   void *initialguessctx;
206:   void *data;

208:   /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
209:    * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user
210:    * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by
211:    * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates
212:    * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes
213:    * to the original level will no longer propagate to that level.
214:    */
215:   DM originaldm;

217:   void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */
218: };
219: PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM,DMKSP*);
220: PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM,DMKSP*);
221: PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM,DM);

223: /*
224:        These allow the various Krylov methods to apply to either the linear system or its transpose.
225: */
226: PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpace(KSP ksp,Vec y)
227: {
230:   if (ksp->pc_side == PC_LEFT) {
231:     Mat          A;
232:     MatNullSpace nullsp;
233:     PCGetOperators(ksp->pc,&A,NULL);
234:     MatGetNullSpace(A,&nullsp);
235:     if (nullsp) {
236:       MatNullSpaceRemove(nullsp,y);
237:     }
238:   }
239:   return(0);
240: }

242: PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp,Vec y)
243: {
246:   if (ksp->pc_side == PC_LEFT) {
247:     Mat          A;
248:     MatNullSpace nullsp;
249:     PCGetOperators(ksp->pc,&A,NULL);
250:     MatGetTransposeNullSpace(A,&nullsp);
251:     if (nullsp) {
252:       MatNullSpaceRemove(nullsp,y);
253:     }
254:   }
255:   return(0);
256: }

258: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y)
259: {
262:   if (!ksp->transpose_solve) {MatMult(A,x,y);}
263:   else                       {MatMultTranspose(A,x,y);}
264:   return(0);
265: }

267: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y)
268: {
271:   if (!ksp->transpose_solve) {MatMultTranspose(A,x,y);}
272:   else                       {MatMult(A,x,y);}
273:   return(0);
274: }

276: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y)
277: {
280:   if (!ksp->transpose_solve) {
281:     PCApply(ksp->pc,x,y);
282:     KSP_RemoveNullSpace(ksp,y);
283:   } else {
284:     PCApplyTranspose(ksp->pc,x,y);
285:     KSP_RemoveNullSpaceTranspose(ksp,y);
286:   }
287:   return(0);
288: }

290: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y)
291: {
294:   if (!ksp->transpose_solve) {
295:     PCApplyTranspose(ksp->pc,x,y);
296:     KSP_RemoveNullSpaceTranspose(ksp,y);
297:   } else {
298:     PCApply(ksp->pc,x,y);
299:     KSP_RemoveNullSpace(ksp,y);
300:   }
301:   return(0);
302: }

304: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w)
305: {
308:   if (!ksp->transpose_solve) {
309:     PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
310:     KSP_RemoveNullSpace(ksp,y);
311:   } else {
312:     PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
313:     KSP_RemoveNullSpaceTranspose(ksp,y);
314:   }
315:   return(0);
316: }

318: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w)
319: {
322:   if (!ksp->transpose_solve) {
323:     PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
324:   } else {
325:     PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
326:   }
327:   return(0);
328: }

330: PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization;
331: PETSC_EXTERN PetscLogEvent KSP_SetUp;
332: PETSC_EXTERN PetscLogEvent KSP_Solve;
333: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0;
334: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1;
335: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2;
336: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3;
337: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4;
338: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S;
339: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L;
340: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U;

342: PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*);
343: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*);

345: /*MC
346:    KSPCheckDot - Checks if the result of a dot product used by the corresponding KSP contains Inf or NaN. These indicate that the previous 
347:       application of the preconditioner generated an error

349:    Collective on KSP

351:    Input Parameter:
352: .  ksp - the linear solver (KSP) context.

354:    Output Parameter:
355: .  beta - the result of the inner product

357:    Level: developer

359:    Developer Note:
360:    this is used to manage returning from KSP solvers whose preconditioners have failed in some way

362: .keywords: KSP, PC, divergence, convergence

364: .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckNorm(), KSPCheckSolve()
365: M*/
366: #define KSPCheckDot(ksp,beta)           \
367:   if (PetscIsInfOrNanScalar(beta)) { \
368:     if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");\
369:     else {\
371:       PCFailedReason pcreason;\
372:       PetscInt       sendbuf,pcreason_max; \
373:       PCGetFailedReason(ksp->pc,&pcreason);\
374:       sendbuf = (PetscInt)pcreason; \
375:       MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp)); \
376:       if (pcreason_max) {\
377:         ksp->reason = KSP_DIVERGED_PC_FAILED;\
378:         VecSetInf(ksp->vec_sol);\
379:       } else {\
380:         ksp->reason = KSP_DIVERGED_NANORINF;\
381:       }\
382:       return(0);\
383:     }\
384:   }

386: /*MC
387:    KSPCheckNorm - Checks if the result of a norm used by the corresponding KSP contains Inf or NaN. These indicate that the previous
388:       application of the preconditioner generated an error

390:    Collective on KSP

392:    Input Parameter:
393: .  ksp - the linear solver (KSP) context.

395:    Output Parameter:
396: .  beta - the result of the norm

398:    Level: developer

400:    Developer Note:
401:    this is used to manage returning from KSP solvers whose preconditioners have failed in some way

403: .keywords: KSP, PC, divergence, convergence

405: .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckDot(), KSPCheckSolve()
406: M*/
407: #define KSPCheckNorm(ksp,beta)           \
408:   if (PetscIsInfOrNanReal(beta)) { \
409:     if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf norm");\
410:     else {\
412:       PCFailedReason pcreason;\
413:       PetscInt       sendbuf,pcreason_max; \
414:       PCGetFailedReason(ksp->pc,&pcreason);\
415:       sendbuf = (PetscInt)pcreason; \
416:       MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp)); \
417:       if (pcreason_max) {\
418:         ksp->reason = KSP_DIVERGED_PC_FAILED;\
419:         VecSetInf(ksp->vec_sol);\
420:       } else {\
421:         ksp->reason = KSP_DIVERGED_NANORINF;\
422:       }\
423:       return(0);\
424:     }\
425:   }

427: #endif