Actual source code: itcreate.c

petsc-master 2020-07-14
Report Typos and Errors
  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: /*@C
 19:   KSPLoad - Loads a KSP that has been stored in binary  with KSPView().

 21:   Collective on viewer

 23:   Input Parameters:
 24: + newdm - the newly loaded KSP, this needs to have been created with KSPCreate() or
 25:            some related function before a call to KSPLoad().
 26: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()

 28:    Level: intermediate

 30:   Notes:
 31:    The type is determined by the data in the file, any type set into the KSP before this call is ignored.

 33:   Notes for advanced users:
 34:   Most users should not need to know the details of the binary storage
 35:   format, since KSPLoad() and KSPView() completely hide these details.
 36:   But for anyone who's interested, the standard binary matrix storage
 37:   format is
 38: .vb
 39:      has not yet been determined
 40: .ve

 42: .seealso: PetscViewerBinaryOpen(), KSPView(), MatLoad(), VecLoad()
 43: @*/
 44: PetscErrorCode  KSPLoad(KSP newdm, PetscViewer viewer)
 45: {
 47:   PetscBool      isbinary;
 48:   PetscInt       classid;
 49:   char           type[256];
 50:   PC             pc;

 55:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
 56:   if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");

 58:   PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
 59:   if (classid != KSP_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not KSP next in file");
 60:   PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
 61:   KSPSetType(newdm, type);
 62:   if (newdm->ops->load) {
 63:     (*newdm->ops->load)(newdm,viewer);
 64:   }
 65:   KSPGetPC(newdm,&pc);
 66:   PCLoad(pc,viewer);
 67:   return(0);
 68: }

 70:  #include <petscdraw.h>
 71: #if defined(PETSC_HAVE_SAWS)
 72:  #include <petscviewersaws.h>
 73: #endif
 74: /*@C
 75:    KSPView - Prints the KSP data structure.

 77:    Collective on ksp

 79:    Input Parameters:
 80: +  ksp - the Krylov space context
 81: -  viewer - visualization context

 83:    Options Database Keys:
 84: .  -ksp_view - print the ksp data structure at the end of a KSPSolve call

 86:    Note:
 87:    The available visualization contexts include
 88: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
 89: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
 90:          output where only the first processor opens
 91:          the file.  All other processors send their
 92:          data to the first processor to print.

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

 97:    Level: beginner

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

111:   if (!viewer) {
112:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
113:   }

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

158:     PetscObjectGetComm((PetscObject)ksp,&comm);
159:     MPI_Comm_rank(comm,&rank);
160:     if (!rank) {
161:       PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT);
162:       PetscStrncpy(type,((PetscObject)ksp)->type_name,256);
163:       PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR);
164:     }
165:     if (ksp->ops->view) {
166:       (*ksp->ops->view)(ksp,viewer);
167:     }
168:   } else if (isstring) {
169:     const char *type;
170:     KSPGetType(ksp,&type);
171:     PetscViewerStringSPrintf(viewer," KSPType: %-7.7s",type);
172:     if (ksp->ops->view) {(*ksp->ops->view)(ksp,viewer);}
173:   } else if (isdraw) {
174:     PetscDraw draw;
175:     char      str[36];
176:     PetscReal x,y,bottom,h;
177:     PetscBool flg;

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

196:     PetscObjectGetName((PetscObject)ksp,&name);
197:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
198:     if (!((PetscObject)ksp)->amsmem && !rank) {
199:       char       dir[1024];

201:       PetscObjectViewSAWs((PetscObject)ksp,viewer);
202:       PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/its",name);
203:       PetscStackCallSAWs(SAWs_Register,(dir,&ksp->its,1,SAWs_READ,SAWs_INT));
204:       if (!ksp->res_hist) {
205:         KSPSetResidualHistory(ksp,NULL,PETSC_DECIDE,PETSC_TRUE);
206:       }
207:       PetscSNPrintf(dir,1024,"/PETSc/Objects/%s/res_hist",name);
208:       PetscStackCallSAWs(SAWs_Register,(dir,ksp->res_hist,10,SAWs_READ,SAWs_DOUBLE));
209:     }
210: #endif
211:   } else if (ksp->ops->view) {
212:     (*ksp->ops->view)(ksp,viewer);
213:   }
214:   if (ksp->pc) {
215:     PCView(ksp->pc,viewer);
216:   }
217:   if (isdraw) {
218:     PetscDraw draw;
219:     PetscViewerDrawGetDraw(viewer,0,&draw);
220:     PetscDrawPopCurrentPoint(draw);
221:   }
222:   return(0);
223: }

225: /*@C
226:    KSPViewFromOptions - View from Options

228:    Collective on KSP

230:    Input Parameters:
231: +  A - Krylov solver context
232: .  obj - Optional object
233: -  name - command line option

235:    Level: intermediate
236: .seealso:  KSP, KSPView, PetscObjectViewFromOptions(), KSPCreate()
237: @*/
238: PetscErrorCode  KSPViewFromOptions(KSP A,PetscObject obj,const char name[])
239: {

244:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
245:   return(0);
246: }

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

251:    Logically Collective on ksp

253:    Input Parameter:
254: +  ksp - Krylov solver context
255: -  normtype - one of
256: $   KSP_NORM_NONE - skips computing the norm, this should generally only be used if you are using
257: $                 the Krylov method as a smoother with a fixed small number of iterations.
258: $                 Implicitly sets KSPConvergedSkip() as KSP convergence test.
259: $                 Note that certain algorithms such as KSPGMRES ALWAYS require the norm calculation,
260: $                 for these methods the norms are still computed, they are just not used in
261: $                 the convergence test. 
262: $   KSP_NORM_PRECONDITIONED - the default for left preconditioned solves, uses the l2 norm
263: $                 of the preconditioned residual P^{-1}(b - A x)
264: $   KSP_NORM_UNPRECONDITIONED - uses the l2 norm of the true b - Ax residual.
265: $   KSP_NORM_NATURAL - supported  by KSPCG, KSPCR, KSPCGNE, KSPCGS


268:    Options Database Key:
269: .   -ksp_norm_type <none,preconditioned,unpreconditioned,natural>

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

276:    Developer Notes:
277:    Supported combinations of norm and preconditioner side are set using KSPSetSupportedNorm().

279:    Level: advanced

281: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPConvergedSkip(), KSPSetCheckNormIteration(), KSPSetPCSide(), KSPGetPCSide(), KSPNormType
282: @*/
283: PetscErrorCode  KSPSetNormType(KSP ksp,KSPNormType normtype)
284: {
288:   ksp->normtype = ksp->normtype_set = normtype;
289:   return(0);
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.

296:    Logically Collective on ksp

298:    Input Parameter:
299: +  ksp - Krylov solver context
300: -  it  - use -1 to check at all iterations

302:    Notes:
303:    Currently only works with KSPCG, KSPBCGS and KSPIBCGS

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

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

311: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPConvergedSkip(), KSPSetNormType()
312: @*/
313: PetscErrorCode  KSPSetCheckNormIteration(KSP ksp,PetscInt it)
314: {
318:   ksp->chknorm = it;
319:   return(0);
320: }

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


328:    Logically Collective on ksp

330:    Input Parameter:
331: +  ksp - Krylov solver context
332: -  flg - PETSC_TRUE or PETSC_FALSE

334:    Options Database Keys:
335: .  -ksp_lag_norm - lag the calculated residual norm

337:    Notes:
338:    Currently only works with KSPIBCGS.

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

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

345: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPConvergedSkip(), KSPSetNormType(), KSPSetCheckNormIteration()
346: @*/
347: PetscErrorCode  KSPSetLagNorm(KSP ksp,PetscBool flg)
348: {
352:   ksp->lagnorm = flg;
353:   return(0);
354: }

356: /*@
357:    KSPSetSupportedNorm - Sets a norm and preconditioner side supported by a KSP

359:    Logically Collective

361:    Input Arguments:
362: +  ksp - Krylov method
363: .  normtype - supported norm type
364: .  pcside - preconditioner side that can be used with this norm
365: -  priority - positive integer preference for this combination; larger values have higher priority

367:    Level: developer

369:    Notes:
370:    This function should be called from the implementation files KSPCreate_XXX() to declare
371:    which norms and preconditioner sides are supported. Users should not need to call this
372:    function.

374: .seealso: KSPSetNormType(), KSPSetPCSide()
375: @*/
376: PetscErrorCode KSPSetSupportedNorm(KSP ksp,KSPNormType normtype,PCSide pcside,PetscInt priority)
377: {

381:   ksp->normsupporttable[normtype][pcside] = priority;
382:   return(0);
383: }

385: PetscErrorCode KSPNormSupportTableReset_Private(KSP ksp)
386: {

390:   PetscMemzero(ksp->normsupporttable,sizeof(ksp->normsupporttable));
391:   ksp->pc_side  = ksp->pc_side_set;
392:   ksp->normtype = ksp->normtype_set;
393:   return(0);
394: }

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

401:   best = 0;
402:   for (i=0; i<KSP_NORM_MAX; i++) {
403:     for (j=0; j<PC_SIDE_MAX; j++) {
404:       if ((ksp->normtype == KSP_NORM_DEFAULT || ksp->normtype == i) && (ksp->pc_side == PC_SIDE_DEFAULT || ksp->pc_side == j) && (ksp->normsupporttable[i][j] > best)) {
405:         best  = ksp->normsupporttable[i][j];
406:         ibest = i;
407:         jbest = j;
408:       }
409:     }
410:   }
411:   if (best < 1 && errorifnotsupported) {
412:     if (ksp->normtype == KSP_NORM_DEFAULT && ksp->pc_side == PC_SIDE_DEFAULT) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"The %s KSP implementation did not call KSPSetSupportedNorm()",((PetscObject)ksp)->type_name);
413:     if (ksp->normtype == KSP_NORM_DEFAULT) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSP %s does not support %s",((PetscObject)ksp)->type_name,PCSides[ksp->pc_side]);
414:     if (ksp->pc_side == PC_SIDE_DEFAULT) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSP %s does not support %s",((PetscObject)ksp)->type_name,KSPNormTypes[ksp->normtype]);
415:     SETERRQ3(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSP %s does not support %s with %s",((PetscObject)ksp)->type_name,KSPNormTypes[ksp->normtype],PCSides[ksp->pc_side]);
416:   }
417:   if (normtype) *normtype = (KSPNormType)ibest;
418:   if (pcside)   *pcside   = (PCSide)jbest;
419:   return(0);
420: }

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

425:    Not Collective

427:    Input Parameter:
428: .  ksp - Krylov solver context

430:    Output Parameter:
431: .  normtype - norm that is used for convergence testing

433:    Level: advanced

435: .seealso: KSPNormType, KSPSetNormType(), KSPConvergedSkip()
436: @*/
437: PetscErrorCode  KSPGetNormType(KSP ksp, KSPNormType *normtype)
438: {

444:   KSPSetUpNorms_Private(ksp,PETSC_TRUE,&ksp->normtype,&ksp->pc_side);
445:   *normtype = ksp->normtype;
446:   return(0);
447: }

449: #if defined(PETSC_HAVE_SAWS)
450:  #include <petscviewersaws.h>
451: #endif

453: /*@
454:    KSPSetOperators - Sets the matrix associated with the linear system
455:    and a (possibly) different one associated with the preconditioner.

457:    Collective on ksp

459:    Input Parameters:
460: +  ksp - the KSP context
461: .  Amat - the matrix that defines the linear system
462: -  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.

464:    Notes:

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

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

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

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

477:     Level: beginner

479:    Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
480:       are created in PC and returned to the user. In this case, if both operators
481:       mat and pmat are requested, two DIFFERENT operators will be returned. If
482:       only one is requested both operators in the PC will be the same (i.e. as
483:       if one had called KSP/PCSetOperators() with the same argument for both Mats).
484:       The user must set the sizes of the returned matrices and their type etc just
485:       as if the user created them with MatCreate(). For example,

487: $         KSP/PCGetOperators(ksp/pc,&mat,NULL); is equivalent to
488: $           set size, type, etc of mat

490: $         MatCreate(comm,&mat);
491: $         KSP/PCSetOperators(ksp/pc,mat,mat);
492: $         PetscObjectDereference((PetscObject)mat);
493: $           set size, type, etc of mat

495:      and

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

500: $         MatCreate(comm,&mat);
501: $         MatCreate(comm,&pmat);
502: $         KSP/PCSetOperators(ksp/pc,mat,pmat);
503: $         PetscObjectDereference((PetscObject)mat);
504: $         PetscObjectDereference((PetscObject)pmat);
505: $           set size, type, etc of mat and pmat

507:     The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
508:     of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
509:     managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
510:     at this is when you create a SNES you do not NEED to create a KSP and attach it to
511:     the SNES object (the SNES object manages it for you). Similarly when you create a KSP
512:     you do not need to attach a PC to it (the KSP object manages the PC object for you).
513:     Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
514:     it can be created for you?

516: .seealso: KSPSolve(), KSPGetPC(), PCGetOperators(), PCSetOperators(), KSPGetOperators(), KSPSetComputeOperators(), KSPSetComputeInitialGuess(), KSPSetComputeRHS()
517: @*/
518: PetscErrorCode  KSPSetOperators(KSP ksp,Mat Amat,Mat Pmat)
519: {

528:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
529:   PCSetOperators(ksp->pc,Amat,Pmat);
530:   if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX;  /* so that next solve call will call PCSetUp() on new matrix */
531:   return(0);
532: }

534: /*@
535:    KSPGetOperators - Gets the matrix associated with the linear system
536:    and a (possibly) different one associated with the preconditioner.

538:    Collective on ksp

540:    Input Parameter:
541: .  ksp - the KSP context

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

547:     Level: intermediate

549:    Notes:
550:     DOES NOT increase the reference counts of the matrix, so you should NOT destroy them.

552: .seealso: KSPSolve(), KSPGetPC(), PCGetOperators(), PCSetOperators(), KSPSetOperators(), KSPGetOperatorsSet()
553: @*/
554: PetscErrorCode  KSPGetOperators(KSP ksp,Mat *Amat,Mat *Pmat)
555: {

560:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
561:   PCGetOperators(ksp->pc,Amat,Pmat);
562:   return(0);
563: }

565: /*@C
566:    KSPGetOperatorsSet - Determines if the matrix associated with the linear system and
567:    possibly a different one associated with the preconditioner have been set in the KSP.

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

571:    Input Parameter:
572: .  pc - the KSP context

574:    Output Parameters:
575: +  mat - the matrix associated with the linear system was set
576: -  pmat - matrix associated with the preconditioner was set, usually the same

578:    Level: intermediate

580: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators(), PCGetOperatorsSet()
581: @*/
582: PetscErrorCode  KSPGetOperatorsSet(KSP ksp,PetscBool  *mat,PetscBool  *pmat)
583: {

588:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
589:   PCGetOperatorsSet(ksp->pc,mat,pmat);
590:   return(0);
591: }

593: /*@C
594:    KSPSetPreSolve - Sets a function that is called before every KSPSolve() is started

596:    Logically Collective on ksp

598:    Input Parameters:
599: +   ksp - the solver object
600: .   presolve - the function to call before the solve
601: -   prectx - any context needed by the function

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

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

611:    Level: developer

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

624: /*@C
625:    KSPSetPostSolve - Sets a function that is called after every KSPSolve() completes (whether it converges or not)

627:    Logically Collective on ksp

629:    Input Parameters:
630: +   ksp - the solver object
631: .   postsolve - the function to call after the solve
632: -   postctx - any context needed by the function

634:    Level: developer

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

639: +  ksp - the KSP context
640: .  rhs - the right-hand side vector
641: .  x - the solution vector
642: -  ctx - optional user-provided context

644: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSP, KSPSetPreSolve()
645: @*/
646: PetscErrorCode  KSPSetPostSolve(KSP ksp,PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*),void *postctx)
647: {
650:   ksp->postsolve = postsolve;
651:   ksp->postctx   = postctx;
652:   return(0);
653: }

655: /*@
656:    KSPCreate - Creates the default KSP context.

658:    Collective

660:    Input Parameter:
661: .  comm - MPI communicator

663:    Output Parameter:
664: .  ksp - location to put the KSP context

666:    Notes:
667:    The default KSP type is GMRES with a restart of 30, using modified Gram-Schmidt
668:    orthogonalization.

670:    Level: beginner

672: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSP
673: @*/
674: PetscErrorCode  KSPCreate(MPI_Comm comm,KSP *inksp)
675: {
676:   KSP            ksp;
678:   void           *ctx;

682:   *inksp = NULL;
683:   KSPInitializePackage();

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

687:   ksp->max_it  = 10000;
688:   ksp->pc_side = ksp->pc_side_set = PC_SIDE_DEFAULT;
689:   ksp->rtol    = 1.e-5;
690: #if defined(PETSC_USE_REAL_SINGLE)
691:   ksp->abstol  = 1.e-25;
692: #else
693:   ksp->abstol  = 1.e-50;
694: #endif
695:   ksp->divtol  = 1.e4;

697:   ksp->chknorm        = -1;
698:   ksp->normtype       = ksp->normtype_set = KSP_NORM_DEFAULT;
699:   ksp->rnorm          = 0.0;
700:   ksp->its            = 0;
701:   ksp->guess_zero     = PETSC_TRUE;
702:   ksp->calc_sings     = PETSC_FALSE;
703:   ksp->res_hist       = NULL;
704:   ksp->res_hist_alloc = NULL;
705:   ksp->res_hist_len   = 0;
706:   ksp->res_hist_max   = 0;
707:   ksp->res_hist_reset = PETSC_TRUE;
708:   ksp->numbermonitors = 0;
709:   ksp->setfromoptionscalled = 0;

711:   KSPConvergedDefaultCreate(&ctx);
712:   KSPSetConvergenceTest(ksp,KSPConvergedDefault,ctx,KSPConvergedDefaultDestroy);
713:   ksp->ops->buildsolution = KSPBuildSolutionDefault;
714:   ksp->ops->buildresidual = KSPBuildResidualDefault;

716:   ksp->vec_sol    = NULL;
717:   ksp->vec_rhs    = NULL;
718:   ksp->pc         = NULL;
719:   ksp->data       = NULL;
720:   ksp->nwork      = 0;
721:   ksp->work       = NULL;
722:   ksp->reason     = KSP_CONVERGED_ITERATING;
723:   ksp->setupstage = KSP_SETUP_NEW;

725:   KSPNormSupportTableReset_Private(ksp);

727:   *inksp = ksp;
728:   return(0);
729: }

731: /*@C
732:    KSPSetType - Builds KSP for a particular solver.

734:    Logically Collective on ksp

736:    Input Parameters:
737: +  ksp      - the Krylov space context
738: -  type - a known method

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

744:    Notes:
745:    See "petsc/include/petscksp.h" for available methods (for instance,
746:    KSPCG or KSPGMRES).

748:   Normally, it is best to use the KSPSetFromOptions() command and
749:   then set the KSP type from the options database rather than by using
750:   this routine.  Using the options database provides the user with
751:   maximum flexibility in evaluating the many different Krylov methods.
752:   The KSPSetType() routine is provided for those situations where it
753:   is necessary to set the iterative solver independently of the command
754:   line or options database.  This might be the case, for example, when
755:   the choice of iterative solver changes during the execution of the
756:   program, and the user's application is taking responsibility for
757:   choosing the appropriate method.  In other words, this routine is
758:   not for beginners.

760:   Level: intermediate

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

765: .seealso: PCSetType(), KSPType, KSPRegister(), KSPCreate()

767: @*/
768: PetscErrorCode  KSPSetType(KSP ksp, KSPType type)
769: {
770:   PetscErrorCode ierr,(*r)(KSP);
771:   PetscBool      match;


777:   PetscObjectTypeCompare((PetscObject)ksp,type,&match);
778:   if (match) return(0);

780:   PetscFunctionListFind(KSPList,type,&r);
781:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested KSP type %s",type);
782:   /* Destroy the previous private KSP context */
783:   if (ksp->ops->destroy) {
784:     (*ksp->ops->destroy)(ksp);
785:     ksp->ops->destroy = NULL;
786:   }
787:   /* Reinitialize function pointers in KSPOps structure */
788:   PetscMemzero(ksp->ops,sizeof(struct _KSPOps));
789:   ksp->ops->buildsolution = KSPBuildSolutionDefault;
790:   ksp->ops->buildresidual = KSPBuildResidualDefault;
791:   KSPNormSupportTableReset_Private(ksp);
792:   ksp->setupnewmatrix     = PETSC_FALSE; // restore default (setup not called in case of new matrix)
793:   /* Call the KSPCreate_XXX routine for this particular Krylov solver */
794:   ksp->setupstage = KSP_SETUP_NEW;
795:   (*r)(ksp);
796:   PetscObjectChangeTypeName((PetscObject)ksp,type);
797:   return(0);
798: }

800: /*@C
801:    KSPGetType - Gets the KSP type as a string from the KSP object.

803:    Not Collective

805:    Input Parameter:
806: .  ksp - Krylov context

808:    Output Parameter:
809: .  name - name of KSP method

811:    Level: intermediate

813: .seealso: KSPSetType()
814: @*/
815: PetscErrorCode  KSPGetType(KSP ksp,KSPType *type)
816: {
820:   *type = ((PetscObject)ksp)->type_name;
821:   return(0);
822: }

824: /*@C
825:   KSPRegister -  Adds a method to the Krylov subspace solver package.

827:    Not Collective

829:    Input Parameters:
830: +  name_solver - name of a new user-defined solver
831: -  routine_create - routine to create method context

833:    Notes:
834:    KSPRegister() may be called multiple times to add several user-defined solvers.

836:    Sample usage:
837: .vb
838:    KSPRegister("my_solver",MySolverCreate);
839: .ve

841:    Then, your solver can be chosen with the procedural interface via
842: $     KSPSetType(ksp,"my_solver")
843:    or at runtime via the option
844: $     -ksp_type my_solver

846:    Level: advanced

848: .seealso: KSPRegisterAll()
849: @*/
850: PetscErrorCode  KSPRegister(const char sname[],PetscErrorCode (*function)(KSP))
851: {

855:   KSPInitializePackage();
856:   PetscFunctionListAdd(&KSPList,sname,function);
857:   return(0);
858: }