xref: /petsc/include/petsc/private/kspimpl.h (revision 2d6c3ceebcb340e82a066ce435cc2ceddc132986)
1 
2 #ifndef _KSPIMPL_H
3   #define _KSPIMPL_H
4 
5   #include <petscksp.h>
6   #include <petscds.h>
7   #include <petsc/private/petscimpl.h>
8 
9 /* SUBMANSEC = KSP */
10 
11 PETSC_EXTERN PetscBool      KSPRegisterAllCalled;
12 PETSC_EXTERN PetscBool      KSPMonitorRegisterAllCalled;
13 PETSC_EXTERN PetscErrorCode KSPRegisterAll(void);
14 PETSC_EXTERN PetscErrorCode KSPMonitorRegisterAll(void);
15 PETSC_EXTERN PetscErrorCode KSPGuessRegisterAll(void);
16 PETSC_EXTERN PetscErrorCode KSPMatRegisterAll(void);
17 
18 typedef struct _KSPOps *KSPOps;
19 
20 struct _KSPOps {
21   PetscErrorCode (*buildsolution)(KSP, Vec, Vec *);      /* Returns a pointer to the solution, or
22                                                           calculates the solution in a
23                                                           user-provided area. */
24   PetscErrorCode (*buildresidual)(KSP, Vec, Vec, Vec *); /* Returns a pointer to the residual, or
25                                                           calculates the residual in a
26                                                           user-provided area.  */
27   PetscErrorCode (*solve)(KSP);                          /* actual solver */
28   PetscErrorCode (*matsolve)(KSP, Mat, Mat);             /* multiple dense RHS solver */
29   PetscErrorCode (*setup)(KSP);
30   PetscErrorCode (*setfromoptions)(KSP, PetscOptionItems *);
31   PetscErrorCode (*publishoptions)(KSP);
32   PetscErrorCode (*computeextremesingularvalues)(KSP, PetscReal *, PetscReal *);
33   PetscErrorCode (*computeeigenvalues)(KSP, PetscInt, PetscReal *, PetscReal *, PetscInt *);
34   PetscErrorCode (*computeritz)(KSP, PetscBool, PetscBool, PetscInt *, Vec[], PetscReal *, PetscReal *);
35   PetscErrorCode (*destroy)(KSP);
36   PetscErrorCode (*view)(KSP, PetscViewer);
37   PetscErrorCode (*reset)(KSP);
38   PetscErrorCode (*load)(KSP, PetscViewer);
39 };
40 
41 typedef struct _KSPGuessOps *KSPGuessOps;
42 
43 struct _KSPGuessOps {
44   PetscErrorCode (*formguess)(KSPGuess, Vec, Vec); /* Form initial guess */
45   PetscErrorCode (*update)(KSPGuess, Vec, Vec);    /* Update database */
46   PetscErrorCode (*setfromoptions)(KSPGuess);
47   PetscErrorCode (*settolerance)(KSPGuess, PetscReal);
48   PetscErrorCode (*setup)(KSPGuess);
49   PetscErrorCode (*destroy)(KSPGuess);
50   PetscErrorCode (*view)(KSPGuess, PetscViewer);
51   PetscErrorCode (*reset)(KSPGuess);
52 };
53 
54 /*
55    Defines the KSPGuess data structure.
56 */
57 struct _p_KSPGuess {
58   PETSCHEADER(struct _KSPGuessOps);
59   KSP              ksp;       /* the parent KSP */
60   Mat              A;         /* the current linear operator */
61   PetscObjectState omatstate; /* previous linear operator state */
62   void            *data;      /* pointer to the specific implementation */
63 };
64 
65 PETSC_EXTERN PetscErrorCode KSPGuessCreate_Fischer(KSPGuess);
66 PETSC_EXTERN PetscErrorCode KSPGuessCreate_POD(KSPGuess);
67 
68   /*
69      Maximum number of monitors you can run with a single KSP
70 */
71   #define MAXKSPMONITORS    5
72   #define MAXKSPREASONVIEWS 5
73 typedef enum {
74   KSP_SETUP_NEW,
75   KSP_SETUP_NEWMATRIX,
76   KSP_SETUP_NEWRHS
77 } KSPSetUpStage;
78 
79 /*
80    Defines the KSP data structure.
81 */
82 struct _p_KSP {
83   PETSCHEADER(struct _KSPOps);
84   DM        dm;
85   PetscBool dmAuto;   /* DM was created automatically by KSP */
86   PetscBool dmActive; /* KSP should use DM for computing operators */
87   /*------------------------- User parameters--------------------------*/
88   PetscInt  max_it; /* maximum number of iterations */
89   PetscInt  min_it; /* minimum number of iterations */
90   KSPGuess  guess;
91   PetscBool guess_zero,                                  /* flag for whether initial guess is 0 */
92     guess_not_read,                                      /* guess is not read, does not need to be zeroed */
93     calc_sings,                                          /* calculate extreme Singular Values */
94     calc_ritz,                                           /* calculate (harmonic) Ritz pairs */
95     guess_knoll;                                         /* use initial guess of PCApply(ksp->B,b */
96   PCSide    pc_side;                                     /* flag for left, right, or symmetric preconditioning */
97   PetscInt  normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */
98   PetscReal rtol,                                        /* relative tolerance */
99     abstol,                                              /* absolute tolerance */
100     ttol,                                                /* (not set by user)  */
101     divtol;                                              /* divergence tolerance */
102   PetscReal          rnorm0;                             /* initial residual norm (used for divergence testing) */
103   PetscReal          rnorm;                              /* current residual norm */
104   KSPConvergedReason reason;
105   PetscBool          errorifnotconverged; /* create an error if the KSPSolve() does not converge */
106 
107   Vec        vec_sol, vec_rhs; /* pointer to where user has stashed
108                                       the solution and rhs, these are
109                                       never touched by the code, only
110                                       passed back to the user */
111   PetscReal *res_hist;         /* If !0 stores residual each at iteration */
112   PetscReal *res_hist_alloc;   /* If !0 means user did not provide buffer, needs deallocation */
113   size_t     res_hist_len;     /* current size of residual history array */
114   size_t     res_hist_max;     /* actual amount of storage in residual history */
115   PetscBool  res_hist_reset;   /* reset history to length zero for each new solve */
116   PetscReal *err_hist;         /* If !0 stores error at each iteration */
117   PetscReal *err_hist_alloc;   /* If !0 means user did not provide buffer, needs deallocation */
118   size_t     err_hist_len;     /* current size of error history array */
119   size_t     err_hist_max;     /* actual amount of storage in error history */
120   PetscBool  err_hist_reset;   /* reset history to length zero for each new solve */
121 
122   PetscInt  chknorm; /* only compute/check norm if iterations is great than this */
123   PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the
124                                         MPI_Allreduce() for computing the inner products for the next iteration. */
125 
126   PetscInt nmax; /* maximum number of right-hand sides to be handled simultaneously */
127 
128   /* --------User (or default) routines (most return -1 on error) --------*/
129   PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP, PetscInt, PetscReal, void *); /* returns control to user after */
130   PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void **);                   /* */
131   void     *monitorcontext[MAXKSPMONITORS];                                    /* residual calculation, allows user */
132   PetscInt  numbermonitors;                                                    /* to, for instance, print residual norm, etc. */
133   PetscBool pauseFinal;                                                        /* Pause all drawing monitor at the final iterate */
134 
135   PetscErrorCode (*reasonview[MAXKSPREASONVIEWS])(KSP, void *);    /* KSP converged reason view */
136   PetscErrorCode (*reasonviewdestroy[MAXKSPREASONVIEWS])(void **); /* Optional destroy routine */
137   void    *reasonviewcontext[MAXKSPREASONVIEWS];                   /* User context */
138   PetscInt numberreasonviews;                                      /* Number if reason viewers */
139 
140   PetscErrorCode (*converged)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *);
141   PetscErrorCode (*convergeddestroy)(void *);
142   void *cnvP;
143 
144   void *user; /* optional user-defined context */
145 
146   PC pc;
147 
148   void *data; /* holder for misc stuff associated
149                                    with a particular iterative solver */
150 
151   PetscBool         view, viewPre, viewRate, viewMat, viewPMat, viewRhs, viewSol, viewMatExp, viewEV, viewSV, viewEVExp, viewFinalRes, viewPOpExp, viewDScale;
152   PetscViewer       viewer, viewerPre, viewerRate, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp, viewerDScale;
153   PetscViewerFormat format, formatPre, formatRate, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp, formatDScale;
154 
155   /* ----------------Default work-area management -------------------- */
156   PetscInt nwork;
157   Vec     *work;
158 
159   KSPSetUpStage setupstage;
160   PetscBool     setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */
161 
162   PetscInt its;      /* number of iterations so far computed in THIS linear solve*/
163   PetscInt totalits; /* number of iterations used by this KSP object since it was created */
164 
165   PetscBool transpose_solve; /* solve transpose system instead */
166   struct {
167     Mat       AT, BT;
168     PetscBool use_explicittranspose; /* transpose the system explicitly in KSPSolveTranspose */
169     PetscBool reuse_transpose;       /* reuse the previous transposed system */
170   } transpose;
171 
172   KSPNormType normtype; /* type of norm used for convergence tests */
173 
174   PCSide      pc_side_set;  /* PC type set explicitly by user */
175   KSPNormType normtype_set; /* Norm type set explicitly by user */
176 
177   /*   Allow diagonally scaling the matrix before computing the preconditioner or using
178        the Krylov method. Note this is NOT just Jacobi preconditioning */
179 
180   PetscBool dscale;     /* diagonal scale system; used with KSPSetDiagonalScale() */
181   PetscBool dscalefix;  /* unscale system after solve */
182   PetscBool dscalefix2; /* system has been unscaled */
183   Vec       diagonal;   /* 1/sqrt(diag of matrix) */
184   Vec       truediagonal;
185 
186   /* Allow declaring convergence when negative curvature is detected */
187   PetscBool converged_neg_curve;
188 
189   PetscInt  setfromoptionscalled;
190   PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */
191 
192   PetscErrorCode (*presolve)(KSP, Vec, Vec, void *);
193   PetscErrorCode (*postsolve)(KSP, Vec, Vec, void *);
194   void *prectx, *postctx;
195 };
196 
197 typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
198   PetscReal coef;
199   PetscReal bnrm;
200 } KSPDynTolCtx;
201 
202 typedef struct {
203   PetscBool initialrtol;    /* default relative residual decrease is computed from initial residual, not rhs */
204   PetscBool mininitialrtol; /* default relative residual decrease is computed from min of initial residual and rhs */
205   PetscBool convmaxits;     /* if true, the convergence test returns KSP_CONVERGED_ITS if the maximum number of iterations is reached */
206   Vec       work;
207 } KSPConvergedDefaultCtx;
208 
209 static inline PetscErrorCode KSPLogResidualHistory(KSP ksp, PetscReal norm)
210 {
211   PetscFunctionBegin;
212   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
213   if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) ksp->res_hist[ksp->res_hist_len++] = norm;
214   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
215   PetscFunctionReturn(PETSC_SUCCESS);
216 }
217 
218 static inline PetscErrorCode KSPLogErrorHistory(KSP ksp)
219 {
220   DM dm;
221 
222   PetscFunctionBegin;
223   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
224   PetscCall(KSPGetDM(ksp, &dm));
225   if (dm && ksp->err_hist && ksp->err_hist_max > ksp->err_hist_len) {
226     PetscSimplePointFunc exactSol;
227     void                *exactCtx;
228     PetscDS              ds;
229     Vec                  u;
230     PetscReal            error;
231     PetscInt             Nf;
232 
233     PetscCall(KSPBuildSolution(ksp, NULL, &u));
234     /* TODO Was needed to correct for Newton solution, but I just need to set a solution */
235     //PetscCall(VecScale(u, -1.0));
236     /* TODO Case when I have a solution */
237     if (0) {
238       PetscCall(DMGetDS(dm, &ds));
239       PetscCall(PetscDSGetNumFields(ds, &Nf));
240       PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle number of fields %" PetscInt_FMT " > 1 right now", Nf);
241       PetscCall(PetscDSGetExactSolution(ds, 0, &exactSol, &exactCtx));
242       PetscCall(DMComputeL2FieldDiff(dm, 0.0, &exactSol, &exactCtx, u, &error));
243     } else {
244       /* The null solution A 0 = 0 */
245       PetscCall(VecNorm(u, NORM_2, &error));
246     }
247     ksp->err_hist[ksp->err_hist_len++] = error;
248   }
249   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
250   PetscFunctionReturn(PETSC_SUCCESS);
251 }
252 
253 static inline PetscScalar KSPNoisyHash_Private(PetscInt xx)
254 {
255   unsigned int x = (unsigned int)xx;
256   x              = ((x >> 16) ^ x) * 0x45d9f3b;
257   x              = ((x >> 16) ^ x) * 0x45d9f3b;
258   x              = ((x >> 16) ^ x);
259   return (PetscScalar)((PetscInt64)x - 2147483648) * 5.e-10; /* center around zero, scaled about -1. to 1.*/
260 }
261 
262 static inline PetscErrorCode KSPSetNoisy_Private(Vec v)
263 {
264   PetscScalar *a;
265   PetscInt     n, istart;
266 
267   PetscFunctionBegin;
268   PetscCall(VecGetOwnershipRange(v, &istart, NULL));
269   PetscCall(VecGetLocalSize(v, &n));
270   PetscCall(VecGetArrayWrite(v, &a));
271   for (PetscInt i = 0; i < n; ++i) a[i] = KSPNoisyHash_Private(i + istart);
272   PetscCall(VecRestoreArrayWrite(v, &a));
273   PetscFunctionReturn(PETSC_SUCCESS);
274 }
275 
276 PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP, PetscBool, KSPNormType *, PCSide *);
277 
278 PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP, PetscInt, const PetscReal *, const PetscReal *);
279 
280 typedef struct _p_DMKSP  *DMKSP;
281 typedef struct _DMKSPOps *DMKSPOps;
282 struct _DMKSPOps {
283   PetscErrorCode (*computeoperators)(KSP, Mat, Mat, void *);
284   PetscErrorCode (*computerhs)(KSP, Vec, void *);
285   PetscErrorCode (*computeinitialguess)(KSP, Vec, void *);
286   PetscErrorCode (*destroy)(DMKSP *);
287   PetscErrorCode (*duplicate)(DMKSP, DMKSP);
288 };
289 
290 struct _p_DMKSP {
291   PETSCHEADER(struct _DMKSPOps);
292   void *operatorsctx;
293   void *rhsctx;
294   void *initialguessctx;
295   void *data;
296 
297   /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
298    * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user
299    * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by
300    * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates
301    * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes
302    * to the original level will no longer propagate to that level.
303    */
304   DM originaldm;
305 
306   void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */
307 };
308 PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM, DMKSP *);
309 PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM, DMKSP *);
310 PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM, DM);
311 
312 /*
313        These allow the various Krylov methods to apply to either the linear system or its transpose.
314 */
315 static inline PetscErrorCode KSP_RemoveNullSpace(KSP ksp, Vec y)
316 {
317   PetscFunctionBegin;
318   if (ksp->pc_side == PC_LEFT) {
319     Mat          A;
320     MatNullSpace nullsp;
321 
322     PetscCall(PCGetOperators(ksp->pc, &A, NULL));
323     PetscCall(MatGetNullSpace(A, &nullsp));
324     if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y));
325   }
326   PetscFunctionReturn(PETSC_SUCCESS);
327 }
328 
329 static inline PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp, Vec y)
330 {
331   PetscFunctionBegin;
332   if (ksp->pc_side == PC_LEFT) {
333     Mat          A;
334     MatNullSpace nullsp;
335 
336     PetscCall(PCGetOperators(ksp->pc, &A, NULL));
337     PetscCall(MatGetTransposeNullSpace(A, &nullsp));
338     if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y));
339   }
340   PetscFunctionReturn(PETSC_SUCCESS);
341 }
342 
343 static inline PetscErrorCode KSP_MatMult(KSP ksp, Mat A, Vec x, Vec y)
344 {
345   PetscFunctionBegin;
346   if (ksp->transpose_solve) PetscCall(MatMultTranspose(A, x, y));
347   else PetscCall(MatMult(A, x, y));
348   PetscFunctionReturn(PETSC_SUCCESS);
349 }
350 
351 static inline PetscErrorCode KSP_MatMultTranspose(KSP ksp, Mat A, Vec x, Vec y)
352 {
353   PetscFunctionBegin;
354   if (ksp->transpose_solve) PetscCall(MatMult(A, x, y));
355   else PetscCall(MatMultTranspose(A, x, y));
356   PetscFunctionReturn(PETSC_SUCCESS);
357 }
358 
359 static inline PetscErrorCode KSP_MatMultHermitianTranspose(KSP ksp, Mat A, Vec x, Vec y)
360 {
361   PetscFunctionBegin;
362   if (!ksp->transpose_solve) PetscCall(MatMultHermitianTranspose(A, x, y));
363   else {
364     Vec w;
365 
366     PetscCall(VecDuplicate(x, &w));
367     PetscCall(VecCopy(x, w));
368     PetscCall(VecConjugate(w));
369     PetscCall(MatMult(A, w, y));
370     PetscCall(VecDestroy(&w));
371     PetscCall(VecConjugate(y));
372   }
373   PetscFunctionReturn(PETSC_SUCCESS);
374 }
375 
376 static inline PetscErrorCode KSP_PCApply(KSP ksp, Vec x, Vec y)
377 {
378   PetscFunctionBegin;
379   if (ksp->transpose_solve) {
380     PetscCall(PCApplyTranspose(ksp->pc, x, y));
381     PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y));
382   } else {
383     PetscCall(PCApply(ksp->pc, x, y));
384     PetscCall(KSP_RemoveNullSpace(ksp, y));
385   }
386   PetscFunctionReturn(PETSC_SUCCESS);
387 }
388 
389 static inline PetscErrorCode KSP_PCApplyTranspose(KSP ksp, Vec x, Vec y)
390 {
391   PetscFunctionBegin;
392   if (ksp->transpose_solve) {
393     PetscCall(PCApply(ksp->pc, x, y));
394     PetscCall(KSP_RemoveNullSpace(ksp, y));
395   } else {
396     PetscCall(PCApplyTranspose(ksp->pc, x, y));
397     PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y));
398   }
399   PetscFunctionReturn(PETSC_SUCCESS);
400 }
401 
402 static inline PetscErrorCode KSP_PCApplyHermitianTranspose(KSP ksp, Vec x, Vec y)
403 {
404   PetscFunctionBegin;
405   PetscCall(VecConjugate(x));
406   PetscCall(KSP_PCApplyTranspose(ksp, x, y));
407   PetscCall(VecConjugate(x));
408   PetscCall(VecConjugate(y));
409   PetscFunctionReturn(PETSC_SUCCESS);
410 }
411 
412 static inline PetscErrorCode KSP_PCMatApply(KSP ksp, Mat X, Mat Y)
413 {
414   PetscFunctionBegin;
415   if (ksp->transpose_solve) {
416     PetscBool flg;
417     PetscCall(PetscObjectTypeCompareAny((PetscObject)ksp->pc, &flg, PCNONE, PCICC, PCCHOLESKY, ""));
418     PetscCheck(flg, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "PCMatApplyTranspose() not yet implemented for nonsymmetric PC");
419   }
420   PetscCall(PCMatApply(ksp->pc, X, Y));
421   PetscFunctionReturn(PETSC_SUCCESS);
422 }
423 
424 static inline PetscErrorCode KSP_PCMatApplyTranspose(KSP ksp, Mat X, Mat Y)
425 {
426   PetscFunctionBegin;
427   if (!ksp->transpose_solve) {
428     PetscBool flg;
429     PetscCall(PetscObjectTypeCompareAny((PetscObject)ksp->pc, &flg, PCNONE, PCICC, PCCHOLESKY, ""));
430     PetscCheck(flg, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "PCMatApplyTranspose() not yet implemented for nonsymmetric PC");
431   }
432   PetscCall(PCMatApply(ksp->pc, X, Y));
433   PetscFunctionReturn(PETSC_SUCCESS);
434 }
435 
436 static inline PetscErrorCode KSP_PCApplyBAorAB(KSP ksp, Vec x, Vec y, Vec w)
437 {
438   PetscFunctionBegin;
439   if (ksp->transpose_solve) {
440     PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w));
441     PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y));
442   } else {
443     PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w));
444     PetscCall(KSP_RemoveNullSpace(ksp, y));
445   }
446   PetscFunctionReturn(PETSC_SUCCESS);
447 }
448 
449 static inline PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp, Vec x, Vec y, Vec w)
450 {
451   PetscFunctionBegin;
452   if (ksp->transpose_solve) PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w));
453   else PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w));
454   PetscFunctionReturn(PETSC_SUCCESS);
455 }
456 
457 PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization;
458 PETSC_EXTERN PetscLogEvent KSP_SetUp;
459 PETSC_EXTERN PetscLogEvent KSP_Solve;
460 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0;
461 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1;
462 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2;
463 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3;
464 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4;
465 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S;
466 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L;
467 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U;
468 PETSC_EXTERN PetscLogEvent KSP_SolveTranspose;
469 PETSC_EXTERN PetscLogEvent KSP_MatSolve;
470 PETSC_EXTERN PetscLogEvent KSP_MatSolveTranspose;
471 
472 PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat, IS, IS, IS, IS, MatReuse, Mat *, MatSchurComplementAinvType, MatReuse, Mat *);
473 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC, PetscBool *);
474 
475   /*MC
476    KSPCheckDot - Checks if the result of a dot product used by the corresponding `KSP` contains Inf or NaN. These indicate that the previous
477       application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`.
478 
479    Collective
480 
481    Input Parameter:
482 .  ksp - the linear solver `KSP` context.
483 
484    Output Parameter:
485 .  beta - the result of the inner product
486 
487    Level: developer
488 
489    Developer Notes:
490    Used to manage returning from `KSP` solvers whose preconditioners have failed, possibly only a subset of MPI ranks, in some way
491 
492    It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products.
493 
494 .seealso: `PCFailedReason`, `KSPConvergedReason`, `PCGetFailedReasonRank()`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckNorm()`, `KSPCheckSolve()`
495 M*/
496   #define KSPCheckDot(ksp, beta) \
497     do { \
498       if (PetscIsInfOrNanScalar(beta)) { \
499         PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to Nan or Inf inner product"); \
500         { \
501           PCFailedReason pcreason; \
502           PetscInt       sendbuf, recvbuf; \
503           PetscCall(PCGetFailedReasonRank(ksp->pc, &pcreason)); \
504           sendbuf = (PetscInt)pcreason; \
505           PetscCall(MPIU_Allreduce(&sendbuf, &recvbuf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ksp))); \
506           if (recvbuf) { \
507             PetscCall(PCSetFailedReason(ksp->pc, (PCFailedReason)recvbuf)); \
508             ksp->reason = KSP_DIVERGED_PC_FAILED; \
509             PetscCall(VecSetInf(ksp->vec_sol)); \
510           } else { \
511             ksp->reason = KSP_DIVERGED_NANORINF; \
512           } \
513           PetscFunctionReturn(PETSC_SUCCESS); \
514         } \
515       } \
516     } while (0)
517 
518   /*MC
519    KSPCheckNorm - Checks if the result of a norm used by the corresponding `KSP` contains `inf` or `NaN`. These indicate that the previous
520       application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`.
521 
522    Collective
523 
524    Input Parameter:
525 .  ksp - the linear solver `KSP` context.
526 
527    Output Parameter:
528 .  beta - the result of the norm
529 
530    Level: developer
531 
532    Developer Notes:
533    Used to manage returning from `KSP` solvers whose preconditioners have failed, possibly only a subset of MPI ranks, in some way.
534 
535    It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products.
536 
537 .seealso: `PCFailedReason`, `KSPConvergedReason`, `PCGetFailedReasonRank()`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckDot()`, `KSPCheckSolve()`
538 M*/
539   #define KSPCheckNorm(ksp, beta) \
540     do { \
541       if (PetscIsInfOrNanReal(beta)) { \
542         PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to Nan or Inf norm"); \
543         { \
544           PCFailedReason pcreason; \
545           PetscInt       sendbuf, recvbuf; \
546           PetscCall(PCGetFailedReasonRank(ksp->pc, &pcreason)); \
547           sendbuf = (PetscInt)pcreason; \
548           PetscCall(MPIU_Allreduce(&sendbuf, &recvbuf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ksp))); \
549           if (recvbuf) { \
550             PetscCall(PCSetFailedReason(ksp->pc, (PCFailedReason)recvbuf)); \
551             ksp->reason = KSP_DIVERGED_PC_FAILED; \
552             PetscCall(VecSetInf(ksp->vec_sol)); \
553             ksp->rnorm = beta; \
554           } else { \
555             PetscCall(PCSetFailedReason(ksp->pc, PC_NOERROR)); \
556             ksp->reason = KSP_DIVERGED_NANORINF; \
557             ksp->rnorm  = beta; \
558           } \
559           PetscFunctionReturn(PETSC_SUCCESS); \
560         } \
561       } \
562     } while (0)
563 
564 #endif
565 
566 PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]);
567 PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP, PetscInt, PetscReal *);
568