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;

 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: /*@C
 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: `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;

 52:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);

 55:   PetscViewerBinaryRead(viewer, &classid, 1, NULL, PETSC_INT);
 57:   PetscViewerBinaryRead(viewer, type, 256, NULL, PETSC_CHAR);
 58:   KSPSetType(newdm, type);
 59:   PetscTryTypeMethod(newdm, load, viewer);
 60:   KSPGetPC(newdm, &pc);
 61:   PCLoad(pc, viewer);
 62:   return 0;
 63: }

 65: #include <petscdraw.h>
 66: #if defined(PETSC_HAVE_SAWS)
 67: #include <petscviewersaws.h>
 68: #endif
 69: /*@C
 70:    KSPView - Prints the `KSP` data structure.

 72:    Collective

 74:    Input Parameters:
 75: +  ksp - the Krylov space context
 76: -  viewer - visualization context

 78:    Options Database Keys:
 79: .  -ksp_view - print the `KSP` data structure at the end of each `KSPSolve()` call

 81:    Notes:
 82:    The available visualization contexts include
 83: +     `PETSC_VIEWER_STDOUT_SELF` - standard output (default)
 84: -     `PETSC_VIEWER_STDOUT_WORLD` - synchronized standard
 85:          output where only the first processor opens
 86:          the file.  All other processors send their
 87:          data to the first processor to print.

 89:    The available formats include
 90: +     `PETSC_VIEWER_DEFAULT` - standard output (default)
 91: -     `PETSC_VIEWER_ASCII_INFO_DETAIL` - more verbose output for PCBJACOBI and PCASM

 93:    The user can open an alternative visualization context with
 94:    `PetscViewerASCIIOpen()` - output to a specified file.

 96:   In the debugger you can do call `KSPView(ksp,0)` to display the `KSP`. (The same holds for any PETSc object viewer).

 98:    Level: beginner

100: .seealso: `KSP`, `PetscViewer`, `PCView()`, `PetscViewerASCIIOpen()`
101: @*/
102: PetscErrorCode KSPView(KSP ksp, PetscViewer viewer)
103: {
104:   PetscBool iascii, isbinary, isdraw, isstring;
105: #if defined(PETSC_HAVE_SAWS)
106:   PetscBool issaws;
107: #endif

110:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp), &viewer);

114:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
115:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary);
116:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
117:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
118: #if defined(PETSC_HAVE_SAWS)
119:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSAWS, &issaws);
120: #endif
121:   if (iascii) {
122:     PetscObjectPrintClassNamePrefixType((PetscObject)ksp, viewer);
123:     PetscViewerASCIIPushTab(viewer);
124:     PetscTryTypeMethod(ksp, view, viewer);
125:     PetscViewerASCIIPopTab(viewer);
126:     if (ksp->guess_zero) {
127:       PetscViewerASCIIPrintf(viewer, "  maximum iterations=%" PetscInt_FMT ", initial guess is zero\n", ksp->max_it);
128:     } else {
129:       PetscViewerASCIIPrintf(viewer, "  maximum iterations=%" PetscInt_FMT ", nonzero initial guess\n", ksp->max_it);
130:     }
131:     if (ksp->guess_knoll) PetscViewerASCIIPrintf(viewer, "  using preconditioner applied to right hand side for initial guess\n");
132:     PetscViewerASCIIPrintf(viewer, "  tolerances:  relative=%g, absolute=%g, divergence=%g\n", (double)ksp->rtol, (double)ksp->abstol, (double)ksp->divtol);
133:     if (ksp->pc_side == PC_RIGHT) {
134:       PetscViewerASCIIPrintf(viewer, "  right preconditioning\n");
135:     } else if (ksp->pc_side == PC_SYMMETRIC) {
136:       PetscViewerASCIIPrintf(viewer, "  symmetric preconditioning\n");
137:     } else {
138:       PetscViewerASCIIPrintf(viewer, "  left preconditioning\n");
139:     }
140:     if (ksp->guess) {
141:       PetscViewerASCIIPushTab(viewer);
142:       KSPGuessView(ksp->guess, viewer);
143:       PetscViewerASCIIPopTab(viewer);
144:     }
145:     if (ksp->dscale) PetscViewerASCIIPrintf(viewer, "  diagonally scaled system\n");
146:     PetscViewerASCIIPrintf(viewer, "  using %s norm type for convergence test\n", KSPNormTypes[ksp->normtype]);
147:   } else if (isbinary) {
148:     PetscInt    classid = KSP_FILE_CLASSID;
149:     MPI_Comm    comm;
150:     PetscMPIInt rank;
151:     char        type[256];

153:     PetscObjectGetComm((PetscObject)ksp, &comm);
154:     MPI_Comm_rank(comm, &rank);
155:     if (rank == 0) {
156:       PetscViewerBinaryWrite(viewer, &classid, 1, PETSC_INT);
157:       PetscStrncpy(type, ((PetscObject)ksp)->type_name, 256);
158:       PetscViewerBinaryWrite(viewer, type, 256, PETSC_CHAR);
159:     }
160:     PetscTryTypeMethod(ksp, view, viewer);
161:   } else if (isstring) {
162:     const char *type;
163:     KSPGetType(ksp, &type);
164:     PetscViewerStringSPrintf(viewer, " KSPType: %-7.7s", type);
165:     PetscTryTypeMethod(ksp, view, viewer);
166:   } else if (isdraw) {
167:     PetscDraw draw;
168:     char      str[36];
169:     PetscReal x, y, bottom, h;
170:     PetscBool flg;

172:     PetscViewerDrawGetDraw(viewer, 0, &draw);
173:     PetscDrawGetCurrentPoint(draw, &x, &y);
174:     PetscObjectTypeCompare((PetscObject)ksp, KSPPREONLY, &flg);
175:     if (!flg) {
176:       PetscStrncpy(str, "KSP: ", sizeof(str));
177:       PetscStrlcat(str, ((PetscObject)ksp)->type_name, sizeof(str));
178:       PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h);
179:       bottom = y - h;
180:     } else {
181:       bottom = y;
182:     }
183:     PetscDrawPushCurrentPoint(draw, x, bottom);
184: #if defined(PETSC_HAVE_SAWS)
185:   } else if (issaws) {
186:     PetscMPIInt rank;
187:     const char *name;

189:     PetscObjectGetName((PetscObject)ksp, &name);
190:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
191:     if (!((PetscObject)ksp)->amsmem && rank == 0) {
192:       char dir[1024];

194:       PetscObjectViewSAWs((PetscObject)ksp, viewer);
195:       PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/its", name);
196:       SAWs_Register, (dir, &ksp->its, 1, SAWs_READ, SAWs_INT);
197:       if (!ksp->res_hist) KSPSetResidualHistory(ksp, NULL, PETSC_DECIDE, PETSC_TRUE);
198:       PetscSNPrintf(dir, 1024, "/PETSc/Objects/%s/res_hist", name);
199:       SAWs_Register, (dir, ksp->res_hist, 10, SAWs_READ, SAWs_DOUBLE);
200:     }
201: #endif
202:   } else PetscTryTypeMethod(ksp, view, viewer);
203:   if (ksp->pc) PCView(ksp->pc, viewer);
204:   if (isdraw) {
205:     PetscDraw draw;
206:     PetscViewerDrawGetDraw(viewer, 0, &draw);
207:     PetscDrawPopCurrentPoint(draw);
208:   }
209:   return 0;
210: }

212: /*@C
213:    KSPViewFromOptions - View a `KSP` object based on values in the options database

215:    Collective

217:    Input Parameters:
218: +  A - Krylov solver context
219: .  obj - Optional object
220: -  name - command line option

222:    Level: intermediate

224: .seealso: `KSP`, `KSPView`, `PetscObjectViewFromOptions()`, `KSPCreate()`
225: @*/
226: PetscErrorCode KSPViewFromOptions(KSP A, PetscObject obj, const char name[])
227: {
229:   PetscObjectViewFromOptions((PetscObject)A, obj, name);
230:   return 0;
231: }

233: /*@
234:    KSPSetNormType - Sets the norm that is used for convergence testing.

236:    Logically Collective

238:    Input Parameters:
239: +  ksp - Krylov solver context
240: -  normtype - one of
241: .vb
242:    KSP_NORM_NONE - skips computing the norm, this should generally only be used if you are using
243:                  the Krylov method as a smoother with a fixed small number of iterations.
244:                  Implicitly sets KSPConvergedSkip() as KSP convergence test.
245:                  Note that certain algorithms such as KSPGMRES ALWAYS require the norm calculation,
246:                  for these methods the norms are still computed, they are just not used in
247:                  the convergence test.
248:    KSP_NORM_PRECONDITIONED - the default for left preconditioned solves, uses the l2 norm
249:                  of the preconditioned residual P^{-1}(b - A x)
250:    KSP_NORM_UNPRECONDITIONED - uses the l2 norm of the true b - Ax residual.
251:    KSP_NORM_NATURAL - supported  by KSPCG, KSPCR, KSPCGNE, KSPCGS
252: .ve

254:    Options Database Key:
255: .   -ksp_norm_type <none,preconditioned,unpreconditioned,natural> - set `KSP` norm type

257:    Level: advanced

259:    Note:
260:    Not all combinations of preconditioner side (see `KSPSetPCSide()`) and norm type are supported by all Krylov methods.
261:    If only one is set, PETSc tries to automatically change the other to find a compatible pair.  If no such combination
262:    is supported, PETSc will generate an error.

264:    Developer Note:
265:    Supported combinations of norm and preconditioner side are set using `KSPSetSupportedNorm()`.

267: .seealso: `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSPConvergedSkip()`, `KSPSetCheckNormIteration()`, `KSPSetPCSide()`, `KSPGetPCSide()`, `KSPNormType`
268: @*/
269: PetscErrorCode KSPSetNormType(KSP ksp, KSPNormType normtype)
270: {
273:   ksp->normtype = ksp->normtype_set = normtype;
274:   return 0;
275: }

277: /*@
278:    KSPSetCheckNormIteration - Sets the first iteration at which the norm of the residual will be
279:      computed and used in the convergence test.

281:    Logically Collective

283:    Input Parameters:
284: +  ksp - Krylov solver context
285: -  it  - use -1 to check at all iterations

287:    Notes:
288:    Currently only works with `KSPCG`, `KSPBCGS` and `KSPIBCGS`

290:    Use `KSPSetNormType`(ksp,`KSP_NORM_NONE`) to never check the norm

292:    On steps where the norm is not computed, the previous norm is still in the variable, so if you run with, for example,
293:     -ksp_monitor the residual norm will appear to be unchanged for several iterations (though it is not really unchanged).
294:    Level: advanced

296: .seealso: `KSP`, `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSPConvergedSkip()`, `KSPSetNormType()`
297: @*/
298: PetscErrorCode KSPSetCheckNormIteration(KSP ksp, PetscInt it)
299: {
302:   ksp->chknorm = it;
303:   return 0;
304: }

306: /*@
307:    KSPSetLagNorm - Lags the residual norm calculation so that it is computed as part of the `MPI_Allreduce()` for
308:    computing the inner products for the next iteration.  This can reduce communication costs at the expense of doing
309:    one additional iteration.

311:    Logically Collective

313:    Input Parameters:
314: +  ksp - Krylov solver context
315: -  flg - `PETSC_TRUE` or `PETSC_FALSE`

317:    Options Database Keys:
318: .  -ksp_lag_norm - lag the calculated residual norm

320:    Level: advanced

322:    Notes:
323:    Currently only works with `KSPIBCGS`.

325:    Use `KSPSetNormType`(ksp,`KSP_NORM_NONE`) to never check the norm

327:    If you lag the norm and run with, for example, -ksp_monitor, the residual norm reported will be the lagged one.

329: .seealso: `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSPConvergedSkip()`, `KSPSetNormType()`, `KSPSetCheckNormIteration()`
330: @*/
331: PetscErrorCode KSPSetLagNorm(KSP ksp, PetscBool flg)
332: {
335:   ksp->lagnorm = flg;
336:   return 0;
337: }

339: /*@
340:    KSPSetSupportedNorm - Sets a norm and preconditioner side supported by a `KSP`

342:    Logically Collective

344:    Input Parameters:
345: +  ksp - Krylov method
346: .  normtype - supported norm type
347: .  pcside - preconditioner side that can be used with this norm
348: -  priority - positive integer preference for this combination; larger values have higher priority

350:    Level: developer

352:    Note:
353:    This function should be called from the implementation files `KSPCreate_XXX()` to declare
354:    which norms and preconditioner sides are supported. Users should not need to call this
355:    function.

357: .seealso: `KSP`, `KSPNormType`, `PCSide`, `KSPSetNormType()`, `KSPSetPCSide()`
358: @*/
359: PetscErrorCode KSPSetSupportedNorm(KSP ksp, KSPNormType normtype, PCSide pcside, PetscInt priority)
360: {
362:   ksp->normsupporttable[normtype][pcside] = priority;
363:   return 0;
364: }

366: PetscErrorCode KSPNormSupportTableReset_Private(KSP ksp)
367: {
368:   PetscMemzero(ksp->normsupporttable, sizeof(ksp->normsupporttable));
369:   ksp->pc_side  = ksp->pc_side_set;
370:   ksp->normtype = ksp->normtype_set;
371:   return 0;
372: }

374: PetscErrorCode KSPSetUpNorms_Private(KSP ksp, PetscBool errorifnotsupported, KSPNormType *normtype, PCSide *pcside)
375: {
376:   PetscInt i, j, best, ibest = 0, jbest = 0;

378:   best = 0;
379:   for (i = 0; i < KSP_NORM_MAX; i++) {
380:     for (j = 0; j < PC_SIDE_MAX; j++) {
381:       if ((ksp->normtype == KSP_NORM_DEFAULT || ksp->normtype == i) && (ksp->pc_side == PC_SIDE_DEFAULT || ksp->pc_side == j) && (ksp->normsupporttable[i][j] > best)) {
382:         best  = ksp->normsupporttable[i][j];
383:         ibest = i;
384:         jbest = j;
385:       }
386:     }
387:   }
388:   if (best < 1 && errorifnotsupported) {
392:     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]);
393:   }
394:   if (normtype) *normtype = (KSPNormType)ibest;
395:   if (pcside) *pcside = (PCSide)jbest;
396:   return 0;
397: }

399: /*@
400:    KSPGetNormType - Gets the norm that is used for convergence testing.

402:    Not Collective

404:    Input Parameter:
405: .  ksp - Krylov solver context

407:    Output Parameter:
408: .  normtype - norm that is used for convergence testing

410:    Level: advanced

412: .seealso: `KSPNormType`, `KSPSetNormType()`, `KSPConvergedSkip()`
413: @*/
414: PetscErrorCode KSPGetNormType(KSP ksp, KSPNormType *normtype)
415: {
418:   KSPSetUpNorms_Private(ksp, PETSC_TRUE, &ksp->normtype, &ksp->pc_side);
419:   *normtype = ksp->normtype;
420:   return 0;
421: }

423: #if defined(PETSC_HAVE_SAWS)
424: #include <petscviewersaws.h>
425: #endif

427: /*@
428:    KSPSetOperators - Sets the matrix associated with the linear system
429:    and a (possibly) different one from which the preconditioner will be built

431:    Collective

433:    Input Parameters:
434: +  ksp - the `KSP` context
435: .  Amat - the matrix that defines the linear system
436: -  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.

438:     Level: beginner

440:    Notes:
441:     If you know the operator Amat has a null space you can use `MatSetNullSpace()` and `MatSetTransposeNullSpace()` to supply the null
442:     space to Amat and the `KSP` solvers will automatically use that null space as needed during the solution process.

444:     All future calls to `KSPSetOperators()` must use the same size matrices!

446:     Passing a NULL for Amat or Pmat removes the matrix that is currently used.

448:     If you wish to replace either Amat or Pmat but leave the other one untouched then
449:     first call KSPGetOperators() to get the one you wish to keep, call `PetscObjectReference()`
450:     on it and then pass it back in in your call to `KSPSetOperators()`.

452:    Developer Notes:
453:    If the operators have NOT been set with `KSPSetOperators()` then the operators
454:       are created in the `PC` and returned to the user. In this case, if both operators
455:       mat and pmat are requested, two DIFFERENT operators will be returned. If
456:       only one is requested both operators in the PC will be the same (i.e. as
457:       if one had called `KSPSetOperators()` with the same argument for both `Mat`s).
458:       The user must set the sizes of the returned matrices and their type etc just
459:       as if the user created them with `MatCreate()`. For example,

461: .vb
462:          KSPGetOperators(ksp/pc,&mat,NULL); is equivalent to
463:            set size, type, etc of mat

465:          MatCreate(comm,&mat);
466:          KSP/PCSetOperators(ksp/pc,mat,mat);
467:          PetscObjectDereference((PetscObject)mat);
468:            set size, type, etc of mat

470:      and

472:          KSP/PCGetOperators(ksp/pc,&mat,&pmat); is equivalent to
473:            set size, type, etc of mat and pmat

475:          MatCreate(comm,&mat);
476:          MatCreate(comm,&pmat);
477:          KSP/PCSetOperators(ksp/pc,mat,pmat);
478:          PetscObjectDereference((PetscObject)mat);
479:          PetscObjectDereference((PetscObject)pmat);
480:            set size, type, etc of mat and pmat
481: .ve

483:     The rationale for this support is so that when creating a `TS`, `SNES`, or `KSP` the hierarchy
484:     of underlying objects (i.e. `SNES`, `KSP`, `PC`, `Mat`) and their livespans can be completely
485:     managed by the top most level object (i.e. the `TS`, `SNES`, or `KSP`). Another way to look
486:     at this is when you create a `SNES` you do not NEED to create a `KSP` and attach it to
487:     the `SNES` object (the `SNES` object manages it for you). Similarly when you create a `KSP`
488:     you do not need to attach a `PC` to it (the `KSP` object manages the `PC` object for you).
489:     Thus, why should YOU have to create the `Mat` and attach it to the `SNES`/`KSP`/`PC`, when
490:     it can be created for you?

492: .seealso: `KSP`, `Mat`, `KSPSolve()`, `KSPGetPC()`, `PCGetOperators()`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetComputeOperators()`, `KSPSetComputeInitialGuess()`, `KSPSetComputeRHS()`
493: @*/
494: PetscErrorCode KSPSetOperators(KSP ksp, Mat Amat, Mat Pmat)
495: {
501:   if (!ksp->pc) KSPGetPC(ksp, &ksp->pc);
502:   PCSetOperators(ksp->pc, Amat, Pmat);
503:   if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX; /* so that next solve call will call PCSetUp() on new matrix */
504:   return 0;
505: }

507: /*@
508:    KSPGetOperators - Gets the matrix associated with the linear system
509:    and a (possibly) different one used to construct the preconditioner.

511:    Collective

513:    Input Parameter:
514: .  ksp - the `KSP` context

516:    Output Parameters:
517: +  Amat - the matrix that defines the linear system
518: -  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.

520:     Level: intermediate

522:    Note:
523:     DOES NOT increase the reference counts of the matrix, so you should NOT destroy them.

525: .seealso: `KSP`, `KSPSolve()`, `KSPGetPC()`, `PCGetOperators()`, `PCSetOperators()`, `KSPSetOperators()`, `KSPGetOperatorsSet()`
526: @*/
527: PetscErrorCode KSPGetOperators(KSP ksp, Mat *Amat, Mat *Pmat)
528: {
530:   if (!ksp->pc) KSPGetPC(ksp, &ksp->pc);
531:   PCGetOperators(ksp->pc, Amat, Pmat);
532:   return 0;
533: }

535: /*@C
536:    KSPGetOperatorsSet - Determines if the matrix associated with the linear system and
537:    possibly a different one associated with the preconditioner have been set in the `KSP`.

539:    Not collective, though the results on all processes should be the same

541:    Input Parameter:
542: .  pc - the `KSP` context

544:    Output Parameters:
545: +  mat - the matrix associated with the linear system was set
546: -  pmat - matrix associated with the preconditioner was set, usually the same

548:    Level: intermediate

550:    Note:
551:    This routine exists because if you call `KSPGetOperators()` on a `KSP` that does not yet have operators they are
552:    automatically created in the call.

554: .seealso: `KSP`, `PCSetOperators()`, `KSPGetOperators()`, `KSPSetOperators()`, `PCGetOperators()`, `PCGetOperatorsSet()`
555: @*/
556: PetscErrorCode KSPGetOperatorsSet(KSP ksp, PetscBool *mat, PetscBool *pmat)
557: {
559:   if (!ksp->pc) KSPGetPC(ksp, &ksp->pc);
560:   PCGetOperatorsSet(ksp->pc, mat, pmat);
561:   return 0;
562: }

564: /*@C
565:    KSPSetPreSolve - Sets a function that is called at the beginning of each `KSPSolve()`

567:    Logically Collective

569:    Input Parameters:
570: +   ksp - the solver object
571: .   presolve - the function to call before the solve
572: -   prectx - any context needed by the function

574:    Calling sequence of presolve:
575: $  func(KSP ksp,Vec rhs,Vec x,void *ctx)

577: +  ksp - the `KSP` context
578: .  rhs - the right-hand side vector
579: .  x - the solution vector
580: -  ctx - optional user-provided context

582:    Level: developer

584: .seealso: `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPSetPostSolve()`, `PCEISENSTAT`
585: @*/
586: PetscErrorCode KSPSetPreSolve(KSP ksp, PetscErrorCode (*presolve)(KSP, Vec, Vec, void *), void *prectx)
587: {
589:   ksp->presolve = presolve;
590:   ksp->prectx   = prectx;
591:   return 0;
592: }

594: /*@C
595:    KSPSetPostSolve - Sets a function that is called at the end of each `KSPSolve()` (whether it converges or not)

597:    Logically Collective

599:    Input Parameters:
600: +   ksp - the solver object
601: .   postsolve - the function to call after the solve
602: -   postctx - any context needed by the function

604:    Calling sequence of postsolve:
605: $  func(KSP ksp,Vec rhs,Vec x,void *ctx)

607: +  ksp - the `KSP` context
608: .  rhs - the right-hand side vector
609: .  x - the solution vector
610: -  ctx - optional user-provided context

612:    Level: developer

614: .seealso: `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPSetPreSolve()`, `PCEISENSTAT`
615: @*/
616: PetscErrorCode KSPSetPostSolve(KSP ksp, PetscErrorCode (*postsolve)(KSP, Vec, Vec, void *), void *postctx)
617: {
619:   ksp->postsolve = postsolve;
620:   ksp->postctx   = postctx;
621:   return 0;
622: }

624: /*@
625:    KSPCreate - Creates the `KSP` context.

627:    Collective

629:    Input Parameter:
630: .  comm - MPI communicator

632:    Output Parameter:
633: .  ksp - location to put the `KSP` context

635:    Note:
636:    The default `KSPType` is `KSPGMRES` with a restart of 30, using modified Gram-Schmidt orthogonalization.

638:    Level: beginner

640: .seealso: [](chapter_ksp), `KSPSetUp()`, `KSPSolve()`, `KSPDestroy()`, `KSP`, `KSPGMRES`, `KSPType`
641: @*/
642: PetscErrorCode KSPCreate(MPI_Comm comm, KSP *inksp)
643: {
644:   KSP   ksp;
645:   void *ctx;

648:   *inksp = NULL;
649:   KSPInitializePackage();

651:   PetscHeaderCreate(ksp, KSP_CLASSID, "KSP", "Krylov Method", "KSP", comm, KSPDestroy, KSPView);

653:   ksp->max_it  = 10000;
654:   ksp->pc_side = ksp->pc_side_set = PC_SIDE_DEFAULT;
655:   ksp->rtol                       = 1.e-5;
656: #if defined(PETSC_USE_REAL_SINGLE)
657:   ksp->abstol = 1.e-25;
658: #else
659:   ksp->abstol = 1.e-50;
660: #endif
661:   ksp->divtol = 1.e4;

663:   ksp->chknorm  = -1;
664:   ksp->normtype = ksp->normtype_set = KSP_NORM_DEFAULT;
665:   ksp->rnorm                        = 0.0;
666:   ksp->its                          = 0;
667:   ksp->guess_zero                   = PETSC_TRUE;
668:   ksp->calc_sings                   = PETSC_FALSE;
669:   ksp->res_hist                     = NULL;
670:   ksp->res_hist_alloc               = NULL;
671:   ksp->res_hist_len                 = 0;
672:   ksp->res_hist_max                 = 0;
673:   ksp->res_hist_reset               = PETSC_TRUE;
674:   ksp->err_hist                     = NULL;
675:   ksp->err_hist_alloc               = NULL;
676:   ksp->err_hist_len                 = 0;
677:   ksp->err_hist_max                 = 0;
678:   ksp->err_hist_reset               = PETSC_TRUE;
679:   ksp->numbermonitors               = 0;
680:   ksp->numberreasonviews            = 0;
681:   ksp->setfromoptionscalled         = 0;
682:   ksp->nmax                         = PETSC_DECIDE;

684:   KSPConvergedDefaultCreate(&ctx);
685:   KSPSetConvergenceTest(ksp, KSPConvergedDefault, ctx, KSPConvergedDefaultDestroy);
686:   ksp->ops->buildsolution = KSPBuildSolutionDefault;
687:   ksp->ops->buildresidual = KSPBuildResidualDefault;

689:   ksp->vec_sol    = NULL;
690:   ksp->vec_rhs    = NULL;
691:   ksp->pc         = NULL;
692:   ksp->data       = NULL;
693:   ksp->nwork      = 0;
694:   ksp->work       = NULL;
695:   ksp->reason     = KSP_CONVERGED_ITERATING;
696:   ksp->setupstage = KSP_SETUP_NEW;

698:   KSPNormSupportTableReset_Private(ksp);

700:   *inksp = ksp;
701:   return 0;
702: }

704: /*@C
705:    KSPSetType - Builds the `KSP` datastructure for a particular `KSPType`

707:    Logically Collective

709:    Input Parameters:
710: +  ksp  - the Krylov space context
711: -  type - a known method

713:    Options Database Key:
714: .  -ksp_type  <method> - Sets the method; use -help for a list  of available methods (for instance, cg or gmres)

716:    Notes:
717:    See "petsc/include/petscksp.h" for available methods (for instance, `KSPCG` or `KSPGMRES`).

719:   Normally, it is best to use the `KSPSetFromOptions()` command and
720:   then set the `KSP` type from the options database rather than by using
721:   this routine.  Using the options database provides the user with
722:   maximum flexibility in evaluating the many different Krylov methods.
723:   The `KSPSetType()` routine is provided for those situations where it
724:   is necessary to set the iterative solver independently of the command
725:   line or options database.  This might be the case, for example, when
726:   the choice of iterative solver changes during the execution of the
727:   program, and the user's application is taking responsibility for
728:   choosing the appropriate method.  In other words, this routine is
729:   not for beginners.

731:   Level: intermediate

733:   Developer Note:
734:   `KSPRegister()` is used to add Krylov types to `KSPList` from which they are accessed by `KSPSetType()`.

736: .seealso: [](chapter_ksp), `PCSetType()`, `KSPType`, `KSPRegister()`, `KSPCreate()`, `KSP`
737: @*/
738: PetscErrorCode KSPSetType(KSP ksp, KSPType type)
739: {
740:   PetscBool match;
741:   PetscErrorCode (*r)(KSP);


746:   PetscObjectTypeCompare((PetscObject)ksp, type, &match);
747:   if (match) return 0;

749:   PetscFunctionListFind(KSPList, type, &r);
751:   /* Destroy the previous private KSP context */
752:   PetscTryTypeMethod(ksp, destroy);
753:   ksp->ops->destroy = NULL;

755:   /* Reinitialize function pointers in KSPOps structure */
756:   PetscMemzero(ksp->ops, sizeof(struct _KSPOps));
757:   ksp->ops->buildsolution = KSPBuildSolutionDefault;
758:   ksp->ops->buildresidual = KSPBuildResidualDefault;
759:   KSPNormSupportTableReset_Private(ksp);
760:   ksp->setupnewmatrix = PETSC_FALSE; // restore default (setup not called in case of new matrix)
761:   /* Call the KSPCreate_XXX routine for this particular Krylov solver */
762:   ksp->setupstage = KSP_SETUP_NEW;
763:   (*r)(ksp);
764:   PetscObjectChangeTypeName((PetscObject)ksp, type);
765:   return 0;
766: }

768: /*@C
769:    KSPGetType - Gets the `KSP` type as a string from the KSP object.

771:    Not Collective

773:    Input Parameter:
774: .  ksp - Krylov context

776:    Output Parameter:
777: .  name - name of the `KSP` method

779:    Level: intermediate

781: .seealso: [](chapter_ksp), `KSPType`, `KSP`, `KSPSetType()`
782: @*/
783: PetscErrorCode KSPGetType(KSP ksp, KSPType *type)
784: {
787:   *type = ((PetscObject)ksp)->type_name;
788:   return 0;
789: }

791: /*@C
792:   KSPRegister -  Adds a method, `KSPType`, to the Krylov subspace solver package.

794:    Not Collective

796:    Input Parameters:
797: +  name_solver - name of a new user-defined solver
798: -  routine_create - routine to create method context

800:    Level: advanced

802:    Note:
803:    `KSPRegister()` may be called multiple times to add several user-defined solvers.

805:    Sample usage:
806: .vb
807:    KSPRegister("my_solver",MySolverCreate);
808: .ve

810:    Then, your solver can be chosen with the procedural interface via
811: $    ` KSPSetType`(ksp,"my_solver")
812:    or at runtime via the option
813: $     -ksp_type my_solver

815: .seealso: [](chapter_ksp), `KSP`, `KSPType`, `KSPSetType`, `KSPRegisterAll()`
816: @*/
817: PetscErrorCode KSPRegister(const char sname[], PetscErrorCode (*function)(KSP))
818: {
819:   KSPInitializePackage();
820:   PetscFunctionListAdd(&KSPList, sname, function);
821:   return 0;
822: }

824: PetscErrorCode KSPMonitorMakeKey_Internal(const char name[], PetscViewerType vtype, PetscViewerFormat format, char key[])
825: {
826:   PetscStrncpy(key, name, PETSC_MAX_PATH_LEN);
827:   PetscStrlcat(key, ":", PETSC_MAX_PATH_LEN);
828:   PetscStrlcat(key, vtype, PETSC_MAX_PATH_LEN);
829:   PetscStrlcat(key, ":", PETSC_MAX_PATH_LEN);
830:   PetscStrlcat(key, PetscViewerFormats[format], PETSC_MAX_PATH_LEN);
831:   return 0;
832: }

834: /*@C
835:   KSPMonitorRegister -  Registers a Krylov subspace solver monitor routine that may be accessed with `KSPMonitorSetFromOptions()`

837:   Not Collective

839:   Input Parameters:
840: + name    - name of a new monitor routine
841: . vtype   - A `PetscViewerType` for the output
842: . format  - A `PetscViewerFormat` for the output
843: . monitor - Monitor routine
844: . create  - Creation routine, or NULL
845: - destroy - Destruction routine, or NULL

847:   Level: advanced

849:   Note:
850:   `KSPMonitorRegister()` may be called multiple times to add several user-defined monitors.

852:   Sample usage:
853: .vb
854:   KSPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
855: .ve

857:   Then, your monitor can be chosen with the procedural interface via
858: $     KSPMonitorSetFromOptions(ksp,"-ksp_monitor_my_monitor","my_monitor",NULL)
859:   or at runtime via the option
860: $     -ksp_monitor_my_monitor

862: .seealso: [](chapter_ksp), `KSP`, `KSPMonitorSet()`, `KSPMonitorRegisterAll()`, `KSPMonitorSetFromOptions()`
863: @*/
864: PetscErrorCode KSPMonitorRegister(const char name[], PetscViewerType vtype, PetscViewerFormat format, PetscErrorCode (*monitor)(KSP, PetscInt, PetscReal, PetscViewerAndFormat *), PetscErrorCode (*create)(PetscViewer, PetscViewerFormat, void *, PetscViewerAndFormat **), PetscErrorCode (*destroy)(PetscViewerAndFormat **))
865: {
866:   char key[PETSC_MAX_PATH_LEN];

868:   KSPInitializePackage();
869:   KSPMonitorMakeKey_Internal(name, vtype, format, key);
870:   PetscFunctionListAdd(&KSPMonitorList, key, monitor);
871:   if (create) PetscFunctionListAdd(&KSPMonitorCreateList, key, create);
872:   if (destroy) PetscFunctionListAdd(&KSPMonitorDestroyList, key, destroy);
873:   return 0;
874: }