Actual source code: iterativ.c

petsc-3.15.0 2021-04-05
Report Typos and Errors
  1: /*
  2:    This file contains some simple default routines.
  3:    These routines should be SHORT, since they will be included in every
  4:    executable image that uses the iterative routines (note that, through
  5:    the registry system, we provide a way to load only the truly necessary
  6:    files)
  7:  */
  8: #include <petsc/private/kspimpl.h>
  9: #include <petscdmshell.h>
 10: #include <petscdraw.h>

 12: /*@
 13:    KSPGetResidualNorm - Gets the last (approximate preconditioned)
 14:    residual norm that has been computed.

 16:    Not Collective

 18:    Input Parameters:
 19: .  ksp - the iterative context

 21:    Output Parameters:
 22: .  rnorm - residual norm

 24:    Level: intermediate

 26: .seealso: KSPBuildResidual()
 27: @*/
 28: PetscErrorCode  KSPGetResidualNorm(KSP ksp,PetscReal *rnorm)
 29: {
 33:   *rnorm = ksp->rnorm;
 34:   return(0);
 35: }

 37: /*@
 38:    KSPGetIterationNumber - Gets the current iteration number; if the
 39:          KSPSolve() is complete, returns the number of iterations
 40:          used.

 42:    Not Collective

 44:    Input Parameters:
 45: .  ksp - the iterative context

 47:    Output Parameters:
 48: .  its - number of iterations

 50:    Level: intermediate

 52:    Notes:
 53:       During the ith iteration this returns i-1
 54: .seealso: KSPBuildResidual(), KSPGetResidualNorm(), KSPGetTotalIterations()
 55: @*/
 56: PetscErrorCode  KSPGetIterationNumber(KSP ksp,PetscInt *its)
 57: {
 61:   *its = ksp->its;
 62:   return(0);
 63: }

 65: /*@
 66:    KSPGetTotalIterations - Gets the total number of iterations this KSP object has performed since was created, counted over all linear solves

 68:    Not Collective

 70:    Input Parameters:
 71: .  ksp - the iterative context

 73:    Output Parameters:
 74: .  its - total number of iterations

 76:    Level: intermediate

 78:    Notes:
 79:     Use KSPGetIterationNumber() to get the count for the most recent solve only
 80:    If this is called within a linear solve (such as in a KSPMonitor routine) then it does not include iterations within that current solve

 82: .seealso: KSPBuildResidual(), KSPGetResidualNorm(), KSPGetIterationNumber()
 83: @*/
 84: PetscErrorCode  KSPGetTotalIterations(KSP ksp,PetscInt *its)
 85: {
 89:   *its = ksp->totalits;
 90:   return(0);
 91: }

 93: /*@C
 94:   KSPMonitorResidual - Print the preconditioned residual norm at each iteration of an iterative solver.

 96:   Collective on ksp

 98:   Input Parameters:
 99: + ksp   - iterative context
100: . n     - iteration number
101: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
102: - vf    - The viewer context

104:   Options Database Key:
105: . -ksp_monitor - Activates KSPMonitorResidual()

107:   Level: intermediate

109: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
110: @*/
111: PetscErrorCode KSPMonitorResidual(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
112: {
113:   PetscViewer       viewer = vf->viewer;
114:   PetscViewerFormat format = vf->format;
115:   PetscInt          tablevel;
116:   const char       *prefix;
117:   PetscErrorCode    ierr;

121:   PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
122:   PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
123:   PetscViewerPushFormat(viewer, format);
124:   PetscViewerASCIIAddTab(viewer, tablevel);
125:   if (n == 0 && prefix) {PetscViewerASCIIPrintf(viewer, "  Residual norms for %s solve.\n", prefix);}
126:   PetscViewerASCIIPrintf(viewer, "%3D KSP Residual norm %14.12e \n", n, (double) rnorm);
127:   PetscViewerASCIISubtractTab(viewer, tablevel);
128:   PetscViewerPopFormat(viewer);
129:   return(0);
130: }

132: /*@C
133:   KSPMonitorResidualDraw - Plots the preconditioned residual at each iteration of an iterative solver.

135:   Collective on ksp

137:   Input Parameters:
138: + ksp   - iterative context
139: . n     - iteration number
140: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
141: - vf    - The viewer context

143:   Options Database Key:
144: . -ksp_monitor draw - Activates KSPMonitorResidualDraw()

146:   Level: intermediate

148: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
149: @*/
150: PetscErrorCode KSPMonitorResidualDraw(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
151: {
152:   PetscViewer       viewer = vf->viewer;
153:   PetscViewerFormat format = vf->format;
154:   Vec               r;
155:   PetscErrorCode    ierr;

159:   PetscViewerPushFormat(viewer, format);
160:   KSPBuildResidual(ksp, NULL, NULL, &r);
161:   PetscObjectSetName((PetscObject) r, "Residual");
162:   PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", (PetscObject) ksp);
163:   VecView(r, viewer);
164:   PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", NULL);
165:   VecDestroy(&r);
166:   PetscViewerPopFormat(viewer);
167:   return(0);
168: }

170: /*@C
171:   KSPMonitorResidualDrawLG - Plots the preconditioned residual norm at each iteration of an iterative solver.

173:   Collective on ksp

175:   Input Parameters:
176: + ksp   - iterative context
177: . n     - iteration number
178: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
179: - vf    - The viewer context

181:   Options Database Key:
182: . -ksp_monitor draw::draw_lg - Activates KSPMonitorResidualDrawLG()

184:   Level: intermediate

186: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
187: @*/
188: PetscErrorCode KSPMonitorResidualDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
189: {
190:   PetscViewer        viewer = vf->viewer;
191:   PetscViewerFormat  format = vf->format;
192:   PetscDrawLG        lg     = vf->lg;
193:   KSPConvergedReason reason;
194:   PetscReal          x, y;
195:   PetscErrorCode     ierr;

200:   PetscViewerPushFormat(viewer, format);
201:   if (!n) {PetscDrawLGReset(lg);}
202:   x = (PetscReal) n;
203:   if (rnorm > 0.0) y = PetscLog10Real(rnorm);
204:   else y = -15.0;
205:   PetscDrawLGAddPoint(lg, &x, &y);
206:   KSPGetConvergedReason(ksp, &reason);
207:   if (n <= 20 || !(n % 5) || reason) {
208:     PetscDrawLGDraw(lg);
209:     PetscDrawLGSave(lg);
210:   }
211:   PetscViewerPopFormat(viewer);
212:   return(0);
213: }

215: /*@C
216:   KSPMonitorResidualDrawLGCreate - Creates the plotter for the preconditioned residual.

218:   Collective on ksp

220:   Input Parameters:
221: + viewer - The PetscViewer
222: . format - The viewer format
223: - ctx    - An optional user context

225:   Output Parameter:
226: . vf    - The viewer context

228:   Level: intermediate

230: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
231: @*/
232: PetscErrorCode KSPMonitorResidualDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
233: {

237:   PetscViewerAndFormatCreate(viewer, format, vf);
238:   (*vf)->data = ctx;
239:   KSPMonitorLGCreate(PetscObjectComm((PetscObject) viewer), NULL, NULL, "Log Residual Norm", 1, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
240:   return(0);
241: }

243: /*
244:   This is the same as KSPMonitorResidual() except it prints fewer digits of the residual as the residual gets smaller.
245:   This is because the later digits are meaningless and are often different on different machines; by using this routine different
246:   machines will usually generate the same output.

248:   Deprecated: Intentionally has no manual page
249: */
250: PetscErrorCode KSPMonitorResidualShort(KSP ksp, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf)
251: {
252:   PetscViewer       viewer = vf->viewer;
253:   PetscViewerFormat format = vf->format;
254:   PetscInt          tablevel;
255:   const char       *prefix;
256:   PetscErrorCode    ierr;

260:   PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
261:   PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
262:   PetscViewerPushFormat(viewer, format);
263:   PetscViewerASCIIAddTab(viewer, tablevel);
264:   if (its == 0 && prefix)  {PetscViewerASCIIPrintf(viewer, "  Residual norms for %s solve.\n", prefix);}
265:   if (fnorm > 1.e-9)       {PetscViewerASCIIPrintf(viewer, "%3D KSP Residual norm %g \n", its, (double) fnorm);}
266:   else if (fnorm > 1.e-11) {PetscViewerASCIIPrintf(viewer, "%3D KSP Residual norm %5.3e \n", its, (double) fnorm);}
267:   else                     {PetscViewerASCIIPrintf(viewer, "%3D KSP Residual norm < 1.e-11\n", its);}
268:   PetscViewerASCIISubtractTab(viewer, tablevel);
269:   PetscViewerPopFormat(viewer);
270:   return(0);
271: }

273: PetscErrorCode KSPMonitorRange_Private(KSP ksp, PetscInt it, PetscReal *per)
274: {
275:   Vec                resid;
276:   const PetscScalar *r;
277:   PetscReal          rmax, pwork;
278:   PetscInt           i, n, N;
279:   PetscErrorCode     ierr;

282:   KSPBuildResidual(ksp, NULL, NULL, &resid);
283:   VecNorm(resid, NORM_INFINITY, &rmax);
284:   VecGetLocalSize(resid, &n);
285:   VecGetSize(resid, &N);
286:   VecGetArrayRead(resid, &r);
287:   pwork = 0.0;
288:   for (i = 0; i < n; ++i) pwork += (PetscAbsScalar(r[i]) > .20*rmax);
289:   VecRestoreArrayRead(resid, &r);
290:   VecDestroy(&resid);
291:   MPIU_Allreduce(&pwork, per, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) ksp));
292:   *per = *per/N;
293:   return(0);
294: }

296: /*@C
297:   KSPMonitorResidualRange - Prints the percentage of residual elements that are more then 10 percent of the maximum value.

299:   Collective on ksp

301:   Input Parameters:
302: + ksp   - iterative context
303: . it    - iteration number
304: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
305: - vf    - The viewer context

307:   Options Database Key:
308: . -ksp_monitor_range - Activates KSPMonitorResidualRange()

310:   Level: intermediate

312: .seealso: KSPMonitorSet(), KSPMonitorResidual()
313: @*/
314: PetscErrorCode KSPMonitorResidualRange(KSP ksp, PetscInt it, PetscReal rnorm, PetscViewerAndFormat *vf)
315: {
316:   static PetscReal  prev;
317:   PetscViewer       viewer = vf->viewer;
318:   PetscViewerFormat format = vf->format;
319:   PetscInt          tablevel;
320:   const char       *prefix;
321:   PetscReal         perc, rel;
322:   PetscErrorCode    ierr;

326:   PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
327:   PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
328:   PetscViewerPushFormat(viewer, format);
329:   PetscViewerASCIIAddTab(viewer, tablevel);
330:   if (!it) prev = rnorm;
331:   if (it == 0 && prefix) {PetscViewerASCIIPrintf(viewer, "  Residual norms for %s solve.\n", prefix);}
332:   KSPMonitorRange_Private(ksp, it, &perc);
333:   rel  = (prev - rnorm)/prev;
334:   prev = rnorm;
335:   PetscViewerASCIIPrintf(viewer, "%3D KSP preconditioned resid norm %14.12e Percent values above 20 percent of maximum %5.2f relative decrease %5.2e ratio %5.2e \n", it, (double) rnorm, (double) (100.0*perc), (double) rel, (double) (rel/perc));
336:   PetscViewerASCIISubtractTab(viewer, tablevel);
337:   PetscViewerPopFormat(viewer);
338:   return(0);
339: }

341: /*@C
342:   KSPMonitorTrueResidual - Prints the true residual norm, as well as the preconditioned residual norm, at each iteration of an iterative solver.

344:   Collective on ksp

346:   Input Parameters:
347: + ksp   - iterative context
348: . n     - iteration number
349: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
350: - vf    - The viewer context

352:   Options Database Key:
353: . -ksp_monitor_true_residual - Activates KSPMonitorTrueResidual()

355:   Notes:
356:   When using right preconditioning, these values are equivalent.

358:   Level: intermediate

360: .seealso: KSPMonitorSet(), KSPMonitorResidual(),KSPMonitorTrueResidualMaxNorm()
361: @*/
362: PetscErrorCode KSPMonitorTrueResidual(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
363: {
364:   PetscViewer       viewer = vf->viewer;
365:   PetscViewerFormat format = vf->format;
366:   Vec               r;
367:   PetscReal         truenorm, bnorm;
368:   char              normtype[256];
369:   PetscInt          tablevel;
370:   const char       *prefix;
371:   PetscErrorCode    ierr;

375:   PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
376:   PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
377:   PetscStrncpy(normtype, KSPNormTypes[ksp->normtype], sizeof(normtype));
378:   PetscStrtolower(normtype);
379:   KSPBuildResidual(ksp, NULL, NULL, &r);
380:   VecNorm(r, NORM_2, &truenorm);
381:   VecNorm(ksp->vec_rhs, NORM_2, &bnorm);
382:   VecDestroy(&r);

384:   PetscViewerPushFormat(viewer, format);
385:   PetscViewerASCIIAddTab(viewer, tablevel);
386:   if (n == 0 && prefix) {PetscViewerASCIIPrintf(viewer, "  Residual norms for %s solve.\n", prefix);}
387:   PetscViewerASCIIPrintf(viewer, "%3D KSP %s resid norm %14.12e true resid norm %14.12e ||r(i)||/||b|| %14.12e\n", n, normtype, (double) rnorm, (double) truenorm, (double) (truenorm/bnorm));
388:   PetscViewerASCIISubtractTab(viewer, tablevel);
389:   PetscViewerPopFormat(viewer);
390:   return(0);
391: }

393: /*@C
394:   KSPMonitorTrueResidualDraw - Plots the true residual at each iteration of an iterative solver.

396:   Collective on ksp

398:   Input Parameters:
399: + ksp   - iterative context
400: . n     - iteration number
401: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
402: - vf    - The viewer context

404:   Options Database Key:
405: . -ksp_monitor_true_residual draw - Activates KSPMonitorResidualDraw()

407:   Level: intermediate

409: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
410: @*/
411: PetscErrorCode KSPMonitorTrueResidualDraw(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
412: {
413:   PetscViewer       viewer = vf->viewer;
414:   PetscViewerFormat format = vf->format;
415:   Vec               r;
416:   PetscErrorCode    ierr;

420:   PetscViewerPushFormat(viewer, format);
421:   KSPBuildResidual(ksp, NULL, NULL, &r);
422:   PetscObjectSetName((PetscObject) r, "Residual");
423:   PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", (PetscObject) ksp);
424:   VecView(r, viewer);
425:   PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", NULL);
426:   VecDestroy(&r);
427:   PetscViewerPopFormat(viewer);
428:   return(0);
429: }

431: /*@C
432:   KSPMonitorTrueResidualDrawLG - Plots the true residual norm at each iteration of an iterative solver.

434:   Collective on ksp

436:   Input Parameters:
437: + ksp   - iterative context
438: . n     - iteration number
439: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
440: - vf    - The viewer context

442:   Options Database Key:
443: . -ksp_monitor_true_residual draw::draw_lg - Activates KSPMonitorTrueResidualDrawLG()

445:   Level: intermediate

447: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
448: @*/
449: PetscErrorCode KSPMonitorTrueResidualDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
450: {
451:   PetscViewer        viewer = vf->viewer;
452:   PetscViewerFormat  format = vf->format;
453:   PetscDrawLG        lg     = vf->lg;
454:   Vec                r;
455:   KSPConvergedReason reason;
456:   PetscReal          truenorm, x[2], y[2];
457:   PetscErrorCode     ierr;

462:   KSPBuildResidual(ksp, NULL, NULL, &r);
463:   VecNorm(r, NORM_2, &truenorm);
464:   VecDestroy(&r);
465:   PetscViewerPushFormat(viewer, format);
466:   if (!n) {PetscDrawLGReset(lg);}
467:   x[0] = (PetscReal) n;
468:   if (rnorm > 0.0) y[0] = PetscLog10Real(rnorm);
469:   else y[0] = -15.0;
470:   x[1] = (PetscReal) n;
471:   if (truenorm > 0.0) y[1] = PetscLog10Real(truenorm);
472:   else y[1] = -15.0;
473:   PetscDrawLGAddPoint(lg, x, y);
474:   KSPGetConvergedReason(ksp, &reason);
475:   if (n <= 20 || !(n % 5) || reason) {
476:     PetscDrawLGDraw(lg);
477:     PetscDrawLGSave(lg);
478:   }
479:   PetscViewerPopFormat(viewer);
480:   return(0);
481: }

483: /*@C
484:   KSPMonitorTrueResidualDrawLGCreate - Creates the plotter for the preconditioned residual.

486:   Collective on ksp

488:   Input Parameters:
489: + viewer - The PetscViewer
490: . format - The viewer format
491: - ctx    - An optional user context

493:   Output Parameter:
494: . vf    - The viewer context

496:   Level: intermediate

498: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
499: @*/
500: PetscErrorCode KSPMonitorTrueResidualDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
501: {
502:   const char    *names[] = {"preconditioned", "true"};

506:   PetscViewerAndFormatCreate(viewer, format, vf);
507:   (*vf)->data = ctx;
508:   KSPMonitorLGCreate(PetscObjectComm((PetscObject) viewer), NULL, NULL, "Log Residual Norm", 2, names, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
509:   return(0);
510: }

512: /*@C
513:   KSPMonitorTrueResidualMax - Prints the true residual max norm at each iteration of an iterative solver.

515:   Collective on ksp

517:   Input Parameters:
518: + ksp   - iterative context
519: . n     - iteration number
520: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
521: - vf    - The viewer context

523:   Options Database Key:
524: . -ksp_monitor_true_residual_max - Activates KSPMonitorTrueResidualMax()

526:   Notes:
527:   When using right preconditioning, these values are equivalent.

529:   Level: intermediate

531: .seealso: KSPMonitorSet(), KSPMonitorResidual(),KSPMonitorTrueResidualMaxNorm()
532: @*/
533: PetscErrorCode KSPMonitorTrueResidualMax(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
534: {
535:   PetscViewer       viewer = vf->viewer;
536:   PetscViewerFormat format = vf->format;
537:   Vec               r;
538:   PetscReal         truenorm, bnorm;
539:   char              normtype[256];
540:   PetscInt          tablevel;
541:   const char       *prefix;
542:   PetscErrorCode    ierr;

546:   PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
547:   PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
548:   PetscStrncpy(normtype, KSPNormTypes[ksp->normtype], sizeof(normtype));
549:   PetscStrtolower(normtype);
550:   KSPBuildResidual(ksp, NULL, NULL, &r);
551:   VecNorm(r, NORM_INFINITY, &truenorm);
552:   VecNorm(ksp->vec_rhs, NORM_INFINITY, &bnorm);
553:   VecDestroy(&r);

555:   PetscViewerPushFormat(viewer, format);
556:   PetscViewerASCIIAddTab(viewer, tablevel);
557:   if (n == 0 && prefix) {PetscViewerASCIIPrintf(viewer, "  Residual norms for %s solve.\n", prefix);}
558:   PetscViewerASCIIPrintf(viewer, "%3D KSP %s true resid max norm %14.12e ||r(i)||/||b|| %14.12e\n", n, normtype, (double) truenorm, (double) (truenorm/bnorm));
559:   PetscViewerASCIISubtractTab(viewer, tablevel);
560:   PetscViewerPopFormat(viewer);
561:   return(0);
562: }

564: /*@C
565:   KSPMonitorError - Prints the error norm, as well as the preconditioned residual norm, at each iteration of an iterative solver.

567:   Collective on ksp

569:   Input Parameters:
570: + ksp   - iterative context
571: . n     - iteration number
572: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
573: - vf    - The viewer context

575:   Options Database Key:
576: . -ksp_monitor_error - Activates KSPMonitorError()

578:   Level: intermediate

580: .seealso: KSPMonitorSet(), KSPMonitorResidual(),KSPMonitorTrueResidualMaxNorm()
581: @*/
582: PetscErrorCode KSPMonitorError(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
583: {
584:   PetscViewer       viewer = vf->viewer;
585:   PetscViewerFormat format = vf->format;
586:   DM                dm;
587:   Vec               sol;
588:   PetscReal        *errors;
589:   PetscInt          Nf, f;
590:   PetscInt          tablevel;
591:   const char       *prefix;
592:   PetscErrorCode    ierr;

596:   PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
597:   PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
598:   KSPGetDM(ksp, &dm);
599:   DMGetNumFields(dm, &Nf);
600:   DMGetGlobalVector(dm, &sol);
601:   KSPBuildSolution(ksp, sol, NULL);
602:   /* TODO: Make a different monitor that flips sign for SNES, Newton system is A dx = -b, so we need to negate the solution */
603:   VecScale(sol, -1.0);
604:   PetscCalloc1(Nf, &errors);
605:   DMComputeError(dm, sol, errors, NULL);

607:   PetscViewerPushFormat(viewer, format);
608:   PetscViewerASCIIAddTab(viewer, tablevel);
609:   if (n == 0 && prefix) {PetscViewerASCIIPrintf(viewer, "  Error norms for %s solve.\n", prefix);}
610:   PetscViewerASCIIPrintf(viewer, "%3D KSP Error norm %s", n, Nf > 1 ? "[" : "");
611:   PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
612:   for (f = 0; f < Nf; ++f) {
613:     if (f > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
614:     PetscViewerASCIIPrintf(viewer, "%14.12e", (double) errors[f]);
615:   }
616:   PetscViewerASCIIPrintf(viewer, "%s resid norm %14.12e\n", Nf > 1 ? "]" : "", (double) rnorm);
617:   PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
618:   PetscViewerASCIISubtractTab(viewer, tablevel);
619:   PetscViewerPopFormat(viewer);
620:   DMRestoreGlobalVector(dm, &sol);
621:   return(0);
622: }

624: /*@C
625:   KSPMonitorErrorDraw - Plots the error at each iteration of an iterative solver.

627:   Collective on ksp

629:   Input Parameters:
630: + ksp   - iterative context
631: . n     - iteration number
632: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
633: - vf    - The viewer context

635:   Options Database Key:
636: . -ksp_monitor_error draw - Activates KSPMonitorErrorDraw()

638:   Level: intermediate

640: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
641: @*/
642: PetscErrorCode KSPMonitorErrorDraw(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
643: {
644:   PetscViewer       viewer = vf->viewer;
645:   PetscViewerFormat format = vf->format;
646:   DM                dm;
647:   Vec               sol, e;
648:   PetscErrorCode    ierr;

652:   PetscViewerPushFormat(viewer, format);
653:   KSPGetDM(ksp, &dm);
654:   DMGetGlobalVector(dm, &sol);
655:   KSPBuildSolution(ksp, sol, NULL);
656:   DMComputeError(dm, sol, NULL, &e);
657:   PetscObjectSetName((PetscObject) e, "Error");
658:   PetscObjectCompose((PetscObject) e, "__Vec_bc_zero__", (PetscObject) ksp);
659:   VecView(e, viewer);
660:   PetscObjectCompose((PetscObject) e, "__Vec_bc_zero__", NULL);
661:   VecDestroy(&e);
662:   DMRestoreGlobalVector(dm, &sol);
663:   PetscViewerPopFormat(viewer);
664:   return(0);
665: }

667: /*@C
668:   KSPMonitorErrorDrawLG - Plots the error and residual norm at each iteration of an iterative solver.

670:   Collective on ksp

672:   Input Parameters:
673: + ksp   - iterative context
674: . n     - iteration number
675: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
676: - vf    - The viewer context

678:   Options Database Key:
679: . -ksp_monitor_error draw::draw_lg - Activates KSPMonitorTrueResidualDrawLG()

681:   Level: intermediate

683: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
684: @*/
685: PetscErrorCode KSPMonitorErrorDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
686: {
687:   PetscViewer        viewer = vf->viewer;
688:   PetscViewerFormat  format = vf->format;
689:   PetscDrawLG        lg     = vf->lg;
690:   DM                 dm;
691:   Vec                sol;
692:   KSPConvergedReason reason;
693:   PetscReal         *x, *errors;
694:   PetscInt           Nf, f;
695:   PetscErrorCode     ierr;

700:   KSPGetDM(ksp, &dm);
701:   DMGetNumFields(dm, &Nf);
702:   DMGetGlobalVector(dm, &sol);
703:   KSPBuildSolution(ksp, sol, NULL);
704:   /* TODO: Make a different monitor that flips sign for SNES, Newton system is A dx = -b, so we need to negate the solution */
705:   VecScale(sol, -1.0);
706:   PetscCalloc2(Nf+1, &x, Nf+1, &errors);
707:   DMComputeError(dm, sol, errors, NULL);

709:   PetscViewerPushFormat(viewer, format);
710:   if (!n) {PetscDrawLGReset(lg);}
711:   for (f = 0; f < Nf; ++f) {
712:     x[f]      = (PetscReal) n;
713:     errors[f] = errors[f] > 0.0 ? PetscLog10Real(errors[f]) : -15.;
714:   }
715:   x[Nf]      = (PetscReal) n;
716:   errors[Nf] = rnorm > 0.0 ? PetscLog10Real(rnorm) : -15.;
717:   PetscDrawLGAddPoint(lg, x, errors);
718:   KSPGetConvergedReason(ksp, &reason);
719:   if (n <= 20 || !(n % 5) || reason) {
720:     PetscDrawLGDraw(lg);
721:     PetscDrawLGSave(lg);
722:   }
723:   PetscViewerPopFormat(viewer);
724:   return(0);
725: }

727: /*@C
728:   KSPMonitorErrorDrawLGCreate - Creates the plotter for the error and preconditioned residual.

730:   Collective on ksp

732:   Input Parameters:
733: + viewer - The PetscViewer
734: . format - The viewer format
735: - ctx    - An optional user context

737:   Output Parameter:
738: . vf    - The viewer context

740:   Level: intermediate

742: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
743: @*/
744: PetscErrorCode KSPMonitorErrorDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
745: {
746:   KSP            ksp = (KSP) ctx;
747:   DM             dm;
748:   char         **names;
749:   PetscInt       Nf, f;

753:   KSPGetDM(ksp, &dm);
754:   DMGetNumFields(dm, &Nf);
755:   PetscMalloc1(Nf+1, &names);
756:   for (f = 0; f < Nf; ++f) {
757:     PetscObject disc;
758:     const char *fname;
759:     char        lname[PETSC_MAX_PATH_LEN];

761:     DMGetField(dm, f, NULL, &disc);
762:     PetscObjectGetName(disc, &fname);
763:     PetscStrncpy(lname, fname, PETSC_MAX_PATH_LEN);
764:     PetscStrlcat(lname, " Error", PETSC_MAX_PATH_LEN);
765:     PetscStrallocpy(lname, &names[f]);
766:   }
767:   PetscStrallocpy("residual", &names[Nf]);
768:   PetscViewerAndFormatCreate(viewer, format, vf);
769:   (*vf)->data = ctx;
770:   KSPMonitorLGCreate(PetscObjectComm((PetscObject) viewer), NULL, NULL, "Log Error Norm", Nf+1, (const char **) names, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
771:   for (f = 0; f <= Nf; ++f) {PetscFree(names[f]);}
772:   PetscFree(names);
773:   return(0);
774: }

776: /*@C
777:   KSPMonitorSolution - Print the solution norm at each iteration of an iterative solver.

779:   Collective on ksp

781:   Input Parameters:
782: + ksp   - iterative context
783: . n     - iteration number
784: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
785: - vf    - The viewer context

787:   Options Database Key:
788: . -ksp_monitor_solution - Activates KSPMonitorSolution()

790:   Level: intermediate

792: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
793: @*/
794: PetscErrorCode KSPMonitorSolution(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
795: {
796:   PetscViewer       viewer = vf->viewer;
797:   PetscViewerFormat format = vf->format;
798:   Vec               x;
799:   PetscReal         snorm;
800:   PetscInt          tablevel;
801:   const char       *prefix;
802:   PetscErrorCode    ierr;

806:   KSPBuildSolution(ksp, NULL, &x);
807:   VecNorm(x, NORM_2, &snorm);
808:   PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
809:   PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
810:   PetscViewerPushFormat(viewer, format);
811:   PetscViewerASCIIAddTab(viewer, tablevel);
812:   if (n == 0 && prefix) {PetscViewerASCIIPrintf(viewer, "  Solution norms for %s solve.\n", prefix);}
813:   PetscViewerASCIIPrintf(viewer, "%3D KSP Solution norm %14.12e \n", n, (double) snorm);
814:   PetscViewerASCIISubtractTab(viewer, tablevel);
815:   PetscViewerPopFormat(viewer);
816:   return(0);
817: }

819: /*@C
820:   KSPMonitorSolutionDraw - Plots the solution at each iteration of an iterative solver.

822:   Collective on ksp

824:   Input Parameters:
825: + ksp   - iterative context
826: . n     - iteration number
827: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
828: - vf    - The viewer context

830:   Options Database Key:
831: . -ksp_monitor_solution draw - Activates KSPMonitorSolutionDraw()

833:   Level: intermediate

835: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
836: @*/
837: PetscErrorCode KSPMonitorSolutionDraw(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
838: {
839:   PetscViewer       viewer = vf->viewer;
840:   PetscViewerFormat format = vf->format;
841:   Vec               x;
842:   PetscErrorCode    ierr;

846:   KSPBuildSolution(ksp, NULL, &x);
847:   PetscViewerPushFormat(viewer, format);
848:   PetscObjectSetName((PetscObject) x, "Solution");
849:   PetscObjectCompose((PetscObject) x, "__Vec_bc_zero__", (PetscObject) ksp);
850:   VecView(x, viewer);
851:   PetscObjectCompose((PetscObject) x, "__Vec_bc_zero__", NULL);
852:   PetscViewerPopFormat(viewer);
853:   return(0);
854: }

856: /*@C
857:   KSPMonitorSolutionDrawLG - Plots the solution norm at each iteration of an iterative solver.

859:   Collective on ksp

861:   Input Parameters:
862: + ksp   - iterative context
863: . n     - iteration number
864: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
865: - vf    - The viewer context

867:   Options Database Key:
868: . -ksp_monitor_solution draw::draw_lg - Activates KSPMonitorSolutionDrawLG()

870:   Level: intermediate

872: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
873: @*/
874: PetscErrorCode KSPMonitorSolutionDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
875: {
876:   PetscViewer        viewer = vf->viewer;
877:   PetscViewerFormat  format = vf->format;
878:   PetscDrawLG        lg     = vf->lg;
879:   Vec                u;
880:   KSPConvergedReason reason;
881:   PetscReal          snorm, x, y;
882:   PetscErrorCode     ierr;

887:   KSPBuildSolution(ksp, NULL, &u);
888:   VecNorm(u, NORM_2, &snorm);
889:   PetscViewerPushFormat(viewer, format);
890:   if (!n) {PetscDrawLGReset(lg);}
891:   x = (PetscReal) n;
892:   if (snorm > 0.0) y = PetscLog10Real(snorm);
893:   else y = -15.0;
894:   PetscDrawLGAddPoint(lg, &x, &y);
895:   KSPGetConvergedReason(ksp, &reason);
896:   if (n <= 20 || !(n % 5) || reason) {
897:     PetscDrawLGDraw(lg);
898:     PetscDrawLGSave(lg);
899:   }
900:   PetscViewerPopFormat(viewer);
901:   return(0);
902: }

904: /*@C
905:   KSPMonitorSolutionDrawLGCreate - Creates the plotter for the solution.

907:   Collective on ksp

909:   Input Parameters:
910: + viewer - The PetscViewer
911: . format - The viewer format
912: - ctx    - An optional user context

914:   Output Parameter:
915: . vf    - The viewer context

917:   Level: intermediate

919: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
920: @*/
921: PetscErrorCode KSPMonitorSolutionDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
922: {

926:   PetscViewerAndFormatCreate(viewer, format, vf);
927:   (*vf)->data = ctx;
928:   KSPMonitorLGCreate(PetscObjectComm((PetscObject) viewer), NULL, NULL, "Log Solution Norm", 1, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
929:   return(0);
930: }

932: /*@C
933:   KSPMonitorSingularValue - Prints the two norm of the true residual and estimation of the extreme singular values of the preconditioned problem at each iteration.

935:   Logically Collective on ksp

937:   Input Parameters:
938: + ksp   - the iterative context
939: . n     - the iteration
940: . rnorm - the two norm of the residual
941: - vf    - The viewer context

943:   Options Database Key:
944: . -ksp_monitor_singular_value - Activates KSPMonitorSingularValue()

946:   Notes:
947:   The CG solver uses the Lanczos technique for eigenvalue computation,
948:   while GMRES uses the Arnoldi technique; other iterative methods do
949:   not currently compute singular values.

951:   Level: intermediate

953: .seealso: KSPComputeExtremeSingularValues()
954: @*/
955: PetscErrorCode KSPMonitorSingularValue(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
956: {
957:   PetscViewer       viewer = vf->viewer;
958:   PetscViewerFormat format = vf->format;
959:   PetscReal         emin, emax;
960:   PetscInt          tablevel;
961:   const char       *prefix;
962:   PetscErrorCode    ierr;

967:   PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
968:   PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
969:   PetscViewerPushFormat(viewer, format);
970:   PetscViewerASCIIAddTab(viewer, tablevel);
971:   if (n == 0 && prefix) {PetscViewerASCIIPrintf(viewer, "  Residual norms for %s solve.\n", prefix);}
972:   if (!ksp->calc_sings) {
973:     PetscViewerASCIIPrintf(viewer, "%3D KSP Residual norm %14.12e \n", n, (double) rnorm);
974:   } else {
975:     KSPComputeExtremeSingularValues(ksp, &emax, &emin);
976:     PetscViewerASCIIPrintf(viewer, "%3D KSP Residual norm %14.12e %% max %14.12e min %14.12e max/min %14.12e\n", n, (double) rnorm, (double) emax, (double) emin, (double) (emax/emin));
977:   }
978:   PetscViewerASCIISubtractTab(viewer, tablevel);
979:   PetscViewerPopFormat(viewer);
980:   return(0);
981: }

983: /*@C
984:   KSPMonitorSingularValueCreate - Creates the singular value monitor.

986:   Collective on ksp

988:   Input Parameters:
989: + viewer - The PetscViewer
990: . format - The viewer format
991: - ctx    - An optional user context

993:   Output Parameter:
994: . vf    - The viewer context

996:   Level: intermediate

998: .seealso: KSPMonitorSet(), KSPMonitorSingularValue()
999: @*/
1000: PetscErrorCode KSPMonitorSingularValueCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
1001: {
1002:   KSP            ksp = (KSP) ctx;

1006:   PetscViewerAndFormatCreate(viewer, format, vf);
1007:   (*vf)->data = ctx;
1008:   KSPSetComputeSingularValues(ksp, PETSC_TRUE);
1009:   return(0);
1010: }

1012: /*@C
1013:    KSPMonitorDynamicTolerance - Recompute the inner tolerance in every outer iteration in an adaptive way.

1015:    Collective on ksp

1017:    Input Parameters:
1018: +  ksp   - iterative context
1019: .  n     - iteration number (not used)
1020: .  fnorm - the current residual norm
1021: -  dummy - some context as a C struct. fields:
1022:              coef: a scaling coefficient. default 1.0. can be passed through
1023:                    -sub_ksp_dynamic_tolerance_param
1024:              bnrm: norm of the right-hand side. store it to avoid repeated calculation

1026:    Notes:
1027:    This may be useful for a flexibly preconditioner Krylov method to
1028:    control the accuracy of the inner solves needed to gaurantee the
1029:    convergence of the outer iterations.

1031:    Level: advanced

1033: .seealso: KSPMonitorDynamicToleranceDestroy()
1034: @*/
1035: PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
1036: {
1038:   PC             pc;
1039:   PetscReal      outer_rtol, outer_abstol, outer_dtol, inner_rtol;
1040:   PetscInt       outer_maxits,nksp,first,i;
1041:   KSPDynTolCtx   *scale   = (KSPDynTolCtx*)dummy;
1042:   KSP            *subksp = NULL;
1043:   KSP            kspinner;
1044:   PetscBool      flg;

1047:   KSPGetPC(ksp, &pc);

1049:   /* compute inner_rtol */
1050:   if (scale->bnrm < 0.0) {
1051:     Vec b;
1052:     KSPGetRhs(ksp, &b);
1053:     VecNorm(b, NORM_2, &(scale->bnrm));
1054:   }
1055:   KSPGetTolerances(ksp, &outer_rtol, &outer_abstol, &outer_dtol, &outer_maxits);
1056:   inner_rtol = PetscMin(scale->coef * scale->bnrm * outer_rtol / fnorm, 0.999);
1057:   /*PetscPrintf(PETSC_COMM_WORLD, "        Inner rtol = %g\n", (double)inner_rtol);*/

1059:   /* if pc is ksp */
1060:   PetscObjectTypeCompare((PetscObject)pc,PCKSP,&flg);
1061:   if (flg) {
1062:     PCKSPGetKSP(pc, &kspinner);
1063:     KSPSetTolerances(kspinner, inner_rtol, outer_abstol, outer_dtol, outer_maxits);
1064:     return(0);
1065:   }

1067:   /* if pc is bjacobi */
1068:   PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);
1069:   if (flg) {
1070:     PCBJacobiGetSubKSP(pc, &nksp, &first, &subksp);
1071:     if (subksp) {
1072:       for (i=0; i<nksp; i++) {
1073:         KSPSetTolerances(subksp[i], inner_rtol, outer_abstol, outer_dtol, outer_maxits);
1074:       }
1075:       return(0);
1076:     }
1077:   }

1079:   /* if pc is deflation*/
1080:   PetscObjectTypeCompare((PetscObject)pc,PCDEFLATION,&flg);
1081:   if (flg) {
1082:     PCDeflationGetCoarseKSP(pc,&kspinner);
1083:     KSPSetTolerances(kspinner,inner_rtol,outer_abstol,outer_dtol,PETSC_DEFAULT);
1084:     return(0);
1085:   }

1087:   /* todo: dynamic tolerance may apply to other types of pc too */
1088:   return(0);
1089: }

1091: /*
1092:   Destroy the dummy context used in KSPMonitorDynamicTolerance()
1093: */
1094: PetscErrorCode KSPMonitorDynamicToleranceDestroy(void **dummy)
1095: {

1099:   PetscFree(*dummy);
1100:   return(0);
1101: }

1103: /*@C
1104:    KSPConvergedSkip - Convergence test that do not return as converged
1105:    until the maximum number of iterations is reached.

1107:    Collective on ksp

1109:    Input Parameters:
1110: +  ksp   - iterative context
1111: .  n     - iteration number
1112: .  rnorm - 2-norm residual value (may be estimated)
1113: -  dummy - unused convergence context

1115:    Returns:
1116: .  reason - KSP_CONVERGED_ITERATING, KSP_CONVERGED_ITS

1118:    Notes:
1119:    This should be used as the convergence test with the option
1120:    KSPSetNormType(ksp,KSP_NORM_NONE), since norms of the residual are
1121:    not computed. Convergence is then declared after the maximum number
1122:    of iterations have been reached. Useful when one is using CG or
1123:    BiCGStab as a smoother.

1125:    Level: advanced

1127: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
1128: @*/
1129: PetscErrorCode  KSPConvergedSkip(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
1130: {
1134:   *reason = KSP_CONVERGED_ITERATING;
1135:   if (n >= ksp->max_it) *reason = KSP_CONVERGED_ITS;
1136:   return(0);
1137: }


1140: /*@C
1141:    KSPConvergedDefaultCreate - Creates and initializes the space used by the KSPConvergedDefault() function context

1143:    Note Collective

1145:    Output Parameter:
1146: .  ctx - convergence context

1148:    Level: intermediate

1150: .seealso: KSPConvergedDefault(), KSPConvergedDefaultDestroy(), KSPSetConvergenceTest(), KSPSetTolerances(),
1151:           KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultSetConvergedMaxits()
1152: @*/
1153: PetscErrorCode  KSPConvergedDefaultCreate(void **ctx)
1154: {
1155:   PetscErrorCode         ierr;
1156:   KSPConvergedDefaultCtx *cctx;

1159:   PetscNew(&cctx);
1160:   *ctx = cctx;
1161:   return(0);
1162: }

1164: /*@
1165:    KSPConvergedDefaultSetUIRNorm - makes the default convergence test use || B*(b - A*(initial guess))||
1166:       instead of || B*b ||. In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
1167:       is used there is no B in the above formula. UIRNorm is short for Use Initial Residual Norm.

1169:    Collective on ksp

1171:    Input Parameters:
1172: .  ksp   - iterative context

1174:    Options Database:
1175: .   -ksp_converged_use_initial_residual_norm

1177:    Notes:
1178:    Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.

1180:    The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
1181:    are defined in petscksp.h.

1183:    If the convergence test is not KSPConvergedDefault() then this is ignored.

1185:    If right preconditioning is being used then B does not appear in the above formula.


1188:    Level: intermediate

1190: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultSetConvergedMaxits()
1191: @*/
1192: PetscErrorCode  KSPConvergedDefaultSetUIRNorm(KSP ksp)
1193: {
1194:   KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;

1198:   if (ksp->converged != KSPConvergedDefault) return(0);
1199:   if (ctx->mininitialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
1200:   ctx->initialrtol = PETSC_TRUE;
1201:   return(0);
1202: }

1204: /*@
1205:    KSPConvergedDefaultSetUMIRNorm - makes the default convergence test use min(|| B*(b - A*(initial guess))||,|| B*b ||)
1206:       In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
1207:       is used there is no B in the above formula. UMIRNorm is short for Use Minimum Initial Residual Norm.

1209:    Collective on ksp

1211:    Input Parameters:
1212: .  ksp   - iterative context

1214:    Options Database:
1215: .   -ksp_converged_use_min_initial_residual_norm

1217:    Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.

1219:    The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
1220:    are defined in petscksp.h.

1222:    Level: intermediate

1224: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetConvergedMaxits()
1225: @*/
1226: PetscErrorCode  KSPConvergedDefaultSetUMIRNorm(KSP ksp)
1227: {
1228:   KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;

1232:   if (ksp->converged != KSPConvergedDefault) return(0);
1233:   if (ctx->initialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
1234:   ctx->mininitialrtol = PETSC_TRUE;
1235:   return(0);
1236: }

1238: /*@
1239:    KSPConvergedDefaultSetConvergedMaxits - allows the default convergence test to declare convergence and return KSP_CONVERGED_ITS if the maximum number of iterations is reached

1241:    Collective on ksp

1243:    Input Parameters:
1244: +  ksp - iterative context
1245: -  flg - boolean flag

1247:    Options Database:
1248: .   -ksp_converged_maxits

1250:    Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.

1252:    The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
1253:    are defined in petscksp.h.

1255:    Level: intermediate

1257: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultSetUIRNorm()
1258: @*/
1259: PetscErrorCode  KSPConvergedDefaultSetConvergedMaxits(KSP ksp, PetscBool flg)
1260: {
1261:   KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;

1266:   if (ksp->converged != KSPConvergedDefault) return(0);
1267:   ctx->convmaxits = flg;
1268:   return(0);
1269: }

1271: /*@C
1272:    KSPConvergedDefault - Determines convergence of the linear iterative solvers by default

1274:    Collective on ksp

1276:    Input Parameters:
1277: +  ksp   - iterative context
1278: .  n     - iteration number
1279: .  rnorm - residual norm (may be estimated, depending on the method may be the preconditioned residual norm)
1280: -  ctx - convergence context which must be created by KSPConvergedDefaultCreate()

1282:    Output Parameter:
1283: +   positive - if the iteration has converged
1284: .   negative - if the iteration has diverged
1285: -   KSP_CONVERGED_ITERATING - otherwise.

1287:    Notes:
1288:    KSPConvergedDefault() reaches convergence when   rnorm < MAX (rtol * rnorm_0, abstol);
1289:    Divergence is detected if rnorm > dtol * rnorm_0, or when failures are detected throughout the iteration.
1290:    By default, reaching the maximum number of iterations is considered divergence (i.e. KSP_DIVERGED_ITS).
1291:    In order to have PETSc declaring convergence in such a case (i.e. KSP_CONVERGED_ITS), users can use KSPConvergedDefaultSetConvergedMaxits()

1293:    where:
1294: +     rtol = relative tolerance,
1295: .     abstol = absolute tolerance.
1296: .     dtol = divergence tolerance,
1297: -     rnorm_0 is the two norm of the right hand side (or the preconditioned norm, depending on what was set with
1298:           KSPSetNormType(). When initial guess is non-zero you
1299:           can call KSPConvergedDefaultSetUIRNorm() to use the norm of (b - A*(initial guess))
1300:           as the starting point for relative norm convergence testing, that is as rnorm_0

1302:    Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.

1304:    Use KSPSetNormType() (or -ksp_norm_type <none,preconditioned,unpreconditioned,natural>) to change the norm used for computing rnorm

1306:    The precise values of reason are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.

1308:    This routine is used by KSP by default so the user generally never needs call it directly.

1310:    Use KSPSetConvergenceTest() to provide your own test instead of using this one.

1312:    Level: intermediate

1314: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
1315:           KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultSetConvergedMaxits(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy()
1316: @*/
1317: PetscErrorCode  KSPConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
1318: {
1319:   PetscErrorCode         ierr;
1320:   KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
1321:   KSPNormType            normtype;

1327:   if (!cctx) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_NULL,"Convergence context must have been created with KSPConvergedDefaultCreate()");
1328:   *reason = KSP_CONVERGED_ITERATING;

1330:   if (cctx->convmaxits && n >= ksp->max_it) {
1331:     *reason = KSP_CONVERGED_ITS;
1332:     PetscInfo1(ksp,"Linear solver has converged. Maximum number of iterations reached %D\n",n);
1333:     return(0);
1334:   }
1335:   KSPGetNormType(ksp,&normtype);
1336:   if (normtype == KSP_NORM_NONE) return(0);

1338:   if (!n) {
1339:     /* if user gives initial guess need to compute norm of b */
1340:     if (!ksp->guess_zero && !cctx->initialrtol) {
1341:       PetscReal snorm = 0.0;
1342:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
1343:         PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");
1344:         VecNorm(ksp->vec_rhs,NORM_2,&snorm);        /*     <- b'*b */
1345:       } else {
1346:         Vec z;
1347:         /* Should avoid allocating the z vector each time but cannot stash it in cctx because if KSPReset() is called the vector size might change */
1348:         VecDuplicate(ksp->vec_rhs,&z);
1349:         KSP_PCApply(ksp,ksp->vec_rhs,z);
1350:         if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
1351:           PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");
1352:           VecNorm(z,NORM_2,&snorm);                 /*    dp <- b'*B'*B*b */
1353:         } else if (ksp->normtype == KSP_NORM_NATURAL) {
1354:           PetscScalar norm;
1355:           PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");
1356:           VecDot(ksp->vec_rhs,z,&norm);
1357:           snorm = PetscSqrtReal(PetscAbsScalar(norm));                            /*    dp <- b'*B*b */
1358:         }
1359:         VecDestroy(&z);
1360:       }
1361:       /* handle special case of zero RHS and nonzero guess */
1362:       if (!snorm) {
1363:         PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");
1364:         snorm = rnorm;
1365:       }
1366:       if (cctx->mininitialrtol) ksp->rnorm0 = PetscMin(snorm,rnorm);
1367:       else ksp->rnorm0 = snorm;
1368:     } else {
1369:       ksp->rnorm0 = rnorm;
1370:     }
1371:     ksp->ttol = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
1372:   }

1374:   if (n <= ksp->chknorm) return(0);

1376:   if (PetscIsInfOrNanReal(rnorm)) {
1377:     PCFailedReason pcreason;
1378:     PetscInt       sendbuf,recvbuf;
1379:     PCGetFailedReasonRank(ksp->pc,&pcreason);
1380:     sendbuf = (PetscInt)pcreason;
1381:     MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPIU_MAX,PetscObjectComm((PetscObject)ksp));
1382:     if (recvbuf) {
1383:       *reason = KSP_DIVERGED_PC_FAILED;
1384:       PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf);
1385:       PetscInfo(ksp,"Linear solver pcsetup fails, declaring divergence \n");
1386:     } else {
1387:       *reason = KSP_DIVERGED_NANORINF;
1388:       PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");
1389:     }
1390:   } else if (rnorm <= ksp->ttol) {
1391:     if (rnorm < ksp->abstol) {
1392:       PetscInfo3(ksp,"Linear solver has converged. Residual norm %14.12e is less than absolute tolerance %14.12e at iteration %D\n",(double)rnorm,(double)ksp->abstol,n);
1393:       *reason = KSP_CONVERGED_ATOL;
1394:     } else {
1395:       if (cctx->initialrtol) {
1396:         PetscInfo4(ksp,"Linear solver has converged. Residual norm %14.12e is less than relative tolerance %14.12e times initial residual norm %14.12e at iteration %D\n",(double)rnorm,(double)ksp->rtol,(double)ksp->rnorm0,n);
1397:       } else {
1398:         PetscInfo4(ksp,"Linear solver has converged. Residual norm %14.12e is less than relative tolerance %14.12e times initial right hand side norm %14.12e at iteration %D\n",(double)rnorm,(double)ksp->rtol,(double)ksp->rnorm0,n);
1399:       }
1400:       *reason = KSP_CONVERGED_RTOL;
1401:     }
1402:   } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
1403:     PetscInfo3(ksp,"Linear solver is diverging. Initial right hand size norm %14.12e, current residual norm %14.12e at iteration %D\n",(double)ksp->rnorm0,(double)rnorm,n);
1404:     *reason = KSP_DIVERGED_DTOL;
1405:   }
1406:   return(0);
1407: }

1409: /*@C
1410:    KSPConvergedDefaultDestroy - Frees the space used by the KSPConvergedDefault() function context

1412:    Not Collective

1414:    Input Parameters:
1415: .  ctx - convergence context

1417:    Level: intermediate

1419: .seealso: KSPConvergedDefault(), KSPConvergedDefaultCreate(), KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(),
1420:           KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
1421: @*/
1422: PetscErrorCode  KSPConvergedDefaultDestroy(void *ctx)
1423: {
1424:   PetscErrorCode         ierr;
1425:   KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;

1428:   VecDestroy(&cctx->work);
1429:   PetscFree(ctx);
1430:   return(0);
1431: }

1433: /*
1434:    KSPBuildSolutionDefault - Default code to create/move the solution.

1436:    Collective on ksp

1438:    Input Parameters:
1439: +  ksp - iterative context
1440: -  v   - pointer to the user's vector

1442:    Output Parameter:
1443: .  V - pointer to a vector containing the solution

1445:    Level: advanced

1447:    Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations

1449: .seealso: KSPGetSolution(), KSPBuildResidualDefault()
1450: */
1451: PetscErrorCode KSPBuildSolutionDefault(KSP ksp,Vec v,Vec *V)
1452: {

1456:   if (ksp->pc_side == PC_RIGHT) {
1457:     if (ksp->pc) {
1458:       if (v) {
1459:         KSP_PCApply(ksp,ksp->vec_sol,v); *V = v;
1460:       } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with right preconditioner");
1461:     } else {
1462:       if (v) {
1463:         VecCopy(ksp->vec_sol,v); *V = v;
1464:       } else *V = ksp->vec_sol;
1465:     }
1466:   } else if (ksp->pc_side == PC_SYMMETRIC) {
1467:     if (ksp->pc) {
1468:       if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
1469:       if (v) {
1470:         PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v);
1471:         *V = v;
1472:       } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner");
1473:     } else {
1474:       if (v) {
1475:         VecCopy(ksp->vec_sol,v); *V = v;
1476:       } else *V = ksp->vec_sol;
1477:     }
1478:   } else {
1479:     if (v) {
1480:       VecCopy(ksp->vec_sol,v); *V = v;
1481:     } else *V = ksp->vec_sol;
1482:   }
1483:   return(0);
1484: }

1486: /*
1487:    KSPBuildResidualDefault - Default code to compute the residual.

1489:    Collecive on ksp

1491:    Input Parameters:
1492: .  ksp - iterative context
1493: .  t   - pointer to temporary vector
1494: .  v   - pointer to user vector

1496:    Output Parameter:
1497: .  V - pointer to a vector containing the residual

1499:    Level: advanced

1501:    Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations

1503: .seealso: KSPBuildSolutionDefault()
1504: */
1505: PetscErrorCode KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec *V)
1506: {
1508:   Mat            Amat,Pmat;

1511:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
1512:   PCGetOperators(ksp->pc,&Amat,&Pmat);
1513:   KSPBuildSolution(ksp,t,NULL);
1514:   KSP_MatMult(ksp,Amat,t,v);
1515:   VecAYPX(v,-1.0,ksp->vec_rhs);
1516:   *V   = v;
1517:   return(0);
1518: }

1520: /*@C
1521:   KSPCreateVecs - Gets a number of work vectors.

1523:   Collective on ksp

1525:   Input Parameters:
1526: + ksp  - iterative context
1527: . rightn  - number of right work vectors
1528: - leftn   - number of left work vectors to allocate

1530:   Output Parameter:
1531: +  right - the array of vectors created
1532: -  left - the array of left vectors

1534:    Note: The right vector has as many elements as the matrix has columns. The left
1535:      vector has as many elements as the matrix has rows.

1537:    The vectors are new vectors that are not owned by the KSP, they should be destroyed with calls to VecDestroyVecs() when no longer needed.

1539:    Developers Note: First tries to duplicate the rhs and solution vectors of the KSP, if they do not exist tries to get them from the matrix, if
1540:                     that does not exist tries to get them from the DM (if it is provided).

1542:    Level: advanced

1544: .seealso:   MatCreateVecs(), VecDestroyVecs()

1546: @*/
1547: PetscErrorCode KSPCreateVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
1548: {
1550:   Vec            vecr = NULL,vecl = NULL;
1551:   PetscBool      matset,pmatset,isshell,preferdm = PETSC_FALSE;
1552:   Mat            mat = NULL;

1555:   if (ksp->dm) {
1556:     PetscObjectTypeCompare((PetscObject) ksp->dm, DMSHELL, &isshell);
1557:     preferdm = isshell ? PETSC_FALSE : PETSC_TRUE;
1558:   }
1559:   if (rightn) {
1560:     if (!right) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for right vectors but did not pass a pointer to hold them");
1561:     if (ksp->vec_sol) vecr = ksp->vec_sol;
1562:     else {
1563:       if (preferdm) {
1564:         DMGetGlobalVector(ksp->dm,&vecr);
1565:       } else if (ksp->pc) {
1566:         PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
1567:         /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
1568:         if (matset) {
1569:           PCGetOperators(ksp->pc,&mat,NULL);
1570:           MatCreateVecs(mat,&vecr,NULL);
1571:         } else if (pmatset) {
1572:           PCGetOperators(ksp->pc,NULL,&mat);
1573:           MatCreateVecs(mat,&vecr,NULL);
1574:         }
1575:       }
1576:       if (!vecr && ksp->dm) {
1577:         DMGetGlobalVector(ksp->dm,&vecr);
1578:       }
1579:       if (!vecr) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
1580:     }
1581:     VecDuplicateVecs(vecr,rightn,right);
1582:     if (!ksp->vec_sol) {
1583:       if (preferdm) {
1584:         DMRestoreGlobalVector(ksp->dm,&vecr);
1585:       } else if (mat) {
1586:         VecDestroy(&vecr);
1587:       } else if (ksp->dm) {
1588:         DMRestoreGlobalVector(ksp->dm,&vecr);
1589:       }
1590:     }
1591:   }
1592:   if (leftn) {
1593:     if (!left) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for left vectors but did not pass a pointer to hold them");
1594:     if (ksp->vec_rhs) vecl = ksp->vec_rhs;
1595:     else {
1596:       if (preferdm) {
1597:         DMGetGlobalVector(ksp->dm,&vecl);
1598:       } else if (ksp->pc) {
1599:         PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
1600:         /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
1601:         if (matset) {
1602:           PCGetOperators(ksp->pc,&mat,NULL);
1603:           MatCreateVecs(mat,NULL,&vecl);
1604:         } else if (pmatset) {
1605:           PCGetOperators(ksp->pc,NULL,&mat);
1606:           MatCreateVecs(mat,NULL,&vecl);
1607:         }
1608:       }
1609:       if (!vecl && ksp->dm) {
1610:         DMGetGlobalVector(ksp->dm,&vecl);
1611:       }
1612:       if (!vecl) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
1613:     }
1614:     VecDuplicateVecs(vecl,leftn,left);
1615:     if (!ksp->vec_rhs) {
1616:       if (preferdm) {
1617:         DMRestoreGlobalVector(ksp->dm,&vecl);
1618:       } else if (mat) {
1619:         VecDestroy(&vecl);
1620:       } else if (ksp->dm) {
1621:         DMRestoreGlobalVector(ksp->dm,&vecl);
1622:       }
1623:     }
1624:   }
1625:   return(0);
1626: }

1628: /*@C
1629:   KSPSetWorkVecs - Sets a number of work vectors into a KSP object

1631:   Collective on ksp

1633:   Input Parameters:
1634: + ksp  - iterative context
1635: - nw   - number of work vectors to allocate

1637:   Level: developer

1639:   Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations

1641: @*/
1642: PetscErrorCode KSPSetWorkVecs(KSP ksp,PetscInt nw)
1643: {

1647:   VecDestroyVecs(ksp->nwork,&ksp->work);
1648:   ksp->nwork = nw;
1649:   KSPCreateVecs(ksp,nw,&ksp->work,0,NULL);
1650:   PetscLogObjectParents(ksp,nw,ksp->work);
1651:   return(0);
1652: }

1654: /*
1655:   KSPDestroyDefault - Destroys a iterative context variable for methods with
1656:   no separate context.  Preferred calling sequence KSPDestroy().

1658:   Input Parameter:
1659: . ksp - the iterative context

1661:    Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations

1663: */
1664: PetscErrorCode KSPDestroyDefault(KSP ksp)
1665: {

1670:   PetscFree(ksp->data);
1671:   return(0);
1672: }

1674: /*@
1675:    KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.

1677:    Not Collective

1679:    Input Parameter:
1680: .  ksp - the KSP context

1682:    Output Parameter:
1683: .  reason - negative value indicates diverged, positive value converged, see KSPConvergedReason

1685:    Possible values for reason: See also manual page for each reason
1686: $  KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
1687: $  KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
1688: $  KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration, or when the KSPConvergedSkip() convergence test routine is set.
1689: $  KSP_CONVERGED_CG_NEG_CURVE (see note below)
1690: $  KSP_CONVERGED_CG_CONSTRAINED (see note below)
1691: $  KSP_CONVERGED_STEP_LENGTH (see note below)
1692: $  KSP_CONVERGED_ITERATING (returned if the solver is not yet finished)
1693: $  KSP_DIVERGED_ITS  (required more than its to reach convergence)
1694: $  KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
1695: $  KSP_DIVERGED_NANORINF (residual norm became Not-a-number or Inf likely due to 0/0)
1696: $  KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
1697: $  KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial residual. Try a different preconditioner, or a different initial Level.)

1699:    Options Database:
1700: .   -ksp_converged_reason - prints the reason to standard out

1702:    Notes:
1703:     If this routine is called before or doing the KSPSolve() the value of KSP_CONVERGED_ITERATING is returned

1705:    The values  KSP_CONVERGED_CG_NEG_CURVE, KSP_CONVERGED_CG_CONSTRAINED, and KSP_CONVERGED_STEP_LENGTH are returned only by the special KSPNASH, KSPSTCG, and KSPGLTR
1706:    solvers which are used by the SNESNEWTONTR (trust region) solver.

1708:    Level: intermediate

1710: .seealso: KSPSetConvergenceTest(), KSPConvergedDefault(), KSPSetTolerances(), KSPConvergedReason,
1711:           KSPConvergedReasonView()
1712: @*/
1713: PetscErrorCode  KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
1714: {
1718:   *reason = ksp->reason;
1719:   return(0);
1720: }

1722: /*@C
1723:    KSPGetConvergedReasonString - Return a human readable string for ksp converged reason

1725:    Not Collective

1727:    Input Parameter:
1728: .  ksp - the KSP context

1730:    Output Parameter:
1731: .  strreason - a human readable string that describes ksp converged reason

1733:    Level: beginner

1735: .seealso: KSPGetConvergedReason()
1736: @*/
1737: PetscErrorCode KSPGetConvergedReasonString(KSP ksp,const char** strreason)
1738: {
1742:   *strreason = KSPConvergedReasons[ksp->reason];
1743:   return(0);
1744: }

1746: #include <petsc/private/dmimpl.h>
1747: /*@
1748:    KSPSetDM - Sets the DM that may be used by some preconditioners

1750:    Logically Collective on ksp

1752:    Input Parameters:
1753: +  ksp - the preconditioner context
1754: -  dm - the dm, cannot be NULL

1756:    Notes:
1757:    If this is used then the KSP will attempt to use the DM to create the matrix and use the routine set with
1758:    DMKSPSetComputeOperators(). Use KSPSetDMActive(ksp,PETSC_FALSE) to instead use the matrix you've provided with
1759:    KSPSetOperators().

1761:    A DM can only be used for solving one problem at a time because information about the problem is stored on the DM,
1762:    even when not using interfaces like DMKSPSetComputeOperators().  Use DMClone() to get a distinct DM when solving
1763:    different problems using the same function space.

1765:    Level: intermediate

1767: .seealso: KSPGetDM(), KSPSetDMActive(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess(), DMKSPSetComputeOperators(), DMKSPSetComputeRHS(), DMKSPSetComputeInitialGuess()
1768: @*/
1769: PetscErrorCode  KSPSetDM(KSP ksp,DM dm)
1770: {
1772:   PC             pc;

1777:   PetscObjectReference((PetscObject)dm);
1778:   if (ksp->dm) {                /* Move the DMSNES context over to the new DM unless the new DM already has one */
1779:     if (ksp->dm->dmksp && !dm->dmksp) {
1780:       DMKSP kdm;
1781:       DMCopyDMKSP(ksp->dm,dm);
1782:       DMGetDMKSP(ksp->dm,&kdm);
1783:       if (kdm->originaldm == ksp->dm) kdm->originaldm = dm; /* Grant write privileges to the replacement DM */
1784:     }
1785:     DMDestroy(&ksp->dm);
1786:   }
1787:   ksp->dm       = dm;
1788:   ksp->dmAuto   = PETSC_FALSE;
1789:   KSPGetPC(ksp,&pc);
1790:   PCSetDM(pc,dm);
1791:   ksp->dmActive = PETSC_TRUE;
1792:   return(0);
1793: }

1795: /*@
1796:    KSPSetDMActive - Indicates the DM should be used to generate the linear system matrix and right hand side

1798:    Logically Collective on ksp

1800:    Input Parameters:
1801: +  ksp - the preconditioner context
1802: -  flg - use the DM

1804:    Notes:
1805:    By default KSPSetDM() sets the DM as active, call KSPSetDMActive(ksp,PETSC_FALSE); after KSPSetDM(ksp,dm) to not have the KSP object use the DM to generate the matrices.

1807:    Level: intermediate

1809: .seealso: KSPGetDM(), KSPSetDM(), SNESSetDM(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess()
1810: @*/
1811: PetscErrorCode  KSPSetDMActive(KSP ksp,PetscBool flg)
1812: {
1816:   ksp->dmActive = flg;
1817:   return(0);
1818: }

1820: /*@
1821:    KSPGetDM - Gets the DM that may be used by some preconditioners

1823:    Not Collective

1825:    Input Parameter:
1826: . ksp - the preconditioner context

1828:    Output Parameter:
1829: .  dm - the dm

1831:    Level: intermediate


1834: .seealso: KSPSetDM(), KSPSetDMActive()
1835: @*/
1836: PetscErrorCode  KSPGetDM(KSP ksp,DM *dm)
1837: {

1842:   if (!ksp->dm) {
1843:     DMShellCreate(PetscObjectComm((PetscObject)ksp),&ksp->dm);
1844:     ksp->dmAuto = PETSC_TRUE;
1845:   }
1846:   *dm = ksp->dm;
1847:   return(0);
1848: }

1850: /*@
1851:    KSPSetApplicationContext - Sets the optional user-defined context for the linear solver.

1853:    Logically Collective on ksp

1855:    Input Parameters:
1856: +  ksp - the KSP context
1857: -  usrP - optional user context

1859:    Fortran Notes:
1860:     To use this from Fortran you must write a Fortran interface definition for this
1861:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

1863:    Level: intermediate

1865: .seealso: KSPGetApplicationContext()
1866: @*/
1867: PetscErrorCode  KSPSetApplicationContext(KSP ksp,void *usrP)
1868: {
1870:   PC             pc;

1874:   ksp->user = usrP;
1875:   KSPGetPC(ksp,&pc);
1876:   PCSetApplicationContext(pc,usrP);
1877:   return(0);
1878: }

1880: /*@
1881:    KSPGetApplicationContext - Gets the user-defined context for the linear solver.

1883:    Not Collective

1885:    Input Parameter:
1886: .  ksp - KSP context

1888:    Output Parameter:
1889: .  usrP - user context

1891:    Fortran Notes:
1892:     To use this from Fortran you must write a Fortran interface definition for this
1893:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

1895:    Level: intermediate

1897: .seealso: KSPSetApplicationContext()
1898: @*/
1899: PetscErrorCode  KSPGetApplicationContext(KSP ksp,void *usrP)
1900: {
1903:   *(void**)usrP = ksp->user;
1904:   return(0);
1905: }

1907: #include <petsc/private/pcimpl.h>

1909: /*@
1910:    KSPCheckSolve - Checks if the PCSetUp() or KSPSolve() failed and set the error flag for the outer PC. A KSP_DIVERGED_ITS is
1911:          not considered a failure in this context

1913:    Collective on ksp

1915:    Input Parameter:
1916: +  ksp - the linear solver (KSP) context.
1917: .  pc - the preconditioner context
1918: -  vec - a vector that will be initialized with Inf to indicate lack of convergence

1920:    Notes: this may be called by a subset of the processes in the PC

1922:    Level: developer

1924:    Developer Note: this is used to manage returning from preconditioners whose inner KSP solvers have failed in some way

1926: .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckNorm(), KSPCheckDot()
1927: @*/
1928: PetscErrorCode KSPCheckSolve(KSP ksp,PC pc,Vec vec)
1929: {
1930:   PetscErrorCode     ierr;
1931:   PCFailedReason     pcreason;
1932:   PC                 subpc;

1935:   KSPGetPC(ksp,&subpc);
1936:   PCGetFailedReason(subpc,&pcreason);
1937:   if (pcreason || (ksp->reason < 0 && ksp->reason != KSP_DIVERGED_ITS)) {
1938:     if (pc->erroriffailure) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_NOT_CONVERGED,"Detected not converged in KSP inner solve: KSP reason %s PC reason %s",KSPConvergedReasons[ksp->reason],PCFailedReasons[pcreason]);
1939:     else {
1940:       PetscInfo2(ksp,"Detected not converged in KSP inner solve: KSP reason %s PC reason %s\n",KSPConvergedReasons[ksp->reason],PCFailedReasons[pcreason]);
1941:       pc->failedreason = PC_SUBPC_ERROR;
1942:       if (vec) {
1943:         VecSetInf(vec);
1944:       }
1945:     }
1946:   }
1947:   return(0);
1948: }