1 #pragma once
2
3 #include <petscksp.h>
4 #include <petscds.h>
5 #include <petsc/private/petscimpl.h>
6
7 /* SUBMANSEC = KSP */
8
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);
15
16 typedef struct _KSPOps *KSPOps;
17
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 };
38
39 typedef struct _KSPGuessOps *KSPGuessOps;
40
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 };
51
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 };
62
63 PETSC_EXTERN PetscErrorCode KSPGuessCreate_Fischer(KSPGuess);
64 PETSC_EXTERN PetscErrorCode KSPGuessCreate_POD(KSPGuess);
65
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;
76
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 */
104
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 */
119
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. */
123
124 PetscInt nmax; /* maximum number of right-hand sides to be handled simultaneously */
125
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 */
132
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 */
139
140 KSPConvergenceTestFn *converged;
141 PetscCtxDestroyFn *convergeddestroy;
142 void *cnvP;
143
144 PetscCtx ctx; /* optional user-defined context */
145
146 PC pc;
147
148 void *data; /* holder for misc stuff associated with a particular iterative solver */
149
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;
153
154 /* ----------------Default work-area management -------------------- */
155 PetscInt nwork;
156 Vec *work;
157
158 KSPSetUpStage setupstage;
159 PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */
160
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 */
163
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;
170
171 KSPNormType normtype; /* type of norm used for convergence tests */
172
173 PCSide pc_side_set; /* PC type set explicitly by user */
174 KSPNormType normtype_set; /* Norm type set explicitly by user */
175
176 /* Allow diagonally scaling the matrix before computing the preconditioner or using
177 the Krylov method. Note this is NOT just Jacobi preconditioning */
178
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;
184
185 /* Allow declaring convergence when negative curvature is detected */
186 PetscBool converged_neg_curve;
187
188 PetscInt setfromoptionscalled;
189 PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */
190
191 PetscErrorCode (*presolve)(KSP, Vec, Vec, void *);
192 PetscErrorCode (*postsolve)(KSP, Vec, Vec, void *);
193 void *prectx, *postctx;
194
195 PetscInt nestlevel; /* how many levels of nesting does the KSP have */
196 };
197
198 typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
199 PetscReal coef;
200 PetscReal bnrm;
201 } KSPDynTolCtx;
202
203 typedef struct {
204 PetscBool initialrtol; /* default relative residual decrease is computed from initial residual, not rhs */
205 PetscBool mininitialrtol; /* default relative residual decrease is computed from min of initial residual and rhs */
206 PetscBool convmaxits; /* if true, the convergence test returns KSP_CONVERGED_ITS if the maximum number of iterations is reached */
207 Vec work;
208 } KSPConvergedDefaultCtx;
209
KSPLogResidualHistory(KSP ksp,PetscReal norm)210 static inline PetscErrorCode KSPLogResidualHistory(KSP ksp, PetscReal norm)
211 {
212 PetscFunctionBegin;
213 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
214 if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) ksp->res_hist[ksp->res_hist_len++] = norm;
215 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
216 PetscFunctionReturn(PETSC_SUCCESS);
217 }
218
KSPLogErrorHistory(KSP ksp)219 static inline PetscErrorCode KSPLogErrorHistory(KSP ksp)
220 {
221 DM dm;
222
223 PetscFunctionBegin;
224 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
225 PetscCall(KSPGetDM(ksp, &dm));
226 if (dm && ksp->err_hist && ksp->err_hist_max > ksp->err_hist_len) {
227 PetscSimplePointFn *exactSol;
228 void *exactCtx;
229 PetscDS ds;
230 Vec u;
231 PetscReal error;
232 PetscInt Nf;
233
234 PetscCall(KSPBuildSolution(ksp, NULL, &u));
235 /* TODO Was needed to correct for Newton solution, but I just need to set a solution */
236 //PetscCall(VecScale(u, -1.0));
237 /* TODO Case when I have a solution */
238 if (0) {
239 PetscCall(DMGetDS(dm, &ds));
240 PetscCall(PetscDSGetNumFields(ds, &Nf));
241 PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle number of fields %" PetscInt_FMT " > 1 right now", Nf);
242 PetscCall(PetscDSGetExactSolution(ds, 0, &exactSol, &exactCtx));
243 PetscCall(DMComputeL2FieldDiff(dm, 0.0, &exactSol, &exactCtx, u, &error));
244 } else {
245 /* The null solution A 0 = 0 */
246 PetscCall(VecNorm(u, NORM_2, &error));
247 }
248 ksp->err_hist[ksp->err_hist_len++] = error;
249 }
250 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
251 PetscFunctionReturn(PETSC_SUCCESS);
252 }
253
KSPNoisyHash_Private(PetscInt xx)254 static inline PetscScalar KSPNoisyHash_Private(PetscInt xx)
255 {
256 unsigned int x = (unsigned int)xx;
257 x = ((x >> 16) ^ x) * 0x45d9f3b;
258 x = ((x >> 16) ^ x) * 0x45d9f3b;
259 x = ((x >> 16) ^ x);
260 return (PetscScalar)(((PetscInt64)x - 2147483648) * 5.e-10); /* center around zero, scaled about -1. to 1.*/
261 }
262
KSPSetNoisy_Private(Mat A,Vec v)263 static inline PetscErrorCode KSPSetNoisy_Private(Mat A, Vec v)
264 {
265 PetscScalar *a;
266 PetscInt n, istart;
267 MatNullSpace nullsp = NULL;
268
269 PetscFunctionBegin;
270 if (A) PetscCall(MatGetNullSpace(A, &nullsp));
271 PetscCall(VecGetOwnershipRange(v, &istart, NULL));
272 PetscCall(VecGetLocalSize(v, &n));
273 PetscCall(VecGetArrayWrite(v, &a));
274 for (PetscInt i = 0; i < n; ++i) a[i] = KSPNoisyHash_Private(i + istart);
275 PetscCall(VecRestoreArrayWrite(v, &a));
276 if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, v));
277 PetscFunctionReturn(PETSC_SUCCESS);
278 }
279
280 PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP, PetscBool, KSPNormType *, PCSide *);
281
282 PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP, PetscInt, const PetscReal *, const PetscReal *);
283
284 typedef struct _p_DMKSP *DMKSP;
285 typedef struct _DMKSPOps *DMKSPOps;
286 struct _DMKSPOps {
287 KSPComputeOperatorsFn *computeoperators;
288 KSPComputeRHSFn *computerhs;
289 KSPComputeInitialGuessFn *computeinitialguess;
290 PetscErrorCode (*destroy)(DMKSP *);
291 PetscErrorCode (*duplicate)(DMKSP, DMKSP);
292 };
293
294 /*S
295 DMKSP - Object held by a `DM` that contains all the callback functions and their contexts needed by a `KSP`
296
297 Level: developer
298
299 Notes:
300 Users provides callback functions and their contexts to `KSP` using, for example, `KSPSetComputeRHS()`. These values are stored
301 in a `DMKSP` that is contained in the `DM` associated with the `KSP`. If no `DM` was provided by
302 the user with `KSPSetDM()` it is automatically created by `KSPGetDM()` with `DMShellCreate()`.
303
304 Users very rarely need to worked directly with the `DMKSP` object, rather they work with the `KSP` and the `DM` they created
305
306 Multiple `DM` can share a single `DMKSP`, often each `DM` is associated with
307 a grid refinement level. `DMGetDMKSP()` returns the `DMKSP` associated with a `DM`. `DMGetDMKSPWrite()` returns a unique
308 `DMKSP` that is only associated with the current `DM`, making a copy of the shared `DMKSP` if needed (copy-on-write).
309
310 Developer Notes:
311 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`.
312 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,
313 but we need to also support different callbacks and contexts on each level. The copy-on-write approach of `DMGetDMKSPWrite()` makes this possible.
314
315 The `originaldm` inside the `DMKSP` is NOT reference counted (to prevent a reference count loop between a `DM` and a `DMKSP`).
316 The `DM` on which this context was first created is cached here to implement one-way
317 copy-on-write. When `DMGetDMKSPWrite()` sees a request using a different `DM`, it makes a copy of the `TSDM`. Thus, if a user
318 only interacts directly with one level, e.g., using `TSSetIFunction()`, then coarse levels of a multilevel item
319 integrator are built, then the user changes the routine with another call to `TSSetIFunction()`, it automatically
320 propagates to all the levels. If instead, they get out a specific level and set the function on that level,
321 subsequent changes to the original level will no longer propagate to that level.
322
323 .seealso: [](ch_ts), `KSP`, `KSPCreate()`, `DM`, `DMGetDMKSPWrite()`, `DMGetDMKSP()`, `DMSNES`, `DMTS`, `DMKSPSetComputeOperators()`, `DMKSPGetComputeOperators()`,
324 `DMKSPSetComputeRHS()`, `DMKSPSetComputeInitialGuess()`
325 S*/
326 struct _p_DMKSP {
327 PETSCHEADER(struct _DMKSPOps);
328 void *operatorsctx;
329 void *rhsctx;
330 void *initialguessctx;
331 void *data;
332
333 /* See developer note for `DMKSP` above */
334 DM originaldm;
335
336 void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */
337 };
338 PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM, DMKSP *);
339 PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM, DMKSP *);
340 PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM, DM);
341
342 /*
343 These allow the various Krylov methods to apply to either the linear system or its transpose.
344 */
KSP_RemoveNullSpace(KSP ksp,Vec y)345 static inline PetscErrorCode KSP_RemoveNullSpace(KSP ksp, Vec y)
346 {
347 PetscFunctionBegin;
348 if (ksp->pc_side == PC_LEFT) {
349 Mat A;
350 MatNullSpace nullsp;
351
352 PetscCall(PCGetOperators(ksp->pc, &A, NULL));
353 PetscCall(MatGetNullSpace(A, &nullsp));
354 if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y));
355 }
356 PetscFunctionReturn(PETSC_SUCCESS);
357 }
358
KSP_RemoveNullSpaceTranspose(KSP ksp,Vec y)359 static inline PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp, Vec y)
360 {
361 PetscFunctionBegin;
362 if (ksp->pc_side == PC_LEFT) {
363 Mat A;
364 MatNullSpace nullsp;
365
366 PetscCall(PCGetOperators(ksp->pc, &A, NULL));
367 PetscCall(MatGetTransposeNullSpace(A, &nullsp));
368 if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y));
369 }
370 PetscFunctionReturn(PETSC_SUCCESS);
371 }
372
KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y)373 static inline PetscErrorCode KSP_MatMult(KSP ksp, Mat A, Vec x, Vec y)
374 {
375 PetscFunctionBegin;
376 if (ksp->transpose_solve) PetscCall(MatMultTranspose(A, x, y));
377 else PetscCall(MatMult(A, x, y));
378 PetscFunctionReturn(PETSC_SUCCESS);
379 }
380
KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y)381 static inline PetscErrorCode KSP_MatMultTranspose(KSP ksp, Mat A, Vec x, Vec y)
382 {
383 PetscFunctionBegin;
384 if (ksp->transpose_solve) PetscCall(MatMult(A, x, y));
385 else PetscCall(MatMultTranspose(A, x, y));
386 PetscFunctionReturn(PETSC_SUCCESS);
387 }
388
KSP_MatMultHermitianTranspose(KSP ksp,Mat A,Vec x,Vec y)389 static inline PetscErrorCode KSP_MatMultHermitianTranspose(KSP ksp, Mat A, Vec x, Vec y)
390 {
391 PetscFunctionBegin;
392 if (!ksp->transpose_solve) PetscCall(MatMultHermitianTranspose(A, x, y));
393 else {
394 Vec w;
395
396 PetscCall(VecDuplicate(x, &w));
397 PetscCall(VecCopy(x, w));
398 PetscCall(VecConjugate(w));
399 PetscCall(MatMult(A, w, y));
400 PetscCall(VecDestroy(&w));
401 PetscCall(VecConjugate(y));
402 }
403 PetscFunctionReturn(PETSC_SUCCESS);
404 }
405
KSP_PCApply(KSP ksp,Vec x,Vec y)406 static inline PetscErrorCode KSP_PCApply(KSP ksp, Vec x, Vec y)
407 {
408 PetscFunctionBegin;
409 if (ksp->transpose_solve) {
410 PetscCall(PCApplyTranspose(ksp->pc, x, y));
411 PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y));
412 } else {
413 PetscCall(PCApply(ksp->pc, x, y));
414 PetscCall(KSP_RemoveNullSpace(ksp, y));
415 }
416 PetscFunctionReturn(PETSC_SUCCESS);
417 }
418
KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y)419 static inline PetscErrorCode KSP_PCApplyTranspose(KSP ksp, Vec x, Vec y)
420 {
421 PetscFunctionBegin;
422 if (ksp->transpose_solve) {
423 PetscCall(PCApply(ksp->pc, x, y));
424 PetscCall(KSP_RemoveNullSpace(ksp, y));
425 } else {
426 PetscCall(PCApplyTranspose(ksp->pc, x, y));
427 PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y));
428 }
429 PetscFunctionReturn(PETSC_SUCCESS);
430 }
431
KSP_PCApplyHermitianTranspose(KSP ksp,Vec x,Vec y)432 static inline PetscErrorCode KSP_PCApplyHermitianTranspose(KSP ksp, Vec x, Vec y)
433 {
434 PetscFunctionBegin;
435 PetscCall(VecConjugate(x));
436 PetscCall(KSP_PCApplyTranspose(ksp, x, y));
437 PetscCall(VecConjugate(x));
438 PetscCall(VecConjugate(y));
439 PetscFunctionReturn(PETSC_SUCCESS);
440 }
441
KSP_PCMatApply(KSP ksp,Mat X,Mat Y)442 static inline PetscErrorCode KSP_PCMatApply(KSP ksp, Mat X, Mat Y)
443 {
444 PetscFunctionBegin;
445 if (ksp->transpose_solve) PetscCall(PCMatApplyTranspose(ksp->pc, X, Y));
446 else PetscCall(PCMatApply(ksp->pc, X, Y));
447 PetscFunctionReturn(PETSC_SUCCESS);
448 }
449
KSP_PCMatApplyTranspose(KSP ksp,Mat X,Mat Y)450 static inline PetscErrorCode KSP_PCMatApplyTranspose(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 }
457
KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w)458 static inline PetscErrorCode KSP_PCApplyBAorAB(KSP ksp, Vec x, Vec y, Vec w)
459 {
460 PetscFunctionBegin;
461 if (ksp->transpose_solve) {
462 PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w));
463 PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y));
464 } else {
465 PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w));
466 PetscCall(KSP_RemoveNullSpace(ksp, y));
467 }
468 PetscFunctionReturn(PETSC_SUCCESS);
469 }
470
KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w)471 static inline PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp, Vec x, Vec y, Vec w)
472 {
473 PetscFunctionBegin;
474 if (ksp->transpose_solve) PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w));
475 else PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w));
476 PetscFunctionReturn(PETSC_SUCCESS);
477 }
478
479 PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization;
480 PETSC_EXTERN PetscLogEvent KSP_SetUp;
481 PETSC_EXTERN PetscLogEvent KSP_Solve;
482 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0;
483 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1;
484 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2;
485 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3;
486 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4;
487 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S;
488 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L;
489 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U;
490 PETSC_EXTERN PetscLogEvent KSP_SolveTranspose;
491 PETSC_EXTERN PetscLogEvent KSP_MatSolve;
492 PETSC_EXTERN PetscLogEvent KSP_MatSolveTranspose;
493
494 PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat, IS, IS, IS, IS, MatReuse, Mat *, MatSchurComplementAinvType, MatReuse, Mat *);
495 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC, PetscBool *);
496
497 /*MC
498 KSPCheckDot - Checks if the result of a dot product used by the corresponding `KSP` contains infinity or NaN. These indicate that the previous
499 application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`.
500
501 Collective
502
503 Input Parameter:
504 . ksp - the linear solver `KSP` context.
505
506 Output Parameter:
507 . beta - the result of the inner product
508
509 Level: developer
510
511 Developer Notes:
512 Used to manage returning from `KSP` solvers collectively whose preconditioners have failed, possibly only a subset of MPI processes, in some way
513
514 It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products.
515
516 .seealso: `PCFailedReason`, `KSPConvergedReason`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckNorm()`, `KSPCheckSolve()`,
517 `KSPSetErrorIfNotConverged()`
518 M*/
519 #define KSPCheckDot(ksp, beta) \
520 do { \
521 if (PetscIsInfOrNanScalar(beta)) { \
522 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to infinity or NaN inner product"); \
523 { \
524 PCFailedReason pcreason; \
525 PetscCall(PCReduceFailedReason(ksp->pc)); \
526 PetscCall(PCGetFailedReason(ksp->pc, &pcreason)); \
527 PetscCall(VecFlag(ksp->vec_sol, pcreason)); \
528 if (pcreason) { \
529 ksp->reason = KSP_DIVERGED_PC_FAILED; \
530 } else { \
531 ksp->reason = KSP_DIVERGED_NANORINF; \
532 } \
533 PetscFunctionReturn(PETSC_SUCCESS); \
534 } \
535 } \
536 } while (0)
537
538 /*MC
539 KSPCheckNorm - Checks if the result of a norm used by the corresponding `KSP` contains `inf` or `NaN`. These indicate that the previous
540 application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`.
541
542 Collective
543
544 Input Parameter:
545 . ksp - the linear solver `KSP` context.
546
547 Output Parameter:
548 . beta - the result of the norm
549
550 Level: developer
551
552 Developer Notes:
553 Used to manage returning from `KSP` solvers collectively whose preconditioners have failed, possibly only a subset of MPI processes, in some way.
554
555 It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products.
556
557 .seealso: `PCFailedReason`, `KSPConvergedReason`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckDot()`, `KSPCheckSolve()`,
558 `KSPSetErrorIfNotConverged()`
559 M*/
560 #define KSPCheckNorm(ksp, beta) \
561 do { \
562 if (PetscIsInfOrNanReal(beta)) { \
563 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to infinity or NaN norm"); \
564 { \
565 PCFailedReason pcreason; \
566 PetscCall(PCReduceFailedReason(ksp->pc)); \
567 PetscCall(PCGetFailedReason(ksp->pc, &pcreason)); \
568 PetscCall(VecFlag(ksp->vec_sol, pcreason)); \
569 if (pcreason) { \
570 ksp->reason = KSP_DIVERGED_PC_FAILED; \
571 } else { \
572 ksp->reason = KSP_DIVERGED_NANORINF; \
573 } \
574 ksp->rnorm = beta; \
575 PetscFunctionReturn(PETSC_SUCCESS); \
576 } \
577 } \
578 } while (0)
579
580 PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]);
581 PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP, PetscInt, PetscReal *);
582