Actual source code: kspimpl.h

  1: #pragma once

  3: #include <petscksp.h>
  4: #include <petscds.h>
  5: #include <petsc/private/petscimpl.h>

  7: /* SUBMANSEC = KSP */

  9: PETSC_EXTERN PetscBool      KSPRegisterAllCalled;
 10: PETSC_EXTERN PetscBool      KSPMonitorRegisterAllCalled;
 11: PETSC_EXTERN PetscErrorCode KSPRegisterAll(void);
 12: PETSC_EXTERN PetscErrorCode KSPMonitorRegisterAll(void);
 13: PETSC_EXTERN PetscErrorCode KSPGuessRegisterAll(void);
 14: PETSC_EXTERN PetscErrorCode KSPMatRegisterAll(void);

 16: typedef struct _KSPOps *KSPOps;

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

 39: typedef struct _KSPGuessOps *KSPGuessOps;

 41: struct _KSPGuessOps {
 42:   PetscErrorCode (*formguess)(KSPGuess, Vec, Vec); /* Form initial guess */
 43:   PetscErrorCode (*update)(KSPGuess, Vec, Vec);    /* Update database */
 44:   PetscErrorCode (*setfromoptions)(KSPGuess);
 45:   PetscErrorCode (*settolerance)(KSPGuess, PetscReal);
 46:   PetscErrorCode (*setup)(KSPGuess);
 47:   PetscErrorCode (*destroy)(KSPGuess);
 48:   PetscErrorCode (*view)(KSPGuess, PetscViewer);
 49:   PetscErrorCode (*reset)(KSPGuess);
 50: };

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

 63: PETSC_EXTERN PetscErrorCode KSPGuessCreate_Fischer(KSPGuess);
 64: PETSC_EXTERN PetscErrorCode KSPGuessCreate_POD(KSPGuess);

 66: /*
 67:      Maximum number of monitors you can run with a single KSP
 68: */
 69: #define MAXKSPMONITORS    5
 70: #define MAXKSPREASONVIEWS 5
 71: typedef enum {
 72:   KSP_SETUP_NEW = 0,
 73:   KSP_SETUP_NEWMATRIX,
 74:   KSP_SETUP_NEWRHS
 75: } KSPSetUpStage;

 77: /*
 78:    Defines the KSP data structure.
 79: */
 80: struct _p_KSP {
 81:   PETSCHEADER(struct _KSPOps);
 82:   DM          dm;
 83:   PetscBool   dmAuto;   /* DM was created automatically by KSP */
 84:   KSPDMActive dmActive; /* KSP should use DM for computing operators */
 85:   /*------------------------- User parameters--------------------------*/
 86:   PetscObjectParameterDeclare(PetscInt, max_it); /* maximum number of iterations */
 87:   PetscInt  min_it;                              /* minimum number of iterations */
 88:   KSPGuess  guess;
 89:   PetscBool guess_zero,                                 /* flag for whether initial guess is 0 */
 90:     guess_not_read,                                     /* guess is not read, does not need to be zeroed */
 91:     calc_sings,                                         /* calculate extreme Singular Values */
 92:     calc_ritz,                                          /* calculate (harmonic) Ritz pairs */
 93:     guess_knoll;                                        /* use initial guess of PCApply(ksp->B,b */
 94:   PCSide   pc_side;                                     /* flag for left, right, or symmetric preconditioning */
 95:   PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */
 96:   PetscObjectParameterDeclare(PetscReal, rtol);         /* relative tolerance */
 97:   PetscObjectParameterDeclare(PetscReal, abstol);       /* absolute tolerance */
 98:   PetscObjectParameterDeclare(PetscReal, ttol);         /* (not set by user)  */
 99:   PetscObjectParameterDeclare(PetscReal, divtol);       /* divergence tolerance */
100:   PetscReal          rnorm0;                            /* initial residual norm (used for divergence testing) */
101:   PetscReal          rnorm;                             /* current residual norm */
102:   KSPConvergedReason reason;
103:   PetscBool          errorifnotconverged; /* create an error if the KSPSolve() does not converge */

105:   Vec        vec_sol, vec_rhs; /* pointer to where user has stashed
106:                                       the solution and rhs, these are
107:                                       never touched by the code, only
108:                                       passed back to the user */
109:   PetscReal *res_hist;         /* If !0 stores residual each at iteration */
110:   PetscReal *res_hist_alloc;   /* If !0 means user did not provide buffer, needs deallocation */
111:   PetscCount res_hist_len;     /* current entry count of residual history array */
112:   PetscCount res_hist_max;     /* total entry count of storage in residual history */
113:   PetscBool  res_hist_reset;   /* reset history to length zero for each new solve */
114:   PetscReal *err_hist;         /* If !0 stores error at each iteration */
115:   PetscReal *err_hist_alloc;   /* If !0 means user did not provide buffer, needs deallocation */
116:   PetscCount err_hist_len;     /* current entry count of error history array */
117:   PetscCount err_hist_max;     /* total entry count of storage in error history */
118:   PetscBool  err_hist_reset;   /* reset history to length zero for each new solve */

120:   PetscInt  chknorm; /* only compute/check norm if iterations is great than this */
121:   PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the
122:                                         MPI_Allreduce() for computing the inner products for the next iteration. */

124:   PetscInt nmax; /* maximum number of right-hand sides to be handled simultaneously */

126:   /* --------User (or default) routines (most return -1 on error) --------*/
127:   KSPMonitorFn      *monitor[MAXKSPMONITORS];
128:   PetscCtxDestroyFn *monitordestroy[MAXKSPMONITORS];
129:   void              *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */
130:   PetscInt           numbermonitors;                 /* to, for instance, print residual norm, etc. */
131:   PetscBool          pauseFinal;                     /* Pause all drawing monitor at the final iterate */

133:   PetscViewer               convergedreasonviewer;
134:   PetscViewerFormat         convergedreasonformat;
135:   KSPConvergedReasonViewFn *reasonview[MAXKSPREASONVIEWS];        /* KSP converged reason view */
136:   PetscCtxDestroyFn        *reasonviewdestroy[MAXKSPREASONVIEWS]; /* optional destroy routine */
137:   void                     *reasonviewcontext[MAXKSPREASONVIEWS]; /* viewer context */
138:   PetscInt                  numberreasonviews;                    /* current number of reason viewers */

140:   KSPConvergenceTestFn *converged;
141:   PetscCtxDestroyFn    *convergeddestroy;
142:   void                 *cnvP;

144:   PetscCtx ctx; /* optional user-defined context */

146:   PC pc;

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

150:   PetscBool         view, viewPre, viewRate, viewMat, viewPMat, viewRhs, viewSol, viewMatExp, viewEV, viewSV, viewEVExp, viewFinalRes, viewPOpExp, viewDScale;
151:   PetscViewer       viewer, viewerPre, viewerRate, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp, viewerDScale;
152:   PetscViewerFormat format, formatPre, formatRate, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp, formatDScale;

154:   /* ----------------Default work-area management -------------------- */
155:   PetscInt nwork;
156:   Vec     *work;

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

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

164:   PetscBool transpose_solve; /* solve transpose system instead */
165:   struct {
166:     Mat       AT, BT;
167:     PetscBool use_explicittranspose; /* transpose the system explicitly in KSP[Mat]SolveTranspose() */
168:     PetscBool reuse_transpose;       /* reuse the previous transposed system */
169:   } transpose;

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

173:   PCSide      pc_side_set;  /* PC type set explicitly by user */
174:   KSPNormType normtype_set; /* Norm type set explicitly by user */

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

179:   PetscBool dscale;     /* diagonal scale system; used with KSPSetDiagonalScale() */
180:   PetscBool dscalefix;  /* unscale system after solve */
181:   PetscBool dscalefix2; /* system has been unscaled */
182:   Vec       diagonal;   /* 1/sqrt(diag of matrix) */
183:   Vec       truediagonal;

185:   /* Allow declaring convergence when negative curvature is detected */
186:   PetscBool converged_neg_curve;

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

191:   /* User-defined pre/post solve callbacks */
192:   PetscErrorCode (*presolve)(KSP, Vec, Vec, void *);
193:   PetscErrorCode (*postsolve)(KSP, Vec, Vec, void *);
194:   void *prectx, *postctx;

196:   /* PETSc internal pre/post solve callbacks for Eisenstat and Walker trick */
197:   PetscErrorCode (*presolve_ew)(KSP, Vec, Vec, void *);
198:   PetscErrorCode (*postsolve_ew)(KSP, Vec, Vec, void *);
199:   void *prectx_ew, *postctx_ew;

201:   PetscInt nestlevel; /* how many levels of nesting does the KSP have */
202: };

204: typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
205:   PetscReal coef;
206:   PetscReal bnrm;
207: } KSPDynTolCtx;

209: typedef struct {
210:   PetscBool initialrtol;    /* default relative residual decrease is computed from initial residual, not rhs */
211:   PetscBool mininitialrtol; /* default relative residual decrease is computed from min of initial residual and rhs */
212:   PetscBool convmaxits;     /* if true, the convergence test returns KSP_CONVERGED_ITS if the maximum number of iterations is reached */
213:   Vec       work;
214: } KSPConvergedDefaultCtx;

216: static inline PetscErrorCode KSPLogResidualHistory(KSP ksp, PetscReal norm)
217: {
218:   PetscFunctionBegin;
219:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
220:   if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) ksp->res_hist[ksp->res_hist_len++] = norm;
221:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
222:   PetscFunctionReturn(PETSC_SUCCESS);
223: }

225: static inline PetscErrorCode KSPLogErrorHistory(KSP ksp)
226: {
227:   DM dm;

229:   PetscFunctionBegin;
230:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
231:   PetscCall(KSPGetDM(ksp, &dm));
232:   if (dm && ksp->err_hist && ksp->err_hist_max > ksp->err_hist_len) {
233:     PetscSimplePointFn *exactSol;
234:     void               *exactCtx;
235:     PetscDS             ds;
236:     Vec                 u;
237:     PetscReal           error;
238:     PetscInt            Nf;

240:     PetscCall(KSPBuildSolution(ksp, NULL, &u));
241:     /* TODO Was needed to correct for Newton solution, but I just need to set a solution */
242:     //PetscCall(VecScale(u, -1.0));
243:     /* TODO Case when I have a solution */
244:     if (0) {
245:       PetscCall(DMGetDS(dm, &ds));
246:       PetscCall(PetscDSGetNumFields(ds, &Nf));
247:       PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle number of fields %" PetscInt_FMT " > 1 right now", Nf);
248:       PetscCall(PetscDSGetExactSolution(ds, 0, &exactSol, &exactCtx));
249:       PetscCall(DMComputeL2FieldDiff(dm, 0.0, &exactSol, &exactCtx, u, &error));
250:     } else {
251:       /* The null solution A 0 = 0 */
252:       PetscCall(VecNorm(u, NORM_2, &error));
253:     }
254:     ksp->err_hist[ksp->err_hist_len++] = error;
255:   }
256:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
257:   PetscFunctionReturn(PETSC_SUCCESS);
258: }

260: static inline PetscScalar KSPNoisyHash_Private(PetscInt xx)
261: {
262:   unsigned int x = (unsigned int)xx;
263:   x              = ((x >> 16) ^ x) * 0x45d9f3b;
264:   x              = ((x >> 16) ^ x) * 0x45d9f3b;
265:   x              = ((x >> 16) ^ x);
266:   return (PetscScalar)(((PetscInt64)x - 2147483648) * 5.e-10); /* center around zero, scaled about -1. to 1.*/
267: }

269: static inline PetscErrorCode KSPSetNoisy_Private(Mat A, Vec v)
270: {
271:   PetscScalar *a;
272:   PetscInt     n, istart;
273:   MatNullSpace nullsp = NULL;

275:   PetscFunctionBegin;
276:   if (A) PetscCall(MatGetNullSpace(A, &nullsp));
277:   PetscCall(VecGetOwnershipRange(v, &istart, NULL));
278:   PetscCall(VecGetLocalSize(v, &n));
279:   PetscCall(VecGetArrayWrite(v, &a));
280:   for (PetscInt i = 0; i < n; ++i) a[i] = KSPNoisyHash_Private(i + istart);
281:   PetscCall(VecRestoreArrayWrite(v, &a));
282:   if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, v));
283:   PetscFunctionReturn(PETSC_SUCCESS);
284: }

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

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

290: typedef struct _p_DMKSP  *DMKSP;
291: typedef struct _DMKSPOps *DMKSPOps;
292: struct _DMKSPOps {
293:   KSPCreateOperatorsFn     *createoperators;
294:   KSPComputeOperatorsFn    *computeoperators;
295:   KSPComputeRHSFn          *computerhs;
296:   KSPComputeInitialGuessFn *computeinitialguess;
297:   PetscErrorCode (*destroy)(DMKSP *);
298:   PetscErrorCode (*duplicate)(DMKSP, DMKSP);
299: };

301: /*S
302:    DMKSP - Object held by a `DM` that contains all the callback functions and their contexts needed by a `KSP`

304:    Level: developer

306:    Notes:
307:    Users provides callback functions and their contexts to `KSP` using, for example, `KSPSetComputeRHS()`. These values are stored
308:    in a `DMKSP` that is contained in the `DM` associated with the `KSP`. If no `DM` was provided by
309:    the user with `KSPSetDM()` it is automatically created by `KSPGetDM()` with `DMShellCreate()`.

311:    Users very rarely need to worked directly with the `DMKSP` object, rather they work with the `KSP` and the `DM` they created

313:    Multiple `DM` can share a single `DMKSP`, often each `DM` is associated with
314:    a grid refinement level. `DMGetDMKSP()` returns the `DMKSP` associated with a `DM`. `DMGetDMKSPWrite()` returns a unique
315:    `DMKSP` that is only associated with the current `DM`, making a copy of the shared `DMKSP` if needed (copy-on-write).

317:    Developer Notes:
318:    It is rather subtle why `DMKSP`, `DMSNES`, and `DMTS` are needed instead of simply storing the user callback functions and contexts in `DM` or `KSP`, `SNES`, or `TS`.
319:    It is to support composable solvers such as geometric multigrid. We want, by default, the same callback functions and contexts for all the levels in the computation,
320:    but we need to also support different callbacks and contexts on each level. The copy-on-write approach of `DMGetDMKSPWrite()` makes this possible.

322:    The `originaldm` inside the `DMKSP` is NOT reference counted (to prevent a reference count loop between a `DM` and a `DMKSP`).
323:    The `DM` on which this context was first created is cached here to implement one-way
324:    copy-on-write. When `DMGetDMKSPWrite()` sees a request using a different `DM`, it makes a copy of the `TSDM`. Thus, if a user
325:    only interacts directly with one level, e.g., using `TSSetIFunction()`, then coarse levels of a multilevel item
326:    integrator are built, then the user changes the routine with another call to `TSSetIFunction()`, it automatically
327:    propagates to all the levels. If instead, they get out a specific level and set the function on that level,
328:    subsequent changes to the original level will no longer propagate to that level.

330: .seealso: [](ch_ts), `KSP`, `KSPCreate()`, `DM`, `DMGetDMKSPWrite()`, `DMGetDMKSP()`,  `DMSNES`, `DMTS`, `DMKSPSetComputeOperators()`, `DMKSPGetComputeOperators()`,
331:           `DMKSPSetComputeRHS()`, `DMKSPSetComputeInitialGuess()`
332: S*/
333: struct _p_DMKSP {
334:   PETSCHEADER(struct _DMKSPOps);
335:   void *createoperatorsctx;
336:   void *operatorsctx;
337:   void *rhsctx;
338:   void *initialguessctx;
339:   void *data;

341:   /* See developer note for `DMKSP` above */
342:   DM originaldm;

344:   void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */
345: };
346: PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM, DMKSP *);
347: PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM, DMKSP *);
348: PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM, DM);

350: /*
351:        These allow the various Krylov methods to apply to either the linear system or its transpose.
352: */
353: static inline PetscErrorCode KSP_RemoveNullSpace(KSP ksp, Vec y)
354: {
355:   PetscFunctionBegin;
356:   if (ksp->pc_side == PC_LEFT) {
357:     Mat          A;
358:     MatNullSpace nullsp;

360:     PetscCall(PCGetOperators(ksp->pc, &A, NULL));
361:     PetscCall(MatGetNullSpace(A, &nullsp));
362:     if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y));
363:   }
364:   PetscFunctionReturn(PETSC_SUCCESS);
365: }

367: static inline PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp, Vec y)
368: {
369:   PetscFunctionBegin;
370:   if (ksp->pc_side == PC_LEFT) {
371:     Mat          A;
372:     MatNullSpace nullsp;

374:     PetscCall(PCGetOperators(ksp->pc, &A, NULL));
375:     PetscCall(MatGetTransposeNullSpace(A, &nullsp));
376:     if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y));
377:   }
378:   PetscFunctionReturn(PETSC_SUCCESS);
379: }

381: static inline PetscErrorCode KSP_MatMult(KSP ksp, Mat A, Vec x, Vec y)
382: {
383:   PetscFunctionBegin;
384:   if (ksp->transpose_solve) PetscCall(MatMultTranspose(A, x, y));
385:   else PetscCall(MatMult(A, x, y));
386:   PetscFunctionReturn(PETSC_SUCCESS);
387: }

389: static inline PetscErrorCode KSP_MatMultTranspose(KSP ksp, Mat A, Vec x, Vec y)
390: {
391:   PetscFunctionBegin;
392:   if (ksp->transpose_solve) PetscCall(MatMult(A, x, y));
393:   else PetscCall(MatMultTranspose(A, x, y));
394:   PetscFunctionReturn(PETSC_SUCCESS);
395: }

397: static inline PetscErrorCode KSP_MatMultHermitianTranspose(KSP ksp, Mat A, Vec x, Vec y)
398: {
399:   PetscFunctionBegin;
400:   if (!ksp->transpose_solve) PetscCall(MatMultHermitianTranspose(A, x, y));
401:   else {
402:     Vec w;

404:     PetscCall(VecDuplicate(x, &w));
405:     PetscCall(VecCopy(x, w));
406:     PetscCall(VecConjugate(w));
407:     PetscCall(MatMult(A, w, y));
408:     PetscCall(VecDestroy(&w));
409:     PetscCall(VecConjugate(y));
410:   }
411:   PetscFunctionReturn(PETSC_SUCCESS);
412: }

414: static inline PetscErrorCode KSP_PCApply(KSP ksp, Vec x, Vec y)
415: {
416:   PetscFunctionBegin;
417:   if (ksp->transpose_solve) {
418:     PetscCall(PCApplyTranspose(ksp->pc, x, y));
419:     PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y));
420:   } else {
421:     PetscCall(PCApply(ksp->pc, x, y));
422:     PetscCall(KSP_RemoveNullSpace(ksp, y));
423:   }
424:   PetscFunctionReturn(PETSC_SUCCESS);
425: }

427: static inline PetscErrorCode KSP_PCApplyTranspose(KSP ksp, Vec x, Vec y)
428: {
429:   PetscFunctionBegin;
430:   if (ksp->transpose_solve) {
431:     PetscCall(PCApply(ksp->pc, x, y));
432:     PetscCall(KSP_RemoveNullSpace(ksp, y));
433:   } else {
434:     PetscCall(PCApplyTranspose(ksp->pc, x, y));
435:     PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y));
436:   }
437:   PetscFunctionReturn(PETSC_SUCCESS);
438: }

440: static inline PetscErrorCode KSP_PCApplyHermitianTranspose(KSP ksp, Vec x, Vec y)
441: {
442:   PetscFunctionBegin;
443:   PetscCall(VecConjugate(x));
444:   PetscCall(KSP_PCApplyTranspose(ksp, x, y));
445:   PetscCall(VecConjugate(x));
446:   PetscCall(VecConjugate(y));
447:   PetscFunctionReturn(PETSC_SUCCESS);
448: }

450: static inline PetscErrorCode KSP_PCMatApply(KSP ksp, Mat X, Mat Y)
451: {
452:   PetscFunctionBegin;
453:   if (ksp->transpose_solve) PetscCall(PCMatApplyTranspose(ksp->pc, X, Y));
454:   else PetscCall(PCMatApply(ksp->pc, X, Y));
455:   PetscFunctionReturn(PETSC_SUCCESS);
456: }

458: static inline PetscErrorCode KSP_PCMatApplyTranspose(KSP ksp, Mat X, Mat Y)
459: {
460:   PetscFunctionBegin;
461:   if (!ksp->transpose_solve) PetscCall(PCMatApplyTranspose(ksp->pc, X, Y));
462:   else PetscCall(PCMatApply(ksp->pc, X, Y));
463:   PetscFunctionReturn(PETSC_SUCCESS);
464: }

466: static inline PetscErrorCode KSP_PCApplyBAorAB(KSP ksp, Vec x, Vec y, Vec w)
467: {
468:   PetscFunctionBegin;
469:   if (ksp->transpose_solve) {
470:     PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w));
471:     PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y));
472:   } else {
473:     PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w));
474:     PetscCall(KSP_RemoveNullSpace(ksp, y));
475:   }
476:   PetscFunctionReturn(PETSC_SUCCESS);
477: }

479: static inline PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp, Vec x, Vec y, Vec w)
480: {
481:   PetscFunctionBegin;
482:   if (ksp->transpose_solve) PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w));
483:   else PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w));
484:   PetscFunctionReturn(PETSC_SUCCESS);
485: }

487: PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization;
488: PETSC_EXTERN PetscLogEvent KSP_SetUp;
489: PETSC_EXTERN PetscLogEvent KSP_Solve;
490: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0;
491: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1;
492: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2;
493: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3;
494: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4;
495: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S;
496: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L;
497: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U;
498: PETSC_EXTERN PetscLogEvent KSP_SolveTranspose;
499: PETSC_EXTERN PetscLogEvent KSP_MatSolve;
500: PETSC_EXTERN PetscLogEvent KSP_MatSolveTranspose;

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

505: /*MC
506:    KSPCheckDot - Checks if the result of a dot product used by the corresponding `KSP` contains infinity or NaN. These indicate that the previous
507:       application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`.

509:    Collective

511:    Input Parameter:
512: .  ksp - the linear solver `KSP` context.

514:    Output Parameter:
515: .  beta - the result of the inner product

517:    Level: developer

519:    Developer Notes:
520:    Used to manage returning from `KSP` solvers collectively whose preconditioners have failed, possibly only a subset of MPI processes, in some way

522:    It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products.

524: .seealso: `PCFailedReason`, `KSPConvergedReason`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckNorm()`, `KSPCheckSolve()`,
525:           `KSPSetErrorIfNotConverged()`
526: M*/
527: #define KSPCheckDot(ksp, beta) \
528:   do { \
529:     if (PetscIsInfOrNanScalar(beta)) { \
530:       PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to infinity or NaN inner product"); \
531:       { \
532:         PCFailedReason pcreason; \
533:         PetscCall(PCReduceFailedReason(ksp->pc)); \
534:         PetscCall(PCGetFailedReason(ksp->pc, &pcreason)); \
535:         PetscCall(VecFlag(ksp->vec_sol, pcreason)); \
536:         if (pcreason) { \
537:           ksp->reason = KSP_DIVERGED_PC_FAILED; \
538:         } else { \
539:           ksp->reason = KSP_DIVERGED_NANORINF; \
540:         } \
541:         PetscFunctionReturn(PETSC_SUCCESS); \
542:       } \
543:     } \
544:   } while (0)

546: /*MC
547:    KSPCheckNorm - Checks if the result of a norm used by the corresponding `KSP` contains `inf` or `NaN`. These indicate that the previous
548:       application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`.

550:    Collective

552:    Input Parameter:
553: .  ksp - the linear solver `KSP` context.

555:    Output Parameter:
556: .  beta - the result of the norm

558:    Level: developer

560:    Developer Notes:
561:    Used to manage returning from `KSP` solvers collectively whose preconditioners have failed, possibly only a subset of MPI processes, in some way.

563:    It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products.

565: .seealso: `PCFailedReason`, `KSPConvergedReason`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckDot()`, `KSPCheckSolve()`,
566:           `KSPSetErrorIfNotConverged()`
567: M*/
568: #define KSPCheckNorm(ksp, beta) \
569:   do { \
570:     if (PetscIsInfOrNanReal(beta)) { \
571:       PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to infinity or NaN norm"); \
572:       { \
573:         PCFailedReason pcreason; \
574:         PetscCall(PCReduceFailedReason(ksp->pc)); \
575:         PetscCall(PCGetFailedReason(ksp->pc, &pcreason)); \
576:         PetscCall(VecFlag(ksp->vec_sol, pcreason)); \
577:         if (pcreason) { \
578:           ksp->reason = KSP_DIVERGED_PC_FAILED; \
579:         } else { \
580:           ksp->reason = KSP_DIVERGED_NANORINF; \
581:         } \
582:         ksp->rnorm = beta; \
583:         PetscFunctionReturn(PETSC_SUCCESS); \
584:       } \
585:     } \
586:   } while (0)

588: PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]);
589: PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP, PetscInt, PetscReal *);