Actual source code: itcreate.c
1: /*
2: The basic KSP routines, Create, View etc. are here.
3: */
4: #include <petsc/private/kspimpl.h>
6: /* Logging support */
7: PetscClassId KSP_CLASSID;
8: PetscClassId DMKSP_CLASSID;
9: PetscClassId KSPGUESS_CLASSID;
10: PetscLogEvent KSP_GMRESOrthogonalization, KSP_SetUp, KSP_Solve, KSP_SolveTranspose, KSP_MatSolve, KSP_MatSolveTranspose;
12: /*
13: Contains the list of registered KSP routines
14: */
15: PetscFunctionList KSPList = NULL;
16: PetscBool KSPRegisterAllCalled = PETSC_FALSE;
18: /*
19: Contains the list of registered KSP monitors
20: */
21: PetscFunctionList KSPMonitorList = NULL;
22: PetscFunctionList KSPMonitorCreateList = NULL;
23: PetscFunctionList KSPMonitorDestroyList = NULL;
24: PetscBool KSPMonitorRegisterAllCalled = PETSC_FALSE;
26: /*@
27: KSPLoad - Loads a `KSP` that has been stored in a `PETSCVIEWERBINARY` with `KSPView()`.
29: Collective
31: Input Parameters:
32: + newdm - the newly loaded `KSP`, this needs to have been created with `KSPCreate()` or
33: some related function before a call to `KSPLoad()`.
34: - viewer - binary file viewer, obtained from `PetscViewerBinaryOpen()`
36: Level: intermediate
38: Note:
39: The type is determined by the data in the file, any type set into the `KSP` before this call is ignored.
41: .seealso: [](ch_ksp), `KSP`, `PetscViewerBinaryOpen()`, `KSPView()`, `MatLoad()`, `VecLoad()`
42: @*/
43: PetscErrorCode KSPLoad(KSP newdm, PetscViewer viewer)
44: {
45: PetscBool isbinary;
46: PetscInt classid;
47: char type[256];
48: PC pc;
50: PetscFunctionBegin;
53: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
54: PetscCheck(isbinary, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerBinaryOpen()");
56: PetscCall(PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT));
57: PetscCheck(classid == KSP_FILE_CLASSID, PetscObjectComm((PetscObject)newdm), PETSC_ERR_ARG_WRONG, "Not KSP next in file");
58: PetscCall(PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR));
59: PetscCall(KSPSetType(newdm, type));
60: PetscTryTypeMethod(newdm, load, viewer);
61: PetscCall(KSPGetPC(newdm, &pc));
62: PetscCall(PCLoad(pc, viewer));
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: #include <petscdraw.h>
67: #if defined(PETSC_HAVE_SAWS)
68: #include <petscviewersaws.h>
69: #endif
70: /*@
71: KSPView - Prints the various parameters currently set in the `KSP` object. For example, the convergence tolerances and `KSPType`.
72: Also views the `PC` and `Mat` contained by the `KSP` with `PCView()` and `MatView()`.
74: Collective
76: Input Parameters:
77: + ksp - the Krylov space context
78: - viewer - visualization context
80: Options Database Key:
81: . -ksp_view - print the `KSP` data structure at the end of each `KSPSolve()` call
83: Level: beginner
85: Notes:
86: The available visualization contexts include
87: + `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
88: - `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
89: output where only the first processor opens
90: the file. All other processors send their
91: data to the first processor to print.
93: The available formats include
94: + `PETSC_VIEWER_DEFAULT` - standard output (default)
95: - `PETSC_VIEWER_ASCII_INFO_DETAIL` - more verbose output for `PCBJACOBI` and `PCASM`
97: The user can open an alternative visualization context with
98: `PetscViewerASCIIOpen()` - output to a specified file.
100: Use `KSPViewFromOptions()` to allow the user to select many different `PetscViewerType` and formats from the options database.
102: In the debugger you can do call `KSPView(ksp,0)` to display the `KSP`. (The same holds for any PETSc object viewer).
104: .seealso: [](ch_ksp), `KSP`, `PetscViewer`, `PCView()`, `PetscViewerASCIIOpen()`, `KSPViewFromOptions()`
105: @*/
106: PetscErrorCode KSPView(KSP ksp, PetscViewer viewer)
107: {
108: PetscBool isascii, isbinary, isdraw, isstring;
109: #if defined(PETSC_HAVE_SAWS)
110: PetscBool issaws;
111: #endif
113: PetscFunctionBegin;
115: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp), &viewer));
117: PetscCheckSameComm(ksp, 1, viewer, 2);
119: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
120: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
121: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
122: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
123: #if defined(PETSC_HAVE_SAWS)
124: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws));
125: #endif
126: if (isascii) {
127: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)ksp, viewer));
128: PetscCall(PetscViewerASCIIPushTab(viewer));
129: PetscTryTypeMethod(ksp, view, viewer);
130: PetscCall(PetscViewerASCIIPopTab(viewer));
131: if (ksp->guess_zero) {
132: PetscCall(PetscViewerASCIIPrintf(viewer, " maximum iterations=%" PetscInt_FMT ", initial guess is zero\n", ksp->max_it));
133: } else {
134: PetscCall(PetscViewerASCIIPrintf(viewer, " maximum iterations=%" PetscInt_FMT ", nonzero initial guess\n", ksp->max_it));
135: }
136: if (ksp->min_it) PetscCall(PetscViewerASCIIPrintf(viewer, " minimum iterations=%" PetscInt_FMT "\n", ksp->min_it));
137: if (ksp->guess_knoll) PetscCall(PetscViewerASCIIPrintf(viewer, " using preconditioner applied to right-hand side for initial guess\n"));
138: PetscCall(PetscViewerASCIIPrintf(viewer, " tolerances: relative=%g, absolute=%g, divergence=%g\n", (double)ksp->rtol, (double)ksp->abstol, (double)ksp->divtol));
139: if (ksp->pc_side == PC_RIGHT) {
140: PetscCall(PetscViewerASCIIPrintf(viewer, " right preconditioning\n"));
141: } else if (ksp->pc_side == PC_SYMMETRIC) {
142: PetscCall(PetscViewerASCIIPrintf(viewer, " symmetric preconditioning\n"));
143: } else {
144: PetscCall(PetscViewerASCIIPrintf(viewer, " left preconditioning\n"));
145: }
146: if (ksp->guess) {
147: PetscCall(PetscViewerASCIIPushTab(viewer));
148: PetscCall(KSPGuessView(ksp->guess, viewer));
149: PetscCall(PetscViewerASCIIPopTab(viewer));
150: }
151: if (ksp->dscale) PetscCall(PetscViewerASCIIPrintf(viewer, " diagonally scaled system\n"));
152: if (ksp->converged == KSPConvergedSkip || ksp->normtype == KSP_NORM_NONE) PetscCall(PetscViewerASCIIPrintf(viewer, " not checking for convergence\n"));
153: else PetscCall(PetscViewerASCIIPrintf(viewer, " using %s norm type for convergence test\n", KSPNormTypes[ksp->normtype]));
154: } else if (isbinary) {
155: PetscInt classid = KSP_FILE_CLASSID;
156: MPI_Comm comm;
157: PetscMPIInt rank;
158: char type[256];
160: PetscCall(PetscObjectGetComm((PetscObject)ksp, &comm));
161: PetscCallMPI(MPI_Comm_rank(comm, &rank));
162: if (rank == 0) {
163: PetscCall(PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT));
164: PetscCall(PetscStrncpy(type, ((PetscObject)ksp)->type_name, 256));
165: PetscCall(PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR));
166: }
167: PetscTryTypeMethod(ksp, view, viewer);
168: } else if (isstring) {
169: const char *type;
170: PetscCall(KSPGetType(ksp, &type));
171: PetscCall(PetscViewerStringSPrintf(viewer, " KSPType: %-7.7s", type));
172: PetscTryTypeMethod(ksp, view, viewer);
173: } else if (isdraw) {
174: PetscDraw draw;
175: char str[36];
176: PetscReal x, y, bottom, h;
177: PetscBool flg;
179: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
180: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
181: PetscCall(PetscObjectTypeCompare((PetscObject)ksp, KSPPREONLY, &flg));
182: if (!flg) {
183: PetscCall(PetscStrncpy(str, "KSP: ", sizeof(str)));
184: PetscCall(PetscStrlcat(str, ((PetscObject)ksp)->type_name, sizeof(str)));
185: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
186: bottom = y - h;
187: } else {
188: bottom = y;
189: }
190: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
191: #if defined(PETSC_HAVE_SAWS)
192: } else if (issaws) {
193: PetscMPIInt rank;
194: const char *name;
196: PetscCall(PetscObjectGetName((PetscObject)ksp, &name));
197: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
198: if (!((PetscObject)ksp)->amsmem && rank == 0) {
199: char dir[1024];
201: PetscCall(PetscObjectViewSAWs((PetscObject)ksp, viewer));
202: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/its", name));
203: PetscCallSAWs(SAWs_Register, (dir, &ksp->its, 1, SAWs_READ, SAWs_INT));
204: if (!ksp->res_hist) PetscCall(KSPSetResidualHistory(ksp, NULL, PETSC_DECIDE, PETSC_TRUE));
205: PetscCall(PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/res_hist", name));
206: PetscCallSAWs(SAWs_Register, (dir, ksp->res_hist, 10, SAWs_READ, SAWs_DOUBLE));
207: }
208: #endif
209: } else PetscTryTypeMethod(ksp, view, viewer);
210: if (ksp->pc) PetscCall(PCView(ksp->pc, viewer));
211: if (isdraw) {
212: PetscDraw draw;
213: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
214: PetscCall(PetscDrawPopCurrentPoint(draw));
215: }
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: /*@
220: KSPViewFromOptions - View (print) a `KSP` object based on values in the options database. Also views the `PC` and `Mat` contained by the `KSP`
221: with `PCView()` and `MatView()`.
223: Collective
225: Input Parameters:
226: + A - Krylov solver context
227: . obj - Optional object that provides the options prefix used to query the options database
228: - name - command line option
230: Options Database Key:
231: . -name [viewertype][:...] - option name and values. See `PetscObjectViewFromOptions()` for the possible arguments
233: Level: intermediate
235: .seealso: [](ch_ksp), `KSP`, `KSPView()`, `PetscObjectViewFromOptions()`, `KSPCreate()`
236: @*/
237: PetscErrorCode KSPViewFromOptions(KSP A, PetscObject obj, const char name[])
238: {
239: PetscFunctionBegin;
241: PetscCall(PetscObjectViewFromOptions((PetscObject)A, obj, name));
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: /*@
246: KSPSetNormType - Sets the type of residual norm that is used for convergence testing in `KSPSolve()` for the given `KSP` context
248: Logically Collective
250: Input Parameters:
251: + ksp - Krylov solver context
252: - normtype - one of
253: .vb
254: KSP_NORM_NONE - skips computing the norm, this should generally only be used if you are using
255: the Krylov method as a smoother with a fixed small number of iterations.
256: Implicitly sets `KSPConvergedSkip()` as the `KSP` convergence test.
257: Note that certain algorithms such as `KSPGMRES` ALWAYS require the norm calculation,
258: for these methods the norms are still computed, they are just not used in
259: the convergence test.
260: KSP_NORM_PRECONDITIONED - the default for left-preconditioned solves, uses the 2-norm
261: of the preconditioned residual $B^{-1}(b - A x)$.
262: KSP_NORM_UNPRECONDITIONED - uses the 2-norm of the true $b - Ax$ residual.
263: KSP_NORM_NATURAL - uses the $A$ norm of the true $b - Ax$ residual; supported by `KSPCG`, `KSPCR`, `KSPCGNE`, `KSPCGS`
264: .ve
266: Options Database Key:
267: . -ksp_norm_type (none|preconditioned|unpreconditioned|natural) - set `KSP` norm type
269: Level: advanced
271: Notes:
272: The norm is always of the equations residual $\| b - A x^n \|$ (or an approximation to that norm), they are never a norm of the error in the equation.
274: Not all combinations of preconditioner side (see `KSPSetPCSide()`) and norm types are supported by all Krylov methods.
275: If only one is set, PETSc tries to automatically change the other to find a compatible pair. If no such combination
276: is supported, PETSc will generate an error.
278: Developer Note:
279: Supported combinations of norm and preconditioner side are set using `KSPSetSupportedNorm()` for each `KSPType`.
281: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSPConvergedSkip()`, `KSPSetCheckNormIteration()`, `KSPSetPCSide()`, `KSPGetPCSide()`, `KSPNormType`
282: @*/
283: PetscErrorCode KSPSetNormType(KSP ksp, KSPNormType normtype)
284: {
285: PetscFunctionBegin;
288: ksp->normtype = ksp->normtype_set = normtype;
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: /*@
293: KSPSetCheckNormIteration - Sets the first iteration at which the norm of the residual will be
294: computed and used in the convergence test of `KSPSolve()` for the given `KSP` context
296: Logically Collective
298: Input Parameters:
299: + ksp - Krylov solver context
300: - it - use -1 to check at all iterations
302: Level: advanced
304: Notes:
305: Currently only works with `KSPCG`, `KSPBCGS` and `KSPIBCGS`
307: Use `KSPSetNormType`(ksp,`KSP_NORM_NONE`) to never check the norm
309: On steps where the norm is not computed, the previous norm is still in the variable, so if you run with, for example,
310: `-ksp_monitor` the residual norm will appear to be unchanged for several iterations (though it is not really unchanged).
312: Certain methods such as `KSPGMRES` always compute the residual norm, this routine will not change that computation, but it will
313: prevent the computed norm from being checked.
315: .seealso: [](ch_ksp), `KSP`, `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSPConvergedSkip()`, `KSPSetNormType()`, `KSPSetLagNorm()`
316: @*/
317: PetscErrorCode KSPSetCheckNormIteration(KSP ksp, PetscInt it)
318: {
319: PetscFunctionBegin;
322: ksp->chknorm = it;
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: /*@
327: KSPSetLagNorm - Lags the residual norm calculation so that it is computed as part of the `MPI_Allreduce()` used for
328: computing the inner products needed for the next iteration.
330: Logically Collective
332: Input Parameters:
333: + ksp - Krylov solver context
334: - flg - `PETSC_TRUE` or `PETSC_FALSE`
336: Options Database Key:
337: . -ksp_lag_norm - lag the calculated residual norm
339: Level: advanced
341: Notes:
342: Currently only works with `KSPIBCGS`.
344: This can reduce communication costs at the expense of doing
345: one additional iteration because the norm used in the convergence test of `KSPSolve()` is one iteration behind the actual
346: current residual norm (which has not yet been computed due to the lag).
348: Use `KSPSetNormType`(ksp,`KSP_NORM_NONE`) to never check the norm
350: If you lag the norm and run with, for example, `-ksp_monitor`, the residual norm reported will be the lagged one.
352: `KSPSetCheckNormIteration()` is an alternative way of avoiding the expense of computing the residual norm at each iteration.
354: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSPConvergedSkip()`, `KSPSetNormType()`, `KSPSetCheckNormIteration()`
355: @*/
356: PetscErrorCode KSPSetLagNorm(KSP ksp, PetscBool flg)
357: {
358: PetscFunctionBegin;
361: ksp->lagnorm = flg;
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: /*@
366: KSPSetSupportedNorm - Sets a norm and preconditioner side supported by a `KSPType`
368: Logically Collective
370: Input Parameters:
371: + ksp - Krylov method
372: . normtype - supported norm type of the type `KSPNormType`
373: . pcside - preconditioner side, of the type `PCSide` that can be used with this `KSPNormType`
374: - priority - positive integer preference for this combination; larger values have higher priority
376: Level: developer
378: Notes:
379: This function should be called from the implementation files `KSPCreate_XXX()` to declare
380: which norms and preconditioner sides are supported. Users should not call this
381: function.
383: This function can be called multiple times for each combination of `KSPNormType` and `PCSide`
384: the `KSPType` supports
386: .seealso: [](ch_ksp), `KSP`, `KSPNormType`, `PCSide`, `KSPSetNormType()`, `KSPSetPCSide()`
387: @*/
388: PetscErrorCode KSPSetSupportedNorm(KSP ksp, KSPNormType normtype, PCSide pcside, PetscInt priority)
389: {
390: PetscFunctionBegin;
392: ksp->normsupporttable[normtype][pcside] = priority;
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: static PetscErrorCode KSPNormSupportTableReset_Private(KSP ksp)
397: {
398: PetscFunctionBegin;
399: PetscCall(PetscMemzero(ksp->normsupporttable, sizeof(ksp->normsupporttable)));
400: ksp->pc_side = ksp->pc_side_set;
401: ksp->normtype = ksp->normtype_set;
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: PetscErrorCode KSPSetUpNorms_Private(KSP ksp, PetscBool errorifnotsupported, KSPNormType *normtype, PCSide *pcside)
406: {
407: PetscInt i, j, best, ibest = 0, jbest = 0;
409: PetscFunctionBegin;
410: best = 0;
411: for (i = 0; i < KSP_NORM_MAX; i++) {
412: for (j = 0; j < PC_SIDE_MAX; j++) {
413: if ((ksp->normtype == KSP_NORM_DEFAULT || ksp->normtype == i) && (ksp->pc_side == PC_SIDE_DEFAULT || ksp->pc_side == j) && ksp->normsupporttable[i][j] > best) {
414: best = ksp->normsupporttable[i][j];
415: ibest = i;
416: jbest = j;
417: }
418: }
419: }
420: if (best < 1 && errorifnotsupported) {
421: PetscCheck(ksp->normtype != KSP_NORM_DEFAULT || ksp->pc_side != PC_SIDE_DEFAULT, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "The %s KSP implementation did not call KSPSetSupportedNorm()", ((PetscObject)ksp)->type_name);
422: PetscCheck(ksp->normtype != KSP_NORM_DEFAULT, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "KSP %s does not support preconditioner side %s", ((PetscObject)ksp)->type_name, PCSides[ksp->pc_side]);
423: PetscCheck(ksp->pc_side != PC_SIDE_DEFAULT, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "KSP %s does not support norm type %s", ((PetscObject)ksp)->type_name, KSPNormTypes[ksp->normtype]);
424: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "KSP %s does not support norm type %s with preconditioner side %s", ((PetscObject)ksp)->type_name, KSPNormTypes[ksp->normtype], PCSides[ksp->pc_side]);
425: }
426: if (normtype) *normtype = (KSPNormType)ibest;
427: if (pcside) *pcside = (PCSide)jbest;
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: /*@
432: KSPGetNormType - Gets the `KSPNormType` that is used for convergence testing during `KSPSolve()` for this `KSP` context
434: Not Collective
436: Input Parameter:
437: . ksp - Krylov solver context
439: Output Parameter:
440: . normtype - the `KSPNormType` that is used for convergence testing
442: Level: advanced
444: .seealso: [](ch_ksp), `KSPNormType`, `KSPSetNormType()`, `KSPConvergedSkip()`
445: @*/
446: PetscErrorCode KSPGetNormType(KSP ksp, KSPNormType *normtype)
447: {
448: PetscFunctionBegin;
450: PetscAssertPointer(normtype, 2);
451: PetscCall(KSPSetUpNorms_Private(ksp, PETSC_TRUE, &ksp->normtype, &ksp->pc_side));
452: *normtype = ksp->normtype;
453: PetscFunctionReturn(PETSC_SUCCESS);
454: }
456: #if defined(PETSC_HAVE_SAWS)
457: #include <petscviewersaws.h>
458: #endif
460: /*@
461: KSPSetOperators - Sets the matrix associated with the linear system
462: and a (possibly) different one from which the preconditioner will be built into the `KSP` context. The matrix will then be used during `KSPSolve()`
464: Collective
466: Input Parameters:
467: + ksp - the `KSP` context
468: . Amat - the matrix that defines the linear system
469: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as `Amat`.
471: Level: beginner
473: Notes:
474: .vb
475: KSPSetOperators(ksp, Amat, Pmat);
476: .ve
477: is the same as
478: .vb
479: KSPGetPC(ksp, &pc);
480: PCSetOperators(pc, Amat, Pmat);
481: .ve
482: and is equivalent to
483: .vb
484: PCCreate(PetscObjectComm((PetscObject)ksp), &pc);
485: PCSetOperators(pc, Amat, Pmat);
486: KSPSetPC(ksp, pc);
487: .ve
489: If you know the operator `Amat` has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
490: space to `Amat` and the `KSP` solvers will automatically use that null space as needed during the solution process.
492: All future calls to `KSPSetOperators()` must use the same size matrices, unless `KSPReset()` is called!
494: Passing a `NULL` for `Amat` or `Pmat` removes the matrix that is currently being used from the `KSP` context.
496: If you wish to replace either `Amat` or `Pmat` but leave the other one untouched then
497: first call `KSPGetOperators()` to get the one you wish to keep, call `PetscObjectReference()`
498: on it and then pass it back in your call to `KSPSetOperators()`.
500: Developer Notes:
501: If the operators have NOT been set with `KSPSetOperators()` then the operators
502: are created in the `PC` and returned to the user. In this case, if both operators
503: mat and pmat are requested, two DIFFERENT operators will be returned. If
504: only one is requested both operators in the `PC` will be the same (i.e. as
505: if one had called `KSPSetOperators()` with the same argument for both `Mat`s).
506: The user must set the sizes of the returned matrices and their type etc just
507: as if the user created them with `MatCreate()`. For example,
509: .vb
510: KSPGetOperators(ksp/pc,&mat,NULL); is equivalent to
511: set size, type, etc of mat
513: MatCreate(comm,&mat);
514: KSP/PCSetOperators(ksp/pc,mat,mat);
515: PetscObjectDereference((PetscObject)mat);
516: set size, type, etc of mat
518: and
520: KSP/PCGetOperators(ksp/pc,&mat,&pmat); is equivalent to
521: set size, type, etc of mat and pmat
523: MatCreate(comm,&mat);
524: MatCreate(comm,&pmat);
525: KSP/PCSetOperators(ksp/pc,mat,pmat);
526: PetscObjectDereference((PetscObject)mat);
527: PetscObjectDereference((PetscObject)pmat);
528: set size, type, etc of mat and pmat
529: .ve
531: The rationale for this support is so that when creating a `TS`, `SNES`, or `KSP` the hierarchy
532: of underlying objects (i.e. `SNES`, `KSP`, `PC`, `Mat`) and their lifespans can be completely
533: managed by the top most level object (i.e. the `TS`, `SNES`, or `KSP`). Another way to look
534: at this is when you create a `SNES` you do not NEED to create a `KSP` and attach it to
535: the `SNES` object (the `SNES` object manages it for you). Similarly when you create a `KSP`
536: you do not need to attach a `PC` to it (the `KSP` object manages the `PC` object for you).
537: Thus, why should YOU have to create the `Mat` and attach it to the `SNES`/`KSP`/`PC`, when
538: it can be created for you?
540: .seealso: [](ch_ksp), `KSP`, `Mat`, `KSPSolve()`, `KSPGetPC()`, `PCGetOperators()`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetComputeOperators()`, `KSPSetComputeInitialGuess()`, `KSPSetComputeRHS()`
541: @*/
542: PetscErrorCode KSPSetOperators(KSP ksp, Mat Amat, Mat Pmat)
543: {
544: PetscFunctionBegin;
548: if (Amat) PetscCheckSameComm(ksp, 1, Amat, 2);
549: if (Pmat) PetscCheckSameComm(ksp, 1, Pmat, 3);
550: if (!ksp->pc) PetscCall(KSPGetPC(ksp, &ksp->pc));
551: PetscCall(PCSetOperators(ksp->pc, Amat, Pmat));
552: if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX; /* so that next solve call will call PCSetUp() on new matrix */
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: /*@
557: KSPGetOperators - Gets the matrix associated with the linear system
558: and a (possibly) different one used to construct the preconditioner from the `KSP` context
560: Collective
562: Input Parameter:
563: . ksp - the `KSP` context
565: Output Parameters:
566: + Amat - the matrix that defines the linear system
567: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as `Amat`.
569: Level: intermediate
571: Notes:
572: If `KSPSetOperators()` has not been called then the `KSP` object will attempt to automatically create the matrix `Amat` and return it
574: Use `KSPGetOperatorsSet()` to determine if matrices have been provided.
576: DOES NOT increase the reference counts of the matrix, so you should NOT destroy them.
578: .seealso: [](ch_ksp), `KSP`, `KSPSolve()`, `KSPGetPC()`, `PCSetOperators()`, `KSPSetOperators()`, `KSPGetOperatorsSet()`
579: @*/
580: PetscErrorCode KSPGetOperators(KSP ksp, Mat *Amat, Mat *Pmat)
581: {
582: PetscFunctionBegin;
584: if (!ksp->pc) PetscCall(KSPGetPC(ksp, &ksp->pc));
585: PetscCall(PCGetOperators(ksp->pc, Amat, Pmat));
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: /*@
590: KSPGetOperatorsSet - Determines if the matrix associated with the linear system and
591: possibly a different one from which the preconditioner will be built have been set in the `KSP` with `KSPSetOperators()`
593: Not Collective, though the results on all processes will be the same
595: Input Parameter:
596: . ksp - the `KSP` context
598: Output Parameters:
599: + mat - the matrix associated with the linear system was set
600: - pmat - matrix from which the preconditioner will be built, usually the same as `mat` was set
602: Level: intermediate
604: Note:
605: This routine exists because if you call `KSPGetOperators()` on a `KSP` that does not yet have operators they are
606: automatically created in the call.
608: .seealso: [](ch_ksp), `KSP`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperators()`, `PCGetOperatorsSet()`
609: @*/
610: PetscErrorCode KSPGetOperatorsSet(KSP ksp, PetscBool *mat, PetscBool *pmat)
611: {
612: PetscFunctionBegin;
614: if (!ksp->pc) PetscCall(KSPGetPC(ksp, &ksp->pc));
615: PetscCall(PCGetOperatorsSet(ksp->pc, mat, pmat));
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: /*@C
620: KSPSetPreSolve - Sets a function that is called at the beginning of each `KSPSolve()`. Used in conjunction with `KSPSetPostSolve()`.
622: Logically Collective
624: Input Parameters:
625: + ksp - the solver object
626: . presolve - the function to call before the solve, see` KSPPSolveFn`
627: - ctx - an optional context needed by the function
629: Level: developer
631: Notes:
632: The function provided here `presolve` is used to modify the right hand side, and possibly the matrix, of the linear system to be solved.
633: The function provided with `KSPSetPostSolve()` then modifies the resulting solution of that linear system to obtain the correct solution
634: to the initial linear system.
636: The functions `PCPreSolve()` and `PCPostSolve()` provide a similar functionality and are used, for example with `PCEISENSTAT`.
638: .seealso: [](ch_ksp), `KSPPSolveFn`, `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPSetPostSolve()`, `PCEISENSTAT`, `PCPreSolve()`, `PCPostSolve()`
639: @*/
640: PetscErrorCode KSPSetPreSolve(KSP ksp, KSPPSolveFn *presolve, PetscCtx ctx)
641: {
642: PetscFunctionBegin;
644: ksp->presolve = presolve;
645: ksp->prectx = ctx;
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: /*@C
650: KSPSetPostSolve - Sets a function that is called at the end of each `KSPSolve()` (whether it converges or not). Used in conjunction with `KSPSetPreSolve()`.
652: Logically Collective
654: Input Parameters:
655: + ksp - the solver object
656: . postsolve - the function to call after the solve, see` KSPPSolveFn`
657: - ctx - an optional context needed by the function
659: Level: developer
661: .seealso: [](ch_ksp), `KSPPSolveFn`, `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPSetPreSolve()`, `PCEISENSTAT`
662: @*/
663: PetscErrorCode KSPSetPostSolve(KSP ksp, KSPPSolveFn *postsolve, PetscCtx ctx)
664: {
665: PetscFunctionBegin;
667: ksp->postsolve = postsolve;
668: ksp->postctx = ctx;
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: /*@
673: KSPSetNestLevel - sets the amount of nesting the `KSP` has. That is the number of levels of `KSP` above this `KSP` in a linear solve.
675: Collective
677: Input Parameters:
678: + ksp - the `KSP`
679: - level - the nest level
681: Level: developer
683: Note:
684: For example, the `KSP` in each block of a `KSPBJACOBI` has a level of 1, while the outer `KSP` has a level of 0.
686: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPGetNestLevel()`, `PCSetKSPNestLevel()`, `PCGetKSPNestLevel()`
687: @*/
688: PetscErrorCode KSPSetNestLevel(KSP ksp, PetscInt level)
689: {
690: PetscFunctionBegin;
693: ksp->nestlevel = level;
694: PetscFunctionReturn(PETSC_SUCCESS);
695: }
697: /*@
698: KSPGetNestLevel - gets the amount of nesting the `KSP` has
700: Not Collective
702: Input Parameter:
703: . ksp - the `KSP`
705: Output Parameter:
706: . level - the nest level
708: Level: developer
710: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPSetNestLevel()`, `PCSetKSPNestLevel()`, `PCGetKSPNestLevel()`
711: @*/
712: PetscErrorCode KSPGetNestLevel(KSP ksp, PetscInt *level)
713: {
714: PetscFunctionBegin;
716: PetscAssertPointer(level, 2);
717: *level = ksp->nestlevel;
718: PetscFunctionReturn(PETSC_SUCCESS);
719: }
721: /*@
722: KSPCreate - Creates the `KSP` context. This `KSP` context is used in PETSc to solve linear systems with `KSPSolve()`
724: Collective
726: Input Parameter:
727: . comm - MPI communicator
729: Output Parameter:
730: . inksp - location to put the `KSP` context
732: Level: beginner
734: Note:
735: The default `KSPType` is `KSPGMRES` with a restart of 30, using modified Gram-Schmidt orthogonalization. The `KSPType` may be
736: changed with `KSPSetType()`
738: .seealso: [](ch_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`, `KSPSetType()`
739: @*/
740: PetscErrorCode KSPCreate(MPI_Comm comm, KSP *inksp)
741: {
742: KSP ksp;
743: PetscCtx ctx;
745: PetscFunctionBegin;
746: PetscAssertPointer(inksp, 2);
747: PetscCall(KSPInitializePackage());
749: PetscCall(PetscHeaderCreate(ksp, KSP_CLASSID, "KSP", "Krylov Method", "KSP", comm, KSPDestroy, KSPView));
750: ksp->default_max_it = ksp->max_it = 10000;
751: ksp->pc_side = ksp->pc_side_set = PC_SIDE_DEFAULT;
753: ksp->default_rtol = ksp->rtol = 1.e-5;
754: ksp->default_abstol = ksp->abstol = PetscDefined(USE_REAL_SINGLE) ? 1.e-25 : 1.e-50;
755: ksp->default_divtol = ksp->divtol = 1.e4;
757: ksp->chknorm = -1;
758: ksp->normtype = ksp->normtype_set = KSP_NORM_DEFAULT;
759: ksp->rnorm = 0.0;
760: ksp->its = 0;
761: ksp->guess_zero = PETSC_TRUE;
762: ksp->calc_sings = PETSC_FALSE;
763: ksp->res_hist = NULL;
764: ksp->res_hist_alloc = NULL;
765: ksp->res_hist_len = 0;
766: ksp->res_hist_max = 0;
767: ksp->res_hist_reset = PETSC_TRUE;
768: ksp->err_hist = NULL;
769: ksp->err_hist_alloc = NULL;
770: ksp->err_hist_len = 0;
771: ksp->err_hist_max = 0;
772: ksp->err_hist_reset = PETSC_TRUE;
773: ksp->numbermonitors = 0;
774: ksp->numberreasonviews = 0;
775: ksp->setfromoptionscalled = 0;
776: ksp->nmax = PETSC_DECIDE;
778: PetscCall(KSPConvergedDefaultCreate(&ctx));
779: PetscCall(KSPSetConvergenceTest(ksp, KSPConvergedDefault, ctx, KSPConvergedDefaultDestroy));
780: ksp->ops->buildsolution = KSPBuildSolutionDefault;
781: ksp->ops->buildresidual = KSPBuildResidualDefault;
783: ksp->vec_sol = NULL;
784: ksp->vec_rhs = NULL;
785: ksp->pc = NULL;
786: ksp->data = NULL;
787: ksp->nwork = 0;
788: ksp->work = NULL;
789: ksp->reason = KSP_CONVERGED_ITERATING;
790: ksp->setupstage = KSP_SETUP_NEW;
792: PetscCall(KSPNormSupportTableReset_Private(ksp));
794: *inksp = ksp;
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
798: /*@
799: KSPSetType - Sets the algorithm/method to be used to solve the linear system with the given `KSP`
801: Logically Collective
803: Input Parameters:
804: + ksp - the Krylov space context
805: - type - a known method
807: Options Database Key:
808: . -ksp_type type - Sets the method; see `KSPType`
810: Level: intermediate
812: Notes:
813: See `KSPType` for available methods (for instance, `KSPCG` or `KSPGMRES`).
815: Normally, it is best to use the `KSPSetFromOptions()` command and
816: then set the `KSP` type from the options database rather than by using
817: this routine. Using the options database provides the user with
818: maximum flexibility in evaluating the many different Krylov methods.
819: The `KSPSetType()` routine is provided for those situations where it
820: is necessary to set the iterative solver independently of the command
821: line or options database. This might be the case, for example, when
822: the choice of iterative solver changes during the execution of the
823: program, and the user's application is taking responsibility for
824: choosing the appropriate method. In other words, this routine is
825: not for beginners.
827: Developer Note:
828: `KSPRegister()` is used to add Krylov types to `KSPList` from which they are accessed by `KSPSetType()`.
830: .seealso: [](ch_ksp), `PCSetType()`, `KSPType`, `KSPRegister()`, `KSPCreate()`, `KSP`
831: @*/
832: PetscErrorCode KSPSetType(KSP ksp, KSPType type)
833: {
834: PetscBool match;
835: PetscErrorCode (*r)(KSP);
837: PetscFunctionBegin;
839: PetscAssertPointer(type, 2);
841: PetscCall(PetscObjectTypeCompare((PetscObject)ksp, type, &match));
842: if (match) PetscFunctionReturn(PETSC_SUCCESS);
844: PetscCall(PetscFunctionListFind(KSPList, type, &r));
845: PetscCheck(r, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested KSP type %s", type);
846: /* Destroy the previous private KSP context */
847: PetscTryTypeMethod(ksp, destroy);
849: /* Reinitialize function pointers in KSPOps structure */
850: PetscCall(PetscMemzero(ksp->ops, sizeof(struct _KSPOps)));
851: ksp->ops->buildsolution = KSPBuildSolutionDefault;
852: ksp->ops->buildresidual = KSPBuildResidualDefault;
853: PetscCall(KSPNormSupportTableReset_Private(ksp));
854: ksp->converged_neg_curve = PETSC_FALSE; // restore default
855: ksp->setupnewmatrix = PETSC_FALSE; // restore default (setup not called in case of new matrix)
856: /* Call the KSPCreate_XXX routine for this particular Krylov solver */
857: ksp->setupstage = KSP_SETUP_NEW;
858: ksp->guess_not_read = PETSC_FALSE; // restore default
859: PetscCall((*r)(ksp));
860: PetscCall(PetscObjectChangeTypeName((PetscObject)ksp, type));
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: /*@
865: KSPGetType - Gets the `KSP` type as a string from the `KSP` object.
867: Not Collective
869: Input Parameter:
870: . ksp - Krylov context
872: Output Parameter:
873: . type - name of the `KSP` method
875: Level: intermediate
877: .seealso: [](ch_ksp), `KSPType`, `KSP`, `KSPSetType()`
878: @*/
879: PetscErrorCode KSPGetType(KSP ksp, KSPType *type)
880: {
881: PetscFunctionBegin;
883: PetscAssertPointer(type, 2);
884: *type = ((PetscObject)ksp)->type_name;
885: PetscFunctionReturn(PETSC_SUCCESS);
886: }
888: /*@C
889: KSPRegister - Adds a method, `KSPType`, to the Krylov subspace solver package.
891: Not Collective, No Fortran Support
893: Input Parameters:
894: + sname - name of a new user-defined solver
895: - function - routine to create method
897: Level: advanced
899: Note:
900: `KSPRegister()` may be called multiple times to add several user-defined solvers.
902: Example Usage:
903: .vb
904: KSPRegister("my_solver", MySolverCreate);
905: .ve
907: Then, your solver can be chosen with the procedural interface via
908: .vb
909: KSPSetType(ksp, "my_solver")
910: .ve
911: or at runtime via the option `-ksp_type my_solver`
913: .seealso: [](ch_ksp), `KSP`, `KSPType`, `KSPSetType`, `KSPRegisterAll()`
914: @*/
915: PetscErrorCode KSPRegister(const char sname[], PetscErrorCode (*function)(KSP))
916: {
917: PetscFunctionBegin;
918: PetscCall(KSPInitializePackage());
919: PetscCall(PetscFunctionListAdd(&KSPList, sname, function));
920: PetscFunctionReturn(PETSC_SUCCESS);
921: }
923: PetscErrorCode KSPMonitorMakeKey_Internal(const char name[], PetscViewerType vtype, PetscViewerFormat format, char key[])
924: {
925: PetscFunctionBegin;
926: PetscCall(PetscStrncpy(key, name, PETSC_MAX_PATH_LEN));
927: PetscCall(PetscStrlcat(key, ":", PETSC_MAX_PATH_LEN));
928: PetscCall(PetscStrlcat(key, vtype, PETSC_MAX_PATH_LEN));
929: PetscCall(PetscStrlcat(key, ":", PETSC_MAX_PATH_LEN));
930: PetscCall(PetscStrlcat(key, PetscViewerFormats[format], PETSC_MAX_PATH_LEN));
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }
934: /*@C
935: KSPMonitorRegister - Registers a Krylov subspace solver monitor routine that may be accessed with `KSPMonitorSetFromOptions()`
937: Not Collective
939: Input Parameters:
940: + name - name of a new monitor type
941: . vtype - A `PetscViewerType` for the output
942: . format - A `PetscViewerFormat` for the output
943: . monitor - Monitor routine, see `KSPMonitorRegisterFn`
944: . create - Creation routine, or `NULL`
945: - destroy - Destruction routine, or `NULL`
947: Level: advanced
949: Notes:
950: `KSPMonitorRegister()` may be called multiple times to add several user-defined monitors.
952: The calling sequence for the given function matches the calling sequence used by `KSPMonitorFn` functions passed to `KSPMonitorSet()` with the additional
953: requirement that its final argument be a `PetscViewerAndFormat`.
955: Example Usage:
956: .vb
957: KSPMonitorRegister("my_monitor", PETSCVIEWERASCII, PETSC_VIEWER_ASCII_INFO_DETAIL, MyMonitor, NULL, NULL);
958: .ve
960: Then, your monitor can be chosen with the procedural interface via
961: .vb
962: KSPMonitorSetFromOptions(ksp, "-ksp_monitor_my_monitor", "my_monitor", NULL)
963: .ve
964: or at runtime via the option `-ksp_monitor_my_monitor`
966: .seealso: [](ch_ksp), `KSP`, `KSPMonitorSet()`, `KSPMonitorRegisterAll()`, `KSPMonitorSetFromOptions()`
967: @*/
968: PetscErrorCode KSPMonitorRegister(const char name[], PetscViewerType vtype, PetscViewerFormat format, KSPMonitorRegisterFn *monitor, KSPMonitorRegisterCreateFn *create, KSPMonitorRegisterDestroyFn *destroy)
969: {
970: char key[PETSC_MAX_PATH_LEN];
972: PetscFunctionBegin;
973: PetscCall(KSPInitializePackage());
974: PetscCall(KSPMonitorMakeKey_Internal(name, vtype, format, key));
975: PetscCall(PetscFunctionListAdd(&KSPMonitorList, key, monitor));
976: if (create) PetscCall(PetscFunctionListAdd(&KSPMonitorCreateList, key, create));
977: if (destroy) PetscCall(PetscFunctionListAdd(&KSPMonitorDestroyList, key, destroy));
978: PetscFunctionReturn(PETSC_SUCCESS);
979: }