Actual source code: itfunc.c

  1: /*
  2:       Interface KSP routines that the user calls.
  3: */

  5: #include <petsc/private/kspimpl.h>
  6: #include <petsc/private/matimpl.h>
  7: #include <petscdm.h>

  9: /* number of nested levels of KSPSetUp/Solve(). This is used to determine if KSP_DIVERGED_ITS should be fatal. */
 10: static PetscInt level = 0;

 12: static inline PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format)
 13: {

 15:   PetscViewerPushFormat(viewer, format);
 16:   PetscObjectView(obj, viewer);
 17:   PetscViewerPopFormat(viewer);
 18:   return(0);
 19: }

 21: /*@
 22:    KSPComputeExtremeSingularValues - Computes the extreme singular values
 23:    for the preconditioned operator. Called after or during KSPSolve().

 25:    Not Collective

 27:    Input Parameter:
 28: .  ksp - iterative context obtained from KSPCreate()

 30:    Output Parameters:
 31: .  emin, emax - extreme singular values

 33:    Options Database Keys:
 34: .  -ksp_view_singularvalues - compute extreme singular values and print when KSPSolve() completes.

 36:    Notes:
 37:    One must call KSPSetComputeSingularValues() before calling KSPSetUp()
 38:    (or use the option -ksp_view_eigenvalues) in order for this routine to work correctly.

 40:    Many users may just want to use the monitoring routine
 41:    KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
 42:    to print the extreme singular values at each iteration of the linear solve.

 44:    Estimates of the smallest singular value may be very inaccurate, especially if the Krylov method has not converged.
 45:    The largest singular value is usually accurate to within a few percent if the method has converged, but is still not
 46:    intended for eigenanalysis.

 48:    Disable restarts if using KSPGMRES, otherwise this estimate will only be using those iterations after the last
 49:    restart. See KSPGMRESSetRestart() for more details.

 51:    Level: advanced

 53: .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeEigenvalues(), KSP
 54: @*/
 55: PetscErrorCode  KSPComputeExtremeSingularValues(KSP ksp,PetscReal *emax,PetscReal *emin)
 56: {

 62:   if (ksp->ops->computeextremesingularvalues) {
 63:     (*ksp->ops->computeextremesingularvalues)(ksp,emax,emin);
 64:   } else {
 65:     *emin = -1.0;
 66:     *emax = -1.0;
 67:   }
 68:   return 0;
 69: }

 71: /*@
 72:    KSPComputeEigenvalues - Computes the extreme eigenvalues for the
 73:    preconditioned operator. Called after or during KSPSolve().

 75:    Not Collective

 77:    Input Parameters:
 78: +  ksp - iterative context obtained from KSPCreate()
 79: -  n - size of arrays r and c. The number of eigenvalues computed (neig) will, in
 80:        general, be less than this.

 82:    Output Parameters:
 83: +  r - real part of computed eigenvalues, provided by user with a dimension of at least n
 84: .  c - complex part of computed eigenvalues, provided by user with a dimension of at least n
 85: -  neig - actual number of eigenvalues computed (will be less than or equal to n)

 87:    Options Database Keys:
 88: .  -ksp_view_eigenvalues - Prints eigenvalues to stdout

 90:    Notes:
 91:    The number of eigenvalues estimated depends on the size of the Krylov space
 92:    generated during the KSPSolve() ; for example, with
 93:    CG it corresponds to the number of CG iterations, for GMRES it is the number
 94:    of GMRES iterations SINCE the last restart. Any extra space in r[] and c[]
 95:    will be ignored.

 97:    KSPComputeEigenvalues() does not usually provide accurate estimates; it is
 98:    intended only for assistance in understanding the convergence of iterative
 99:    methods, not for eigenanalysis. For accurate computation of eigenvalues we recommend using
100:    the excellent package SLEPc.

102:    One must call KSPSetComputeEigenvalues() before calling KSPSetUp()
103:    in order for this routine to work correctly.

105:    Many users may just want to use the monitoring routine
106:    KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
107:    to print the singular values at each iteration of the linear solve.

109:    Level: advanced

111: .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), KSP
112: @*/
113: PetscErrorCode  KSPComputeEigenvalues(KSP ksp,PetscInt n,PetscReal r[],PetscReal c[],PetscInt *neig)
114: {

122:   if (n && ksp->ops->computeeigenvalues) {
123:     (*ksp->ops->computeeigenvalues)(ksp,n,r,c,neig);
124:   } else {
125:     *neig = 0;
126:   }
127:   return 0;
128: }

130: /*@
131:    KSPComputeRitz - Computes the Ritz or harmonic Ritz pairs associated to the
132:    smallest or largest in modulus, for the preconditioned operator.
133:    Called after KSPSolve().

135:    Not Collective

137:    Input Parameters:
138: +  ksp   - iterative context obtained from KSPCreate()
139: .  ritz  - PETSC_TRUE or PETSC_FALSE for Ritz pairs or harmonic Ritz pairs, respectively
140: -  small - PETSC_TRUE or PETSC_FALSE for smallest or largest (harmonic) Ritz values, respectively

142:    Output Parameters:
143: +  nrit  - On input number of (harmonic) Ritz pairs to compute; on output, actual number of computed (harmonic) Ritz pairs
144: .  S     - an array of the Ritz vectors
145: .  tetar - real part of the Ritz values
146: -  tetai - imaginary part of the Ritz values

148:    Notes:
149:    -For GMRES, the (harmonic) Ritz pairs are computed from the Hessenberg matrix obtained during
150:    the last complete cycle, or obtained at the end of the solution if the method is stopped before
151:    a restart. Then, the number of actual (harmonic) Ritz pairs computed is less or equal to the restart
152:    parameter for GMRES if a complete cycle has been performed or less or equal to the number of GMRES
153:    iterations.
154:    -Moreover, for real matrices, the (harmonic) Ritz pairs are possibly complex-valued. In such a case,
155:    the routine selects the complex (harmonic) Ritz value and its conjugate, and two successive entries of S
156:    are equal to the real and the imaginary parts of the associated vectors.
157:    -the (harmonic) Ritz pairs are given in order of increasing (harmonic) Ritz values in modulus
158:    -this is currently not implemented when PETSc is built with complex numbers

160:    One must call KSPSetComputeRitz() before calling KSPSetUp()
161:    in order for this routine to work correctly.

163:    Level: advanced

165: .seealso: KSPSetComputeRitz(), KSP
166: @*/
167: PetscErrorCode  KSPComputeRitz(KSP ksp,PetscBool ritz,PetscBool small,PetscInt *nrit,Vec S[],PetscReal tetar[],PetscReal tetai[])
168: {
171:   if (ksp->ops->computeritz) (*ksp->ops->computeritz)(ksp,ritz,small,nrit,S,tetar,tetai);
172:   return 0;
173: }
174: /*@
175:    KSPSetUpOnBlocks - Sets up the preconditioner for each block in
176:    the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
177:    methods.

179:    Collective on ksp

181:    Input Parameter:
182: .  ksp - the KSP context

184:    Notes:
185:    KSPSetUpOnBlocks() is a routine that the user can optinally call for
186:    more precise profiling (via -log_view) of the setup phase for these
187:    block preconditioners.  If the user does not call KSPSetUpOnBlocks(),
188:    it will automatically be called from within KSPSolve().

190:    Calling KSPSetUpOnBlocks() is the same as calling PCSetUpOnBlocks()
191:    on the PC context within the KSP context.

193:    Level: advanced

195: .seealso: PCSetUpOnBlocks(), KSPSetUp(), PCSetUp(), KSP
196: @*/
197: PetscErrorCode  KSPSetUpOnBlocks(KSP ksp)
198: {
199:   PC             pc;
200:   PCFailedReason pcreason;

203:   level++;
204:   KSPGetPC(ksp,&pc);
205:   PCSetUpOnBlocks(pc);
206:   PCGetFailedReasonRank(pc,&pcreason);
207:   level--;
208:   /*
209:      This is tricky since only a subset of MPI ranks may set this; each KSPSolve_*() is responsible for checking
210:      this flag and initializing an appropriate vector with VecSetInf() so that the first norm computation can
211:      produce a result at KSPCheckNorm() thus communicating the known problem to all MPI ranks so they may
212:      terminate the Krylov solve. For many KSP implementations this is handled within KSPInitialResidual()
213:   */
214:   if (pcreason) {
215:     ksp->reason = KSP_DIVERGED_PC_FAILED;
216:   }
217:   return 0;
218: }

220: /*@
221:    KSPSetReusePreconditioner - reuse the current preconditioner, do not construct a new one even if the operator changes

223:    Collective on ksp

225:    Input Parameters:
226: +  ksp   - iterative context obtained from KSPCreate()
227: -  flag - PETSC_TRUE to reuse the current preconditioner

229:    Level: intermediate

231: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), PCSetReusePreconditioner(), KSP
232: @*/
233: PetscErrorCode  KSPSetReusePreconditioner(KSP ksp,PetscBool flag)
234: {
235:   PC             pc;

238:   KSPGetPC(ksp,&pc);
239:   PCSetReusePreconditioner(pc,flag);
240:   return 0;
241: }

243: /*@
244:    KSPGetReusePreconditioner - Determines if the KSP reuses the current preconditioner even if the operator in the preconditioner has changed.

246:    Collective on ksp

248:    Input Parameters:
249: .  ksp   - iterative context obtained from KSPCreate()

251:    Output Parameters:
252: .  flag - the boolean flag

254:    Level: intermediate

256: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), KSPSetReusePreconditioner(), KSP
257: @*/
258: PetscErrorCode  KSPGetReusePreconditioner(KSP ksp,PetscBool *flag)
259: {
262:   *flag = PETSC_FALSE;
263:   if (ksp->pc) {
264:     PCGetReusePreconditioner(ksp->pc,flag);
265:   }
266:   return 0;
267: }

269: /*@
270:    KSPSetSkipPCSetFromOptions - prevents KSPSetFromOptions() from call PCSetFromOptions(). This is used if the same PC is shared by more than one KSP so its options are not resetable for each KSP

272:    Collective on ksp

274:    Input Parameters:
275: +  ksp   - iterative context obtained from KSPCreate()
276: -  flag - PETSC_TRUE to skip calling the PCSetFromOptions()

278:    Level: intermediate

280: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), PCSetReusePreconditioner(), KSP
281: @*/
282: PetscErrorCode  KSPSetSkipPCSetFromOptions(KSP ksp,PetscBool flag)
283: {
285:   ksp->skippcsetfromoptions = flag;
286:   return 0;
287: }

289: /*@
290:    KSPSetUp - Sets up the internal data structures for the
291:    later use of an iterative solver.

293:    Collective on ksp

295:    Input Parameter:
296: .  ksp   - iterative context obtained from KSPCreate()

298:    Level: developer

300: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), KSP
301: @*/
302: PetscErrorCode KSPSetUp(KSP ksp)
303: {
304:   Mat            A,B;
305:   Mat            mat,pmat;
306:   MatNullSpace   nullsp;
307:   PCFailedReason pcreason;

310:   level++;

312:   /* reset the convergence flag from the previous solves */
313:   ksp->reason = KSP_CONVERGED_ITERATING;

315:   if (!((PetscObject)ksp)->type_name) {
316:     KSPSetType(ksp,KSPGMRES);
317:   }
318:   KSPSetUpNorms_Private(ksp,PETSC_TRUE,&ksp->normtype,&ksp->pc_side);

320:   if (ksp->dmActive && !ksp->setupstage) {
321:     /* first time in so build matrix and vector data structures using DM */
322:     if (!ksp->vec_rhs) DMCreateGlobalVector(ksp->dm,&ksp->vec_rhs);
323:     if (!ksp->vec_sol) DMCreateGlobalVector(ksp->dm,&ksp->vec_sol);
324:     DMCreateMatrix(ksp->dm,&A);
325:     KSPSetOperators(ksp,A,A);
326:     PetscObjectDereference((PetscObject)A);
327:   }

329:   if (ksp->dmActive) {
330:     DMKSP kdm;
331:     DMGetDMKSP(ksp->dm,&kdm);

333:     if (kdm->ops->computeinitialguess && ksp->setupstage != KSP_SETUP_NEWRHS) {
334:       /* only computes initial guess the first time through */
335:       (*kdm->ops->computeinitialguess)(ksp,ksp->vec_sol,kdm->initialguessctx);
336:       KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
337:     }
338:     if (kdm->ops->computerhs) {
339:       (*kdm->ops->computerhs)(ksp,ksp->vec_rhs,kdm->rhsctx);
340:     }

342:     if (ksp->setupstage != KSP_SETUP_NEWRHS) {
343:       if (kdm->ops->computeoperators) {
344:         KSPGetOperators(ksp,&A,&B);
345:         (*kdm->ops->computeoperators)(ksp,A,B,kdm->operatorsctx);
346:       } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You called KSPSetDM() but did not use DMKSPSetComputeOperators() or KSPSetDMActive(ksp,PETSC_FALSE);");
347:     }
348:   }

350:   if (ksp->setupstage == KSP_SETUP_NEWRHS) {
351:     level--;
352:     return 0;
353:   }
354:   PetscLogEventBegin(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);

356:   switch (ksp->setupstage) {
357:   case KSP_SETUP_NEW:
358:     (*ksp->ops->setup)(ksp);
359:     break;
360:   case KSP_SETUP_NEWMATRIX: {   /* This should be replaced with a more general mechanism */
361:     if (ksp->setupnewmatrix) {
362:       (*ksp->ops->setup)(ksp);
363:     }
364:   } break;
365:   default: break;
366:   }

368:   if (!ksp->pc) KSPGetPC(ksp,&ksp->pc);
369:   PCGetOperators(ksp->pc,&mat,&pmat);
370:   /* scale the matrix if requested */
371:   if (ksp->dscale) {
372:     PetscScalar *xx;
373:     PetscInt    i,n;
374:     PetscBool   zeroflag = PETSC_FALSE;
375:     if (!ksp->pc) KSPGetPC(ksp,&ksp->pc);
376:     if (!ksp->diagonal) { /* allocate vector to hold diagonal */
377:       MatCreateVecs(pmat,&ksp->diagonal,NULL);
378:     }
379:     MatGetDiagonal(pmat,ksp->diagonal);
380:     VecGetLocalSize(ksp->diagonal,&n);
381:     VecGetArray(ksp->diagonal,&xx);
382:     for (i=0; i<n; i++) {
383:       if (xx[i] != 0.0) xx[i] = 1.0/PetscSqrtReal(PetscAbsScalar(xx[i]));
384:       else {
385:         xx[i]    = 1.0;
386:         zeroflag = PETSC_TRUE;
387:       }
388:     }
389:     VecRestoreArray(ksp->diagonal,&xx);
390:     if (zeroflag) {
391:       PetscInfo(ksp,"Zero detected in diagonal of matrix, using 1 at those locations\n");
392:     }
393:     MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
394:     if (mat != pmat) MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);
395:     ksp->dscalefix2 = PETSC_FALSE;
396:   }
397:   PetscLogEventEnd(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);
398:   PCSetErrorIfFailure(ksp->pc,ksp->errorifnotconverged);
399:   PCSetUp(ksp->pc);
400:   PCGetFailedReasonRank(ksp->pc,&pcreason);
401:   /* TODO: this code was wrong and is still wrong, there is no way to propagate the failure to all processes; their is no code to handle a ksp->reason on only some ranks */
402:   if (pcreason) {
403:     ksp->reason = KSP_DIVERGED_PC_FAILED;
404:   }

406:   MatGetNullSpace(mat,&nullsp);
407:   if (nullsp) {
408:     PetscBool test = PETSC_FALSE;
409:     PetscOptionsGetBool(((PetscObject)ksp)->options,((PetscObject)ksp)->prefix,"-ksp_test_null_space",&test,NULL);
410:     if (test) {
411:       MatNullSpaceTest(nullsp,mat,NULL);
412:     }
413:   }
414:   ksp->setupstage = KSP_SETUP_NEWRHS;
415:   level--;
416:   return 0;
417: }

419: /*@C
420:    KSPConvergedReasonView - Displays the reason a KSP solve converged or diverged to a viewer

422:    Collective on ksp

424:    Parameter:
425: +  ksp - iterative context obtained from KSPCreate()
426: -  viewer - the viewer to display the reason

428:    Options Database Keys:
429: +  -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations
430: -  -ksp_converged_reason ::failed - only print reason and number of iterations when diverged

432:    Notes:
433:      To change the format of the output call PetscViewerPushFormat(viewer,format) before this call. Use PETSC_VIEWER_DEFAULT for the default,
434:      use PETSC_VIEWER_FAILED to only display a reason if it fails.

436:    Level: beginner

438: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
439:           KSPSolveTranspose(), KSPGetIterationNumber(), KSP, KSPGetConvergedReason(), PetscViewerPushFormat(), PetscViewerPopFormat()
440: @*/
441: PetscErrorCode KSPConvergedReasonView(KSP ksp, PetscViewer viewer)
442: {
443:   PetscBool         isAscii;
444:   PetscViewerFormat format;

446:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));
447:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
448:   if (isAscii) {
449:     PetscViewerGetFormat(viewer, &format);
450:     PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
451:     if (ksp->reason > 0 && format != PETSC_VIEWER_FAILED) {
452:       if (((PetscObject) ksp)->prefix) {
453:         PetscViewerASCIIPrintf(viewer,"Linear %s solve converged due to %s iterations %D\n",((PetscObject) ksp)->prefix,KSPConvergedReasons[ksp->reason],ksp->its);
454:       } else {
455:         PetscViewerASCIIPrintf(viewer,"Linear solve converged due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
456:       }
457:     } else if (ksp->reason <= 0) {
458:       if (((PetscObject) ksp)->prefix) {
459:         PetscViewerASCIIPrintf(viewer,"Linear %s solve did not converge due to %s iterations %D\n",((PetscObject) ksp)->prefix,KSPConvergedReasons[ksp->reason],ksp->its);
460:       } else {
461:         PetscViewerASCIIPrintf(viewer,"Linear solve did not converge due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
462:       }
463:       if (ksp->reason == KSP_DIVERGED_PC_FAILED) {
464:         PCFailedReason reason;
465:         PCGetFailedReason(ksp->pc,&reason);
466:         PetscViewerASCIIPrintf(viewer,"               PC failed due to %s \n",PCFailedReasons[reason]);
467:       }
468:     }
469:     PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
470:   }
471:   return 0;
472: }

474: /*@C
475:    KSPConvergedReasonViewSet - Sets an ADDITIONAL function that is to be used at the
476:     end of the linear solver to display the convergence reason of the linear solver.

478:    Logically Collective on KSP

480:    Input Parameters:
481: +  ksp - the KSP context
482: .  f - the ksp converged reason view function
483: .  vctx - [optional] user-defined context for private data for the
484:           ksp converged reason view routine (use NULL if no context is desired)
485: -  reasonviewdestroy - [optional] routine that frees reasonview context
486:           (may be NULL)

488:    Options Database Keys:
489: +    -ksp_converged_reason        - sets a default KSPConvergedReasonView()
490: -    -ksp_converged_reason_view_cancel - cancels all converged reason viewers that have
491:                             been hardwired into a code by
492:                             calls to KSPConvergedReasonViewSet(), but
493:                             does not cancel those set via
494:                             the options database.

496:    Notes:
497:    Several different converged reason view routines may be set by calling
498:    KSPConvergedReasonViewSet() multiple times; all will be called in the
499:    order in which they were set.

501:    Level: intermediate

503: .seealso: KSPConvergedReasonView(), KSPConvergedReasonViewCancel()
504: @*/
505: PetscErrorCode  KSPConvergedReasonViewSet(KSP ksp,PetscErrorCode (*f)(KSP,void*),void *vctx,PetscErrorCode (*reasonviewdestroy)(void**))
506: {
507:   PetscInt       i;
508:   PetscBool      identical;

511:   for (i=0; i<ksp->numberreasonviews;i++) {
512:     PetscMonitorCompare((PetscErrorCode (*)(void))f,vctx,reasonviewdestroy,(PetscErrorCode (*)(void))ksp->reasonview[i],ksp->reasonviewcontext[i],ksp->reasonviewdestroy[i],&identical);
513:     if (identical) return 0;
514:   }
516:   ksp->reasonview[ksp->numberreasonviews]          = f;
517:   ksp->reasonviewdestroy[ksp->numberreasonviews]   = reasonviewdestroy;
518:   ksp->reasonviewcontext[ksp->numberreasonviews++] = (void*)vctx;
519:   return 0;
520: }

522: /*@
523:    KSPConvergedReasonViewCancel - Clears all the reasonview functions for a KSP object.

525:    Collective on KSP

527:    Input Parameter:
528: .  ksp - iterative context obtained from KSPCreate()

530:    Level: intermediate

532: .seealso: KSPCreate(), KSPDestroy(), KSPReset()
533: @*/
534: PetscErrorCode  KSPConvergedReasonViewCancel(KSP ksp)
535: {
536:   PetscInt       i;

539:   for (i=0; i<ksp->numberreasonviews; i++) {
540:     if (ksp->reasonviewdestroy[i]) {
541:       (*ksp->reasonviewdestroy[i])(&ksp->reasonviewcontext[i]);
542:     }
543:   }
544:   ksp->numberreasonviews = 0;
545:   return 0;
546: }

548: /*@
549:   KSPConvergedReasonViewFromOptions - Processes command line options to determine if/how a KSPReason is to be viewed.

551:   Collective on ksp

553:   Input Parameters:
554: . ksp   - the KSP object

556:   Level: intermediate

558: .seealso: KSPConvergedReasonView()
559: @*/
560: PetscErrorCode KSPConvergedReasonViewFromOptions(KSP ksp)
561: {
562:   PetscViewer       viewer;
563:   PetscBool         flg;
564:   PetscViewerFormat format;
565:   PetscInt          i;


568:   /* Call all user-provided reason review routines */
569:   for (i=0; i<ksp->numberreasonviews; i++) {
570:     (*ksp->reasonview[i])(ksp,ksp->reasonviewcontext[i]);
571:   }

573:   /* Call the default PETSc routine */
574:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->options,((PetscObject)ksp)->prefix,"-ksp_converged_reason",&viewer,&format,&flg);
575:   if (flg) {
576:     PetscViewerPushFormat(viewer,format);
577:     KSPConvergedReasonView(ksp, viewer);
578:     PetscViewerPopFormat(viewer);
579:     PetscViewerDestroy(&viewer);
580:   }
581:   return 0;
582: }

584: /*@C
585:   KSPConvergedRateView - Displays the reason a KSP solve converged or diverged to a viewer

587:   Collective on ksp

589:   Input Parameters:
590: +  ksp    - iterative context obtained from KSPCreate()
591: -  viewer - the viewer to display the reason

593:   Options Database Keys:
594: . -ksp_converged_rate - print reason for convergence or divergence and the convergence rate (or 0.0 for divergence)

596:   Notes:
597:   To change the format of the output, call PetscViewerPushFormat(viewer,format) before this call.

599:   Suppose that the residual is reduced linearly, $r_k = c^k r_0$, which means $log r_k = log r_0 + k log c$. After linear regression,
600:   the slope is $\log c$. The coefficient of determination is given by $1 - \frac{\sum_i (y_i - f(x_i))^2}{\sum_i (y_i - \bar y)}$,
601:   see also https://en.wikipedia.org/wiki/Coefficient_of_determination

603:   Level: intermediate

605: .seealso: KSPConvergedReasonView(), KSPGetConvergedRate(), KSPSetTolerances(), KSPConvergedDefault()
606: @*/
607: PetscErrorCode KSPConvergedRateView(KSP ksp, PetscViewer viewer)
608: {
609:   PetscViewerFormat format;
610:   PetscBool         isAscii;
611:   PetscReal         rrate, rRsq, erate = 0.0, eRsq = 0.0;
612:   PetscInt          its;
613:   const char       *prefix, *reason = KSPConvergedReasons[ksp->reason];

615:   KSPGetOptionsPrefix(ksp, &prefix);
616:   KSPGetIterationNumber(ksp, &its);
617:   KSPComputeConvergenceRate(ksp, &rrate, &rRsq, &erate, &eRsq);
618:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));
619:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
620:   if (isAscii) {
621:     PetscViewerGetFormat(viewer, &format);
622:     PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
623:     if (ksp->reason > 0) {
624:       if (prefix) PetscViewerASCIIPrintf(viewer, "Linear %s solve converged due to %s iterations %D", prefix, reason, its);
625:       else        PetscViewerASCIIPrintf(viewer, "Linear solve converged due to %s iterations %D", reason, its);
626:       PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
627:       if (rRsq >= 0.0) PetscViewerASCIIPrintf(viewer, " res rate %g R^2 %g", rrate, rRsq);
628:       if (eRsq >= 0.0) PetscViewerASCIIPrintf(viewer, " error rate %g R^2 %g", erate, eRsq);
629:       PetscViewerASCIIPrintf(viewer, "\n");
630:       PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
631:     } else if (ksp->reason <= 0) {
632:       if (prefix) PetscViewerASCIIPrintf(viewer, "Linear %s solve did not converge due to %s iterations %D", prefix, reason, its);
633:       else        PetscViewerASCIIPrintf(viewer, "Linear solve did not converge due to %s iterations %D", reason, its);
634:       PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
635:       if (rRsq >= 0.0) PetscViewerASCIIPrintf(viewer, " res rate %g R^2 %g", rrate, rRsq);
636:       if (eRsq >= 0.0) PetscViewerASCIIPrintf(viewer, " error rate %g R^2 %g", erate, eRsq);
637:       PetscViewerASCIIPrintf(viewer, "\n");
638:       PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
639:       if (ksp->reason == KSP_DIVERGED_PC_FAILED) {
640:         PCFailedReason reason;
641:         PCGetFailedReason(ksp->pc,&reason);
642:         PetscViewerASCIIPrintf(viewer,"               PC failed due to %s \n",PCFailedReasons[reason]);
643:       }
644:     }
645:     PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
646:   }
647:   return 0;
648: }

650: #include <petscdraw.h>

652: static PetscErrorCode KSPViewEigenvalues_Internal(KSP ksp, PetscBool isExplicit, PetscViewer viewer, PetscViewerFormat format)
653: {
654:   PetscReal     *r, *c;
655:   PetscInt       n, i, neig;
656:   PetscBool      isascii, isdraw;
657:   PetscMPIInt    rank;

659:   MPI_Comm_rank(PetscObjectComm((PetscObject) ksp), &rank);
660:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
661:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERDRAW,  &isdraw);
662:   if (isExplicit) {
663:     VecGetSize(ksp->vec_sol,&n);
664:     PetscMalloc2(n, &r, n, &c);
665:     KSPComputeEigenvaluesExplicitly(ksp, n, r, c);
666:     neig = n;
667:   } else {
668:     PetscInt nits;

670:     KSPGetIterationNumber(ksp, &nits);
671:     n    = nits+2;
672:     if (!nits) {PetscViewerASCIIPrintf(viewer, "Zero iterations in solver, cannot approximate any eigenvalues\n"));PetscFunctionReturn(0;}
673:     PetscMalloc2(n, &r, n, &c);
674:     KSPComputeEigenvalues(ksp, n, r, c, &neig);
675:   }
676:   if (isascii) {
677:     PetscViewerASCIIPrintf(viewer, "%s computed eigenvalues\n", isExplicit ? "Explicitly" : "Iteratively");
678:     for (i = 0; i < neig; ++i) {
679:       if (c[i] >= 0.0) PetscViewerASCIIPrintf(viewer, "%g + %gi\n", (double) r[i],  (double) c[i]);
680:       else             PetscViewerASCIIPrintf(viewer, "%g - %gi\n", (double) r[i], -(double) c[i]);
681:     }
682:   } else if (isdraw && rank == 0) {
683:     PetscDraw   draw;
684:     PetscDrawSP drawsp;

686:     if (format == PETSC_VIEWER_DRAW_CONTOUR) {
687:       KSPPlotEigenContours_Private(ksp,neig,r,c);
688:     } else {
689:       if (!ksp->eigviewer) PetscViewerDrawOpen(PETSC_COMM_SELF,NULL,isExplicit ? "Explicitly Computed Eigenvalues" : "Iteratively Computed Eigenvalues",PETSC_DECIDE,PETSC_DECIDE,400,400,&ksp->eigviewer);
690:       PetscViewerDrawGetDraw(ksp->eigviewer,0,&draw);
691:       PetscDrawSPCreate(draw,1,&drawsp);
692:       PetscDrawSPReset(drawsp);
693:       for (i = 0; i < neig; ++i) PetscDrawSPAddPoint(drawsp,r+i,c+i);
694:       PetscDrawSPDraw(drawsp,PETSC_TRUE);
695:       PetscDrawSPSave(drawsp);
696:       PetscDrawSPDestroy(&drawsp);
697:     }
698:   }
699:   PetscFree2(r, c);
700:   return 0;
701: }

703: static PetscErrorCode KSPViewSingularvalues_Internal(KSP ksp, PetscViewer viewer, PetscViewerFormat format)
704: {
705:   PetscReal      smax, smin;
706:   PetscInt       nits;
707:   PetscBool      isascii;

709:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
710:   KSPGetIterationNumber(ksp, &nits);
711:   if (!nits) {PetscViewerASCIIPrintf(viewer, "Zero iterations in solver, cannot approximate any singular values\n"));PetscFunctionReturn(0;}
712:   KSPComputeExtremeSingularValues(ksp, &smax, &smin);
713:   if (isascii) PetscViewerASCIIPrintf(viewer, "Iteratively computed extreme singular values: max %g min %g max/min %g\n",(double)smax,(double)smin,(double)(smax/smin));
714:   return 0;
715: }

717: static PetscErrorCode KSPViewFinalResidual_Internal(KSP ksp, PetscViewer viewer, PetscViewerFormat format)
718: {
719:   PetscBool      isascii;

721:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
723:   if (isascii) {
724:     Mat       A;
725:     Vec       t;
726:     PetscReal norm;

728:     PCGetOperators(ksp->pc, &A, NULL);
729:     VecDuplicate(ksp->vec_rhs, &t);
730:     KSP_MatMult(ksp, A, ksp->vec_sol, t);
731:     VecAYPX(t, -1.0, ksp->vec_rhs);
732:     VecNorm(t, NORM_2, &norm);
733:     VecDestroy(&t);
734:     PetscViewerASCIIPrintf(viewer, "KSP final norm of residual %g\n", (double) norm);
735:   }
736:   return 0;
737: }

739: static PetscErrorCode KSPMonitorPauseFinal_Internal(KSP ksp)
740: {
741:   PetscInt       i;

743:   if (!ksp->pauseFinal) return 0;
744:   for (i = 0; i < ksp->numbermonitors; ++i) {
745:     PetscViewerAndFormat *vf = (PetscViewerAndFormat *) ksp->monitorcontext[i];
746:     PetscDraw             draw;
747:     PetscReal             lpause;

749:     if (!vf) continue;
750:     if (vf->lg) {
752:       if (((PetscObject) vf->lg)->classid != PETSC_DRAWLG_CLASSID) continue;
753:       PetscDrawLGGetDraw(vf->lg, &draw);
754:       PetscDrawGetPause(draw, &lpause);
755:       PetscDrawSetPause(draw, -1.0);
756:       PetscDrawPause(draw);
757:       PetscDrawSetPause(draw, lpause);
758:     } else {
759:       PetscBool isdraw;

762:       if (((PetscObject) vf->viewer)->classid != PETSC_VIEWER_CLASSID) continue;
763:       PetscObjectTypeCompare((PetscObject) vf->viewer, PETSCVIEWERDRAW, &isdraw);
764:       if (!isdraw) continue;
765:       PetscViewerDrawGetDraw(vf->viewer, 0, &draw);
766:       PetscDrawGetPause(draw, &lpause);
767:       PetscDrawSetPause(draw, -1.0);
768:       PetscDrawPause(draw);
769:       PetscDrawSetPause(draw, lpause);
770:     }
771:   }
772:   return 0;
773: }

775: static PetscErrorCode KSPSolve_Private(KSP ksp,Vec b,Vec x)
776: {
777:   PetscBool       flg = PETSC_FALSE,inXisinB = PETSC_FALSE,guess_zero;
778:   Mat             mat,pmat;
779:   MPI_Comm        comm;
780:   MatNullSpace    nullsp;
781:   Vec             btmp,vec_rhs = NULL;

783:   level++;
784:   comm = PetscObjectComm((PetscObject)ksp);
785:   if (x && x == b) {
787:     VecDuplicate(b,&x);
788:     inXisinB = PETSC_TRUE;
789:   }
790:   if (b) {
791:     PetscObjectReference((PetscObject)b);
792:     VecDestroy(&ksp->vec_rhs);
793:     ksp->vec_rhs = b;
794:   }
795:   if (x) {
796:     PetscObjectReference((PetscObject)x);
797:     VecDestroy(&ksp->vec_sol);
798:     ksp->vec_sol = x;
799:   }

801:   if (ksp->viewPre) ObjectView((PetscObject) ksp, ksp->viewerPre, ksp->formatPre);

803:   if (ksp->presolve) (*ksp->presolve)(ksp,ksp->vec_rhs,ksp->vec_sol,ksp->prectx);

805:   /* reset the residual history list if requested */
806:   if (ksp->res_hist_reset) ksp->res_hist_len = 0;
807:   if (ksp->err_hist_reset) ksp->err_hist_len = 0;

809:   if (ksp->guess) {
810:     PetscObjectState ostate,state;

812:     KSPGuessSetUp(ksp->guess);
813:     PetscObjectStateGet((PetscObject)ksp->vec_sol,&ostate);
814:     KSPGuessFormGuess(ksp->guess,ksp->vec_rhs,ksp->vec_sol);
815:     PetscObjectStateGet((PetscObject)ksp->vec_sol,&state);
816:     if (state != ostate) {
817:       ksp->guess_zero = PETSC_FALSE;
818:     } else {
819:       PetscInfo(ksp,"Using zero initial guess since the KSPGuess object did not change the vector\n");
820:       ksp->guess_zero = PETSC_TRUE;
821:     }
822:   }

824:   /* KSPSetUp() scales the matrix if needed */
825:   KSPSetUp(ksp);
826:   KSPSetUpOnBlocks(ksp);

828:   VecSetErrorIfLocked(ksp->vec_sol,3);

830:   PetscLogEventBegin(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
831:   PCGetOperators(ksp->pc,&mat,&pmat);
832:   /* diagonal scale RHS if called for */
833:   if (ksp->dscale) {
834:     VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
835:     /* second time in, but matrix was scaled back to original */
836:     if (ksp->dscalefix && ksp->dscalefix2) {
837:       Mat mat,pmat;

839:       PCGetOperators(ksp->pc,&mat,&pmat);
840:       MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
841:       if (mat != pmat) MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);
842:     }

844:     /* scale initial guess */
845:     if (!ksp->guess_zero) {
846:       if (!ksp->truediagonal) {
847:         VecDuplicate(ksp->diagonal,&ksp->truediagonal);
848:         VecCopy(ksp->diagonal,ksp->truediagonal);
849:         VecReciprocal(ksp->truediagonal);
850:       }
851:       VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->truediagonal);
852:     }
853:   }
854:   PCPreSolve(ksp->pc,ksp);

856:   if (ksp->guess_zero) VecSet(ksp->vec_sol,0.0);
857:   if (ksp->guess_knoll) { /* The Knoll trick is independent on the KSPGuess specified */
858:     PCApply(ksp->pc,ksp->vec_rhs,ksp->vec_sol);
859:     KSP_RemoveNullSpace(ksp,ksp->vec_sol);
860:     ksp->guess_zero = PETSC_FALSE;
861:   }

863:   /* can we mark the initial guess as zero for this solve? */
864:   guess_zero = ksp->guess_zero;
865:   if (!ksp->guess_zero) {
866:     PetscReal norm;

868:     VecNormAvailable(ksp->vec_sol,NORM_2,&flg,&norm);
869:     if (flg && !norm) ksp->guess_zero = PETSC_TRUE;
870:   }
871:   if (ksp->transpose_solve) {
872:     MatGetNullSpace(pmat,&nullsp);
873:   } else {
874:     MatGetTransposeNullSpace(pmat,&nullsp);
875:   }
876:   if (nullsp) {
877:     VecDuplicate(ksp->vec_rhs,&btmp);
878:     VecCopy(ksp->vec_rhs,btmp);
879:     MatNullSpaceRemove(nullsp,btmp);
880:     vec_rhs      = ksp->vec_rhs;
881:     ksp->vec_rhs = btmp;
882:   }
883:   VecLockReadPush(ksp->vec_rhs);
884:   (*ksp->ops->solve)(ksp);
885:   KSPMonitorPauseFinal_Internal(ksp);

887:   VecLockReadPop(ksp->vec_rhs);
888:   if (nullsp) {
889:     ksp->vec_rhs = vec_rhs;
890:     VecDestroy(&btmp);
891:   }

893:   ksp->guess_zero = guess_zero;

896:   ksp->totalits += ksp->its;

898:   KSPConvergedReasonViewFromOptions(ksp);

900:   if (ksp->viewRate) {
901:     PetscViewerPushFormat(ksp->viewerRate,ksp->formatRate);
902:     KSPConvergedRateView(ksp, ksp->viewerRate);
903:     PetscViewerPopFormat(ksp->viewerRate);
904:   }
905:   PCPostSolve(ksp->pc,ksp);

907:   /* diagonal scale solution if called for */
908:   if (ksp->dscale) {
909:     VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->diagonal);
910:     /* unscale right hand side and matrix */
911:     if (ksp->dscalefix) {
912:       Mat mat,pmat;

914:       VecReciprocal(ksp->diagonal);
915:       VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
916:       PCGetOperators(ksp->pc,&mat,&pmat);
917:       MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
918:       if (mat != pmat) MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);
919:       VecReciprocal(ksp->diagonal);
920:       ksp->dscalefix2 = PETSC_TRUE;
921:     }
922:   }
923:   PetscLogEventEnd(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
924:   if (ksp->guess) {
925:     KSPGuessUpdate(ksp->guess,ksp->vec_rhs,ksp->vec_sol);
926:   }
927:   if (ksp->postsolve) {
928:     (*ksp->postsolve)(ksp,ksp->vec_rhs,ksp->vec_sol,ksp->postctx);
929:   }

931:   PCGetOperators(ksp->pc,&mat,&pmat);
932:   if (ksp->viewEV)       KSPViewEigenvalues_Internal(ksp, PETSC_FALSE, ksp->viewerEV,    ksp->formatEV);
933:   if (ksp->viewEVExp)    KSPViewEigenvalues_Internal(ksp, PETSC_TRUE,  ksp->viewerEVExp, ksp->formatEVExp);
934:   if (ksp->viewSV)       KSPViewSingularvalues_Internal(ksp, ksp->viewerSV, ksp->formatSV);
935:   if (ksp->viewFinalRes) KSPViewFinalResidual_Internal(ksp, ksp->viewerFinalRes, ksp->formatFinalRes);
936:   if (ksp->viewMat)      ObjectView((PetscObject) mat,           ksp->viewerMat,    ksp->formatMat);
937:   if (ksp->viewPMat)     ObjectView((PetscObject) pmat,          ksp->viewerPMat,   ksp->formatPMat);
938:   if (ksp->viewRhs)      ObjectView((PetscObject) ksp->vec_rhs,  ksp->viewerRhs,    ksp->formatRhs);
939:   if (ksp->viewSol)      ObjectView((PetscObject) ksp->vec_sol,  ksp->viewerSol,    ksp->formatSol);
940:   if (ksp->view)         ObjectView((PetscObject) ksp,           ksp->viewer,       ksp->format);
941:   if (ksp->viewDScale)   ObjectView((PetscObject) ksp->diagonal, ksp->viewerDScale, ksp->formatDScale);
942:   if (ksp->viewMatExp)   {
943:     Mat A, B;

945:     PCGetOperators(ksp->pc, &A, NULL);
946:     if (ksp->transpose_solve) {
947:       Mat AT;

949:       MatCreateTranspose(A, &AT);
950:       MatComputeOperator(AT, MATAIJ, &B);
951:       MatDestroy(&AT);
952:     } else {
953:       MatComputeOperator(A, MATAIJ, &B);
954:     }
955:     ObjectView((PetscObject) B, ksp->viewerMatExp, ksp->formatMatExp);
956:     MatDestroy(&B);
957:   }
958:   if (ksp->viewPOpExp)   {
959:     Mat B;

961:     KSPComputeOperator(ksp, MATAIJ, &B);
962:     ObjectView((PetscObject) B, ksp->viewerPOpExp, ksp->formatPOpExp);
963:     MatDestroy(&B);
964:   }

966:   if (inXisinB) {
967:     VecCopy(x,b);
968:     VecDestroy(&x);
969:   }
970:   PetscObjectSAWsBlock((PetscObject)ksp);
971:   if (ksp->errorifnotconverged && ksp->reason < 0 && ((level == 1) || (ksp->reason != KSP_DIVERGED_ITS))) {
972:     if (ksp->reason == KSP_DIVERGED_PC_FAILED) {
973:       PCFailedReason reason;
974:       PCGetFailedReason(ksp->pc,&reason);
975:       SETERRQ(comm,PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged, reason %s PC failed due to %s",KSPConvergedReasons[ksp->reason],PCFailedReasons[reason]);
976:     } else SETERRQ(comm,PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged, reason %s",KSPConvergedReasons[ksp->reason]);
977:   }
978:   level--;
979:   return 0;
980: }

982: /*@
983:    KSPSolve - Solves linear system.

985:    Collective on ksp

987:    Parameters:
988: +  ksp - iterative context obtained from KSPCreate()
989: .  b - the right hand side vector
990: -  x - the solution (this may be the same vector as b, then b will be overwritten with answer)

992:    Options Database Keys:
993: +  -ksp_view_eigenvalues - compute preconditioned operators eigenvalues
994: .  -ksp_view_eigenvalues_explicit - compute the eigenvalues by forming the dense operator and using LAPACK
995: .  -ksp_view_mat binary - save matrix to the default binary viewer
996: .  -ksp_view_pmat binary - save matrix used to build preconditioner to the default binary viewer
997: .  -ksp_view_rhs binary - save right hand side vector to the default binary viewer
998: .  -ksp_view_solution binary - save computed solution vector to the default binary viewer
999:            (can be read later with src/ksp/tutorials/ex10.c for testing solvers)
1000: .  -ksp_view_mat_explicit - for matrix-free operators, computes the matrix entries and views them
1001: .  -ksp_view_preconditioned_operator_explicit - computes the product of the preconditioner and matrix as an explicit matrix and views it
1002: .  -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations
1003: .  -ksp_view_final_residual - print 2-norm of true linear system residual at the end of the solution process
1004: .  -ksp_error_if_not_converged - stop the program as soon as an error is detected in a KSPSolve()
1005: -  -ksp_view - print the ksp data structure at the end of the system solution

1007:    Notes:

1009:    If one uses KSPSetDM() then x or b need not be passed. Use KSPGetSolution() to access the solution in this case.

1011:    The operator is specified with KSPSetOperators().

1013:    KSPSolve() will normally return without generating an error regardless of whether the linear system was solved or if constructing the preconditioner failed.
1014:    Call KSPGetConvergedReason() to determine if the solver converged or failed and why. The option -ksp_error_if_not_converged or function KSPSetErrorIfNotConverged()
1015:    will cause KSPSolve() to error as soon as an error occurs in the linear solver.  In inner KSPSolves() KSP_DIVERGED_ITS is not treated as an error because when using nested solvers
1016:    it may be fine that inner solvers in the preconditioner do not converge during the solution process.

1018:    The number of iterations can be obtained from KSPGetIterationNumber().

1020:    If you provide a matrix that has a MatSetNullSpace() and MatSetTransposeNullSpace() this will use that information to solve singular systems
1021:    in the least squares sense with a norm minimizing solution.
1022: $
1023: $                   A x = b   where b = b_p + b_t where b_t is not in the range of A (and hence by the fundamental theorem of linear algebra is in the nullspace(A') see MatSetNullSpace()
1024: $
1025: $    KSP first removes b_t producing the linear system  A x = b_p (which has multiple solutions) and solves this to find the ||x|| minimizing solution (and hence
1026: $    it finds the solution x orthogonal to the nullspace(A). The algorithm is simply in each iteration of the Krylov method we remove the nullspace(A) from the search
1027: $    direction thus the solution which is a linear combination of the search directions has no component in the nullspace(A).
1028: $
1029: $    We recommend always using GMRES for such singular systems.
1030: $    If nullspace(A) = nullspace(A') (note symmetric matrices always satisfy this property) then both left and right preconditioning will work
1031: $    If nullspace(A) != nullspace(A') then left preconditioning will work but right preconditioning may not work (or it may).

1033:    Developer Note: The reason we cannot always solve  nullspace(A) != nullspace(A') systems with right preconditioning is because we need to remove at each iteration
1034:        the nullspace(AB) from the search direction. While we know the nullspace(A) the nullspace(AB) equals B^-1 times the nullspace(A) but except for trivial preconditioners
1035:        such as diagonal scaling we cannot apply the inverse of the preconditioner to a vector and thus cannot compute the nullspace(AB).

1037:    If using a direct method (e.g., via the KSP solver
1038:    KSPPREONLY and a preconditioner such as PCLU/PCILU),
1039:    then its=1.  See KSPSetTolerances() and KSPConvergedDefault()
1040:    for more details.

1042:    Understanding Convergence:
1043:    The routines KSPMonitorSet(), KSPComputeEigenvalues(), and
1044:    KSPComputeEigenvaluesExplicitly() provide information on additional
1045:    options to monitor convergence and print eigenvalue information.

1047:    Level: beginner

1049: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
1050:           KSPSolveTranspose(), KSPGetIterationNumber(), MatNullSpaceCreate(), MatSetNullSpace(), MatSetTransposeNullSpace(), KSP,
1051:           KSPConvergedReasonView(), KSPCheckSolve(), KSPSetErrorIfNotConverged()
1052: @*/
1053: PetscErrorCode KSPSolve(KSP ksp,Vec b,Vec x)
1054: {
1058:   ksp->transpose_solve = PETSC_FALSE;
1059:   KSPSolve_Private(ksp,b,x);
1060:   return 0;
1061: }

1063: /*@
1064:    KSPSolveTranspose - Solves the transpose of a linear system.

1066:    Collective on ksp

1068:    Input Parameters:
1069: +  ksp - iterative context obtained from KSPCreate()
1070: .  b - right hand side vector
1071: -  x - solution vector

1073:    Notes:
1074:     For complex numbers this solve the non-Hermitian transpose system.

1076:    Developer Notes:
1077:     We need to implement a KSPSolveHermitianTranspose()

1079:    Level: developer

1081: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
1082:           KSPSolve(), KSP
1083: @*/
1084: PetscErrorCode KSPSolveTranspose(KSP ksp,Vec b,Vec x)
1085: {
1089:   if (ksp->transpose.use_explicittranspose) {
1090:     Mat J,Jpre;
1091:     KSPGetOperators(ksp,&J,&Jpre);
1092:     if (!ksp->transpose.reuse_transpose) {
1093:       MatTranspose(J,MAT_INITIAL_MATRIX,&ksp->transpose.AT);
1094:       if (J != Jpre) {
1095:         MatTranspose(Jpre,MAT_INITIAL_MATRIX,&ksp->transpose.BT);
1096:       }
1097:       ksp->transpose.reuse_transpose = PETSC_TRUE;
1098:     } else {
1099:       MatTranspose(J,MAT_REUSE_MATRIX,&ksp->transpose.AT);
1100:       if (J != Jpre) {
1101:         MatTranspose(Jpre,MAT_REUSE_MATRIX,&ksp->transpose.BT);
1102:       }
1103:     }
1104:     if (J == Jpre && ksp->transpose.BT != ksp->transpose.AT) {
1105:       PetscObjectReference((PetscObject)ksp->transpose.AT);
1106:       ksp->transpose.BT = ksp->transpose.AT;
1107:     }
1108:     KSPSetOperators(ksp,ksp->transpose.AT,ksp->transpose.BT);
1109:   } else {
1110:     ksp->transpose_solve = PETSC_TRUE;
1111:   }
1112:   KSPSolve_Private(ksp,b,x);
1113:   return 0;
1114: }

1116: static PetscErrorCode KSPViewFinalMatResidual_Internal(KSP ksp, Mat B, Mat X, PetscViewer viewer, PetscViewerFormat format, PetscInt shift)
1117: {
1118:   Mat            A, R;
1119:   PetscReal      *norms;
1120:   PetscInt       i, N;
1121:   PetscBool      flg;

1123:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &flg);
1124:   if (flg) {
1125:     PCGetOperators(ksp->pc, &A, NULL);
1126:     MatMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &R);
1127:     MatAYPX(R, -1.0, B, SAME_NONZERO_PATTERN);
1128:     MatGetSize(R, NULL, &N);
1129:     PetscMalloc1(N, &norms);
1130:     MatGetColumnNorms(R, NORM_2, norms);
1131:     MatDestroy(&R);
1132:     for (i = 0; i < N; ++i) {
1133:       PetscViewerASCIIPrintf(viewer, "%s #%D %g\n", i == 0 ? "KSP final norm of residual" : "                          ", shift + i, (double)norms[i]);
1134:     }
1135:     PetscFree(norms);
1136:   }
1137:   return 0;
1138: }

1140: /*@
1141:      KSPMatSolve - Solves a linear system with multiple right-hand sides stored as a MATDENSE. Unlike KSPSolve(), B and X must be different matrices.

1143:    Input Parameters:
1144: +     ksp - iterative context
1145: -     B - block of right-hand sides

1147:    Output Parameter:
1148: .     X - block of solutions

1150:    Notes:
1151:      This is a stripped-down version of KSPSolve(), which only handles -ksp_view, -ksp_converged_reason, and -ksp_view_final_residual.

1153:    Level: intermediate

1155: .seealso:  KSPSolve(), MatMatSolve(), MATDENSE, KSPHPDDM, PCBJACOBI, PCASM
1156: @*/
1157: PetscErrorCode KSPMatSolve(KSP ksp, Mat B, Mat X)
1158: {
1159:   Mat            A, P, vB, vX;
1160:   Vec            cb, cx;
1161:   PetscInt       n1, N1, n2, N2, Bbn = PETSC_DECIDE;
1162:   PetscBool      match;

1170:   MatCheckPreallocated(X, 3);
1171:   if (!X->assembled) {
1172:     MatSetOption(X, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
1173:     MatAssemblyBegin(X, MAT_FINAL_ASSEMBLY);
1174:     MatAssemblyEnd(X, MAT_FINAL_ASSEMBLY);
1175:   }
1177:   KSPGetOperators(ksp, &A, &P);
1178:   MatGetLocalSize(B, NULL, &n2);
1179:   MatGetLocalSize(X, NULL, &n1);
1180:   MatGetSize(B, NULL, &N2);
1181:   MatGetSize(X, NULL, &N1);
1183:   PetscObjectBaseTypeCompareAny((PetscObject)B, &match, MATSEQDENSE, MATMPIDENSE, "");
1185:   PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, "");
1187:   KSPSetUp(ksp);
1188:   KSPSetUpOnBlocks(ksp);
1189:   if (ksp->ops->matsolve) {
1190:     if (ksp->guess_zero) {
1191:       MatZeroEntries(X);
1192:     }
1193:     PetscLogEventBegin(KSP_MatSolve, ksp, B, X, 0);
1194:     KSPGetMatSolveBatchSize(ksp, &Bbn);
1195:     /* by default, do a single solve with all columns */
1196:     if (Bbn == PETSC_DECIDE) Bbn = N2;
1198:     PetscInfo(ksp, "KSP type %s solving using batches of width at most %D\n", ((PetscObject)ksp)->type_name, Bbn);
1199:     /* if -ksp_matsolve_batch_size is greater than the actual number of columns, do a single solve with all columns */
1200:     if (Bbn >= N2) {
1201:       (*ksp->ops->matsolve)(ksp, B, X);
1202:       if (ksp->viewFinalRes) {
1203:         KSPViewFinalMatResidual_Internal(ksp, B, X, ksp->viewerFinalRes, ksp->formatFinalRes, 0);
1204:       }

1206:       KSPConvergedReasonViewFromOptions(ksp);

1208:       if (ksp->viewRate) {
1209:         PetscViewerPushFormat(ksp->viewerRate,PETSC_VIEWER_DEFAULT);
1210:         KSPConvergedRateView(ksp, ksp->viewerRate);
1211:         PetscViewerPopFormat(ksp->viewerRate);
1212:       }
1213:     } else {
1214:       for (n2 = 0; n2 < N2; n2 += Bbn) {
1215:         MatDenseGetSubMatrix(B, n2, PetscMin(n2+Bbn, N2), &vB);
1216:         MatDenseGetSubMatrix(X, n2, PetscMin(n2+Bbn, N2), &vX);
1217:         (*ksp->ops->matsolve)(ksp, vB, vX);
1218:         if (ksp->viewFinalRes) {
1219:           KSPViewFinalMatResidual_Internal(ksp, vB, vX, ksp->viewerFinalRes, ksp->formatFinalRes, n2);
1220:         }

1222:         KSPConvergedReasonViewFromOptions(ksp);

1224:         if (ksp->viewRate) {
1225:           PetscViewerPushFormat(ksp->viewerRate,PETSC_VIEWER_DEFAULT);
1226:           KSPConvergedRateView(ksp, ksp->viewerRate);
1227:           PetscViewerPopFormat(ksp->viewerRate);
1228:         }
1229:         MatDenseRestoreSubMatrix(B, &vB);
1230:         MatDenseRestoreSubMatrix(X, &vX);
1231:       }
1232:     }
1233:     if (ksp->viewMat)  ObjectView((PetscObject) A, ksp->viewerMat, ksp->formatMat);
1234:     if (ksp->viewPMat) ObjectView((PetscObject) P, ksp->viewerPMat,ksp->formatPMat);
1235:     if (ksp->viewRhs)  ObjectView((PetscObject) B, ksp->viewerRhs, ksp->formatRhs);
1236:     if (ksp->viewSol)  ObjectView((PetscObject) X, ksp->viewerSol, ksp->formatSol);
1237:     if (ksp->view) {
1238:       KSPView(ksp, ksp->viewer);
1239:     }
1240:     PetscLogEventEnd(KSP_MatSolve, ksp, B, X, 0);
1241:   } else {
1242:     PetscInfo(ksp, "KSP type %s solving column by column\n", ((PetscObject)ksp)->type_name);
1243:     for (n2 = 0; n2 < N2; ++n2) {
1244:       MatDenseGetColumnVecRead(B, n2, &cb);
1245:       MatDenseGetColumnVecWrite(X, n2, &cx);
1246:       KSPSolve(ksp, cb, cx);
1247:       MatDenseRestoreColumnVecWrite(X, n2, &cx);
1248:       MatDenseRestoreColumnVecRead(B, n2, &cb);
1249:     }
1250:   }
1251:   return 0;
1252: }

1254: /*@
1255:      KSPSetMatSolveBatchSize - Sets the maximum number of columns treated simultaneously in KSPMatSolve().

1257:     Logically collective

1259:    Input Parameters:
1260: +     ksp - iterative context
1261: -     bs - batch size

1263:    Level: advanced

1265: .seealso:  KSPMatSolve(), KSPGetMatSolveBatchSize(), -mat_mumps_icntl_27, -matmatmult_Bbn
1266: @*/
1267: PetscErrorCode KSPSetMatSolveBatchSize(KSP ksp, PetscInt bs)
1268: {
1271:   ksp->nmax = bs;
1272:   return 0;
1273: }

1275: /*@
1276:      KSPGetMatSolveBatchSize - Gets the maximum number of columns treated simultaneously in KSPMatSolve().

1278:    Input Parameter:
1279: .     ksp - iterative context

1281:    Output Parameter:
1282: .     bs - batch size

1284:    Level: advanced

1286: .seealso:  KSPMatSolve(), KSPSetMatSolveBatchSize(), -mat_mumps_icntl_27, -matmatmult_Bbn
1287: @*/
1288: PetscErrorCode KSPGetMatSolveBatchSize(KSP ksp, PetscInt *bs)
1289: {
1292:   *bs = ksp->nmax;
1293:   return 0;
1294: }

1296: /*@
1297:    KSPResetViewers - Resets all the viewers set from the options database during KSPSetFromOptions()

1299:    Collective on ksp

1301:    Input Parameter:
1302: .  ksp - iterative context obtained from KSPCreate()

1304:    Level: beginner

1306: .seealso: KSPCreate(), KSPSetUp(), KSPSolve(), KSPSetFromOptions(), KSP
1307: @*/
1308: PetscErrorCode  KSPResetViewers(KSP ksp)
1309: {
1311:   if (!ksp) return 0;
1312:   PetscViewerDestroy(&ksp->viewer);
1313:   PetscViewerDestroy(&ksp->viewerPre);
1314:   PetscViewerDestroy(&ksp->viewerRate);
1315:   PetscViewerDestroy(&ksp->viewerMat);
1316:   PetscViewerDestroy(&ksp->viewerPMat);
1317:   PetscViewerDestroy(&ksp->viewerRhs);
1318:   PetscViewerDestroy(&ksp->viewerSol);
1319:   PetscViewerDestroy(&ksp->viewerMatExp);
1320:   PetscViewerDestroy(&ksp->viewerEV);
1321:   PetscViewerDestroy(&ksp->viewerSV);
1322:   PetscViewerDestroy(&ksp->viewerEVExp);
1323:   PetscViewerDestroy(&ksp->viewerFinalRes);
1324:   PetscViewerDestroy(&ksp->viewerPOpExp);
1325:   PetscViewerDestroy(&ksp->viewerDScale);
1326:   ksp->view         = PETSC_FALSE;
1327:   ksp->viewPre      = PETSC_FALSE;
1328:   ksp->viewMat      = PETSC_FALSE;
1329:   ksp->viewPMat     = PETSC_FALSE;
1330:   ksp->viewRhs      = PETSC_FALSE;
1331:   ksp->viewSol      = PETSC_FALSE;
1332:   ksp->viewMatExp   = PETSC_FALSE;
1333:   ksp->viewEV       = PETSC_FALSE;
1334:   ksp->viewSV       = PETSC_FALSE;
1335:   ksp->viewEVExp    = PETSC_FALSE;
1336:   ksp->viewFinalRes = PETSC_FALSE;
1337:   ksp->viewPOpExp   = PETSC_FALSE;
1338:   ksp->viewDScale   = PETSC_FALSE;
1339:   return 0;
1340: }

1342: /*@
1343:    KSPReset - Resets a KSP context to the kspsetupcalled = 0 state and removes any allocated Vecs and Mats

1345:    Collective on ksp

1347:    Input Parameter:
1348: .  ksp - iterative context obtained from KSPCreate()

1350:    Level: beginner

1352: .seealso: KSPCreate(), KSPSetUp(), KSPSolve(), KSP
1353: @*/
1354: PetscErrorCode  KSPReset(KSP ksp)
1355: {
1357:   if (!ksp) return 0;
1358:   if (ksp->ops->reset) {
1359:     (*ksp->ops->reset)(ksp);
1360:   }
1361:   if (ksp->pc) PCReset(ksp->pc);
1362:   if (ksp->guess) {
1363:     KSPGuess guess = ksp->guess;
1364:     if (guess->ops->reset) (*guess->ops->reset)(guess);
1365:   }
1366:   VecDestroyVecs(ksp->nwork,&ksp->work);
1367:   VecDestroy(&ksp->vec_rhs);
1368:   VecDestroy(&ksp->vec_sol);
1369:   VecDestroy(&ksp->diagonal);
1370:   VecDestroy(&ksp->truediagonal);

1372:   KSPResetViewers(ksp);

1374:   ksp->setupstage = KSP_SETUP_NEW;
1375:   ksp->nmax = PETSC_DECIDE;
1376:   return 0;
1377: }

1379: /*@C
1380:    KSPDestroy - Destroys KSP context.

1382:    Collective on ksp

1384:    Input Parameter:
1385: .  ksp - iterative context obtained from KSPCreate()

1387:    Level: beginner

1389: .seealso: KSPCreate(), KSPSetUp(), KSPSolve(), KSP
1390: @*/
1391: PetscErrorCode  KSPDestroy(KSP *ksp)
1392: {
1393:   PC             pc;

1395:   if (!*ksp) return 0;
1397:   if (--((PetscObject)(*ksp))->refct > 0) {*ksp = NULL; return 0;}

1399:   PetscObjectSAWsViewOff((PetscObject)*ksp);

1401:   /*
1402:    Avoid a cascading call to PCReset(ksp->pc) from the following call:
1403:    PCReset() shouldn't be called from KSPDestroy() as it is unprotected by pc's
1404:    refcount (and may be shared, e.g., by other ksps).
1405:    */
1406:   pc         = (*ksp)->pc;
1407:   (*ksp)->pc = NULL;
1408:   KSPReset((*ksp));
1409:   (*ksp)->pc = pc;
1410:   if ((*ksp)->ops->destroy) (*(*ksp)->ops->destroy)(*ksp);

1412:   if ((*ksp)->transpose.use_explicittranspose) {
1413:     MatDestroy(&(*ksp)->transpose.AT);
1414:     MatDestroy(&(*ksp)->transpose.BT);
1415:     (*ksp)->transpose.reuse_transpose = PETSC_FALSE;
1416:   }

1418:   KSPGuessDestroy(&(*ksp)->guess);
1419:   DMDestroy(&(*ksp)->dm);
1420:   PCDestroy(&(*ksp)->pc);
1421:   PetscFree((*ksp)->res_hist_alloc);
1422:   PetscFree((*ksp)->err_hist_alloc);
1423:   if ((*ksp)->convergeddestroy) {
1424:     (*(*ksp)->convergeddestroy)((*ksp)->cnvP);
1425:   }
1426:   KSPMonitorCancel((*ksp));
1427:   KSPConvergedReasonViewCancel((*ksp));
1428:   PetscViewerDestroy(&(*ksp)->eigviewer);
1429:   PetscHeaderDestroy(ksp);
1430:   return 0;
1431: }

1433: /*@
1434:     KSPSetPCSide - Sets the preconditioning side.

1436:     Logically Collective on ksp

1438:     Input Parameter:
1439: .   ksp - iterative context obtained from KSPCreate()

1441:     Output Parameter:
1442: .   side - the preconditioning side, where side is one of
1443: .vb
1444:       PC_LEFT - left preconditioning (default)
1445:       PC_RIGHT - right preconditioning
1446:       PC_SYMMETRIC - symmetric preconditioning
1447: .ve

1449:     Options Database Keys:
1450: .   -ksp_pc_side <right,left,symmetric> - KSP preconditioner side

1452:     Notes:
1453:     Left preconditioning is used by default for most Krylov methods except KSPFGMRES which only supports right preconditioning.

1455:     For methods changing the side of the preconditioner changes the norm type that is used, see KSPSetNormType().

1457:     Symmetric preconditioning is currently available only for the KSPQCG method. Note, however, that
1458:     symmetric preconditioning can be emulated by using either right or left
1459:     preconditioning and a pre or post processing step.

1461:     Setting the PC side often affects the default norm type.  See KSPSetNormType() for details.

1463:     Level: intermediate

1465: .seealso: KSPGetPCSide(), KSPSetNormType(), KSPGetNormType(), KSP
1466: @*/
1467: PetscErrorCode  KSPSetPCSide(KSP ksp,PCSide side)
1468: {
1471:   ksp->pc_side = ksp->pc_side_set = side;
1472:   return 0;
1473: }

1475: /*@
1476:     KSPGetPCSide - Gets the preconditioning side.

1478:     Not Collective

1480:     Input Parameter:
1481: .   ksp - iterative context obtained from KSPCreate()

1483:     Output Parameter:
1484: .   side - the preconditioning side, where side is one of
1485: .vb
1486:       PC_LEFT - left preconditioning (default)
1487:       PC_RIGHT - right preconditioning
1488:       PC_SYMMETRIC - symmetric preconditioning
1489: .ve

1491:     Level: intermediate

1493: .seealso: KSPSetPCSide(), KSP
1494: @*/
1495: PetscErrorCode  KSPGetPCSide(KSP ksp,PCSide *side)
1496: {
1499:   KSPSetUpNorms_Private(ksp,PETSC_TRUE,&ksp->normtype,&ksp->pc_side);
1500:   *side = ksp->pc_side;
1501:   return 0;
1502: }

1504: /*@
1505:    KSPGetTolerances - Gets the relative, absolute, divergence, and maximum
1506:    iteration tolerances used by the default KSP convergence tests.

1508:    Not Collective

1510:    Input Parameter:
1511: .  ksp - the Krylov subspace context

1513:    Output Parameters:
1514: +  rtol - the relative convergence tolerance
1515: .  abstol - the absolute convergence tolerance
1516: .  dtol - the divergence tolerance
1517: -  maxits - maximum number of iterations

1519:    Notes:
1520:    The user can specify NULL for any parameter that is not needed.

1522:    Level: intermediate

1524:            maximum, iterations

1526: .seealso: KSPSetTolerances(), KSP
1527: @*/
1528: PetscErrorCode  KSPGetTolerances(KSP ksp,PetscReal *rtol,PetscReal *abstol,PetscReal *dtol,PetscInt *maxits)
1529: {
1531:   if (abstol) *abstol = ksp->abstol;
1532:   if (rtol) *rtol = ksp->rtol;
1533:   if (dtol) *dtol = ksp->divtol;
1534:   if (maxits) *maxits = ksp->max_it;
1535:   return 0;
1536: }

1538: /*@
1539:    KSPSetTolerances - Sets the relative, absolute, divergence, and maximum
1540:    iteration tolerances used by the default KSP convergence testers.

1542:    Logically Collective on ksp

1544:    Input Parameters:
1545: +  ksp - the Krylov subspace context
1546: .  rtol - the relative convergence tolerance, relative decrease in the (possibly preconditioned) residual norm
1547: .  abstol - the absolute convergence tolerance   absolute size of the (possibly preconditioned) residual norm
1548: .  dtol - the divergence tolerance,   amount (possibly preconditioned) residual norm can increase before KSPConvergedDefault() concludes that the method is diverging
1549: -  maxits - maximum number of iterations to use

1551:    Options Database Keys:
1552: +  -ksp_atol <abstol> - Sets abstol
1553: .  -ksp_rtol <rtol> - Sets rtol
1554: .  -ksp_divtol <dtol> - Sets dtol
1555: -  -ksp_max_it <maxits> - Sets maxits

1557:    Notes:
1558:    Use PETSC_DEFAULT to retain the default value of any of the tolerances.

1560:    See KSPConvergedDefault() for details how these parameters are used in the default convergence test.  See also KSPSetConvergenceTest()
1561:    for setting user-defined stopping criteria.

1563:    Level: intermediate

1565:            convergence, maximum, iterations

1567: .seealso: KSPGetTolerances(), KSPConvergedDefault(), KSPSetConvergenceTest(), KSP
1568: @*/
1569: PetscErrorCode  KSPSetTolerances(KSP ksp,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt maxits)
1570: {

1577:   if (rtol != PETSC_DEFAULT) {
1579:     ksp->rtol = rtol;
1580:   }
1581:   if (abstol != PETSC_DEFAULT) {
1583:     ksp->abstol = abstol;
1584:   }
1585:   if (dtol != PETSC_DEFAULT) {
1587:     ksp->divtol = dtol;
1588:   }
1589:   if (maxits != PETSC_DEFAULT) {
1591:     ksp->max_it = maxits;
1592:   }
1593:   return 0;
1594: }

1596: /*@
1597:    KSPSetInitialGuessNonzero - Tells the iterative solver that the
1598:    initial guess is nonzero; otherwise KSP assumes the initial guess
1599:    is to be zero (and thus zeros it out before solving).

1601:    Logically Collective on ksp

1603:    Input Parameters:
1604: +  ksp - iterative context obtained from KSPCreate()
1605: -  flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero

1607:    Options database keys:
1608: .  -ksp_initial_guess_nonzero <true,false> - use nonzero initial guess

1610:    Level: beginner

1612:    Notes:
1613:     If this is not called the X vector is zeroed in the call to KSPSolve().

1615: .seealso: KSPGetInitialGuessNonzero(), KSPSetGuessType(), KSPGuessType, KSP
1616: @*/
1617: PetscErrorCode  KSPSetInitialGuessNonzero(KSP ksp,PetscBool flg)
1618: {
1621:   ksp->guess_zero = (PetscBool) !(int)flg;
1622:   return 0;
1623: }

1625: /*@
1626:    KSPGetInitialGuessNonzero - Determines whether the KSP solver is using
1627:    a zero initial guess.

1629:    Not Collective

1631:    Input Parameter:
1632: .  ksp - iterative context obtained from KSPCreate()

1634:    Output Parameter:
1635: .  flag - PETSC_TRUE if guess is nonzero, else PETSC_FALSE

1637:    Level: intermediate

1639: .seealso: KSPSetInitialGuessNonzero(), KSP
1640: @*/
1641: PetscErrorCode  KSPGetInitialGuessNonzero(KSP ksp,PetscBool  *flag)
1642: {
1645:   if (ksp->guess_zero) *flag = PETSC_FALSE;
1646:   else *flag = PETSC_TRUE;
1647:   return 0;
1648: }

1650: /*@
1651:    KSPSetErrorIfNotConverged - Causes KSPSolve() to generate an error if the solver has not converged.

1653:    Logically Collective on ksp

1655:    Input Parameters:
1656: +  ksp - iterative context obtained from KSPCreate()
1657: -  flg - PETSC_TRUE indicates you want the error generated

1659:    Options database keys:
1660: .  -ksp_error_if_not_converged <true,false> - generate an error and stop the program

1662:    Level: intermediate

1664:    Notes:
1665:     Normally PETSc continues if a linear solver fails to converge, you can call KSPGetConvergedReason() after a KSPSolve()
1666:     to determine if it has converged.

1668:    A KSP_DIVERGED_ITS will not generate an error in a KSPSolve() inside a nested linear solver

1670: .seealso: KSPGetErrorIfNotConverged(), KSP
1671: @*/
1672: PetscErrorCode  KSPSetErrorIfNotConverged(KSP ksp,PetscBool flg)
1673: {
1676:   ksp->errorifnotconverged = flg;
1677:   return 0;
1678: }

1680: /*@
1681:    KSPGetErrorIfNotConverged - Will KSPSolve() generate an error if the solver does not converge?

1683:    Not Collective

1685:    Input Parameter:
1686: .  ksp - iterative context obtained from KSPCreate()

1688:    Output Parameter:
1689: .  flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE

1691:    Level: intermediate

1693: .seealso: KSPSetErrorIfNotConverged(), KSP
1694: @*/
1695: PetscErrorCode  KSPGetErrorIfNotConverged(KSP ksp,PetscBool  *flag)
1696: {
1699:   *flag = ksp->errorifnotconverged;
1700:   return 0;
1701: }

1703: /*@
1704:    KSPSetInitialGuessKnoll - Tells the iterative solver to use PCApply(pc,b,..) to compute the initial guess (The Knoll trick)

1706:    Logically Collective on ksp

1708:    Input Parameters:
1709: +  ksp - iterative context obtained from KSPCreate()
1710: -  flg - PETSC_TRUE or PETSC_FALSE

1712:    Level: advanced

1714:    Developer Note: the Knoll trick is not currently implemented using the KSPGuess class

1716: .seealso: KSPGetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero(), KSP
1717: @*/
1718: PetscErrorCode  KSPSetInitialGuessKnoll(KSP ksp,PetscBool flg)
1719: {
1722:   ksp->guess_knoll = flg;
1723:   return 0;
1724: }

1726: /*@
1727:    KSPGetInitialGuessKnoll - Determines whether the KSP solver is using the Knoll trick (using PCApply(pc,b,...) to compute
1728:      the initial guess

1730:    Not Collective

1732:    Input Parameter:
1733: .  ksp - iterative context obtained from KSPCreate()

1735:    Output Parameter:
1736: .  flag - PETSC_TRUE if using Knoll trick, else PETSC_FALSE

1738:    Level: advanced

1740: .seealso: KSPSetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero(), KSP
1741: @*/
1742: PetscErrorCode  KSPGetInitialGuessKnoll(KSP ksp,PetscBool  *flag)
1743: {
1746:   *flag = ksp->guess_knoll;
1747:   return 0;
1748: }

1750: /*@
1751:    KSPGetComputeSingularValues - Gets the flag indicating whether the extreme singular
1752:    values will be calculated via a Lanczos or Arnoldi process as the linear
1753:    system is solved.

1755:    Not Collective

1757:    Input Parameter:
1758: .  ksp - iterative context obtained from KSPCreate()

1760:    Output Parameter:
1761: .  flg - PETSC_TRUE or PETSC_FALSE

1763:    Options Database Key:
1764: .  -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()

1766:    Notes:
1767:    Currently this option is not valid for all iterative methods.

1769:    Many users may just want to use the monitoring routine
1770:    KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1771:    to print the singular values at each iteration of the linear solve.

1773:    Level: advanced

1775: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue(), KSP
1776: @*/
1777: PetscErrorCode  KSPGetComputeSingularValues(KSP ksp,PetscBool  *flg)
1778: {
1781:   *flg = ksp->calc_sings;
1782:   return 0;
1783: }

1785: /*@
1786:    KSPSetComputeSingularValues - Sets a flag so that the extreme singular
1787:    values will be calculated via a Lanczos or Arnoldi process as the linear
1788:    system is solved.

1790:    Logically Collective on ksp

1792:    Input Parameters:
1793: +  ksp - iterative context obtained from KSPCreate()
1794: -  flg - PETSC_TRUE or PETSC_FALSE

1796:    Options Database Key:
1797: .  -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()

1799:    Notes:
1800:    Currently this option is not valid for all iterative methods.

1802:    Many users may just want to use the monitoring routine
1803:    KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1804:    to print the singular values at each iteration of the linear solve.

1806:    Level: advanced

1808: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue(), KSP
1809: @*/
1810: PetscErrorCode  KSPSetComputeSingularValues(KSP ksp,PetscBool flg)
1811: {
1814:   ksp->calc_sings = flg;
1815:   return 0;
1816: }

1818: /*@
1819:    KSPGetComputeEigenvalues - Gets the flag indicating that the extreme eigenvalues
1820:    values will be calculated via a Lanczos or Arnoldi process as the linear
1821:    system is solved.

1823:    Not Collective

1825:    Input Parameter:
1826: .  ksp - iterative context obtained from KSPCreate()

1828:    Output Parameter:
1829: .  flg - PETSC_TRUE or PETSC_FALSE

1831:    Notes:
1832:    Currently this option is not valid for all iterative methods.

1834:    Level: advanced

1836: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly(), KSP
1837: @*/
1838: PetscErrorCode  KSPGetComputeEigenvalues(KSP ksp,PetscBool  *flg)
1839: {
1842:   *flg = ksp->calc_sings;
1843:   return 0;
1844: }

1846: /*@
1847:    KSPSetComputeEigenvalues - Sets a flag so that the extreme eigenvalues
1848:    values will be calculated via a Lanczos or Arnoldi process as the linear
1849:    system is solved.

1851:    Logically Collective on ksp

1853:    Input Parameters:
1854: +  ksp - iterative context obtained from KSPCreate()
1855: -  flg - PETSC_TRUE or PETSC_FALSE

1857:    Notes:
1858:    Currently this option is not valid for all iterative methods.

1860:    Level: advanced

1862: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly(), KSP
1863: @*/
1864: PetscErrorCode  KSPSetComputeEigenvalues(KSP ksp,PetscBool flg)
1865: {
1868:   ksp->calc_sings = flg;
1869:   return 0;
1870: }

1872: /*@
1873:    KSPSetComputeRitz - Sets a flag so that the Ritz or harmonic Ritz pairs
1874:    will be calculated via a Lanczos or Arnoldi process as the linear
1875:    system is solved.

1877:    Logically Collective on ksp

1879:    Input Parameters:
1880: +  ksp - iterative context obtained from KSPCreate()
1881: -  flg - PETSC_TRUE or PETSC_FALSE

1883:    Notes:
1884:    Currently this option is only valid for the GMRES method.

1886:    Level: advanced

1888: .seealso: KSPComputeRitz(), KSP
1889: @*/
1890: PetscErrorCode  KSPSetComputeRitz(KSP ksp, PetscBool flg)
1891: {
1894:   ksp->calc_ritz = flg;
1895:   return 0;
1896: }

1898: /*@
1899:    KSPGetRhs - Gets the right-hand-side vector for the linear system to
1900:    be solved.

1902:    Not Collective

1904:    Input Parameter:
1905: .  ksp - iterative context obtained from KSPCreate()

1907:    Output Parameter:
1908: .  r - right-hand-side vector

1910:    Level: developer

1912: .seealso: KSPGetSolution(), KSPSolve(), KSP
1913: @*/
1914: PetscErrorCode  KSPGetRhs(KSP ksp,Vec *r)
1915: {
1918:   *r = ksp->vec_rhs;
1919:   return 0;
1920: }

1922: /*@
1923:    KSPGetSolution - Gets the location of the solution for the
1924:    linear system to be solved.  Note that this may not be where the solution
1925:    is stored during the iterative process; see KSPBuildSolution().

1927:    Not Collective

1929:    Input Parameters:
1930: .  ksp - iterative context obtained from KSPCreate()

1932:    Output Parameters:
1933: .  v - solution vector

1935:    Level: developer

1937: .seealso: KSPGetRhs(),  KSPBuildSolution(), KSPSolve(), KSP
1938: @*/
1939: PetscErrorCode  KSPGetSolution(KSP ksp,Vec *v)
1940: {
1943:   *v = ksp->vec_sol;
1944:   return 0;
1945: }

1947: /*@
1948:    KSPSetPC - Sets the preconditioner to be used to calculate the
1949:    application of the preconditioner on a vector.

1951:    Collective on ksp

1953:    Input Parameters:
1954: +  ksp - iterative context obtained from KSPCreate()
1955: -  pc   - the preconditioner object (can be NULL)

1957:    Notes:
1958:    Use KSPGetPC() to retrieve the preconditioner context.

1960:    Level: developer

1962: .seealso: KSPGetPC(), KSP
1963: @*/
1964: PetscErrorCode  KSPSetPC(KSP ksp,PC pc)
1965: {
1967:   if (pc) {
1970:   }
1971:   PetscObjectReference((PetscObject)pc);
1972:   PCDestroy(&ksp->pc);
1973:   ksp->pc = pc;
1974:   PetscLogObjectParent((PetscObject)ksp,(PetscObject)ksp->pc);
1975:   return 0;
1976: }

1978: /*@
1979:    KSPGetPC - Returns a pointer to the preconditioner context
1980:    set with KSPSetPC().

1982:    Not Collective

1984:    Input Parameters:
1985: .  ksp - iterative context obtained from KSPCreate()

1987:    Output Parameter:
1988: .  pc - preconditioner context

1990:    Level: developer

1992: .seealso: KSPSetPC(), KSP
1993: @*/
1994: PetscErrorCode  KSPGetPC(KSP ksp,PC *pc)
1995: {
1998:   if (!ksp->pc) {
1999:     PCCreate(PetscObjectComm((PetscObject)ksp),&ksp->pc);
2000:     PetscObjectIncrementTabLevel((PetscObject)ksp->pc,(PetscObject)ksp,0);
2001:     PetscLogObjectParent((PetscObject)ksp,(PetscObject)ksp->pc);
2002:     PetscObjectSetOptions((PetscObject)ksp->pc,((PetscObject)ksp)->options);
2003:   }
2004:   *pc = ksp->pc;
2005:   return 0;
2006: }

2008: /*@
2009:    KSPMonitor - runs the user provided monitor routines, if they exist

2011:    Collective on ksp

2013:    Input Parameters:
2014: +  ksp - iterative context obtained from KSPCreate()
2015: .  it - iteration number
2016: -  rnorm - relative norm of the residual

2018:    Notes:
2019:    This routine is called by the KSP implementations.
2020:    It does not typically need to be called by the user.

2022:    Level: developer

2024: .seealso: KSPMonitorSet()
2025: @*/
2026: PetscErrorCode KSPMonitor(KSP ksp,PetscInt it,PetscReal rnorm)
2027: {
2028:   PetscInt       i, n = ksp->numbermonitors;

2030:   for (i=0; i<n; i++) {
2031:     (*ksp->monitor[i])(ksp,it,rnorm,ksp->monitorcontext[i]);
2032:   }
2033:   return 0;
2034: }

2036: /*@C
2037:    KSPMonitorSet - Sets an ADDITIONAL function to be called at every iteration to monitor
2038:    the residual/error etc.

2040:    Logically Collective on ksp

2042:    Input Parameters:
2043: +  ksp - iterative context obtained from KSPCreate()
2044: .  monitor - pointer to function (if this is NULL, it turns off monitoring
2045: .  mctx    - [optional] context for private data for the
2046:              monitor routine (use NULL if no context is desired)
2047: -  monitordestroy - [optional] routine that frees monitor context
2048:           (may be NULL)

2050:    Calling Sequence of monitor:
2051: $     monitor (KSP ksp, PetscInt it, PetscReal rnorm, void *mctx)

2053: +  ksp - iterative context obtained from KSPCreate()
2054: .  it - iteration number
2055: .  rnorm - (estimated) 2-norm of (preconditioned) residual
2056: -  mctx  - optional monitoring context, as set by KSPMonitorSet()

2058:    Options Database Keys:
2059: +    -ksp_monitor               - sets KSPMonitorResidual()
2060: .    -ksp_monitor draw          - sets KSPMonitorResidualDraw() and plots residual
2061: .    -ksp_monitor draw::draw_lg - sets KSPMonitorResidualDrawLG() and plots residual
2062: .    -ksp_monitor_pause_final   - Pauses any graphics when the solve finishes (only works for internal monitors)
2063: .    -ksp_monitor_true_residual - sets KSPMonitorTrueResidual()
2064: .    -ksp_monitor_true_residual draw::draw_lg - sets KSPMonitorTrueResidualDrawLG() and plots residual
2065: .    -ksp_monitor_max           - sets KSPMonitorTrueResidualMax()
2066: .    -ksp_monitor_singular_value - sets KSPMonitorSingularValue()
2067: -    -ksp_monitor_cancel - cancels all monitors that have
2068:                           been hardwired into a code by
2069:                           calls to KSPMonitorSet(), but
2070:                           does not cancel those set via
2071:                           the options database.

2073:    Notes:
2074:    The default is to do nothing.  To print the residual, or preconditioned
2075:    residual if KSPSetNormType(ksp,KSP_NORM_PRECONDITIONED) was called, use
2076:    KSPMonitorResidual() as the monitoring routine, with a ASCII viewer as the
2077:    context.

2079:    Several different monitoring routines may be set by calling
2080:    KSPMonitorSet() multiple times; all will be called in the
2081:    order in which they were set.

2083:    Fortran Notes:
2084:     Only a single monitor function can be set for each KSP object

2086:    Level: beginner

2088: .seealso: KSPMonitorResidual(), KSPMonitorCancel(), KSP
2089: @*/
2090: PetscErrorCode  KSPMonitorSet(KSP ksp,PetscErrorCode (*monitor)(KSP,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
2091: {
2092:   PetscInt       i;
2093:   PetscBool      identical;

2096:   for (i=0; i<ksp->numbermonitors;i++) {
2097:     PetscMonitorCompare((PetscErrorCode (*)(void))monitor,mctx,monitordestroy,(PetscErrorCode (*)(void))ksp->monitor[i],ksp->monitorcontext[i],ksp->monitordestroy[i],&identical);
2098:     if (identical) return 0;
2099:   }
2101:   ksp->monitor[ksp->numbermonitors]          = monitor;
2102:   ksp->monitordestroy[ksp->numbermonitors]   = monitordestroy;
2103:   ksp->monitorcontext[ksp->numbermonitors++] = (void*)mctx;
2104:   return 0;
2105: }

2107: /*@
2108:    KSPMonitorCancel - Clears all monitors for a KSP object.

2110:    Logically Collective on ksp

2112:    Input Parameters:
2113: .  ksp - iterative context obtained from KSPCreate()

2115:    Options Database Key:
2116: .  -ksp_monitor_cancel - Cancels all monitors that have
2117:     been hardwired into a code by calls to KSPMonitorSet(),
2118:     but does not cancel those set via the options database.

2120:    Level: intermediate

2122: .seealso: KSPMonitorResidual(), KSPMonitorSet(), KSP
2123: @*/
2124: PetscErrorCode  KSPMonitorCancel(KSP ksp)
2125: {
2126:   PetscInt       i;

2129:   for (i=0; i<ksp->numbermonitors; i++) {
2130:     if (ksp->monitordestroy[i]) {
2131:       (*ksp->monitordestroy[i])(&ksp->monitorcontext[i]);
2132:     }
2133:   }
2134:   ksp->numbermonitors = 0;
2135:   return 0;
2136: }

2138: /*@C
2139:    KSPGetMonitorContext - Gets the monitoring context, as set by
2140:    KSPMonitorSet() for the FIRST monitor only.

2142:    Not Collective

2144:    Input Parameter:
2145: .  ksp - iterative context obtained from KSPCreate()

2147:    Output Parameter:
2148: .  ctx - monitoring context

2150:    Level: intermediate

2152: .seealso: KSPMonitorResidual(), KSP
2153: @*/
2154: PetscErrorCode  KSPGetMonitorContext(KSP ksp,void *ctx)
2155: {
2157:   *(void**)ctx = ksp->monitorcontext[0];
2158:   return 0;
2159: }

2161: /*@
2162:    KSPSetResidualHistory - Sets the array used to hold the residual history.
2163:    If set, this array will contain the residual norms computed at each
2164:    iteration of the solver.

2166:    Not Collective

2168:    Input Parameters:
2169: +  ksp - iterative context obtained from KSPCreate()
2170: .  a   - array to hold history
2171: .  na  - size of a
2172: -  reset - PETSC_TRUE indicates the history counter is reset to zero
2173:            for each new linear solve

2175:    Level: advanced

2177:    Notes:
2178:    If provided, he array is NOT freed by PETSc so the user needs to keep track of it and destroy once the KSP object is destroyed.
2179:    If 'a' is NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
2180:    default array of length 10000 is allocated.

2182: .seealso: KSPGetResidualHistory(), KSP

2184: @*/
2185: PetscErrorCode KSPSetResidualHistory(KSP ksp,PetscReal a[],PetscInt na,PetscBool reset)
2186: {

2189:   PetscFree(ksp->res_hist_alloc);
2190:   if (na != PETSC_DECIDE && na != PETSC_DEFAULT && a) {
2191:     ksp->res_hist     = a;
2192:     ksp->res_hist_max = (size_t) na;
2193:   } else {
2194:     if (na != PETSC_DECIDE && na != PETSC_DEFAULT) ksp->res_hist_max = (size_t) na;
2195:     else                                           ksp->res_hist_max = 10000; /* like default ksp->max_it */
2196:     PetscCalloc1(ksp->res_hist_max,&ksp->res_hist_alloc);

2198:     ksp->res_hist = ksp->res_hist_alloc;
2199:   }
2200:   ksp->res_hist_len   = 0;
2201:   ksp->res_hist_reset = reset;
2202:   return 0;
2203: }

2205: /*@C
2206:    KSPGetResidualHistory - Gets the array used to hold the residual history
2207:    and the number of residuals it contains.

2209:    Not Collective

2211:    Input Parameter:
2212: .  ksp - iterative context obtained from KSPCreate()

2214:    Output Parameters:
2215: +  a   - pointer to array to hold history (or NULL)
2216: -  na  - number of used entries in a (or NULL)

2218:    Level: advanced

2220:    Notes:
2221:      This array is borrowed and should not be freed by the caller.
2222:      Can only be called after a KSPSetResidualHistory() otherwise a and na are set to zero

2224:      The Fortran version of this routine has a calling sequence
2225: $   call KSPGetResidualHistory(KSP ksp, integer na, integer ierr)
2226:     note that you have passed a Fortran array into KSPSetResidualHistory() and you need
2227:     to access the residual values from this Fortran array you provided. Only the na (number of
2228:     residual norms currently held) is set.

2230: .seealso: KSPSetResidualHistory(), KSP

2232: @*/
2233: PetscErrorCode KSPGetResidualHistory(KSP ksp, const PetscReal *a[],PetscInt *na)
2234: {
2236:   if (a) *a = ksp->res_hist;
2237:   if (na) *na = (PetscInt) ksp->res_hist_len;
2238:   return 0;
2239: }

2241: /*@
2242:   KSPSetErrorHistory - Sets the array used to hold the error history. If set, this array will contain the error norms computed at each iteration of the solver.

2244:   Not Collective

2246:   Input Parameters:
2247: + ksp   - iterative context obtained from KSPCreate()
2248: . a     - array to hold history
2249: . na    - size of a
2250: - reset - PETSC_TRUE indicates the history counter is reset to zero for each new linear solve

2252:   Level: advanced

2254:   Notes:
2255:   If provided, the array is NOT freed by PETSc so the user needs to keep track of it and destroy once the KSP object is destroyed.
2256:   If 'a' is NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a default array of length 10000 is allocated.

2258: .seealso: KSPGetErrorHistory(), KSPSetResidualHistory(), KSP
2259: @*/
2260: PetscErrorCode KSPSetErrorHistory(KSP ksp, PetscReal a[], PetscInt na, PetscBool reset)
2261: {

2264:   PetscFree(ksp->err_hist_alloc);
2265:   if (na != PETSC_DECIDE && na != PETSC_DEFAULT && a) {
2266:     ksp->err_hist     = a;
2267:     ksp->err_hist_max = (size_t) na;
2268:   } else {
2269:     if (na != PETSC_DECIDE && na != PETSC_DEFAULT) ksp->err_hist_max = (size_t) na;
2270:     else                                           ksp->err_hist_max = 10000; /* like default ksp->max_it */
2271:     PetscCalloc1(ksp->err_hist_max, &ksp->err_hist_alloc);

2273:     ksp->err_hist = ksp->err_hist_alloc;
2274:   }
2275:   ksp->err_hist_len   = 0;
2276:   ksp->err_hist_reset = reset;
2277:   return 0;
2278: }

2280: /*@C
2281:   KSPGetErrorHistory - Gets the array used to hold the error history and the number of residuals it contains.

2283:   Not Collective

2285:   Input Parameter:
2286: . ksp - iterative context obtained from KSPCreate()

2288:   Output Parameters:
2289: + a  - pointer to array to hold history (or NULL)
2290: - na - number of used entries in a (or NULL)

2292:   Level: advanced

2294:   Notes:
2295:   This array is borrowed and should not be freed by the caller.
2296:   Can only be called after a KSPSetErrorHistory() otherwise a and na are set to zero
2297:   The Fortran version of this routine has a calling sequence
2298: $   call KSPGetErrorHistory(KSP ksp, integer na, integer ierr)
2299:   note that you have passed a Fortran array into KSPSetErrorHistory() and you need
2300:   to access the residual values from this Fortran array you provided. Only the na (number of
2301:   residual norms currently held) is set.

2303: .seealso: KSPSetErrorHistory(), KSPGetResidualHistory(), KSP
2304: @*/
2305: PetscErrorCode KSPGetErrorHistory(KSP ksp, const PetscReal *a[], PetscInt *na)
2306: {
2308:   if (a)  *a  = ksp->err_hist;
2309:   if (na) *na = (PetscInt) ksp->err_hist_len;
2310:   return 0;
2311: }

2313: /*
2314:   KSPComputeConvergenceRate - Compute the convergence rate for the iteration

2316:   Not collective

2318:   Input Parameter:
2319: . ksp - The KSP

2321:   Output Parameters:
2322: + cr   - The residual contraction rate
2323: . rRsq - The coefficient of determination, R^2, indicating the linearity of the data
2324: . ce   - The error contraction rate
2325: - eRsq - The coefficient of determination, R^2, indicating the linearity of the data

2327:   Note:
2328:   Suppose that the residual is reduced linearly, $r_k = c^k r_0$, which means $log r_k = log r_0 + k log c$. After linear regression,
2329:   the slope is $\log c$. The coefficient of determination is given by $1 - \frac{\sum_i (y_i - f(x_i))^2}{\sum_i (y_i - \bar y)}$,
2330:   see also https://en.wikipedia.org/wiki/Coefficient_of_determination

2332:   Level: advanced

2334: .seealso: KSPConvergedRateView()
2335: */
2336: PetscErrorCode KSPComputeConvergenceRate(KSP ksp, PetscReal *cr, PetscReal *rRsq, PetscReal *ce, PetscReal *eRsq)
2337: {
2338:   PetscReal      const *hist;
2339:   PetscReal      *x, *y, slope, intercept, mean = 0.0, var = 0.0, res = 0.0;
2340:   PetscInt       n, k;

2342:   if (cr || rRsq) {
2343:     KSPGetResidualHistory(ksp, &hist, &n);
2344:     if (!n) {
2345:       if (cr)   *cr   =  0.0;
2346:       if (rRsq) *rRsq = -1.0;
2347:     } else {
2348:       PetscMalloc2(n, &x, n, &y);
2349:       for (k = 0; k < n; ++k) {
2350:         x[k] = k;
2351:         y[k] = PetscLogReal(hist[k]);
2352:         mean += y[k];
2353:       }
2354:       mean /= n;
2355:       PetscLinearRegression(n, x, y, &slope, &intercept);
2356:       for (k = 0; k < n; ++k) {
2357:         res += PetscSqr(y[k] - (slope*x[k] + intercept));
2358:         var += PetscSqr(y[k] - mean);
2359:       }
2360:       PetscFree2(x, y);
2361:       if (cr)   *cr   = PetscExpReal(slope);
2362:       if (rRsq) *rRsq = var < PETSC_MACHINE_EPSILON ? 0.0 : 1.0 - (res / var);
2363:     }
2364:   }
2365:   if (ce || eRsq) {
2366:     KSPGetErrorHistory(ksp, &hist, &n);
2367:     if (!n) {
2368:       if (ce)   *ce   =  0.0;
2369:       if (eRsq) *eRsq = -1.0;
2370:     } else {
2371:       PetscMalloc2(n, &x, n, &y);
2372:       for (k = 0; k < n; ++k) {
2373:         x[k] = k;
2374:         y[k] = PetscLogReal(hist[k]);
2375:         mean += y[k];
2376:       }
2377:       mean /= n;
2378:       PetscLinearRegression(n, x, y, &slope, &intercept);
2379:       for (k = 0; k < n; ++k) {
2380:         res += PetscSqr(y[k] - (slope*x[k] + intercept));
2381:         var += PetscSqr(y[k] - mean);
2382:       }
2383:       PetscFree2(x, y);
2384:       if (ce)   *ce   = PetscExpReal(slope);
2385:       if (eRsq) *eRsq = var < PETSC_MACHINE_EPSILON ? 0.0 : 1.0 - (res / var);
2386:     }
2387:   }
2388:   return 0;
2389: }

2391: /*@C
2392:    KSPSetConvergenceTest - Sets the function to be used to determine
2393:    convergence.

2395:    Logically Collective on ksp

2397:    Input Parameters:
2398: +  ksp - iterative context obtained from KSPCreate()
2399: .  converge - pointer to the function
2400: .  cctx    - context for private data for the convergence routine (may be null)
2401: -  destroy - a routine for destroying the context (may be null)

2403:    Calling sequence of converge:
2404: $     converge (KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)

2406: +  ksp - iterative context obtained from KSPCreate()
2407: .  it - iteration number
2408: .  rnorm - (estimated) 2-norm of (preconditioned) residual
2409: .  reason - the reason why it has converged or diverged
2410: -  cctx  - optional convergence context, as set by KSPSetConvergenceTest()

2412:    Notes:
2413:    Must be called after the KSP type has been set so put this after
2414:    a call to KSPSetType(), or KSPSetFromOptions().

2416:    The default convergence test, KSPConvergedDefault(), aborts if the
2417:    residual grows to more than 10000 times the initial residual.

2419:    The default is a combination of relative and absolute tolerances.
2420:    The residual value that is tested may be an approximation; routines
2421:    that need exact values should compute them.

2423:    In the default PETSc convergence test, the precise values of reason
2424:    are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.

2426:    Level: advanced

2428: .seealso: KSPConvergedDefault(), KSPGetConvergenceContext(), KSPSetTolerances(), KSP, KSPGetConvergenceTest(), KSPGetAndClearConvergenceTest()
2429: @*/
2430: PetscErrorCode  KSPSetConvergenceTest(KSP ksp,PetscErrorCode (*converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
2431: {
2433:   if (ksp->convergeddestroy) {
2434:     (*ksp->convergeddestroy)(ksp->cnvP);
2435:   }
2436:   ksp->converged        = converge;
2437:   ksp->convergeddestroy = destroy;
2438:   ksp->cnvP             = (void*)cctx;
2439:   return 0;
2440: }

2442: /*@C
2443:    KSPGetConvergenceTest - Gets the function to be used to determine
2444:    convergence.

2446:    Logically Collective on ksp

2448:    Input Parameter:
2449: .   ksp - iterative context obtained from KSPCreate()

2451:    Output Parameters:
2452: +  converge - pointer to convergence test function
2453: .  cctx    - context for private data for the convergence routine (may be null)
2454: -  destroy - a routine for destroying the context (may be null)

2456:    Calling sequence of converge:
2457: $     converge (KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)

2459: +  ksp - iterative context obtained from KSPCreate()
2460: .  it - iteration number
2461: .  rnorm - (estimated) 2-norm of (preconditioned) residual
2462: .  reason - the reason why it has converged or diverged
2463: -  cctx  - optional convergence context, as set by KSPSetConvergenceTest()

2465:    Level: advanced

2467: .seealso: KSPConvergedDefault(), KSPGetConvergenceContext(), KSPSetTolerances(), KSP, KSPSetConvergenceTest(), KSPGetAndClearConvergenceTest()
2468: @*/
2469: PetscErrorCode  KSPGetConvergenceTest(KSP ksp,PetscErrorCode (**converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void **cctx,PetscErrorCode (**destroy)(void*))
2470: {
2472:   if (converge) *converge = ksp->converged;
2473:   if (destroy)  *destroy  = ksp->convergeddestroy;
2474:   if (cctx)     *cctx     = ksp->cnvP;
2475:   return 0;
2476: }

2478: /*@C
2479:    KSPGetAndClearConvergenceTest - Gets the function to be used to determine convergence. Removes the current test without calling destroy on the test context

2481:    Logically Collective on ksp

2483:    Input Parameter:
2484: .   ksp - iterative context obtained from KSPCreate()

2486:    Output Parameters:
2487: +  converge - pointer to convergence test function
2488: .  cctx    - context for private data for the convergence routine
2489: -  destroy - a routine for destroying the context

2491:    Calling sequence of converge:
2492: $     converge (KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)

2494: +  ksp - iterative context obtained from KSPCreate()
2495: .  it - iteration number
2496: .  rnorm - (estimated) 2-norm of (preconditioned) residual
2497: .  reason - the reason why it has converged or diverged
2498: -  cctx  - optional convergence context, as set by KSPSetConvergenceTest()

2500:    Level: advanced

2502:    Notes: This is intended to be used to allow transferring the convergence test (and its context) to another testing object (for example another KSP) and then calling
2503:           KSPSetConvergenceTest() on this original KSP. If you just called KSPGetConvergenceTest() followed by KSPSetConvergenceTest() the original context information
2504:           would be destroyed and hence the transferred context would be invalid and trigger a crash on use

2506: .seealso: KSPConvergedDefault(), KSPGetConvergenceContext(), KSPSetTolerances(), KSP, KSPSetConvergenceTest(), KSPGetConvergenceTest()
2507: @*/
2508: PetscErrorCode  KSPGetAndClearConvergenceTest(KSP ksp,PetscErrorCode (**converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void **cctx,PetscErrorCode (**destroy)(void*))
2509: {
2511:   *converge             = ksp->converged;
2512:   *destroy              = ksp->convergeddestroy;
2513:   *cctx                 = ksp->cnvP;
2514:   ksp->converged        = NULL;
2515:   ksp->cnvP             = NULL;
2516:   ksp->convergeddestroy = NULL;
2517:   return 0;
2518: }

2520: /*@C
2521:    KSPGetConvergenceContext - Gets the convergence context set with
2522:    KSPSetConvergenceTest().

2524:    Not Collective

2526:    Input Parameter:
2527: .  ksp - iterative context obtained from KSPCreate()

2529:    Output Parameter:
2530: .  ctx - monitoring context

2532:    Level: advanced

2534: .seealso: KSPConvergedDefault(), KSPSetConvergenceTest(), KSP
2535: @*/
2536: PetscErrorCode  KSPGetConvergenceContext(KSP ksp,void *ctx)
2537: {
2539:   *(void**)ctx = ksp->cnvP;
2540:   return 0;
2541: }

2543: /*@C
2544:    KSPBuildSolution - Builds the approximate solution in a vector provided.
2545:    This routine is NOT commonly needed (see KSPSolve()).

2547:    Collective on ksp

2549:    Input Parameter:
2550: .  ctx - iterative context obtained from KSPCreate()

2552:    Output Parameter:
2553:    Provide exactly one of
2554: +  v - location to stash solution.
2555: -  V - the solution is returned in this location. This vector is created
2556:        internally. This vector should NOT be destroyed by the user with
2557:        VecDestroy().

2559:    Notes:
2560:    This routine can be used in one of two ways
2561: .vb
2562:       KSPBuildSolution(ksp,NULL,&V);
2563:    or
2564:       KSPBuildSolution(ksp,v,NULL); or KSPBuildSolution(ksp,v,&v);
2565: .ve
2566:    In the first case an internal vector is allocated to store the solution
2567:    (the user cannot destroy this vector). In the second case the solution
2568:    is generated in the vector that the user provides. Note that for certain
2569:    methods, such as KSPCG, the second case requires a copy of the solution,
2570:    while in the first case the call is essentially free since it simply
2571:    returns the vector where the solution already is stored. For some methods
2572:    like GMRES this is a reasonably expensive operation and should only be
2573:    used in truly needed.

2575:    Level: advanced

2577: .seealso: KSPGetSolution(), KSPBuildResidual(), KSP
2578: @*/
2579: PetscErrorCode  KSPBuildSolution(KSP ksp,Vec v,Vec *V)
2580: {
2583:   if (!V) V = &v;
2584:   (*ksp->ops->buildsolution)(ksp,v,V);
2585:   return 0;
2586: }

2588: /*@C
2589:    KSPBuildResidual - Builds the residual in a vector provided.

2591:    Collective on ksp

2593:    Input Parameter:
2594: .  ksp - iterative context obtained from KSPCreate()

2596:    Output Parameters:
2597: +  v - optional location to stash residual.  If v is not provided,
2598:        then a location is generated.
2599: .  t - work vector.  If not provided then one is generated.
2600: -  V - the residual

2602:    Notes:
2603:    Regardless of whether or not v is provided, the residual is
2604:    returned in V.

2606:    Level: advanced

2608: .seealso: KSPBuildSolution()
2609: @*/
2610: PetscErrorCode  KSPBuildResidual(KSP ksp,Vec t,Vec v,Vec *V)
2611: {
2612:   PetscBool      flag = PETSC_FALSE;
2613:   Vec            w    = v,tt = t;

2616:   if (!w) {
2617:     VecDuplicate(ksp->vec_rhs,&w);
2618:     PetscLogObjectParent((PetscObject)ksp,(PetscObject)w);
2619:   }
2620:   if (!tt) {
2621:     VecDuplicate(ksp->vec_sol,&tt); flag = PETSC_TRUE;
2622:     PetscLogObjectParent((PetscObject)ksp,(PetscObject)tt);
2623:   }
2624:   (*ksp->ops->buildresidual)(ksp,tt,w,V);
2625:   if (flag) VecDestroy(&tt);
2626:   return 0;
2627: }

2629: /*@
2630:    KSPSetDiagonalScale - Tells KSP to symmetrically diagonally scale the system
2631:      before solving. This actually CHANGES the matrix (and right hand side).

2633:    Logically Collective on ksp

2635:    Input Parameters:
2636: +  ksp - the KSP context
2637: -  scale - PETSC_TRUE or PETSC_FALSE

2639:    Options Database Key:
2640: +   -ksp_diagonal_scale -
2641: -   -ksp_diagonal_scale_fix - scale the matrix back AFTER the solve

2643:     Notes:
2644:     Scales the matrix by  D^(-1/2)  A  D^(-1/2)  [D^(1/2) x ] = D^(-1/2) b
2645:        where D_{ii} is 1/abs(A_{ii}) unless A_{ii} is zero and then it is 1.

2647:     BE CAREFUL with this routine: it actually scales the matrix and right
2648:     hand side that define the system. After the system is solved the matrix
2649:     and right hand side remain scaled unless you use KSPSetDiagonalScaleFix()

2651:     This should NOT be used within the SNES solves if you are using a line
2652:     search.

2654:     If you use this with the PCType Eisenstat preconditioner than you can
2655:     use the PCEisenstatSetNoDiagonalScaling() option, or -pc_eisenstat_no_diagonal_scaling
2656:     to save some unneeded, redundant flops.

2658:    Level: intermediate

2660: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScaleFix(), KSP
2661: @*/
2662: PetscErrorCode  KSPSetDiagonalScale(KSP ksp,PetscBool scale)
2663: {
2666:   ksp->dscale = scale;
2667:   return 0;
2668: }

2670: /*@
2671:    KSPGetDiagonalScale - Checks if KSP solver scales the matrix and
2672:                           right hand side

2674:    Not Collective

2676:    Input Parameter:
2677: .  ksp - the KSP context

2679:    Output Parameter:
2680: .  scale - PETSC_TRUE or PETSC_FALSE

2682:    Notes:
2683:     BE CAREFUL with this routine: it actually scales the matrix and right
2684:     hand side that define the system. After the system is solved the matrix
2685:     and right hand side remain scaled  unless you use KSPSetDiagonalScaleFix()

2687:    Level: intermediate

2689: .seealso: KSPSetDiagonalScale(), KSPSetDiagonalScaleFix(), KSP
2690: @*/
2691: PetscErrorCode  KSPGetDiagonalScale(KSP ksp,PetscBool  *scale)
2692: {
2695:   *scale = ksp->dscale;
2696:   return 0;
2697: }

2699: /*@
2700:    KSPSetDiagonalScaleFix - Tells KSP to diagonally scale the system
2701:      back after solving.

2703:    Logically Collective on ksp

2705:    Input Parameters:
2706: +  ksp - the KSP context
2707: -  fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
2708:          rescale (default)

2710:    Notes:
2711:      Must be called after KSPSetDiagonalScale()

2713:      Using this will slow things down, because it rescales the matrix before and
2714:      after each linear solve. This is intended mainly for testing to allow one
2715:      to easily get back the original system to make sure the solution computed is
2716:      accurate enough.

2718:    Level: intermediate

2720: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPGetDiagonalScaleFix(), KSP
2721: @*/
2722: PetscErrorCode  KSPSetDiagonalScaleFix(KSP ksp,PetscBool fix)
2723: {
2726:   ksp->dscalefix = fix;
2727:   return 0;
2728: }

2730: /*@
2731:    KSPGetDiagonalScaleFix - Determines if KSP diagonally scales the system
2732:      back after solving.

2734:    Not Collective

2736:    Input Parameter:
2737: .  ksp - the KSP context

2739:    Output Parameter:
2740: .  fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
2741:          rescale (default)

2743:    Notes:
2744:      Must be called after KSPSetDiagonalScale()

2746:      If PETSC_TRUE will slow things down, because it rescales the matrix before and
2747:      after each linear solve. This is intended mainly for testing to allow one
2748:      to easily get back the original system to make sure the solution computed is
2749:      accurate enough.

2751:    Level: intermediate

2753: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPSetDiagonalScaleFix(), KSP
2754: @*/
2755: PetscErrorCode  KSPGetDiagonalScaleFix(KSP ksp,PetscBool  *fix)
2756: {
2759:   *fix = ksp->dscalefix;
2760:   return 0;
2761: }

2763: /*@C
2764:    KSPSetComputeOperators - set routine to compute the linear operators

2766:    Logically Collective

2768:    Input Parameters:
2769: +  ksp - the KSP context
2770: .  func - function to compute the operators
2771: -  ctx - optional context

2773:    Calling sequence of func:
2774: $  func(KSP ksp,Mat A,Mat B,void *ctx)

2776: +  ksp - the KSP context
2777: .  A - the linear operator
2778: .  B - preconditioning matrix
2779: -  ctx - optional user-provided context

2781:    Notes:
2782:     The user provided func() will be called automatically at the very next call to KSPSolve(). It will not be called at future KSPSolve() calls
2783:           unless either KSPSetComputeOperators() or KSPSetOperators() is called before that KSPSolve() is called.

2785:           To reuse the same preconditioner for the next KSPSolve() and not compute a new one based on the most recently computed matrix call KSPSetReusePreconditioner()

2787:    Level: beginner

2789: .seealso: KSPSetOperators(), KSPSetComputeRHS(), DMKSPSetComputeOperators(), KSPSetComputeInitialGuess()
2790: @*/
2791: PetscErrorCode KSPSetComputeOperators(KSP ksp,PetscErrorCode (*func)(KSP,Mat,Mat,void*),void *ctx)
2792: {
2793:   DM             dm;

2796:   KSPGetDM(ksp,&dm);
2797:   DMKSPSetComputeOperators(dm,func,ctx);
2798:   if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX;
2799:   return 0;
2800: }

2802: /*@C
2803:    KSPSetComputeRHS - set routine to compute the right hand side of the linear system

2805:    Logically Collective

2807:    Input Parameters:
2808: +  ksp - the KSP context
2809: .  func - function to compute the right hand side
2810: -  ctx - optional context

2812:    Calling sequence of func:
2813: $  func(KSP ksp,Vec b,void *ctx)

2815: +  ksp - the KSP context
2816: .  b - right hand side of linear system
2817: -  ctx - optional user-provided context

2819:    Notes:
2820:     The routine you provide will be called EACH you call KSPSolve() to prepare the new right hand side for that solve

2822:    Level: beginner

2824: .seealso: KSPSolve(), DMKSPSetComputeRHS(), KSPSetComputeOperators()
2825: @*/
2826: PetscErrorCode KSPSetComputeRHS(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2827: {
2828:   DM             dm;

2831:   KSPGetDM(ksp,&dm);
2832:   DMKSPSetComputeRHS(dm,func,ctx);
2833:   return 0;
2834: }

2836: /*@C
2837:    KSPSetComputeInitialGuess - set routine to compute the initial guess of the linear system

2839:    Logically Collective

2841:    Input Parameters:
2842: +  ksp - the KSP context
2843: .  func - function to compute the initial guess
2844: -  ctx - optional context

2846:    Calling sequence of func:
2847: $  func(KSP ksp,Vec x,void *ctx)

2849: +  ksp - the KSP context
2850: .  x - solution vector
2851: -  ctx - optional user-provided context

2853:    Notes: This should only be used in conjunction with KSPSetComputeRHS(), KSPSetComputeOperators(), otherwise
2854:    call KSPSetInitialGuessNonzero() and set the initial guess values in the solution vector passed to KSPSolve().

2856:    Level: beginner

2858: .seealso: KSPSolve(), KSPSetComputeRHS(), KSPSetComputeOperators(), DMKSPSetComputeInitialGuess()
2859: @*/
2860: PetscErrorCode KSPSetComputeInitialGuess(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2861: {
2862:   DM             dm;

2865:   KSPGetDM(ksp,&dm);
2866:   DMKSPSetComputeInitialGuess(dm,func,ctx);
2867:   return 0;
2868: }

2870: /*@
2871:    KSPSetUseExplicitTranspose - Determines if transpose the system explicitly
2872:    in KSPSolveTranspose.

2874:    Logically Collective on ksp

2876:    Input Parameter:
2877: .  ksp - the KSP context

2879:    Output Parameter:
2880: .  flg - PETSC_TRUE to transpose the system in KSPSolveTranspose, PETSC_FALSE to not
2881:          transpose (default)

2883:    Level: advanced

2885: .seealso: KSPSolveTranspose(), KSP
2886: @*/
2887: PetscErrorCode KSPSetUseExplicitTranspose(KSP ksp,PetscBool flg)
2888: {
2891:   ksp->transpose.use_explicittranspose = flg;
2892:   return 0;
2893: }