Actual source code: iterativ.c

  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: }

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

1142:    Note Collective

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

1147:    Level: intermediate

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

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

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

1168:    Collective on ksp

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

1173:    Options Database:
1174: .   -ksp_converged_use_initial_residual_norm

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

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

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

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

1186:    Level: intermediate

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

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

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

1207:    Collective on ksp

1209:    Input Parameters:
1210: .  ksp   - iterative context

1212:    Options Database:
1213: .   -ksp_converged_use_min_initial_residual_norm

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

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

1220:    Level: intermediate

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

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

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

1239:    Collective on ksp

1241:    Input Parameters:
1242: +  ksp - iterative context
1243: -  flg - boolean flag

1245:    Options Database:
1246: .   -ksp_converged_maxits

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

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

1253:    Level: intermediate

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

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

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

1272:    Collective on ksp

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

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

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

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

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

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

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

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

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

1310:    Level: intermediate

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

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

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

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

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

1374:   if (PetscIsInfOrNanReal(rnorm)) {
1375:     PCFailedReason pcreason;
1376:     PetscInt       sendbuf,recvbuf;
1377:     PCGetFailedReasonRank(ksp->pc,&pcreason);
1378:     sendbuf = (PetscInt)pcreason;
1379:     MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPIU_MAX,PetscObjectComm((PetscObject)ksp));
1380:     if (recvbuf) {
1381:       *reason = KSP_DIVERGED_PC_FAILED;
1382:       PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf);
1383:       PetscInfo(ksp,"Linear solver pcsetup fails, declaring divergence \n");
1384:     } else {
1385:       *reason = KSP_DIVERGED_NANORINF;
1386:       PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");
1387:     }
1388:   } else if (rnorm <= ksp->ttol) {
1389:     if (rnorm < ksp->abstol) {
1390:       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);
1391:       *reason = KSP_CONVERGED_ATOL;
1392:     } else {
1393:       if (cctx->initialrtol) {
1394:         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);
1395:       } else {
1396:         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);
1397:       }
1398:       *reason = KSP_CONVERGED_RTOL;
1399:     }
1400:   } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
1401:     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);
1402:     *reason = KSP_DIVERGED_DTOL;
1403:   }
1404:   return(0);
1405: }

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

1410:    Not Collective

1412:    Input Parameters:
1413: .  ctx - convergence context

1415:    Level: intermediate

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

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

1431: /*
1432:    KSPBuildSolutionDefault - Default code to create/move the solution.

1434:    Collective on ksp

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

1440:    Output Parameter:
1441: .  V - pointer to a vector containing the solution

1443:    Level: advanced

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

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

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

1484: /*
1485:    KSPBuildResidualDefault - Default code to compute the residual.

1487:    Collecive on ksp

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

1494:    Output Parameter:
1495: .  V - pointer to a vector containing the residual

1497:    Level: advanced

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

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

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

1518: /*@C
1519:   KSPCreateVecs - Gets a number of work vectors.

1521:   Collective on ksp

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

1528:   Output Parameter:
1529: +  right - the array of vectors created
1530: -  left - the array of left vectors

1532:    Note: The right vector has as many elements as the matrix has columns. The left
1533:      vector has as many elements as the matrix has rows.

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

1537:    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
1538:                     that does not exist tries to get them from the DM (if it is provided).

1540:    Level: advanced

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

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

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

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

1629:   Collective on ksp

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

1635:   Level: developer

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

1639: @*/
1640: PetscErrorCode KSPSetWorkVecs(KSP ksp,PetscInt nw)
1641: {

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

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

1656:   Input Parameter:
1657: . ksp - the iterative context

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

1661: */
1662: PetscErrorCode KSPDestroyDefault(KSP ksp)
1663: {

1668:   PetscFree(ksp->data);
1669:   return(0);
1670: }

1672: /*@
1673:    KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.

1675:    Not Collective

1677:    Input Parameter:
1678: .  ksp - the KSP context

1680:    Output Parameter:
1681: .  reason - negative value indicates diverged, positive value converged, see KSPConvergedReason

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

1697:    Options Database:
1698: .   -ksp_converged_reason - prints the reason to standard out

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

1703:    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
1704:    solvers which are used by the SNESNEWTONTR (trust region) solver.

1706:    Level: intermediate

1708: .seealso: KSPSetConvergenceTest(), KSPConvergedDefault(), KSPSetTolerances(), KSPConvergedReason,
1709:           KSPConvergedReasonView()
1710: @*/
1711: PetscErrorCode  KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
1712: {
1716:   *reason = ksp->reason;
1717:   return(0);
1718: }

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

1723:    Not Collective

1725:    Input Parameter:
1726: .  ksp - the KSP context

1728:    Output Parameter:
1729: .  strreason - a human readable string that describes ksp converged reason

1731:    Level: beginner

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

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

1748:    Logically Collective on ksp

1750:    Input Parameters:
1751: +  ksp - the preconditioner context
1752: -  dm - the dm, cannot be NULL

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

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

1763:    Level: intermediate

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

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

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

1796:    Logically Collective on ksp

1798:    Input Parameters:
1799: +  ksp - the preconditioner context
1800: -  flg - use the DM

1802:    Notes:
1803:    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.

1805:    Level: intermediate

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

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

1821:    Not Collective

1823:    Input Parameter:
1824: . ksp - the preconditioner context

1826:    Output Parameter:
1827: .  dm - the dm

1829:    Level: intermediate

1831: .seealso: KSPSetDM(), KSPSetDMActive()
1832: @*/
1833: PetscErrorCode  KSPGetDM(KSP ksp,DM *dm)
1834: {

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

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

1850:    Logically Collective on ksp

1852:    Input Parameters:
1853: +  ksp - the KSP context
1854: -  usrP - optional user context

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

1860:    Level: intermediate

1862: .seealso: KSPGetApplicationContext()
1863: @*/
1864: PetscErrorCode  KSPSetApplicationContext(KSP ksp,void *usrP)
1865: {
1867:   PC             pc;

1871:   ksp->user = usrP;
1872:   KSPGetPC(ksp,&pc);
1873:   PCSetApplicationContext(pc,usrP);
1874:   return(0);
1875: }

1877: /*@
1878:    KSPGetApplicationContext - Gets the user-defined context for the linear solver.

1880:    Not Collective

1882:    Input Parameter:
1883: .  ksp - KSP context

1885:    Output Parameter:
1886: .  usrP - user context

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

1892:    Level: intermediate

1894: .seealso: KSPSetApplicationContext()
1895: @*/
1896: PetscErrorCode  KSPGetApplicationContext(KSP ksp,void *usrP)
1897: {
1900:   *(void**)usrP = ksp->user;
1901:   return(0);
1902: }

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

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

1910:    Collective on ksp

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

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

1919:    Level: developer

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

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

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