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: {
32: *rnorm = ksp->rnorm;
33: return 0;
34: }
36: /*@
37: KSPGetIterationNumber - Gets the current iteration number; if the
38: KSPSolve() is complete, returns the number of iterations
39: used.
41: Not Collective
43: Input Parameters:
44: . ksp - the iterative context
46: Output Parameters:
47: . its - number of iterations
49: Level: intermediate
51: Notes:
52: During the ith iteration this returns i-1
53: .seealso: KSPBuildResidual(), KSPGetResidualNorm(), KSPGetTotalIterations()
54: @*/
55: PetscErrorCode KSPGetIterationNumber(KSP ksp,PetscInt *its)
56: {
59: *its = ksp->its;
60: return 0;
61: }
63: /*@
64: KSPGetTotalIterations - Gets the total number of iterations this KSP object has performed since was created, counted over all linear solves
66: Not Collective
68: Input Parameters:
69: . ksp - the iterative context
71: Output Parameters:
72: . its - total number of iterations
74: Level: intermediate
76: Notes:
77: Use KSPGetIterationNumber() to get the count for the most recent solve only
78: If this is called within a linear solve (such as in a KSPMonitor routine) then it does not include iterations within that current solve
80: .seealso: KSPBuildResidual(), KSPGetResidualNorm(), KSPGetIterationNumber()
81: @*/
82: PetscErrorCode KSPGetTotalIterations(KSP ksp,PetscInt *its)
83: {
86: *its = ksp->totalits;
87: return 0;
88: }
90: /*@C
91: KSPMonitorResidual - Print the preconditioned residual norm at each iteration of an iterative solver.
93: Collective on ksp
95: Input Parameters:
96: + ksp - iterative context
97: . n - iteration number
98: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
99: - vf - The viewer context
101: Options Database Key:
102: . -ksp_monitor - Activates KSPMonitorResidual()
104: Level: intermediate
106: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
107: @*/
108: PetscErrorCode KSPMonitorResidual(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
109: {
110: PetscViewer viewer = vf->viewer;
111: PetscViewerFormat format = vf->format;
112: PetscInt tablevel;
113: const char *prefix;
116: PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
117: PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
118: PetscViewerPushFormat(viewer, format);
119: PetscViewerASCIIAddTab(viewer, tablevel);
120: if (n == 0 && prefix) PetscViewerASCIIPrintf(viewer, " Residual norms for %s solve.\n", prefix);
121: PetscViewerASCIIPrintf(viewer, "%3D KSP Residual norm %14.12e \n", n, (double) rnorm);
122: PetscViewerASCIISubtractTab(viewer, tablevel);
123: PetscViewerPopFormat(viewer);
124: return 0;
125: }
127: /*@C
128: KSPMonitorResidualDraw - Plots the preconditioned residual at each iteration of an iterative solver.
130: Collective on ksp
132: Input Parameters:
133: + ksp - iterative context
134: . n - iteration number
135: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
136: - vf - The viewer context
138: Options Database Key:
139: . -ksp_monitor draw - Activates KSPMonitorResidualDraw()
141: Level: intermediate
143: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
144: @*/
145: PetscErrorCode KSPMonitorResidualDraw(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
146: {
147: PetscViewer viewer = vf->viewer;
148: PetscViewerFormat format = vf->format;
149: Vec r;
152: PetscViewerPushFormat(viewer, format);
153: KSPBuildResidual(ksp, NULL, NULL, &r);
154: PetscObjectSetName((PetscObject) r, "Residual");
155: PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", (PetscObject) ksp);
156: VecView(r, viewer);
157: PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", NULL);
158: VecDestroy(&r);
159: PetscViewerPopFormat(viewer);
160: return 0;
161: }
163: /*@C
164: KSPMonitorResidualDrawLG - Plots the preconditioned residual norm at each iteration of an iterative solver.
166: Collective on ksp
168: Input Parameters:
169: + ksp - iterative context
170: . n - iteration number
171: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
172: - vf - The viewer context
174: Options Database Key:
175: . -ksp_monitor draw::draw_lg - Activates KSPMonitorResidualDrawLG()
177: Level: intermediate
179: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
180: @*/
181: PetscErrorCode KSPMonitorResidualDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
182: {
183: PetscViewer viewer = vf->viewer;
184: PetscViewerFormat format = vf->format;
185: PetscDrawLG lg = vf->lg;
186: KSPConvergedReason reason;
187: PetscReal x, y;
191: PetscViewerPushFormat(viewer, format);
192: if (!n) PetscDrawLGReset(lg);
193: x = (PetscReal) n;
194: if (rnorm > 0.0) y = PetscLog10Real(rnorm);
195: else y = -15.0;
196: PetscDrawLGAddPoint(lg, &x, &y);
197: KSPGetConvergedReason(ksp, &reason);
198: if (n <= 20 || !(n % 5) || reason) {
199: PetscDrawLGDraw(lg);
200: PetscDrawLGSave(lg);
201: }
202: PetscViewerPopFormat(viewer);
203: return 0;
204: }
206: /*@C
207: KSPMonitorResidualDrawLGCreate - Creates the plotter for the preconditioned residual.
209: Collective on ksp
211: Input Parameters:
212: + viewer - The PetscViewer
213: . format - The viewer format
214: - ctx - An optional user context
216: Output Parameter:
217: . vf - The viewer context
219: Level: intermediate
221: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
222: @*/
223: PetscErrorCode KSPMonitorResidualDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
224: {
225: PetscViewerAndFormatCreate(viewer, format, vf);
226: (*vf)->data = ctx;
227: KSPMonitorLGCreate(PetscObjectComm((PetscObject) viewer), NULL, NULL, "Log Residual Norm", 1, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
228: return 0;
229: }
231: /*
232: This is the same as KSPMonitorResidual() except it prints fewer digits of the residual as the residual gets smaller.
233: This is because the later digits are meaningless and are often different on different machines; by using this routine different
234: machines will usually generate the same output.
236: Deprecated: Intentionally has no manual page
237: */
238: PetscErrorCode KSPMonitorResidualShort(KSP ksp, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf)
239: {
240: PetscViewer viewer = vf->viewer;
241: PetscViewerFormat format = vf->format;
242: PetscInt tablevel;
243: const char *prefix;
246: PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
247: PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
248: PetscViewerPushFormat(viewer, format);
249: PetscViewerASCIIAddTab(viewer, tablevel);
250: if (its == 0 && prefix) PetscViewerASCIIPrintf(viewer, " Residual norms for %s solve.\n", prefix);
251: if (fnorm > 1.e-9) PetscViewerASCIIPrintf(viewer, "%3D KSP Residual norm %g \n", its, (double) fnorm);
252: else if (fnorm > 1.e-11) PetscViewerASCIIPrintf(viewer, "%3D KSP Residual norm %5.3e \n", its, (double) fnorm);
253: else PetscViewerASCIIPrintf(viewer, "%3D KSP Residual norm < 1.e-11\n", its);
254: PetscViewerASCIISubtractTab(viewer, tablevel);
255: PetscViewerPopFormat(viewer);
256: return 0;
257: }
259: PetscErrorCode KSPMonitorRange_Private(KSP ksp, PetscInt it, PetscReal *per)
260: {
261: Vec resid;
262: const PetscScalar *r;
263: PetscReal rmax, pwork;
264: PetscInt i, n, N;
266: KSPBuildResidual(ksp, NULL, NULL, &resid);
267: VecNorm(resid, NORM_INFINITY, &rmax);
268: VecGetLocalSize(resid, &n);
269: VecGetSize(resid, &N);
270: VecGetArrayRead(resid, &r);
271: pwork = 0.0;
272: for (i = 0; i < n; ++i) pwork += (PetscAbsScalar(r[i]) > .20*rmax);
273: VecRestoreArrayRead(resid, &r);
274: VecDestroy(&resid);
275: MPIU_Allreduce(&pwork, per, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) ksp));
276: *per = *per/N;
277: return 0;
278: }
280: /*@C
281: KSPMonitorResidualRange - Prints the percentage of residual elements that are more then 10 percent of the maximum value.
283: Collective on ksp
285: Input Parameters:
286: + ksp - iterative context
287: . it - iteration number
288: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
289: - vf - The viewer context
291: Options Database Key:
292: . -ksp_monitor_range - Activates KSPMonitorResidualRange()
294: Level: intermediate
296: .seealso: KSPMonitorSet(), KSPMonitorResidual()
297: @*/
298: PetscErrorCode KSPMonitorResidualRange(KSP ksp, PetscInt it, PetscReal rnorm, PetscViewerAndFormat *vf)
299: {
300: static PetscReal prev;
301: PetscViewer viewer = vf->viewer;
302: PetscViewerFormat format = vf->format;
303: PetscInt tablevel;
304: const char *prefix;
305: PetscReal perc, rel;
308: PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
309: PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
310: PetscViewerPushFormat(viewer, format);
311: PetscViewerASCIIAddTab(viewer, tablevel);
312: if (!it) prev = rnorm;
313: if (it == 0 && prefix) PetscViewerASCIIPrintf(viewer, " Residual norms for %s solve.\n", prefix);
314: KSPMonitorRange_Private(ksp, it, &perc);
315: rel = (prev - rnorm)/prev;
316: prev = rnorm;
317: 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));
318: PetscViewerASCIISubtractTab(viewer, tablevel);
319: PetscViewerPopFormat(viewer);
320: return 0;
321: }
323: /*@C
324: KSPMonitorTrueResidual - Prints the true residual norm, as well as the preconditioned residual norm, at each iteration of an iterative solver.
326: Collective on ksp
328: Input Parameters:
329: + ksp - iterative context
330: . n - iteration number
331: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
332: - vf - The viewer context
334: Options Database Key:
335: . -ksp_monitor_true_residual - Activates KSPMonitorTrueResidual()
337: Notes:
338: When using right preconditioning, these values are equivalent.
340: Level: intermediate
342: .seealso: KSPMonitorSet(), KSPMonitorResidual(),KSPMonitorTrueResidualMaxNorm()
343: @*/
344: PetscErrorCode KSPMonitorTrueResidual(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
345: {
346: PetscViewer viewer = vf->viewer;
347: PetscViewerFormat format = vf->format;
348: Vec r;
349: PetscReal truenorm, bnorm;
350: char normtype[256];
351: PetscInt tablevel;
352: const char *prefix;
355: PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
356: PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
357: PetscStrncpy(normtype, KSPNormTypes[ksp->normtype], sizeof(normtype));
358: PetscStrtolower(normtype);
359: KSPBuildResidual(ksp, NULL, NULL, &r);
360: VecNorm(r, NORM_2, &truenorm);
361: VecNorm(ksp->vec_rhs, NORM_2, &bnorm);
362: VecDestroy(&r);
364: PetscViewerPushFormat(viewer, format);
365: PetscViewerASCIIAddTab(viewer, tablevel);
366: if (n == 0 && prefix) PetscViewerASCIIPrintf(viewer, " Residual norms for %s solve.\n", prefix);
367: 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));
368: PetscViewerASCIISubtractTab(viewer, tablevel);
369: PetscViewerPopFormat(viewer);
370: return 0;
371: }
373: /*@C
374: KSPMonitorTrueResidualDraw - Plots the true residual at each iteration of an iterative solver.
376: Collective on ksp
378: Input Parameters:
379: + ksp - iterative context
380: . n - iteration number
381: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
382: - vf - The viewer context
384: Options Database Key:
385: . -ksp_monitor_true_residual draw - Activates KSPMonitorResidualDraw()
387: Level: intermediate
389: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
390: @*/
391: PetscErrorCode KSPMonitorTrueResidualDraw(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
392: {
393: PetscViewer viewer = vf->viewer;
394: PetscViewerFormat format = vf->format;
395: Vec r;
398: PetscViewerPushFormat(viewer, format);
399: KSPBuildResidual(ksp, NULL, NULL, &r);
400: PetscObjectSetName((PetscObject) r, "Residual");
401: PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", (PetscObject) ksp);
402: VecView(r, viewer);
403: PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", NULL);
404: VecDestroy(&r);
405: PetscViewerPopFormat(viewer);
406: return 0;
407: }
409: /*@C
410: KSPMonitorTrueResidualDrawLG - Plots the true residual norm at each iteration of an iterative solver.
412: Collective on ksp
414: Input Parameters:
415: + ksp - iterative context
416: . n - iteration number
417: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
418: - vf - The viewer context
420: Options Database Key:
421: . -ksp_monitor_true_residual draw::draw_lg - Activates KSPMonitorTrueResidualDrawLG()
423: Level: intermediate
425: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
426: @*/
427: PetscErrorCode KSPMonitorTrueResidualDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
428: {
429: PetscViewer viewer = vf->viewer;
430: PetscViewerFormat format = vf->format;
431: PetscDrawLG lg = vf->lg;
432: Vec r;
433: KSPConvergedReason reason;
434: PetscReal truenorm, x[2], y[2];
438: KSPBuildResidual(ksp, NULL, NULL, &r);
439: VecNorm(r, NORM_2, &truenorm);
440: VecDestroy(&r);
441: PetscViewerPushFormat(viewer, format);
442: if (!n) PetscDrawLGReset(lg);
443: x[0] = (PetscReal) n;
444: if (rnorm > 0.0) y[0] = PetscLog10Real(rnorm);
445: else y[0] = -15.0;
446: x[1] = (PetscReal) n;
447: if (truenorm > 0.0) y[1] = PetscLog10Real(truenorm);
448: else y[1] = -15.0;
449: PetscDrawLGAddPoint(lg, x, y);
450: KSPGetConvergedReason(ksp, &reason);
451: if (n <= 20 || !(n % 5) || reason) {
452: PetscDrawLGDraw(lg);
453: PetscDrawLGSave(lg);
454: }
455: PetscViewerPopFormat(viewer);
456: return 0;
457: }
459: /*@C
460: KSPMonitorTrueResidualDrawLGCreate - Creates the plotter for the preconditioned residual.
462: Collective on ksp
464: Input Parameters:
465: + viewer - The PetscViewer
466: . format - The viewer format
467: - ctx - An optional user context
469: Output Parameter:
470: . vf - The viewer context
472: Level: intermediate
474: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
475: @*/
476: PetscErrorCode KSPMonitorTrueResidualDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
477: {
478: const char *names[] = {"preconditioned", "true"};
480: PetscViewerAndFormatCreate(viewer, format, vf);
481: (*vf)->data = ctx;
482: KSPMonitorLGCreate(PetscObjectComm((PetscObject) viewer), NULL, NULL, "Log Residual Norm", 2, names, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
483: return 0;
484: }
486: /*@C
487: KSPMonitorTrueResidualMax - Prints the true residual max norm at each iteration of an iterative solver.
489: Collective on ksp
491: Input Parameters:
492: + ksp - iterative context
493: . n - iteration number
494: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
495: - vf - The viewer context
497: Options Database Key:
498: . -ksp_monitor_true_residual_max - Activates KSPMonitorTrueResidualMax()
500: Notes:
501: When using right preconditioning, these values are equivalent.
503: Level: intermediate
505: .seealso: KSPMonitorSet(), KSPMonitorResidual(),KSPMonitorTrueResidualMaxNorm()
506: @*/
507: PetscErrorCode KSPMonitorTrueResidualMax(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
508: {
509: PetscViewer viewer = vf->viewer;
510: PetscViewerFormat format = vf->format;
511: Vec r;
512: PetscReal truenorm, bnorm;
513: char normtype[256];
514: PetscInt tablevel;
515: const char *prefix;
518: PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
519: PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
520: PetscStrncpy(normtype, KSPNormTypes[ksp->normtype], sizeof(normtype));
521: PetscStrtolower(normtype);
522: KSPBuildResidual(ksp, NULL, NULL, &r);
523: VecNorm(r, NORM_INFINITY, &truenorm);
524: VecNorm(ksp->vec_rhs, NORM_INFINITY, &bnorm);
525: VecDestroy(&r);
527: PetscViewerPushFormat(viewer, format);
528: PetscViewerASCIIAddTab(viewer, tablevel);
529: if (n == 0 && prefix) PetscViewerASCIIPrintf(viewer, " Residual norms for %s solve.\n", prefix);
530: PetscViewerASCIIPrintf(viewer, "%3D KSP %s true resid max norm %14.12e ||r(i)||/||b|| %14.12e\n", n, normtype, (double) truenorm, (double) (truenorm/bnorm));
531: PetscViewerASCIISubtractTab(viewer, tablevel);
532: PetscViewerPopFormat(viewer);
533: return 0;
534: }
536: /*@C
537: KSPMonitorError - Prints the error norm, as well as the preconditioned residual norm, at each iteration of an iterative solver.
539: Collective on ksp
541: Input Parameters:
542: + ksp - iterative context
543: . n - iteration number
544: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
545: - vf - The viewer context
547: Options Database Key:
548: . -ksp_monitor_error - Activates KSPMonitorError()
550: Level: intermediate
552: .seealso: KSPMonitorSet(), KSPMonitorResidual(),KSPMonitorTrueResidualMaxNorm()
553: @*/
554: PetscErrorCode KSPMonitorError(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
555: {
556: PetscViewer viewer = vf->viewer;
557: PetscViewerFormat format = vf->format;
558: DM dm;
559: Vec sol;
560: PetscReal *errors;
561: PetscInt Nf, f;
562: PetscInt tablevel;
563: const char *prefix;
566: PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
567: PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
568: KSPGetDM(ksp, &dm);
569: DMGetNumFields(dm, &Nf);
570: DMGetGlobalVector(dm, &sol);
571: KSPBuildSolution(ksp, sol, NULL);
572: /* TODO: Make a different monitor that flips sign for SNES, Newton system is A dx = -b, so we need to negate the solution */
573: VecScale(sol, -1.0);
574: PetscCalloc1(Nf, &errors);
575: DMComputeError(dm, sol, errors, NULL);
577: PetscViewerPushFormat(viewer, format);
578: PetscViewerASCIIAddTab(viewer, tablevel);
579: if (n == 0 && prefix) PetscViewerASCIIPrintf(viewer, " Error norms for %s solve.\n", prefix);
580: PetscViewerASCIIPrintf(viewer, "%3D KSP Error norm %s", n, Nf > 1 ? "[" : "");
581: PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
582: for (f = 0; f < Nf; ++f) {
583: if (f > 0) PetscViewerASCIIPrintf(viewer, ", ");
584: PetscViewerASCIIPrintf(viewer, "%14.12e", (double) errors[f]);
585: }
586: PetscViewerASCIIPrintf(viewer, "%s resid norm %14.12e\n", Nf > 1 ? "]" : "", (double) rnorm);
587: PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
588: PetscViewerASCIISubtractTab(viewer, tablevel);
589: PetscViewerPopFormat(viewer);
590: DMRestoreGlobalVector(dm, &sol);
591: return 0;
592: }
594: /*@C
595: KSPMonitorErrorDraw - Plots the error at each iteration of an iterative solver.
597: Collective on ksp
599: Input Parameters:
600: + ksp - iterative context
601: . n - iteration number
602: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
603: - vf - The viewer context
605: Options Database Key:
606: . -ksp_monitor_error draw - Activates KSPMonitorErrorDraw()
608: Level: intermediate
610: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
611: @*/
612: PetscErrorCode KSPMonitorErrorDraw(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
613: {
614: PetscViewer viewer = vf->viewer;
615: PetscViewerFormat format = vf->format;
616: DM dm;
617: Vec sol, e;
620: PetscViewerPushFormat(viewer, format);
621: KSPGetDM(ksp, &dm);
622: DMGetGlobalVector(dm, &sol);
623: KSPBuildSolution(ksp, sol, NULL);
624: DMComputeError(dm, sol, NULL, &e);
625: PetscObjectSetName((PetscObject) e, "Error");
626: PetscObjectCompose((PetscObject) e, "__Vec_bc_zero__", (PetscObject) ksp);
627: VecView(e, viewer);
628: PetscObjectCompose((PetscObject) e, "__Vec_bc_zero__", NULL);
629: VecDestroy(&e);
630: DMRestoreGlobalVector(dm, &sol);
631: PetscViewerPopFormat(viewer);
632: return 0;
633: }
635: /*@C
636: KSPMonitorErrorDrawLG - Plots the error and residual norm at each iteration of an iterative solver.
638: Collective on ksp
640: Input Parameters:
641: + ksp - iterative context
642: . n - iteration number
643: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
644: - vf - The viewer context
646: Options Database Key:
647: . -ksp_monitor_error draw::draw_lg - Activates KSPMonitorTrueResidualDrawLG()
649: Level: intermediate
651: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
652: @*/
653: PetscErrorCode KSPMonitorErrorDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
654: {
655: PetscViewer viewer = vf->viewer;
656: PetscViewerFormat format = vf->format;
657: PetscDrawLG lg = vf->lg;
658: DM dm;
659: Vec sol;
660: KSPConvergedReason reason;
661: PetscReal *x, *errors;
662: PetscInt Nf, f;
666: KSPGetDM(ksp, &dm);
667: DMGetNumFields(dm, &Nf);
668: DMGetGlobalVector(dm, &sol);
669: KSPBuildSolution(ksp, sol, NULL);
670: /* TODO: Make a different monitor that flips sign for SNES, Newton system is A dx = -b, so we need to negate the solution */
671: VecScale(sol, -1.0);
672: PetscCalloc2(Nf+1, &x, Nf+1, &errors);
673: DMComputeError(dm, sol, errors, NULL);
675: PetscViewerPushFormat(viewer, format);
676: if (!n) PetscDrawLGReset(lg);
677: for (f = 0; f < Nf; ++f) {
678: x[f] = (PetscReal) n;
679: errors[f] = errors[f] > 0.0 ? PetscLog10Real(errors[f]) : -15.;
680: }
681: x[Nf] = (PetscReal) n;
682: errors[Nf] = rnorm > 0.0 ? PetscLog10Real(rnorm) : -15.;
683: PetscDrawLGAddPoint(lg, x, errors);
684: KSPGetConvergedReason(ksp, &reason);
685: if (n <= 20 || !(n % 5) || reason) {
686: PetscDrawLGDraw(lg);
687: PetscDrawLGSave(lg);
688: }
689: PetscViewerPopFormat(viewer);
690: return 0;
691: }
693: /*@C
694: KSPMonitorErrorDrawLGCreate - Creates the plotter for the error and preconditioned residual.
696: Collective on ksp
698: Input Parameters:
699: + viewer - The PetscViewer
700: . format - The viewer format
701: - ctx - An optional user context
703: Output Parameter:
704: . vf - The viewer context
706: Level: intermediate
708: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
709: @*/
710: PetscErrorCode KSPMonitorErrorDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
711: {
712: KSP ksp = (KSP) ctx;
713: DM dm;
714: char **names;
715: PetscInt Nf, f;
717: KSPGetDM(ksp, &dm);
718: DMGetNumFields(dm, &Nf);
719: PetscMalloc1(Nf+1, &names);
720: for (f = 0; f < Nf; ++f) {
721: PetscObject disc;
722: const char *fname;
723: char lname[PETSC_MAX_PATH_LEN];
725: DMGetField(dm, f, NULL, &disc);
726: PetscObjectGetName(disc, &fname);
727: PetscStrncpy(lname, fname, PETSC_MAX_PATH_LEN);
728: PetscStrlcat(lname, " Error", PETSC_MAX_PATH_LEN);
729: PetscStrallocpy(lname, &names[f]);
730: }
731: PetscStrallocpy("residual", &names[Nf]);
732: PetscViewerAndFormatCreate(viewer, format, vf);
733: (*vf)->data = ctx;
734: KSPMonitorLGCreate(PetscObjectComm((PetscObject) viewer), NULL, NULL, "Log Error Norm", Nf+1, (const char **) names, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
735: for (f = 0; f <= Nf; ++f) PetscFree(names[f]);
736: PetscFree(names);
737: return 0;
738: }
740: /*@C
741: KSPMonitorSolution - Print the solution norm at each iteration of an iterative solver.
743: Collective on ksp
745: Input Parameters:
746: + ksp - iterative context
747: . n - iteration number
748: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
749: - vf - The viewer context
751: Options Database Key:
752: . -ksp_monitor_solution - Activates KSPMonitorSolution()
754: Level: intermediate
756: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
757: @*/
758: PetscErrorCode KSPMonitorSolution(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
759: {
760: PetscViewer viewer = vf->viewer;
761: PetscViewerFormat format = vf->format;
762: Vec x;
763: PetscReal snorm;
764: PetscInt tablevel;
765: const char *prefix;
768: KSPBuildSolution(ksp, NULL, &x);
769: VecNorm(x, NORM_2, &snorm);
770: PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
771: PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
772: PetscViewerPushFormat(viewer, format);
773: PetscViewerASCIIAddTab(viewer, tablevel);
774: if (n == 0 && prefix) PetscViewerASCIIPrintf(viewer, " Solution norms for %s solve.\n", prefix);
775: PetscViewerASCIIPrintf(viewer, "%3D KSP Solution norm %14.12e \n", n, (double) snorm);
776: PetscViewerASCIISubtractTab(viewer, tablevel);
777: PetscViewerPopFormat(viewer);
778: return 0;
779: }
781: /*@C
782: KSPMonitorSolutionDraw - Plots the solution at each iteration of an iterative solver.
784: Collective on ksp
786: Input Parameters:
787: + ksp - iterative context
788: . n - iteration number
789: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
790: - vf - The viewer context
792: Options Database Key:
793: . -ksp_monitor_solution draw - Activates KSPMonitorSolutionDraw()
795: Level: intermediate
797: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
798: @*/
799: PetscErrorCode KSPMonitorSolutionDraw(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
800: {
801: PetscViewer viewer = vf->viewer;
802: PetscViewerFormat format = vf->format;
803: Vec x;
806: KSPBuildSolution(ksp, NULL, &x);
807: PetscViewerPushFormat(viewer, format);
808: PetscObjectSetName((PetscObject) x, "Solution");
809: PetscObjectCompose((PetscObject) x, "__Vec_bc_zero__", (PetscObject) ksp);
810: VecView(x, viewer);
811: PetscObjectCompose((PetscObject) x, "__Vec_bc_zero__", NULL);
812: PetscViewerPopFormat(viewer);
813: return 0;
814: }
816: /*@C
817: KSPMonitorSolutionDrawLG - Plots the solution norm at each iteration of an iterative solver.
819: Collective on ksp
821: Input Parameters:
822: + ksp - iterative context
823: . n - iteration number
824: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
825: - vf - The viewer context
827: Options Database Key:
828: . -ksp_monitor_solution draw::draw_lg - Activates KSPMonitorSolutionDrawLG()
830: Level: intermediate
832: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
833: @*/
834: PetscErrorCode KSPMonitorSolutionDrawLG(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
835: {
836: PetscViewer viewer = vf->viewer;
837: PetscViewerFormat format = vf->format;
838: PetscDrawLG lg = vf->lg;
839: Vec u;
840: KSPConvergedReason reason;
841: PetscReal snorm, x, y;
845: KSPBuildSolution(ksp, NULL, &u);
846: VecNorm(u, NORM_2, &snorm);
847: PetscViewerPushFormat(viewer, format);
848: if (!n) PetscDrawLGReset(lg);
849: x = (PetscReal) n;
850: if (snorm > 0.0) y = PetscLog10Real(snorm);
851: else y = -15.0;
852: PetscDrawLGAddPoint(lg, &x, &y);
853: KSPGetConvergedReason(ksp, &reason);
854: if (n <= 20 || !(n % 5) || reason) {
855: PetscDrawLGDraw(lg);
856: PetscDrawLGSave(lg);
857: }
858: PetscViewerPopFormat(viewer);
859: return 0;
860: }
862: /*@C
863: KSPMonitorSolutionDrawLGCreate - Creates the plotter for the solution.
865: Collective on ksp
867: Input Parameters:
868: + viewer - The PetscViewer
869: . format - The viewer format
870: - ctx - An optional user context
872: Output Parameter:
873: . vf - The viewer context
875: Level: intermediate
877: .seealso: KSPMonitorSet(), KSPMonitorTrueResidual()
878: @*/
879: PetscErrorCode KSPMonitorSolutionDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
880: {
881: PetscViewerAndFormatCreate(viewer, format, vf);
882: (*vf)->data = ctx;
883: KSPMonitorLGCreate(PetscObjectComm((PetscObject) viewer), NULL, NULL, "Log Solution Norm", 1, NULL, PETSC_DECIDE, PETSC_DECIDE, 400, 300, &(*vf)->lg);
884: return 0;
885: }
887: /*@C
888: KSPMonitorSingularValue - Prints the two norm of the true residual and estimation of the extreme singular values of the preconditioned problem at each iteration.
890: Logically Collective on ksp
892: Input Parameters:
893: + ksp - the iterative context
894: . n - the iteration
895: . rnorm - the two norm of the residual
896: - vf - The viewer context
898: Options Database Key:
899: . -ksp_monitor_singular_value - Activates KSPMonitorSingularValue()
901: Notes:
902: The CG solver uses the Lanczos technique for eigenvalue computation,
903: while GMRES uses the Arnoldi technique; other iterative methods do
904: not currently compute singular values.
906: Level: intermediate
908: .seealso: KSPComputeExtremeSingularValues()
909: @*/
910: PetscErrorCode KSPMonitorSingularValue(KSP ksp, PetscInt n, PetscReal rnorm, PetscViewerAndFormat *vf)
911: {
912: PetscViewer viewer = vf->viewer;
913: PetscViewerFormat format = vf->format;
914: PetscReal emin, emax;
915: PetscInt tablevel;
916: const char *prefix;
920: PetscObjectGetTabLevel((PetscObject) ksp, &tablevel);
921: PetscObjectGetOptionsPrefix((PetscObject) ksp, &prefix);
922: PetscViewerPushFormat(viewer, format);
923: PetscViewerASCIIAddTab(viewer, tablevel);
924: if (n == 0 && prefix) PetscViewerASCIIPrintf(viewer, " Residual norms for %s solve.\n", prefix);
925: if (!ksp->calc_sings) {
926: PetscViewerASCIIPrintf(viewer, "%3D KSP Residual norm %14.12e \n", n, (double) rnorm);
927: } else {
928: KSPComputeExtremeSingularValues(ksp, &emax, &emin);
929: 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));
930: }
931: PetscViewerASCIISubtractTab(viewer, tablevel);
932: PetscViewerPopFormat(viewer);
933: return 0;
934: }
936: /*@C
937: KSPMonitorSingularValueCreate - Creates the singular value monitor.
939: Collective on ksp
941: Input Parameters:
942: + viewer - The PetscViewer
943: . format - The viewer format
944: - ctx - An optional user context
946: Output Parameter:
947: . vf - The viewer context
949: Level: intermediate
951: .seealso: KSPMonitorSet(), KSPMonitorSingularValue()
952: @*/
953: PetscErrorCode KSPMonitorSingularValueCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
954: {
955: KSP ksp = (KSP) ctx;
957: PetscViewerAndFormatCreate(viewer, format, vf);
958: (*vf)->data = ctx;
959: KSPSetComputeSingularValues(ksp, PETSC_TRUE);
960: return 0;
961: }
963: /*@C
964: KSPMonitorDynamicTolerance - Recompute the inner tolerance in every outer iteration in an adaptive way.
966: Collective on ksp
968: Input Parameters:
969: + ksp - iterative context
970: . n - iteration number (not used)
971: . fnorm - the current residual norm
972: - dummy - some context as a C struct. fields:
973: coef: a scaling coefficient. default 1.0. can be passed through
974: -sub_ksp_dynamic_tolerance_param
975: bnrm: norm of the right-hand side. store it to avoid repeated calculation
977: Notes:
978: This may be useful for a flexibly preconditioner Krylov method to
979: control the accuracy of the inner solves needed to guarantee the
980: convergence of the outer iterations.
982: Level: advanced
984: .seealso: KSPMonitorDynamicToleranceDestroy()
985: @*/
986: PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
987: {
988: PC pc;
989: PetscReal outer_rtol, outer_abstol, outer_dtol, inner_rtol;
990: PetscInt outer_maxits,nksp,first,i;
991: KSPDynTolCtx *scale = (KSPDynTolCtx*)dummy;
992: KSP *subksp = NULL;
993: KSP kspinner;
994: PetscBool flg;
996: KSPGetPC(ksp, &pc);
998: /* compute inner_rtol */
999: if (scale->bnrm < 0.0) {
1000: Vec b;
1001: KSPGetRhs(ksp, &b);
1002: VecNorm(b, NORM_2, &(scale->bnrm));
1003: }
1004: KSPGetTolerances(ksp, &outer_rtol, &outer_abstol, &outer_dtol, &outer_maxits);
1005: inner_rtol = PetscMin(scale->coef * scale->bnrm * outer_rtol / fnorm, 0.999);
1006: /* PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Inner rtol = %g\n",
1007: (double)inner_rtol)); */
1009: /* if pc is ksp */
1010: PetscObjectTypeCompare((PetscObject)pc,PCKSP,&flg);
1011: if (flg) {
1012: PCKSPGetKSP(pc, &kspinner);
1013: KSPSetTolerances(kspinner, inner_rtol, outer_abstol, outer_dtol, outer_maxits);
1014: return 0;
1015: }
1017: /* if pc is bjacobi */
1018: PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&flg);
1019: if (flg) {
1020: PCBJacobiGetSubKSP(pc, &nksp, &first, &subksp);
1021: if (subksp) {
1022: for (i=0; i<nksp; i++) {
1023: KSPSetTolerances(subksp[i], inner_rtol, outer_abstol, outer_dtol, outer_maxits);
1024: }
1025: return 0;
1026: }
1027: }
1029: /* if pc is deflation*/
1030: PetscObjectTypeCompare((PetscObject)pc,PCDEFLATION,&flg);
1031: if (flg) {
1032: PCDeflationGetCoarseKSP(pc,&kspinner);
1033: KSPSetTolerances(kspinner,inner_rtol,outer_abstol,outer_dtol,PETSC_DEFAULT);
1034: return 0;
1035: }
1037: /* todo: dynamic tolerance may apply to other types of pc too */
1038: return 0;
1039: }
1041: /*
1042: Destroy the dummy context used in KSPMonitorDynamicTolerance()
1043: */
1044: PetscErrorCode KSPMonitorDynamicToleranceDestroy(void **dummy)
1045: {
1046: PetscFree(*dummy);
1047: return 0;
1048: }
1050: /*@C
1051: KSPConvergedSkip - Convergence test that do not return as converged
1052: until the maximum number of iterations is reached.
1054: Collective on ksp
1056: Input Parameters:
1057: + ksp - iterative context
1058: . n - iteration number
1059: . rnorm - 2-norm residual value (may be estimated)
1060: - dummy - unused convergence context
1062: Returns:
1063: . reason - KSP_CONVERGED_ITERATING, KSP_CONVERGED_ITS
1065: Notes:
1066: This should be used as the convergence test with the option
1067: KSPSetNormType(ksp,KSP_NORM_NONE), since norms of the residual are
1068: not computed. Convergence is then declared after the maximum number
1069: of iterations have been reached. Useful when one is using CG or
1070: BiCGStab as a smoother.
1072: Level: advanced
1074: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
1075: @*/
1076: PetscErrorCode KSPConvergedSkip(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
1077: {
1080: *reason = KSP_CONVERGED_ITERATING;
1081: if (n >= ksp->max_it) *reason = KSP_CONVERGED_ITS;
1082: return 0;
1083: }
1085: /*@C
1086: KSPConvergedDefaultCreate - Creates and initializes the space used by the KSPConvergedDefault() function context
1088: Note Collective
1090: Output Parameter:
1091: . ctx - convergence context
1093: Level: intermediate
1095: .seealso: KSPConvergedDefault(), KSPConvergedDefaultDestroy(), KSPSetConvergenceTest(), KSPSetTolerances(),
1096: KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultSetConvergedMaxits()
1097: @*/
1098: PetscErrorCode KSPConvergedDefaultCreate(void **ctx)
1099: {
1100: KSPConvergedDefaultCtx *cctx;
1102: PetscNew(&cctx);
1103: *ctx = cctx;
1104: return 0;
1105: }
1107: /*@
1108: KSPConvergedDefaultSetUIRNorm - makes the default convergence test use || B*(b - A*(initial guess))||
1109: instead of || B*b ||. In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDITIONED)
1110: is used there is no B in the above formula. UIRNorm is short for Use Initial Residual Norm.
1112: Collective on ksp
1114: Input Parameters:
1115: . ksp - iterative context
1117: Options Database:
1118: . -ksp_converged_use_initial_residual_norm <bool> - Use initial residual norm for computing relative convergence
1120: Notes:
1121: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
1123: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
1124: are defined in petscksp.h.
1126: If the convergence test is not KSPConvergedDefault() then this is ignored.
1128: If right preconditioning is being used then B does not appear in the above formula.
1130: Level: intermediate
1132: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultSetConvergedMaxits()
1133: @*/
1134: PetscErrorCode KSPConvergedDefaultSetUIRNorm(KSP ksp)
1135: {
1136: KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
1139: if (ksp->converged != KSPConvergedDefault) return 0;
1141: ctx->initialrtol = PETSC_TRUE;
1142: return 0;
1143: }
1145: /*@
1146: KSPConvergedDefaultSetUMIRNorm - makes the default convergence test use min(|| B*(b - A*(initial guess))||,|| B*b ||)
1147: In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDITIONED)
1148: is used there is no B in the above formula. UMIRNorm is short for Use Minimum Initial Residual Norm.
1150: Collective on ksp
1152: Input Parameters:
1153: . ksp - iterative context
1155: Options Database:
1156: . -ksp_converged_use_min_initial_residual_norm <bool> - Use minimum of initial residual norm and b for computing relative convergence
1158: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
1160: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
1161: are defined in petscksp.h.
1163: Level: intermediate
1165: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetConvergedMaxits()
1166: @*/
1167: PetscErrorCode KSPConvergedDefaultSetUMIRNorm(KSP ksp)
1168: {
1169: KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
1172: if (ksp->converged != KSPConvergedDefault) return 0;
1174: ctx->mininitialrtol = PETSC_TRUE;
1175: return 0;
1176: }
1178: /*@
1179: KSPConvergedDefaultSetConvergedMaxits - allows the default convergence test to declare convergence and return KSP_CONVERGED_ITS if the maximum number of iterations is reached
1181: Collective on ksp
1183: Input Parameters:
1184: + ksp - iterative context
1185: - flg - boolean flag
1187: Options Database:
1188: . -ksp_converged_maxits <bool> - Declare convergence if the maximum number of iterations is reached
1190: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
1192: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
1193: are defined in petscksp.h.
1195: Level: intermediate
1197: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultSetUIRNorm()
1198: @*/
1199: PetscErrorCode KSPConvergedDefaultSetConvergedMaxits(KSP ksp, PetscBool flg)
1200: {
1201: KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
1205: if (ksp->converged != KSPConvergedDefault) return 0;
1206: ctx->convmaxits = flg;
1207: return 0;
1208: }
1210: /*@C
1211: KSPConvergedDefault - Determines convergence of the linear iterative solvers by default
1213: Collective on ksp
1215: Input Parameters:
1216: + ksp - iterative context
1217: . n - iteration number
1218: . rnorm - residual norm (may be estimated, depending on the method may be the preconditioned residual norm)
1219: - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
1221: Output Parameter:
1222: . reason - the convergence reason; it is positive if the iteration has converged,
1223: negative if the iteration has diverged, and KSP_CONVERGED_ITERATING otherwise
1225: Notes:
1226: KSPConvergedDefault() reaches convergence when rnorm < MAX (rtol * rnorm_0, abstol);
1227: Divergence is detected if rnorm > dtol * rnorm_0, or when failures are detected throughout the iteration.
1228: By default, reaching the maximum number of iterations is considered divergence (i.e. KSP_DIVERGED_ITS).
1229: In order to have PETSc declaring convergence in such a case (i.e. KSP_CONVERGED_ITS), users can use KSPConvergedDefaultSetConvergedMaxits()
1231: where:
1232: + rtol - relative tolerance,
1233: . abstol - absolute tolerance.
1234: . dtol - divergence tolerance,
1235: - rnorm_0 - the two norm of the right hand side (or the preconditioned norm, depending on what was set with
1236: KSPSetNormType(). When initial guess is non-zero you
1237: can call KSPConvergedDefaultSetUIRNorm() to use the norm of (b - A*(initial guess))
1238: as the starting point for relative norm convergence testing, that is as rnorm_0
1240: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
1242: Use KSPSetNormType() (or -ksp_norm_type <none,preconditioned,unpreconditioned,natural>) to change the norm used for computing rnorm
1244: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
1246: This routine is used by KSP by default so the user generally never needs call it directly.
1248: Use KSPSetConvergenceTest() to provide your own test instead of using this one.
1250: Level: intermediate
1252: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
1253: KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultSetConvergedMaxits(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy()
1254: @*/
1255: PetscErrorCode KSPConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
1256: {
1257: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
1258: KSPNormType normtype;
1264: *reason = KSP_CONVERGED_ITERATING;
1266: if (cctx->convmaxits && n >= ksp->max_it) {
1267: *reason = KSP_CONVERGED_ITS;
1268: PetscInfo(ksp,"Linear solver has converged. Maximum number of iterations reached %D\n",n);
1269: return 0;
1270: }
1271: KSPGetNormType(ksp,&normtype);
1272: if (normtype == KSP_NORM_NONE) return 0;
1274: if (!n) {
1275: /* if user gives initial guess need to compute norm of b */
1276: if (!ksp->guess_zero && !cctx->initialrtol) {
1277: PetscReal snorm = 0.0;
1278: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
1279: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");
1280: VecNorm(ksp->vec_rhs,NORM_2,&snorm); /* <- b'*b */
1281: } else {
1282: Vec z;
1283: /* Should avoid allocating the z vector each time but cannot stash it in cctx because if KSPReset() is called the vector size might change */
1284: VecDuplicate(ksp->vec_rhs,&z);
1285: KSP_PCApply(ksp,ksp->vec_rhs,z);
1286: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
1287: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");
1288: VecNorm(z,NORM_2,&snorm); /* dp <- b'*B'*B*b */
1289: } else if (ksp->normtype == KSP_NORM_NATURAL) {
1290: PetscScalar norm;
1291: PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");
1292: VecDot(ksp->vec_rhs,z,&norm);
1293: snorm = PetscSqrtReal(PetscAbsScalar(norm)); /* dp <- b'*B*b */
1294: }
1295: VecDestroy(&z);
1296: }
1297: /* handle special case of zero RHS and nonzero guess */
1298: if (!snorm) {
1299: PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");
1300: snorm = rnorm;
1301: }
1302: if (cctx->mininitialrtol) ksp->rnorm0 = PetscMin(snorm,rnorm);
1303: else ksp->rnorm0 = snorm;
1304: } else {
1305: ksp->rnorm0 = rnorm;
1306: }
1307: ksp->ttol = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
1308: }
1310: if (n <= ksp->chknorm) return 0;
1312: if (PetscIsInfOrNanReal(rnorm)) {
1313: PCFailedReason pcreason;
1314: PetscInt sendbuf,recvbuf;
1315: PCGetFailedReasonRank(ksp->pc,&pcreason);
1316: sendbuf = (PetscInt)pcreason;
1317: MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPIU_MAX,PetscObjectComm((PetscObject)ksp));
1318: if (recvbuf) {
1319: *reason = KSP_DIVERGED_PC_FAILED;
1320: PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf);
1321: PetscInfo(ksp,"Linear solver pcsetup fails, declaring divergence \n");
1322: } else {
1323: *reason = KSP_DIVERGED_NANORINF;
1324: PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");
1325: }
1326: } else if (rnorm <= ksp->ttol) {
1327: if (rnorm < ksp->abstol) {
1328: PetscInfo(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);
1329: *reason = KSP_CONVERGED_ATOL;
1330: } else {
1331: if (cctx->initialrtol) {
1332: PetscInfo(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);
1333: } else {
1334: PetscInfo(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);
1335: }
1336: *reason = KSP_CONVERGED_RTOL;
1337: }
1338: } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
1339: PetscInfo(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);
1340: *reason = KSP_DIVERGED_DTOL;
1341: }
1342: return 0;
1343: }
1345: /*@C
1346: KSPConvergedDefaultDestroy - Frees the space used by the KSPConvergedDefault() function context
1348: Not Collective
1350: Input Parameters:
1351: . ctx - convergence context
1353: Level: intermediate
1355: .seealso: KSPConvergedDefault(), KSPConvergedDefaultCreate(), KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(),
1356: KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
1357: @*/
1358: PetscErrorCode KSPConvergedDefaultDestroy(void *ctx)
1359: {
1360: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
1362: VecDestroy(&cctx->work);
1363: PetscFree(ctx);
1364: return 0;
1365: }
1367: /*
1368: KSPBuildSolutionDefault - Default code to create/move the solution.
1370: Collective on ksp
1372: Input Parameters:
1373: + ksp - iterative context
1374: - v - pointer to the user's vector
1376: Output Parameter:
1377: . V - pointer to a vector containing the solution
1379: Level: advanced
1381: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1383: .seealso: KSPGetSolution(), KSPBuildResidualDefault()
1384: */
1385: PetscErrorCode KSPBuildSolutionDefault(KSP ksp,Vec v,Vec *V)
1386: {
1387: if (ksp->pc_side == PC_RIGHT) {
1388: if (ksp->pc) {
1389: if (v) {
1390: KSP_PCApply(ksp,ksp->vec_sol,v); *V = v;
1391: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with right preconditioner");
1392: } else {
1393: if (v) {
1394: VecCopy(ksp->vec_sol,v); *V = v;
1395: } else *V = ksp->vec_sol;
1396: }
1397: } else if (ksp->pc_side == PC_SYMMETRIC) {
1398: if (ksp->pc) {
1400: if (v) {
1401: PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v);
1402: *V = v;
1403: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner");
1404: } else {
1405: if (v) {
1406: VecCopy(ksp->vec_sol,v); *V = v;
1407: } else *V = ksp->vec_sol;
1408: }
1409: } else {
1410: if (v) {
1411: VecCopy(ksp->vec_sol,v); *V = v;
1412: } else *V = ksp->vec_sol;
1413: }
1414: return 0;
1415: }
1417: /*
1418: KSPBuildResidualDefault - Default code to compute the residual.
1420: Collecive on ksp
1422: Input Parameters:
1423: . ksp - iterative context
1424: . t - pointer to temporary vector
1425: . v - pointer to user vector
1427: Output Parameter:
1428: . V - pointer to a vector containing the residual
1430: Level: advanced
1432: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1434: .seealso: KSPBuildSolutionDefault()
1435: */
1436: PetscErrorCode KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec *V)
1437: {
1438: Mat Amat,Pmat;
1440: if (!ksp->pc) KSPGetPC(ksp,&ksp->pc);
1441: PCGetOperators(ksp->pc,&Amat,&Pmat);
1442: KSPBuildSolution(ksp,t,NULL);
1443: KSP_MatMult(ksp,Amat,t,v);
1444: VecAYPX(v,-1.0,ksp->vec_rhs);
1445: *V = v;
1446: return 0;
1447: }
1449: /*@C
1450: KSPCreateVecs - Gets a number of work vectors.
1452: Collective on ksp
1454: Input Parameters:
1455: + ksp - iterative context
1456: . rightn - number of right work vectors
1457: - leftn - number of left work vectors to allocate
1459: Output Parameters:
1460: + right - the array of vectors created
1461: - left - the array of left vectors
1463: Note: The right vector has as many elements as the matrix has columns. The left
1464: vector has as many elements as the matrix has rows.
1466: The vectors are new vectors that are not owned by the KSP, they should be destroyed with calls to VecDestroyVecs() when no longer needed.
1468: 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
1469: that does not exist tries to get them from the DM (if it is provided).
1471: Level: advanced
1473: .seealso: MatCreateVecs(), VecDestroyVecs()
1475: @*/
1476: PetscErrorCode KSPCreateVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
1477: {
1478: Vec vecr = NULL,vecl = NULL;
1479: PetscBool matset,pmatset,isshell,preferdm = PETSC_FALSE;
1480: Mat mat = NULL;
1482: if (ksp->dm) {
1483: PetscObjectTypeCompare((PetscObject) ksp->dm, DMSHELL, &isshell);
1484: preferdm = isshell ? PETSC_FALSE : PETSC_TRUE;
1485: }
1486: if (rightn) {
1488: if (ksp->vec_sol) vecr = ksp->vec_sol;
1489: else {
1490: if (preferdm) {
1491: DMGetGlobalVector(ksp->dm,&vecr);
1492: } else if (ksp->pc) {
1493: PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
1494: /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
1495: if (matset) {
1496: PCGetOperators(ksp->pc,&mat,NULL);
1497: MatCreateVecs(mat,&vecr,NULL);
1498: } else if (pmatset) {
1499: PCGetOperators(ksp->pc,NULL,&mat);
1500: MatCreateVecs(mat,&vecr,NULL);
1501: }
1502: }
1503: if (!vecr && ksp->dm) {
1504: DMGetGlobalVector(ksp->dm,&vecr);
1505: }
1507: }
1508: VecDuplicateVecs(vecr,rightn,right);
1509: if (!ksp->vec_sol) {
1510: if (preferdm) {
1511: DMRestoreGlobalVector(ksp->dm,&vecr);
1512: } else if (mat) {
1513: VecDestroy(&vecr);
1514: } else if (ksp->dm) {
1515: DMRestoreGlobalVector(ksp->dm,&vecr);
1516: }
1517: }
1518: }
1519: if (leftn) {
1521: if (ksp->vec_rhs) vecl = ksp->vec_rhs;
1522: else {
1523: if (preferdm) {
1524: DMGetGlobalVector(ksp->dm,&vecl);
1525: } else if (ksp->pc) {
1526: PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
1527: /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
1528: if (matset) {
1529: PCGetOperators(ksp->pc,&mat,NULL);
1530: MatCreateVecs(mat,NULL,&vecl);
1531: } else if (pmatset) {
1532: PCGetOperators(ksp->pc,NULL,&mat);
1533: MatCreateVecs(mat,NULL,&vecl);
1534: }
1535: }
1536: if (!vecl && ksp->dm) DMGetGlobalVector(ksp->dm,&vecl);
1538: }
1539: VecDuplicateVecs(vecl,leftn,left);
1540: if (!ksp->vec_rhs) {
1541: if (preferdm) {
1542: DMRestoreGlobalVector(ksp->dm,&vecl);
1543: } else if (mat) {
1544: VecDestroy(&vecl);
1545: } else if (ksp->dm) {
1546: DMRestoreGlobalVector(ksp->dm,&vecl);
1547: }
1548: }
1549: }
1550: return 0;
1551: }
1553: /*@C
1554: KSPSetWorkVecs - Sets a number of work vectors into a KSP object
1556: Collective on ksp
1558: Input Parameters:
1559: + ksp - iterative context
1560: - nw - number of work vectors to allocate
1562: Level: developer
1564: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1566: @*/
1567: PetscErrorCode KSPSetWorkVecs(KSP ksp,PetscInt nw)
1568: {
1569: VecDestroyVecs(ksp->nwork,&ksp->work);
1570: ksp->nwork = nw;
1571: KSPCreateVecs(ksp,nw,&ksp->work,0,NULL);
1572: PetscLogObjectParents(ksp,nw,ksp->work);
1573: return 0;
1574: }
1576: /*
1577: KSPDestroyDefault - Destroys a iterative context variable for methods with
1578: no separate context. Preferred calling sequence KSPDestroy().
1580: Input Parameter:
1581: . ksp - the iterative context
1583: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1585: */
1586: PetscErrorCode KSPDestroyDefault(KSP ksp)
1587: {
1589: PetscFree(ksp->data);
1590: return 0;
1591: }
1593: /*@
1594: KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.
1596: Not Collective
1598: Input Parameter:
1599: . ksp - the KSP context
1601: Output Parameter:
1602: . reason - negative value indicates diverged, positive value converged, see KSPConvergedReason
1604: Possible values for reason: See also manual page for each reason
1605: $ KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
1606: $ KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
1607: $ KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration, or when the KSPConvergedSkip() convergence test routine is set.
1608: $ KSP_CONVERGED_CG_NEG_CURVE (see note below)
1609: $ KSP_CONVERGED_CG_CONSTRAINED (see note below)
1610: $ KSP_CONVERGED_STEP_LENGTH (see note below)
1611: $ KSP_CONVERGED_ITERATING (returned if the solver is not yet finished)
1612: $ KSP_DIVERGED_ITS (required more than its to reach convergence)
1613: $ KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
1614: $ KSP_DIVERGED_NANORINF (residual norm became Not-a-number or Inf likely due to 0/0)
1615: $ KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
1616: $ KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial residual. Try a different preconditioner, or a different initial Level.)
1618: Options Database:
1619: . -ksp_converged_reason - prints the reason to standard out
1621: Notes:
1622: If this routine is called before or doing the KSPSolve() the value of KSP_CONVERGED_ITERATING is returned
1624: 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
1625: solvers which are used by the SNESNEWTONTR (trust region) solver.
1627: Level: intermediate
1629: .seealso: KSPSetConvergenceTest(), KSPConvergedDefault(), KSPSetTolerances(), KSPConvergedReason,
1630: KSPConvergedReasonView()
1631: @*/
1632: PetscErrorCode KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
1633: {
1636: *reason = ksp->reason;
1637: return 0;
1638: }
1640: /*@C
1641: KSPGetConvergedReasonString - Return a human readable string for ksp converged reason
1643: Not Collective
1645: Input Parameter:
1646: . ksp - the KSP context
1648: Output Parameter:
1649: . strreason - a human readable string that describes ksp converged reason
1651: Level: beginner
1653: .seealso: KSPGetConvergedReason()
1654: @*/
1655: PetscErrorCode KSPGetConvergedReasonString(KSP ksp,const char** strreason)
1656: {
1659: *strreason = KSPConvergedReasons[ksp->reason];
1660: return 0;
1661: }
1663: #include <petsc/private/dmimpl.h>
1664: /*@
1665: KSPSetDM - Sets the DM that may be used by some preconditioners
1667: Logically Collective on ksp
1669: Input Parameters:
1670: + ksp - the preconditioner context
1671: - dm - the dm, cannot be NULL
1673: Notes:
1674: If this is used then the KSP will attempt to use the DM to create the matrix and use the routine set with
1675: DMKSPSetComputeOperators(). Use KSPSetDMActive(ksp,PETSC_FALSE) to instead use the matrix you've provided with
1676: KSPSetOperators().
1678: A DM can only be used for solving one problem at a time because information about the problem is stored on the DM,
1679: even when not using interfaces like DMKSPSetComputeOperators(). Use DMClone() to get a distinct DM when solving
1680: different problems using the same function space.
1682: Level: intermediate
1684: .seealso: KSPGetDM(), KSPSetDMActive(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess(), DMKSPSetComputeOperators(), DMKSPSetComputeRHS(), DMKSPSetComputeInitialGuess()
1685: @*/
1686: PetscErrorCode KSPSetDM(KSP ksp,DM dm)
1687: {
1688: PC pc;
1692: PetscObjectReference((PetscObject)dm);
1693: if (ksp->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */
1694: if (ksp->dm->dmksp && !dm->dmksp) {
1695: DMKSP kdm;
1696: DMCopyDMKSP(ksp->dm,dm);
1697: DMGetDMKSP(ksp->dm,&kdm);
1698: if (kdm->originaldm == ksp->dm) kdm->originaldm = dm; /* Grant write privileges to the replacement DM */
1699: }
1700: DMDestroy(&ksp->dm);
1701: }
1702: ksp->dm = dm;
1703: ksp->dmAuto = PETSC_FALSE;
1704: KSPGetPC(ksp,&pc);
1705: PCSetDM(pc,dm);
1706: ksp->dmActive = PETSC_TRUE;
1707: return 0;
1708: }
1710: /*@
1711: KSPSetDMActive - Indicates the DM should be used to generate the linear system matrix and right hand side
1713: Logically Collective on ksp
1715: Input Parameters:
1716: + ksp - the preconditioner context
1717: - flg - use the DM
1719: Notes:
1720: 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.
1722: Level: intermediate
1724: .seealso: KSPGetDM(), KSPSetDM(), SNESSetDM(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess()
1725: @*/
1726: PetscErrorCode KSPSetDMActive(KSP ksp,PetscBool flg)
1727: {
1730: ksp->dmActive = flg;
1731: return 0;
1732: }
1734: /*@
1735: KSPGetDM - Gets the DM that may be used by some preconditioners
1737: Not Collective
1739: Input Parameter:
1740: . ksp - the preconditioner context
1742: Output Parameter:
1743: . dm - the dm
1745: Level: intermediate
1747: .seealso: KSPSetDM(), KSPSetDMActive()
1748: @*/
1749: PetscErrorCode KSPGetDM(KSP ksp,DM *dm)
1750: {
1752: if (!ksp->dm) {
1753: DMShellCreate(PetscObjectComm((PetscObject)ksp),&ksp->dm);
1754: ksp->dmAuto = PETSC_TRUE;
1755: }
1756: *dm = ksp->dm;
1757: return 0;
1758: }
1760: /*@
1761: KSPSetApplicationContext - Sets the optional user-defined context for the linear solver.
1763: Logically Collective on ksp
1765: Input Parameters:
1766: + ksp - the KSP context
1767: - usrP - optional user context
1769: Fortran Notes:
1770: To use this from Fortran you must write a Fortran interface definition for this
1771: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1773: Level: intermediate
1775: .seealso: KSPGetApplicationContext()
1776: @*/
1777: PetscErrorCode KSPSetApplicationContext(KSP ksp,void *usrP)
1778: {
1779: PC pc;
1782: ksp->user = usrP;
1783: KSPGetPC(ksp,&pc);
1784: PCSetApplicationContext(pc,usrP);
1785: return 0;
1786: }
1788: /*@
1789: KSPGetApplicationContext - Gets the user-defined context for the linear solver.
1791: Not Collective
1793: Input Parameter:
1794: . ksp - KSP context
1796: Output Parameter:
1797: . usrP - user context
1799: Fortran Notes:
1800: To use this from Fortran you must write a Fortran interface definition for this
1801: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1803: Level: intermediate
1805: .seealso: KSPSetApplicationContext()
1806: @*/
1807: PetscErrorCode KSPGetApplicationContext(KSP ksp,void *usrP)
1808: {
1810: *(void**)usrP = ksp->user;
1811: return 0;
1812: }
1814: #include <petsc/private/pcimpl.h>
1816: /*@
1817: KSPCheckSolve - Checks if the PCSetUp() or KSPSolve() failed and set the error flag for the outer PC. A KSP_DIVERGED_ITS is
1818: not considered a failure in this context
1820: Collective on ksp
1822: Input Parameters:
1823: + ksp - the linear solver (KSP) context.
1824: . pc - the preconditioner context
1825: - vec - a vector that will be initialized with Inf to indicate lack of convergence
1827: Notes: this may be called by a subset of the processes in the PC
1829: Level: developer
1831: Developer Note: this is used to manage returning from preconditioners whose inner KSP solvers have failed in some way
1833: .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckNorm(), KSPCheckDot()
1834: @*/
1835: PetscErrorCode KSPCheckSolve(KSP ksp,PC pc,Vec vec)
1836: {
1837: PCFailedReason pcreason;
1838: PC subpc;
1840: KSPGetPC(ksp,&subpc);
1841: PCGetFailedReason(subpc,&pcreason);
1842: if (pcreason || (ksp->reason < 0 && ksp->reason != KSP_DIVERGED_ITS)) {
1844: else {
1845: PetscInfo(ksp,"Detected not converged in KSP inner solve: KSP reason %s PC reason %s\n",KSPConvergedReasons[ksp->reason],PCFailedReasons[pcreason]);
1846: pc->failedreason = PC_SUBPC_ERROR;
1847: if (vec) VecSetInf(vec);
1848: }
1849: }
1850: return 0;
1851: }