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: MatState 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: Mat mat_rhs; /* borrowed reference to the batch of
110: right-hand sides being solved, set
111: by KSPMatSolve_Private() only while
112: ksp->ops->matsolve runs, never
113: destroyed by the KSP */
114: PetscReal *res_hist; /* If !0 stores residual each at iteration */
115: PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
116: PetscCount res_hist_len; /* current entry count of residual history array */
117: PetscCount res_hist_max; /* total entry count of storage in residual history */
118: PetscBool res_hist_reset; /* reset history to length zero for each new solve */
119: PetscReal *err_hist; /* If !0 stores error at each iteration */
120: PetscReal *err_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
121: PetscCount err_hist_len; /* current entry count of error history array */
122: PetscCount err_hist_max; /* total entry count of storage in error history */
123: PetscBool err_hist_reset; /* reset history to length zero for each new solve */
125: PetscInt chknorm; /* only compute/check norm if iterations is great than this */
126: PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the
127: MPI_Allreduce() for computing the inner products for the next iteration. */
129: PetscInt nmax; /* maximum number of right-hand sides to be handled simultaneously */
131: /* --------User (or default) routines (most return -1 on error) --------*/
132: KSPMonitorFn *monitor[MAXKSPMONITORS];
133: PetscCtxDestroyFn *monitordestroy[MAXKSPMONITORS];
134: void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */
135: PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */
136: PetscBool pauseFinal; /* Pause all drawing monitor at the final iterate */
138: PetscViewer convergedreasonviewer;
139: PetscViewerFormat convergedreasonformat;
140: KSPConvergedReasonViewFn *reasonview[MAXKSPREASONVIEWS]; /* KSP converged reason view */
141: PetscCtxDestroyFn *reasonviewdestroy[MAXKSPREASONVIEWS]; /* optional destroy routine */
142: void *reasonviewcontext[MAXKSPREASONVIEWS]; /* viewer context */
143: PetscInt numberreasonviews; /* current number of reason viewers */
145: KSPConvergenceTestFn *converged;
146: PetscCtxDestroyFn *convergeddestroy;
147: void *cnvP;
149: PetscCtx ctx; /* optional user-defined context */
151: PC pc;
153: void *data; /* holder for misc stuff associated with a particular iterative solver */
155: PetscBool view, viewPre, viewRate, viewMat, viewPMat, viewRhs, viewSol, viewMatExp, viewEV, viewSV, viewEVExp, viewFinalRes, viewPOpExp;
156: PetscViewer viewer, viewerPre, viewerRate, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp;
157: PetscViewerFormat format, formatPre, formatRate, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp;
159: /* ----------------Default work-area management -------------------- */
160: PetscInt nwork;
161: Vec *work;
163: /* state machine to decide whether type specific setup is needed */
164: KSPSetUpStage setupstage;
165: PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */
166: MatState amatstate;
168: PetscInt its; /* number of iterations so far computed in THIS linear solve*/
169: PetscInt totalits; /* number of iterations used by this KSP object since it was created */
171: PetscBool transpose_solve; /* PETSC_TRUE applies the set operators transposed; PETSC_FALSE applies them as-is, including for an explicit transpose solve */
172: struct {
173: Mat A, B, AT, BT; /* parent operators, referenced by the KSP, and their explicit transposes */
174: PetscObjectId Aid, Bid; /* IDs detect MatHeaderReplace() on A and B */
175: PetscObjectState Anonzerostate, Bnonzerostate; /* nonzero states of A and B when the transposes were formed */
176: PetscBool solve_requested; /* PETSC_TRUE means KSP[Mat]SolveTranspose() was requested; PETSC_FALSE means KSP[Mat]Solve() was requested */
177: PetscBool use_explicittranspose; /* transpose the system explicitly in KSP[Mat]SolveTranspose() */
178: PetscBool reuse_transpose; /* reuse the previous transposed system */
179: } transpose;
181: KSPNormType normtype; /* type of norm used for convergence tests */
183: PCSide pc_side_set; /* PC type set explicitly by user */
184: KSPNormType normtype_set; /* Norm type set explicitly by user */
186: /* Allow declaring convergence when negative curvature is detected */
187: PetscBool converged_neg_curve;
189: PetscInt setfromoptionscalled;
190: PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */
192: /* User-defined pre/post solve callbacks */
193: PetscErrorCode (*presolve)(KSP, Vec, Vec, void *);
194: PetscErrorCode (*postsolve)(KSP, Vec, Vec, void *);
195: void *prectx, *postctx;
197: /* PETSc internal pre/post solve callbacks for Eisenstat and Walker trick */
198: PetscErrorCode (*presolve_ew)(KSP, Vec, Vec, void *);
199: PetscErrorCode (*postsolve_ew)(KSP, Vec, Vec, void *);
200: void *prectx_ew, *postctx_ew;
202: PetscInt nestlevel; /* how many levels of nesting does the KSP have */
204: /* orthogonalization */
205: KSPOrthogonalizationFn *orthog; /* orthogonalization function */
206: KSPOrthogonalizationCGSRefinementType cgstype; /* refinement in case of CGS */
207: PetscScalar *orthogwork; /* holds dot products computed in orthogonalization */
208: PetscInt lorthogwork; /* length of orthogwork */
209: };
211: typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
212: PetscReal coef;
213: PetscReal bnrm;
214: } KSPDynTolCtx;
216: typedef struct {
217: PetscBool initialrtol; /* default relative residual decrease is computed from initial residual, not rhs */
218: PetscBool mininitialrtol; /* default relative residual decrease is computed from min of initial residual and rhs */
219: PetscBool convmaxits; /* if true, the convergence test returns KSP_CONVERGED_ITS if the maximum number of iterations is reached */
220: Vec work;
221: } KSPConvergedDefaultCtx;
223: static inline PetscErrorCode KSPLogResidualHistory(KSP ksp, PetscReal norm)
224: {
225: PetscFunctionBegin;
226: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
227: if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) ksp->res_hist[ksp->res_hist_len++] = norm;
228: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: static inline PetscErrorCode KSPLogErrorHistory(KSP ksp)
233: {
234: DM dm;
236: PetscFunctionBegin;
237: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
238: PetscCall(KSPGetDM(ksp, &dm));
239: if (dm && ksp->err_hist && ksp->err_hist_max > ksp->err_hist_len) {
240: PetscSimplePointFn *exactSol;
241: void *exactCtx;
242: PetscDS ds;
243: Vec u;
244: PetscReal error;
245: PetscInt Nf;
247: PetscCall(KSPBuildSolution(ksp, NULL, &u));
248: /* TODO Was needed to correct for Newton solution, but I just need to set a solution */
249: //PetscCall(VecScale(u, -1.0));
250: /* TODO Case when I have a solution */
251: if (0) {
252: PetscCall(DMGetDS(dm, &ds));
253: PetscCall(PetscDSGetNumFields(ds, &Nf));
254: PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle number of fields %" PetscInt_FMT " > 1 right now", Nf);
255: PetscCall(PetscDSGetExactSolution(ds, 0, &exactSol, &exactCtx));
256: PetscCall(DMComputeL2FieldDiff(dm, 0.0, &exactSol, &exactCtx, u, &error));
257: } else {
258: /* The null solution A 0 = 0 */
259: PetscCall(VecNorm(u, NORM_2, &error));
260: }
261: ksp->err_hist[ksp->err_hist_len++] = error;
262: }
263: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: static inline PetscScalar KSPNoisyHash_Private(PetscInt xx)
268: {
269: unsigned int x = (unsigned int)xx;
270: x = ((x >> 16) ^ x) * 0x45d9f3b;
271: x = ((x >> 16) ^ x) * 0x45d9f3b;
272: x = ((x >> 16) ^ x);
273: return (PetscScalar)(((PetscInt64)x - 2147483648) * 5.e-10); /* center around zero, scaled about -1. to 1.*/
274: }
276: static inline PetscErrorCode KSPSetNoisy_Private(Mat A, Vec v)
277: {
278: PetscScalar *a;
279: PetscInt n, istart;
280: MatNullSpace nullsp = NULL;
282: PetscFunctionBegin;
283: if (A) PetscCall(MatGetNullSpace(A, &nullsp));
284: PetscCall(VecGetOwnershipRange(v, &istart, NULL));
285: PetscCall(VecGetLocalSize(v, &n));
286: PetscCall(VecGetArrayWrite(v, &a));
287: for (PetscInt i = 0; i < n; ++i) a[i] = KSPNoisyHash_Private(i + istart);
288: PetscCall(VecRestoreArrayWrite(v, &a));
289: if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, v));
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP, PetscBool, KSPNormType *, PCSide *);
295: PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP, PetscInt, const PetscReal *, const PetscReal *);
297: typedef struct _p_DMKSP *DMKSP;
298: typedef struct _DMKSPOps *DMKSPOps;
299: struct _DMKSPOps {
300: KSPCreateOperatorsFn *createoperators;
301: KSPComputeOperatorsFn *computeoperators;
302: KSPComputeRHSFn *computerhs;
303: KSPComputeInitialGuessFn *computeinitialguess;
304: PetscErrorCode (*destroy)(DMKSP *);
305: PetscErrorCode (*duplicate)(DMKSP, DMKSP);
306: };
308: /*S
309: DMKSP - Object held by a `DM` that contains all the callback functions and their contexts needed by a `KSP`
311: Level: developer
313: Notes:
314: Users provides callback functions and their contexts to `KSP` using, for example, `KSPSetComputeRHS()`. These values are stored
315: in a `DMKSP` that is contained in the `DM` associated with the `KSP`. If no `DM` was provided by
316: the user with `KSPSetDM()` it is automatically created by `KSPGetDM()` with `DMShellCreate()`.
318: Users very rarely need to worked directly with the `DMKSP` object, rather they work with the `KSP` and the `DM` they created
320: Multiple `DM` can share a single `DMKSP`, often each `DM` is associated with
321: a grid refinement level. `DMGetDMKSP()` returns the `DMKSP` associated with a `DM`. `DMGetDMKSPWrite()` returns a unique
322: `DMKSP` that is only associated with the current `DM`, making a copy of the shared `DMKSP` if needed (copy-on-write).
324: Developer Notes:
325: 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`.
326: 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,
327: but we need to also support different callbacks and contexts on each level. The copy-on-write approach of `DMGetDMKSPWrite()` makes this possible.
329: The `originaldm` inside the `DMKSP` is NOT reference counted (to prevent a reference count loop between a `DM` and a `DMKSP`).
330: The `DM` on which this context was first created is cached here to implement one-way
331: copy-on-write. When `DMGetDMKSPWrite()` sees a request using a different `DM`, it makes a copy of the `TSDM`. Thus, if a user
332: only interacts directly with one level, e.g., using `TSSetIFunction()`, then coarse levels of a multilevel item
333: integrator are built, then the user changes the routine with another call to `TSSetIFunction()`, it automatically
334: propagates to all the levels. If instead, they get out a specific level and set the function on that level,
335: subsequent changes to the original level will no longer propagate to that level.
337: .seealso: [](ch_ts), `KSP`, `KSPCreate()`, `DM`, `DMGetDMKSPWrite()`, `DMGetDMKSP()`, `DMSNES`, `DMTS`, `DMKSPSetComputeOperators()`, `DMKSPGetComputeOperators()`,
338: `DMKSPSetComputeRHS()`, `DMKSPSetComputeInitialGuess()`
339: S*/
340: struct _p_DMKSP {
341: PETSCHEADER(struct _DMKSPOps);
342: void *createoperatorsctx;
343: void *operatorsctx;
344: void *rhsctx;
345: void *initialguessctx;
346: void *data;
348: /* See developer note for `DMKSP` above */
349: DM originaldm;
351: void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */
352: };
353: PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM, DMKSP *);
354: PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM, DMKSP *);
355: PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM, DM);
357: /*
358: These allow the various Krylov methods to apply to either the linear system or its transpose.
359: */
360: static inline PetscErrorCode KSP_RemoveNullSpace(KSP ksp, Vec y)
361: {
362: PetscFunctionBegin;
363: if (ksp->pc_side == PC_LEFT) {
364: Mat A;
365: MatNullSpace nullsp;
367: PetscCall(PCGetOperators(ksp->pc, &A, NULL));
368: PetscCall(MatGetNullSpace(A, &nullsp));
369: if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y));
370: }
371: PetscFunctionReturn(PETSC_SUCCESS);
372: }
374: static inline PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp, Vec y)
375: {
376: PetscFunctionBegin;
377: if (ksp->pc_side == PC_LEFT) {
378: Mat A;
379: MatNullSpace nullsp;
381: PetscCall(PCGetOperators(ksp->pc, &A, NULL));
382: PetscCall(MatGetTransposeNullSpace(A, &nullsp));
383: if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y));
384: }
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: /*
389: Block analog of KSP_RemoveNullSpace() and KSP_RemoveNullSpaceTranspose(), the null space is removed from each column of Y. The removal is deliberately
390: kept out of KSP_PCMatApply() so that the behavior of the ksp->ops->matsolve implementations which do not need it, KSPPREONLY and KSPHPDDM, is left unchanged
391: */
392: static inline PetscErrorCode KSP_RemoveNullSpaceMat(KSP ksp, Mat Y)
393: {
394: PetscFunctionBegin;
395: if (ksp->pc_side == PC_LEFT) {
396: Mat A;
397: Vec cy;
398: PetscInt N;
399: MatNullSpace nullsp;
401: PetscCall(PCGetOperators(ksp->pc, &A, NULL));
402: if (ksp->transpose_solve) PetscCall(MatGetTransposeNullSpace(A, &nullsp));
403: else PetscCall(MatGetNullSpace(A, &nullsp));
404: if (nullsp) {
405: /* a matrix-based null space removal with dense kernels could replace this column loop for BLAS-3 performance */
406: PetscCall(MatGetSize(Y, NULL, &N));
407: for (PetscInt i = 0; i < N; i++) {
408: PetscCall(MatDenseGetColumnVec(Y, i, &cy));
409: PetscCall(MatNullSpaceRemove(nullsp, cy));
410: PetscCall(MatDenseRestoreColumnVec(Y, i, &cy));
411: }
412: }
413: }
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: static inline PetscErrorCode KSP_MatMult(KSP ksp, Mat A, Vec x, Vec y)
418: {
419: PetscFunctionBegin;
420: if (ksp->transpose_solve) PetscCall(MatMultTranspose(A, x, y));
421: else PetscCall(MatMult(A, x, y));
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: static inline PetscErrorCode KSP_MatMultTranspose(KSP ksp, Mat A, Vec x, Vec y)
426: {
427: PetscFunctionBegin;
428: if (ksp->transpose_solve) PetscCall(MatMult(A, x, y));
429: else PetscCall(MatMultTranspose(A, x, y));
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: static inline PetscErrorCode KSP_MatMultHermitianTranspose(KSP ksp, Mat A, Vec x, Vec y)
434: {
435: PetscFunctionBegin;
436: if (!ksp->transpose_solve) PetscCall(MatMultHermitianTranspose(A, x, y));
437: else {
438: Vec w;
440: PetscCall(VecDuplicate(x, &w));
441: PetscCall(VecCopy(x, w));
442: PetscCall(VecConjugate(w));
443: PetscCall(MatMult(A, w, y));
444: PetscCall(VecDestroy(&w));
445: PetscCall(VecConjugate(y));
446: }
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: static inline PetscErrorCode KSP_PCApply(KSP ksp, Vec x, Vec y)
451: {
452: PetscFunctionBegin;
453: if (ksp->transpose_solve) {
454: PetscCall(PCApplyTranspose(ksp->pc, x, y));
455: PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y));
456: } else {
457: PetscCall(PCApply(ksp->pc, x, y));
458: PetscCall(KSP_RemoveNullSpace(ksp, y));
459: }
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: static inline PetscErrorCode KSP_PCApplyTranspose(KSP ksp, Vec x, Vec y)
464: {
465: PetscFunctionBegin;
466: if (ksp->transpose_solve) {
467: PetscCall(PCApply(ksp->pc, x, y));
468: PetscCall(KSP_RemoveNullSpace(ksp, y));
469: } else {
470: PetscCall(PCApplyTranspose(ksp->pc, x, y));
471: PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y));
472: }
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: static inline PetscErrorCode KSP_PCApplyHermitianTranspose(KSP ksp, Vec x, Vec y)
477: {
478: PetscFunctionBegin;
479: PetscCall(VecConjugate(x));
480: PetscCall(KSP_PCApplyTranspose(ksp, x, y));
481: PetscCall(VecConjugate(x));
482: PetscCall(VecConjugate(y));
483: PetscFunctionReturn(PETSC_SUCCESS);
484: }
486: static inline PetscErrorCode KSP_PCMatApply(KSP ksp, Mat X, Mat Y)
487: {
488: PetscFunctionBegin;
489: if (ksp->transpose_solve) PetscCall(PCMatApplyTranspose(ksp->pc, X, Y));
490: else PetscCall(PCMatApply(ksp->pc, X, Y));
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: static inline PetscErrorCode KSP_PCMatApplyTranspose(KSP ksp, Mat X, Mat Y)
495: {
496: PetscFunctionBegin;
497: if (!ksp->transpose_solve) PetscCall(PCMatApplyTranspose(ksp->pc, X, Y));
498: else PetscCall(PCMatApply(ksp->pc, X, Y));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: static inline PetscErrorCode KSP_PCApplyBAorAB(KSP ksp, Vec x, Vec y, Vec w)
503: {
504: PetscFunctionBegin;
505: if (ksp->transpose_solve) {
506: PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w));
507: PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y));
508: } else {
509: PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w));
510: PetscCall(KSP_RemoveNullSpace(ksp, y));
511: }
512: PetscFunctionReturn(PETSC_SUCCESS);
513: }
515: static inline PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp, Vec x, Vec y, Vec w)
516: {
517: PetscFunctionBegin;
518: if (ksp->transpose_solve) PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w));
519: else PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: PETSC_EXTERN PetscLogEvent KSP_Orthogonalization;
524: PETSC_EXTERN PetscLogEvent KSP_SetUp;
525: PETSC_EXTERN PetscLogEvent KSP_Solve;
526: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0;
527: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1;
528: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2;
529: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3;
530: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4;
531: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S;
532: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L;
533: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U;
534: PETSC_EXTERN PetscLogEvent KSP_SolveTranspose;
535: PETSC_EXTERN PetscLogEvent KSP_MatSolve;
536: PETSC_EXTERN PetscLogEvent KSP_MatSolveTranspose;
538: PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat, IS, IS, IS, IS, MatReuse, Mat *, MatSchurComplementAinvType, MatReuse, Mat *);
539: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC, PetscBool *);
541: /*MC
542: KSPCheckDot - Checks if the result of a dot product used by the corresponding `KSP` contains infinity or NaN. These indicate that the previous
543: application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`.
545: Collective
547: Input Parameter:
548: . ksp - the linear solver `KSP` context.
550: Output Parameter:
551: . beta - the result of the inner product
553: Level: developer
555: Developer Notes:
556: Used to manage returning from `KSP` solvers collectively whose preconditioners have failed, possibly only a subset of MPI processes, in some way
558: It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products.
560: .seealso: `PCFailedReason`, `KSPConvergedReason`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckNorm()`, `KSPCheckSolve()`,
561: `KSPSetErrorIfNotConverged()`
562: M*/
563: #define KSPCheckDot(ksp, beta) \
564: do { \
565: if (PetscIsInfOrNanScalar(beta)) { \
566: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to infinity or NaN inner product"); \
567: { \
568: PCFailedReason pcreason; \
569: PetscCall(PCReduceFailedReason(ksp->pc)); \
570: PetscCall(PCGetFailedReason(ksp->pc, &pcreason)); \
571: PetscCall(VecFlag(ksp->vec_sol, pcreason)); \
572: if (pcreason) { \
573: ksp->reason = KSP_DIVERGED_PC_FAILED; \
574: } else { \
575: ksp->reason = KSP_DIVERGED_NANORINF; \
576: } \
577: PetscFunctionReturn(PETSC_SUCCESS); \
578: } \
579: } \
580: } while (0)
582: /*MC
583: KSPCheckNorm - Checks if the result of a norm used by the corresponding `KSP` contains `inf` or `NaN`. These indicate that the previous
584: application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`.
586: Collective
588: Input Parameter:
589: . ksp - the linear solver `KSP` context.
591: Output Parameter:
592: . beta - the result of the norm
594: Level: developer
596: Developer Notes:
597: Used to manage returning from `KSP` solvers collectively whose preconditioners have failed, possibly only a subset of MPI processes, in some way.
599: It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products.
601: .seealso: `PCFailedReason`, `KSPConvergedReason`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckDot()`, `KSPCheckSolve()`,
602: `KSPSetErrorIfNotConverged()`
603: M*/
604: #define KSPCheckNorm(ksp, beta) \
605: do { \
606: if (PetscIsInfOrNanReal(beta)) { \
607: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to infinity or NaN norm"); \
608: { \
609: PCFailedReason pcreason; \
610: PetscCall(PCReduceFailedReason(ksp->pc)); \
611: PetscCall(PCGetFailedReason(ksp->pc, &pcreason)); \
612: PetscCall(VecFlag(ksp->vec_sol, pcreason)); \
613: if (pcreason) { \
614: ksp->reason = KSP_DIVERGED_PC_FAILED; \
615: } else { \
616: ksp->reason = KSP_DIVERGED_NANORINF; \
617: } \
618: ksp->rnorm = beta; \
619: PetscFunctionReturn(PETSC_SUCCESS); \
620: } \
621: } \
622: } while (0)
624: PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]);
625: PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP, PetscInt, PetscReal *);