Actual source code: itfunc.c
petsc-3.14.6 2021-03-30
1: /*
2: Interface KSP routines that the user calls.
3: */
5: #include <petsc/private/kspimpl.h>
6: #include <petsc/private/matimpl.h>
7: #include <petscdm.h>
9: PETSC_STATIC_INLINE PetscErrorCode ObjectView(PetscObject obj, PetscViewer viewer, PetscViewerFormat format)
10: {
13: PetscViewerPushFormat(viewer, format);
14: PetscObjectView(obj, viewer);
15: PetscViewerPopFormat(viewer);
16: return(0);
17: }
19: /*@
20: KSPComputeExtremeSingularValues - Computes the extreme singular values
21: for the preconditioned operator. Called after or during KSPSolve().
23: Not Collective
25: Input Parameter:
26: . ksp - iterative context obtained from KSPCreate()
28: Output Parameters:
29: . emin, emax - extreme singular values
31: Options Database Keys:
32: . -ksp_view_singularvalues - compute extreme singular values and print when KSPSolve completes.
34: Notes:
35: One must call KSPSetComputeSingularValues() before calling KSPSetUp()
36: (or use the option -ksp_view_eigenvalues) in order for this routine to work correctly.
38: Many users may just want to use the monitoring routine
39: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
40: to print the extreme singular values at each iteration of the linear solve.
42: Estimates of the smallest singular value may be very inaccurate, especially if the Krylov method has not converged.
43: The largest singular value is usually accurate to within a few percent if the method has converged, but is still not
44: intended for eigenanalysis.
46: Disable restarts if using KSPGMRES, otherwise this estimate will only be using those iterations after the last
47: restart. See KSPGMRESSetRestart() for more details.
49: Level: advanced
51: .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeEigenvalues(), KSP
52: @*/
53: PetscErrorCode KSPComputeExtremeSingularValues(KSP ksp,PetscReal *emax,PetscReal *emin)
54: {
61: if (!ksp->calc_sings) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Singular values not requested before KSPSetUp()");
63: if (ksp->ops->computeextremesingularvalues) {
64: (*ksp->ops->computeextremesingularvalues)(ksp,emax,emin);
65: } else {
66: *emin = -1.0;
67: *emax = -1.0;
68: }
69: return(0);
70: }
72: /*@
73: KSPComputeEigenvalues - Computes the extreme eigenvalues for the
74: preconditioned operator. Called after or during KSPSolve().
76: Not Collective
78: Input Parameters:
79: + ksp - iterative context obtained from KSPCreate()
80: - n - size of arrays r and c. The number of eigenvalues computed (neig) will, in
81: general, be less than this.
83: Output Parameters:
84: + r - real part of computed eigenvalues, provided by user with a dimension of at least n
85: . c - complex part of computed eigenvalues, provided by user with a dimension of at least n
86: - neig - actual number of eigenvalues computed (will be less than or equal to n)
88: Options Database Keys:
89: . -ksp_view_eigenvalues - Prints eigenvalues to stdout
91: Notes:
92: The number of eigenvalues estimated depends on the size of the Krylov space
93: generated during the KSPSolve() ; for example, with
94: CG it corresponds to the number of CG iterations, for GMRES it is the number
95: of GMRES iterations SINCE the last restart. Any extra space in r[] and c[]
96: will be ignored.
98: KSPComputeEigenvalues() does not usually provide accurate estimates; it is
99: intended only for assistance in understanding the convergence of iterative
100: methods, not for eigenanalysis. For accurate computation of eigenvalues we recommend using
101: the excellent package SLEPc.
103: One must call KSPSetComputeEigenvalues() before calling KSPSetUp()
104: in order for this routine to work correctly.
106: Many users may just want to use the monitoring routine
107: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
108: to print the singular values at each iteration of the linear solve.
110: Level: advanced
112: .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), KSP
113: @*/
114: PetscErrorCode KSPComputeEigenvalues(KSP ksp,PetscInt n,PetscReal r[],PetscReal c[],PetscInt *neig)
115: {
122: if (n<0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Requested < 0 Eigenvalues");
124: if (!ksp->calc_sings) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Eigenvalues not requested before KSPSetUp()");
126: if (n && ksp->ops->computeeigenvalues) {
127: (*ksp->ops->computeeigenvalues)(ksp,n,r,c,neig);
128: } else {
129: *neig = 0;
130: }
131: return(0);
132: }
134: /*@
135: KSPComputeRitz - Computes the Ritz or harmonic Ritz pairs associated to the
136: smallest or largest in modulus, for the preconditioned operator.
137: Called after KSPSolve().
139: Not Collective
141: Input Parameters:
142: + ksp - iterative context obtained from KSPCreate()
143: . ritz - PETSC_TRUE or PETSC_FALSE for ritz pairs or harmonic Ritz pairs, respectively
144: . small - PETSC_TRUE or PETSC_FALSE for smallest or largest (harmonic) Ritz values, respectively
145: - nrit - number of (harmonic) Ritz pairs to compute
147: Output Parameters:
148: + nrit - actual number of computed (harmonic) Ritz pairs
149: . S - multidimensional vector with Ritz vectors
150: . tetar - real part of the Ritz values
151: - tetai - imaginary part of the Ritz values
153: Notes:
154: -For GMRES, the (harmonic) Ritz pairs are computed from the Hessenberg matrix obtained during
155: the last complete cycle, or obtained at the end of the solution if the method is stopped before
156: a restart. Then, the number of actual (harmonic) Ritz pairs computed is less or equal to the restart
157: parameter for GMRES if a complete cycle has been performed or less or equal to the number of GMRES
158: iterations.
159: -Moreover, for real matrices, the (harmonic) Ritz pairs are possibly complex-valued. In such a case,
160: the routine selects the complex (harmonic) Ritz value and its conjugate, and two successive columns of S
161: are equal to the real and the imaginary parts of the associated vectors.
162: -the (harmonic) Ritz pairs are given in order of increasing (harmonic) Ritz values in modulus
163: -this is currently not implemented when PETSc is built with complex numbers
165: One must call KSPSetComputeRitz() before calling KSPSetUp()
166: in order for this routine to work correctly.
168: Level: advanced
170: .seealso: KSPSetComputeRitz(), KSP
171: @*/
172: PetscErrorCode KSPComputeRitz(KSP ksp,PetscBool ritz,PetscBool small,PetscInt *nrit,Vec S[],PetscReal tetar[],PetscReal tetai[])
173: {
178: if (!ksp->calc_ritz) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Ritz pairs not requested before KSPSetUp()");
179: if (ksp->ops->computeritz) {(*ksp->ops->computeritz)(ksp,ritz,small,nrit,S,tetar,tetai);}
180: return(0);
181: }
182: /*@
183: KSPSetUpOnBlocks - Sets up the preconditioner for each block in
184: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
185: methods.
187: Collective on ksp
189: Input Parameter:
190: . ksp - the KSP context
192: Notes:
193: KSPSetUpOnBlocks() is a routine that the user can optinally call for
194: more precise profiling (via -log_view) of the setup phase for these
195: block preconditioners. If the user does not call KSPSetUpOnBlocks(),
196: it will automatically be called from within KSPSolve().
198: Calling KSPSetUpOnBlocks() is the same as calling PCSetUpOnBlocks()
199: on the PC context within the KSP context.
201: Level: advanced
203: .seealso: PCSetUpOnBlocks(), KSPSetUp(), PCSetUp(), KSP
204: @*/
205: PetscErrorCode KSPSetUpOnBlocks(KSP ksp)
206: {
207: PC pc;
209: PCFailedReason pcreason;
213: KSPGetPC(ksp,&pc);
214: PCSetUpOnBlocks(pc);
215: PCGetFailedReasonRank(pc,&pcreason);
216: /* TODO: this code was wrong and is still wrong, there is no way to propogate the failure to all processes; their is no code to handle a ksp->reason on only some ranks */
217: if (pcreason) {
218: ksp->reason = KSP_DIVERGED_PC_FAILED;
219: }
220: return(0);
221: }
223: /*@
224: KSPSetReusePreconditioner - reuse the current preconditioner, do not construct a new one even if the operator changes
226: Collective on ksp
228: Input Parameters:
229: + ksp - iterative context obtained from KSPCreate()
230: - flag - PETSC_TRUE to reuse the current preconditioner
232: Level: intermediate
234: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), PCSetReusePreconditioner(), KSP
235: @*/
236: PetscErrorCode KSPSetReusePreconditioner(KSP ksp,PetscBool flag)
237: {
238: PC pc;
243: KSPGetPC(ksp,&pc);
244: PCSetReusePreconditioner(pc,flag);
245: return(0);
246: }
248: /*@
249: KSPGetReusePreconditioner - Determines if the KSP reuses the current preconditioner even if the operator in the preconditioner has changed.
251: Collective on ksp
253: Input Parameters:
254: . ksp - iterative context obtained from KSPCreate()
256: Output Parameters:
257: . flag - the boolean flag
259: Level: intermediate
261: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), KSPSetReusePreconditioner(), KSP
262: @*/
263: PetscErrorCode KSPGetReusePreconditioner(KSP ksp,PetscBool *flag)
264: {
270: *flag = PETSC_FALSE;
271: if (ksp->pc) {
272: PCGetReusePreconditioner(ksp->pc,flag);
273: }
274: return(0);
275: }
277: /*@
278: KSPSetSkipPCSetFromOptions - prevents KSPSetFromOptions() from call PCSetFromOptions(). This is used if the same PC is shared by more than one KSP so its options are not resetable for each KSP
280: Collective on ksp
282: Input Parameters:
283: + ksp - iterative context obtained from KSPCreate()
284: - flag - PETSC_TRUE to skip calling the PCSetFromOptions()
286: Level: intermediate
288: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), PCSetReusePreconditioner(), KSP
289: @*/
290: PetscErrorCode KSPSetSkipPCSetFromOptions(KSP ksp,PetscBool flag)
291: {
294: ksp->skippcsetfromoptions = flag;
295: return(0);
296: }
298: /*@
299: KSPSetUp - Sets up the internal data structures for the
300: later use of an iterative solver.
302: Collective on ksp
304: Input Parameter:
305: . ksp - iterative context obtained from KSPCreate()
307: Level: developer
309: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), KSP
310: @*/
311: PetscErrorCode KSPSetUp(KSP ksp)
312: {
314: Mat A,B;
315: Mat mat,pmat;
316: MatNullSpace nullsp;
317: PCFailedReason pcreason;
322: /* reset the convergence flag from the previous solves */
323: ksp->reason = KSP_CONVERGED_ITERATING;
325: if (!((PetscObject)ksp)->type_name) {
326: KSPSetType(ksp,KSPGMRES);
327: }
328: KSPSetUpNorms_Private(ksp,PETSC_TRUE,&ksp->normtype,&ksp->pc_side);
330: if (ksp->dmActive && !ksp->setupstage) {
331: /* first time in so build matrix and vector data structures using DM */
332: if (!ksp->vec_rhs) {DMCreateGlobalVector(ksp->dm,&ksp->vec_rhs);}
333: if (!ksp->vec_sol) {DMCreateGlobalVector(ksp->dm,&ksp->vec_sol);}
334: DMCreateMatrix(ksp->dm,&A);
335: KSPSetOperators(ksp,A,A);
336: PetscObjectDereference((PetscObject)A);
337: }
339: if (ksp->dmActive) {
340: DMKSP kdm;
341: DMGetDMKSP(ksp->dm,&kdm);
343: if (kdm->ops->computeinitialguess && ksp->setupstage != KSP_SETUP_NEWRHS) {
344: /* only computes initial guess the first time through */
345: (*kdm->ops->computeinitialguess)(ksp,ksp->vec_sol,kdm->initialguessctx);
346: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
347: }
348: if (kdm->ops->computerhs) {
349: (*kdm->ops->computerhs)(ksp,ksp->vec_rhs,kdm->rhsctx);
350: }
352: if (ksp->setupstage != KSP_SETUP_NEWRHS) {
353: if (kdm->ops->computeoperators) {
354: KSPGetOperators(ksp,&A,&B);
355: (*kdm->ops->computeoperators)(ksp,A,B,kdm->operatorsctx);
356: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You called KSPSetDM() but did not use DMKSPSetComputeOperators() or KSPSetDMActive(ksp,PETSC_FALSE);");
357: }
358: }
360: if (ksp->setupstage == KSP_SETUP_NEWRHS) return(0);
361: PetscLogEventBegin(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);
363: switch (ksp->setupstage) {
364: case KSP_SETUP_NEW:
365: (*ksp->ops->setup)(ksp);
366: break;
367: case KSP_SETUP_NEWMATRIX: { /* This should be replaced with a more general mechanism */
368: if (ksp->setupnewmatrix) {
369: (*ksp->ops->setup)(ksp);
370: }
371: } break;
372: default: break;
373: }
375: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
376: PCGetOperators(ksp->pc,&mat,&pmat);
377: /* scale the matrix if requested */
378: if (ksp->dscale) {
379: PetscScalar *xx;
380: PetscInt i,n;
381: PetscBool zeroflag = PETSC_FALSE;
382: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
383: if (!ksp->diagonal) { /* allocate vector to hold diagonal */
384: MatCreateVecs(pmat,&ksp->diagonal,NULL);
385: }
386: MatGetDiagonal(pmat,ksp->diagonal);
387: VecGetLocalSize(ksp->diagonal,&n);
388: VecGetArray(ksp->diagonal,&xx);
389: for (i=0; i<n; i++) {
390: if (xx[i] != 0.0) xx[i] = 1.0/PetscSqrtReal(PetscAbsScalar(xx[i]));
391: else {
392: xx[i] = 1.0;
393: zeroflag = PETSC_TRUE;
394: }
395: }
396: VecRestoreArray(ksp->diagonal,&xx);
397: if (zeroflag) {
398: PetscInfo(ksp,"Zero detected in diagonal of matrix, using 1 at those locations\n");
399: }
400: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
401: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
402: ksp->dscalefix2 = PETSC_FALSE;
403: }
404: PetscLogEventEnd(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);
405: PCSetErrorIfFailure(ksp->pc,ksp->errorifnotconverged);
406: PCSetUp(ksp->pc);
407: PCGetFailedReasonRank(ksp->pc,&pcreason);
408: /* TODO: this code was wrong and is still wrong, there is no way to propogate the failure to all processes; their is no code to handle a ksp->reason on only some ranks */
409: if (pcreason) {
410: ksp->reason = KSP_DIVERGED_PC_FAILED;
411: }
413: MatGetNullSpace(mat,&nullsp);
414: if (nullsp) {
415: PetscBool test = PETSC_FALSE;
416: PetscOptionsGetBool(((PetscObject)ksp)->options,((PetscObject)ksp)->prefix,"-ksp_test_null_space",&test,NULL);
417: if (test) {
418: MatNullSpaceTest(nullsp,mat,NULL);
419: }
420: }
421: ksp->setupstage = KSP_SETUP_NEWRHS;
422: return(0);
423: }
425: /*@C
426: KSPConvergedReasonView - Displays the reason a KSP solve converged or diverged to a viewer
428: Collective on ksp
430: Parameter:
431: + ksp - iterative context obtained from KSPCreate()
432: - viewer - the viewer to display the reason
434: Options Database Keys:
435: + -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations
436: - -ksp_converged_reason ::failed - only print reason and number of iterations when diverged
438: Notes:
439: To change the format of the output call PetscViewerPushFormat(viewer,format) before this call. Use PETSC_VIEWER_DEFAULT for the default,
440: use PETSC_VIEWER_FAILED to only display a reason if it fails.
442: Level: beginner
444: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
445: KSPSolveTranspose(), KSPGetIterationNumber(), KSP, KSPGetConvergedReason(), PetscViewerPushFormat(), PetscViewerPopFormat()
446: @*/
447: PetscErrorCode KSPConvergedReasonView(KSP ksp, PetscViewer viewer)
448: {
449: PetscErrorCode ierr;
450: PetscBool isAscii;
451: PetscViewerFormat format;
454: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));
455: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
456: if (isAscii) {
457: PetscViewerGetFormat(viewer, &format);
458: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
459: if (ksp->reason > 0 && format != PETSC_VIEWER_FAILED) {
460: if (((PetscObject) ksp)->prefix) {
461: PetscViewerASCIIPrintf(viewer,"Linear %s solve converged due to %s iterations %D\n",((PetscObject) ksp)->prefix,KSPConvergedReasons[ksp->reason],ksp->its);
462: } else {
463: PetscViewerASCIIPrintf(viewer,"Linear solve converged due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
464: }
465: } else if (ksp->reason <= 0) {
466: if (((PetscObject) ksp)->prefix) {
467: PetscViewerASCIIPrintf(viewer,"Linear %s solve did not converge due to %s iterations %D\n",((PetscObject) ksp)->prefix,KSPConvergedReasons[ksp->reason],ksp->its);
468: } else {
469: PetscViewerASCIIPrintf(viewer,"Linear solve did not converge due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
470: }
471: if (ksp->reason == KSP_DIVERGED_PC_FAILED) {
472: PCFailedReason reason;
473: PCGetFailedReason(ksp->pc,&reason);
474: PetscViewerASCIIPrintf(viewer," PC failed due to %s \n",PCFailedReasons[reason]);
475: }
476: }
477: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
478: }
479: return(0);
480: }
482: /*@
483: KSPConvergedReasonViewFromOptions - Processes command line options to determine if/how a KSPReason is to be viewed.
485: Collective on ksp
487: Input Parameters:
488: . ksp - the KSP object
490: Level: intermediate
492: @*/
493: PetscErrorCode KSPConvergedReasonViewFromOptions(KSP ksp)
494: {
495: PetscViewer viewer;
496: PetscBool flg;
497: PetscViewerFormat format;
498: PetscErrorCode ierr;
501: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->options,((PetscObject)ksp)->prefix,"-ksp_converged_reason",&viewer,&format,&flg);
502: if (flg) {
503: PetscViewerPushFormat(viewer,format);
504: KSPConvergedReasonView(ksp, viewer);
505: PetscViewerPopFormat(viewer);
506: PetscViewerDestroy(&viewer);
507: }
508: return(0);
509: }
511: #include <petscdraw.h>
513: static PetscErrorCode KSPViewEigenvalues_Internal(KSP ksp, PetscBool isExplicit, PetscViewer viewer, PetscViewerFormat format)
514: {
515: PetscReal *r, *c;
516: PetscInt n, i, neig;
517: PetscBool isascii, isdraw;
518: PetscMPIInt rank;
522: MPI_Comm_rank(PetscObjectComm((PetscObject) ksp), &rank);
523: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
524: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERDRAW, &isdraw);
525: if (isExplicit) {
526: VecGetSize(ksp->vec_sol,&n);
527: PetscMalloc2(n, &r, n, &c);
528: KSPComputeEigenvaluesExplicitly(ksp, n, r, c);
529: neig = n;
530: } else {
531: PetscInt nits;
533: KSPGetIterationNumber(ksp, &nits);
534: n = nits+2;
535: if (!nits) {PetscViewerASCIIPrintf(viewer, "Zero iterations in solver, cannot approximate any eigenvalues\n");return(0);}
536: PetscMalloc2(n, &r, n, &c);
537: KSPComputeEigenvalues(ksp, n, r, c, &neig);
538: }
539: if (isascii) {
540: PetscViewerASCIIPrintf(viewer, "%s computed eigenvalues\n", isExplicit ? "Explicitly" : "Iteratively");
541: for (i = 0; i < neig; ++i) {
542: if (c[i] >= 0.0) {PetscViewerASCIIPrintf(viewer, "%g + %gi\n", (double) r[i], (double) c[i]);}
543: else {PetscViewerASCIIPrintf(viewer, "%g - %gi\n", (double) r[i], -(double) c[i]);}
544: }
545: } else if (isdraw && !rank) {
546: PetscDraw draw;
547: PetscDrawSP drawsp;
549: if (format == PETSC_VIEWER_DRAW_CONTOUR) {
550: KSPPlotEigenContours_Private(ksp,neig,r,c);
551: } else {
552: if (!ksp->eigviewer) {PetscViewerDrawOpen(PETSC_COMM_SELF,NULL,isExplicit ? "Explicitly Computed Eigenvalues" : "Iteratively Computed Eigenvalues",PETSC_DECIDE,PETSC_DECIDE,400,400,&ksp->eigviewer);}
553: PetscViewerDrawGetDraw(ksp->eigviewer,0,&draw);
554: PetscDrawSPCreate(draw,1,&drawsp);
555: PetscDrawSPReset(drawsp);
556: for (i = 0; i < neig; ++i) {PetscDrawSPAddPoint(drawsp,r+i,c+i);}
557: PetscDrawSPDraw(drawsp,PETSC_TRUE);
558: PetscDrawSPSave(drawsp);
559: PetscDrawSPDestroy(&drawsp);
560: }
561: }
562: PetscFree2(r, c);
563: return(0);
564: }
566: static PetscErrorCode KSPViewSingularvalues_Internal(KSP ksp, PetscViewer viewer, PetscViewerFormat format)
567: {
568: PetscReal smax, smin;
569: PetscInt nits;
570: PetscBool isascii;
574: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
575: KSPGetIterationNumber(ksp, &nits);
576: if (!nits) {PetscViewerASCIIPrintf(viewer, "Zero iterations in solver, cannot approximate any singular values\n");return(0);}
577: KSPComputeExtremeSingularValues(ksp, &smax, &smin);
578: if (isascii) {PetscViewerASCIIPrintf(viewer, "Iteratively computed extreme singular values: max %g min %g max/min %g\n",(double)smax,(double)smin,(double)(smax/smin));}
579: return(0);
580: }
582: static PetscErrorCode KSPViewFinalResidual_Internal(KSP ksp, PetscViewer viewer, PetscViewerFormat format)
583: {
584: PetscBool isascii;
588: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &isascii);
589: if (ksp->dscale && !ksp->dscalefix) SETERRQ(PetscObjectComm((PetscObject) ksp), PETSC_ERR_ARG_WRONGSTATE, "Cannot compute final scale with -ksp_diagonal_scale except also with -ksp_diagonal_scale_fix");
590: if (isascii) {
591: Mat A;
592: Vec t;
593: PetscReal norm;
595: PCGetOperators(ksp->pc, &A, NULL);
596: VecDuplicate(ksp->vec_rhs, &t);
597: KSP_MatMult(ksp, A, ksp->vec_sol, t);
598: VecAYPX(t, -1.0, ksp->vec_rhs);
599: VecNorm(t, NORM_2, &norm);
600: VecDestroy(&t);
601: PetscViewerASCIIPrintf(viewer, "KSP final norm of residual %g\n", (double) norm);
602: }
603: return(0);
604: }
606: static PetscErrorCode KSPSolve_Private(KSP ksp,Vec b,Vec x)
607: {
609: PetscBool flg = PETSC_FALSE,inXisinB=PETSC_FALSE,guess_zero;
610: Mat mat,pmat;
611: MPI_Comm comm;
612: MatNullSpace nullsp;
613: Vec btmp,vec_rhs=NULL;
616: comm = PetscObjectComm((PetscObject)ksp);
617: if (x && x == b) {
618: if (!ksp->guess_zero) SETERRQ(comm,PETSC_ERR_ARG_INCOMP,"Cannot use x == b with nonzero initial guess");
619: VecDuplicate(b,&x);
620: inXisinB = PETSC_TRUE;
621: }
622: if (b) {
623: PetscObjectReference((PetscObject)b);
624: VecDestroy(&ksp->vec_rhs);
625: ksp->vec_rhs = b;
626: }
627: if (x) {
628: PetscObjectReference((PetscObject)x);
629: VecDestroy(&ksp->vec_sol);
630: ksp->vec_sol = x;
631: }
633: if (ksp->viewPre) {ObjectView((PetscObject) ksp, ksp->viewerPre, ksp->formatPre);}
635: if (ksp->presolve) {(*ksp->presolve)(ksp,ksp->vec_rhs,ksp->vec_sol,ksp->prectx);}
637: /* reset the residual history list if requested */
638: if (ksp->res_hist_reset) ksp->res_hist_len = 0;
640: PetscLogEventBegin(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
642: if (ksp->guess) {
643: PetscObjectState ostate,state;
645: KSPGuessSetUp(ksp->guess);
646: PetscObjectStateGet((PetscObject)ksp->vec_sol,&ostate);
647: KSPGuessFormGuess(ksp->guess,ksp->vec_rhs,ksp->vec_sol);
648: PetscObjectStateGet((PetscObject)ksp->vec_sol,&state);
649: if (state != ostate) {
650: ksp->guess_zero = PETSC_FALSE;
651: } else {
652: PetscInfo(ksp,"Using zero initial guess since the KSPGuess object did not change the vector\n");
653: ksp->guess_zero = PETSC_TRUE;
654: }
655: }
657: /* KSPSetUp() scales the matrix if needed */
658: KSPSetUp(ksp);
659: KSPSetUpOnBlocks(ksp);
661: VecSetErrorIfLocked(ksp->vec_sol,3);
663: PCGetOperators(ksp->pc,&mat,&pmat);
664: /* diagonal scale RHS if called for */
665: if (ksp->dscale) {
666: VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
667: /* second time in, but matrix was scaled back to original */
668: if (ksp->dscalefix && ksp->dscalefix2) {
669: Mat mat,pmat;
671: PCGetOperators(ksp->pc,&mat,&pmat);
672: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
673: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
674: }
676: /* scale initial guess */
677: if (!ksp->guess_zero) {
678: if (!ksp->truediagonal) {
679: VecDuplicate(ksp->diagonal,&ksp->truediagonal);
680: VecCopy(ksp->diagonal,ksp->truediagonal);
681: VecReciprocal(ksp->truediagonal);
682: }
683: VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->truediagonal);
684: }
685: }
686: PCPreSolve(ksp->pc,ksp);
688: if (ksp->guess_zero) { VecSet(ksp->vec_sol,0.0);}
689: if (ksp->guess_knoll) { /* The Knoll trick is independent on the KSPGuess specified */
690: PCApply(ksp->pc,ksp->vec_rhs,ksp->vec_sol);
691: KSP_RemoveNullSpace(ksp,ksp->vec_sol);
692: ksp->guess_zero = PETSC_FALSE;
693: }
695: /* can we mark the initial guess as zero for this solve? */
696: guess_zero = ksp->guess_zero;
697: if (!ksp->guess_zero) {
698: PetscReal norm;
700: VecNormAvailable(ksp->vec_sol,NORM_2,&flg,&norm);
701: if (flg && !norm) ksp->guess_zero = PETSC_TRUE;
702: }
703: if (ksp->transpose_solve) {
704: MatGetNullSpace(pmat,&nullsp);
705: } else {
706: MatGetTransposeNullSpace(pmat,&nullsp);
707: }
708: if (nullsp) {
709: VecDuplicate(ksp->vec_rhs,&btmp);
710: VecCopy(ksp->vec_rhs,btmp);
711: MatNullSpaceRemove(nullsp,btmp);
712: vec_rhs = ksp->vec_rhs;
713: ksp->vec_rhs = btmp;
714: }
715: VecLockReadPush(ksp->vec_rhs);
716: if (ksp->reason == KSP_DIVERGED_PC_FAILED) {
717: VecSetInf(ksp->vec_sol);
718: }
719: (*ksp->ops->solve)(ksp);
721: VecLockReadPop(ksp->vec_rhs);
722: if (nullsp) {
723: ksp->vec_rhs = vec_rhs;
724: VecDestroy(&btmp);
725: }
727: ksp->guess_zero = guess_zero;
729: if (!ksp->reason) SETERRQ(comm,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
730: ksp->totalits += ksp->its;
732: if (ksp->viewReason) {
733: PetscViewerPushFormat(ksp->viewerReason,ksp->formatReason);
734: KSPConvergedReasonView(ksp, ksp->viewerReason);
735: PetscViewerPopFormat(ksp->viewerReason);
736: }
737: PCPostSolve(ksp->pc,ksp);
739: /* diagonal scale solution if called for */
740: if (ksp->dscale) {
741: VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->diagonal);
742: /* unscale right hand side and matrix */
743: if (ksp->dscalefix) {
744: Mat mat,pmat;
746: VecReciprocal(ksp->diagonal);
747: VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
748: PCGetOperators(ksp->pc,&mat,&pmat);
749: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
750: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
751: VecReciprocal(ksp->diagonal);
752: ksp->dscalefix2 = PETSC_TRUE;
753: }
754: }
755: PetscLogEventEnd(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
756: if (ksp->guess) {
757: KSPGuessUpdate(ksp->guess,ksp->vec_rhs,ksp->vec_sol);
758: }
759: if (ksp->postsolve) {
760: (*ksp->postsolve)(ksp,ksp->vec_rhs,ksp->vec_sol,ksp->postctx);
761: }
763: PCGetOperators(ksp->pc,&mat,&pmat);
764: if (ksp->viewEV) {KSPViewEigenvalues_Internal(ksp, PETSC_FALSE, ksp->viewerEV, ksp->formatEV);}
765: if (ksp->viewEVExp) {KSPViewEigenvalues_Internal(ksp, PETSC_TRUE, ksp->viewerEVExp, ksp->formatEVExp);}
766: if (ksp->viewSV) {KSPViewSingularvalues_Internal(ksp, ksp->viewerSV, ksp->formatSV);}
767: if (ksp->viewFinalRes) {KSPViewFinalResidual_Internal(ksp, ksp->viewerFinalRes, ksp->formatFinalRes);}
768: if (ksp->viewMat) {ObjectView((PetscObject) mat, ksp->viewerMat, ksp->formatMat);}
769: if (ksp->viewPMat) {ObjectView((PetscObject) pmat, ksp->viewerPMat, ksp->formatPMat);}
770: if (ksp->viewRhs) {ObjectView((PetscObject) ksp->vec_rhs, ksp->viewerRhs, ksp->formatRhs);}
771: if (ksp->viewSol) {ObjectView((PetscObject) ksp->vec_sol, ksp->viewerSol, ksp->formatSol);}
772: if (ksp->view) {ObjectView((PetscObject) ksp, ksp->viewer, ksp->format);}
773: if (ksp->viewDScale) {ObjectView((PetscObject) ksp->diagonal, ksp->viewerDScale, ksp->formatDScale);}
774: if (ksp->viewMatExp) {
775: Mat A, B;
777: PCGetOperators(ksp->pc, &A, NULL);
778: if (ksp->transpose_solve) {
779: Mat AT;
781: MatCreateTranspose(A, &AT);
782: MatComputeOperator(AT, MATAIJ, &B);
783: MatDestroy(&AT);
784: } else {
785: MatComputeOperator(A, MATAIJ, &B);
786: }
787: ObjectView((PetscObject) B, ksp->viewerMatExp, ksp->formatMatExp);
788: MatDestroy(&B);
789: }
790: if (ksp->viewPOpExp) {
791: Mat B;
793: KSPComputeOperator(ksp, MATAIJ, &B);
794: ObjectView((PetscObject) B, ksp->viewerPOpExp, ksp->formatPOpExp);
795: MatDestroy(&B);
796: }
798: if (inXisinB) {
799: VecCopy(x,b);
800: VecDestroy(&x);
801: }
802: PetscObjectSAWsBlock((PetscObject)ksp);
803: if (ksp->errorifnotconverged && ksp->reason < 0 && ksp->reason != KSP_DIVERGED_ITS) {
804: if (ksp->reason == KSP_DIVERGED_PC_FAILED) {
805: PCFailedReason reason;
806: PCGetFailedReason(ksp->pc,&reason);
807: SETERRQ2(comm,PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged, reason %s PC failed due to %s",KSPConvergedReasons[ksp->reason],PCFailedReasons[reason]);
808: } else SETERRQ1(comm,PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged, reason %s",KSPConvergedReasons[ksp->reason]);
809: }
810: return(0);
811: }
813: /*@
814: KSPSolve - Solves linear system.
816: Collective on ksp
818: Parameters:
819: + ksp - iterative context obtained from KSPCreate()
820: . b - the right hand side vector
821: - x - the solution (this may be the same vector as b, then b will be overwritten with answer)
823: Options Database Keys:
824: + -ksp_view_eigenvalues - compute preconditioned operators eigenvalues
825: . -ksp_view_eigenvalues_explicit - compute the eigenvalues by forming the dense operator and using LAPACK
826: . -ksp_view_mat binary - save matrix to the default binary viewer
827: . -ksp_view_pmat binary - save matrix used to build preconditioner to the default binary viewer
828: . -ksp_view_rhs binary - save right hand side vector to the default binary viewer
829: . -ksp_view_solution binary - save computed solution vector to the default binary viewer
830: (can be read later with src/ksp/tutorials/ex10.c for testing solvers)
831: . -ksp_view_mat_explicit - for matrix-free operators, computes the matrix entries and views them
832: . -ksp_view_preconditioned_operator_explicit - computes the product of the preconditioner and matrix as an explicit matrix and views it
833: . -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations
834: . -ksp_view_final_residual - print 2-norm of true linear system residual at the end of the solution process
835: - -ksp_view - print the ksp data structure at the end of the system solution
837: Notes:
839: If one uses KSPSetDM() then x or b need not be passed. Use KSPGetSolution() to access the solution in this case.
841: The operator is specified with KSPSetOperators().
843: Call KSPGetConvergedReason() to determine if the solver converged or failed and
844: why. The number of iterations can be obtained from KSPGetIterationNumber().
846: If you provide a matrix that has a MatSetNullSpace() and MatSetTransposeNullSpace() this will use that information to solve singular systems
847: in the least squares sense with a norm minimizing solution.
848: $
849: $ A x = b where b = b_p + b_t where b_t is not in the range of A (and hence by the fundamental theorem of linear algebra is in the nullspace(A') see MatSetNullSpace()
850: $
851: $ KSP first removes b_t producing the linear system A x = b_p (which has multiple solutions) and solves this to find the ||x|| minimizing solution (and hence
852: $ it finds the solution x orthogonal to the nullspace(A). The algorithm is simply in each iteration of the Krylov method we remove the nullspace(A) from the search
853: $ direction thus the solution which is a linear combination of the search directions has no component in the nullspace(A).
854: $
855: $ We recommend always using GMRES for such singular systems.
856: $ If nullspace(A) = nullspace(A') (note symmetric matrices always satisfy this property) then both left and right preconditioning will work
857: $ If nullspace(A) != nullspace(A') then left preconditioning will work but right preconditioning may not work (or it may).
859: Developer Note: The reason we cannot always solve nullspace(A) != nullspace(A') systems with right preconditioning is because we need to remove at each iteration
860: the nullspace(AB) from the search direction. While we know the nullspace(A) the nullspace(AB) equals B^-1 times the nullspace(A) but except for trivial preconditioners
861: such as diagonal scaling we cannot apply the inverse of the preconditioner to a vector and thus cannot compute the nullspace(AB).
864: If using a direct method (e.g., via the KSP solver
865: KSPPREONLY and a preconditioner such as PCLU/PCILU),
866: then its=1. See KSPSetTolerances() and KSPConvergedDefault()
867: for more details.
869: Understanding Convergence:
870: The routines KSPMonitorSet(), KSPComputeEigenvalues(), and
871: KSPComputeEigenvaluesExplicitly() provide information on additional
872: options to monitor convergence and print eigenvalue information.
874: Level: beginner
876: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
877: KSPSolveTranspose(), KSPGetIterationNumber(), MatNullSpaceCreate(), MatSetNullSpace(), MatSetTransposeNullSpace(), KSP,
878: KSPConvergedReasonView()
879: @*/
880: PetscErrorCode KSPSolve(KSP ksp,Vec b,Vec x)
881: {
888: ksp->transpose_solve = PETSC_FALSE;
889: KSPSolve_Private(ksp,b,x);
890: return(0);
891: }
893: /*@
894: KSPSolveTranspose - Solves the transpose of a linear system.
896: Collective on ksp
898: Input Parameters:
899: + ksp - iterative context obtained from KSPCreate()
900: . b - right hand side vector
901: - x - solution vector
903: Notes:
904: For complex numbers this solve the non-Hermitian transpose system.
906: Developer Notes:
907: We need to implement a KSPSolveHermitianTranspose()
909: Level: developer
911: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
912: KSPSolve(), KSP
913: @*/
914: PetscErrorCode KSPSolveTranspose(KSP ksp,Vec b,Vec x)
915: {
922: ksp->transpose_solve = PETSC_TRUE;
923: KSPSolve_Private(ksp,b,x);
924: return(0);
925: }
927: static PetscErrorCode KSPViewFinalMatResidual_Internal(KSP ksp, Mat B, Mat X, PetscViewer viewer, PetscViewerFormat format, PetscInt shift)
928: {
929: Mat A, R;
930: PetscReal *norms;
931: PetscInt i, N;
932: PetscBool flg;
936: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &flg);
937: if (flg) {
938: PCGetOperators(ksp->pc, &A, NULL);
939: MatMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &R);
940: MatAYPX(R, -1.0, B, SAME_NONZERO_PATTERN);
941: MatGetSize(R, NULL, &N);
942: PetscMalloc1(N, &norms);
943: MatGetColumnNorms(R, NORM_2, norms);
944: MatDestroy(&R);
945: for (i = 0; i < N; ++i) {
946: PetscViewerASCIIPrintf(viewer, "%s #%D %g\n", i == 0 ? "KSP final norm of residual" : " ", shift + i, (double)norms[i]);
947: }
948: PetscFree(norms);
949: }
950: return(0);
951: }
953: /*@
954: KSPMatSolve - Solves a linear system with multiple right-hand sides stored as a MATDENSE. Unlike KSPSolve(), B and X must be different matrices.
956: Input Parameters:
957: + ksp - iterative context
958: - B - block of right-hand sides
960: Output Parameter:
961: . X - block of solutions
963: Notes:
964: This is a stripped-down version of KSPSolve(), which only handles -ksp_view, -ksp_converged_reason, and -ksp_view_final_residual.
966: Level: intermediate
968: .seealso: KSPSolve(), MatMatSolve(), MATDENSE, KSPHPDDM, PCBJACOBI, PCASM
969: @*/
970: PetscErrorCode KSPMatSolve(KSP ksp, Mat B, Mat X)
971: {
972: Mat A, vB, vX;
973: Vec cb, cx;
974: PetscInt m1, M1, m2, M2, n1, N1, n2, N2, Bbn = PETSC_DECIDE;
975: PetscBool match;
984: if (!B->assembled) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
985: MatCheckPreallocated(X, 3);
986: if (!X->assembled) {
987: MatSetOption(X, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
988: MatAssemblyBegin(X, MAT_FINAL_ASSEMBLY);
989: MatAssemblyEnd(X, MAT_FINAL_ASSEMBLY);
990: }
991: if (B == X) SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_IDN, "B and X must be different matrices");
992: KSPGetOperators(ksp, &A, NULL);
993: MatGetLocalSize(A, &m1, NULL);
994: MatGetLocalSize(B, &m2, &n2);
995: MatGetSize(A, &M1, NULL);
996: MatGetSize(B, &M2, &N2);
997: if (m1 != m2 || M1 != M2) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot use a block of right-hand sides with (m2,M2) = (%D,%D) for a linear system with (m1,M1) = (%D,%D)", m2, M2, m1, M1);
998: MatGetLocalSize(X, &m1, &n1);
999: MatGetSize(X, &M1, &N1);
1000: if (m1 != m2 || M1 != M2 || n1 != n2 || N1 != N2) SETERRQ8(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible block of right-hand sides (m2,M2)x(n2,N2) = (%D,%D)x(%D,%D) and solutions (m1,M1)x(n1,N1) = (%D,%D)x(%D,%D)", m2, M2, n2, N2, m1, M1, n1, N1);
1001: PetscObjectBaseTypeCompareAny((PetscObject)B, &match, MATSEQDENSE, MATMPIDENSE, "");
1002: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of right-hand sides not stored in a dense Mat");
1003: PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, "");
1004: if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of solutions not stored in a dense Mat");
1005: KSPSetUp(ksp);
1006: KSPSetUpOnBlocks(ksp);
1007: if (ksp->ops->matsolve) {
1008: if (ksp->guess_zero) {
1009: MatZeroEntries(X);
1010: }
1011: PetscLogEventBegin(KSP_MatSolve, ksp, B, X, 0);
1012: KSPGetMatSolveBlockSize(ksp, &Bbn);
1013: /* by default, do a single solve with all columns */
1014: if (Bbn == PETSC_DECIDE) Bbn = N2;
1015: else if (Bbn < 1) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "KSPMatSolve() block size %D must be positive", Bbn);
1016: PetscInfo2(ksp, "KSP type %s solving using blocks of width at most %D\n", ((PetscObject)ksp)->type_name, Bbn);
1017: /* if -ksp_matsolve_block_size is greater than the actual number of columns, do a single solve with all columns */
1018: if (Bbn >= N2) {
1019: (*ksp->ops->matsolve)(ksp, B, X);
1020: if (ksp->viewFinalRes) {
1021: KSPViewFinalMatResidual_Internal(ksp, B, X, ksp->viewerFinalRes, ksp->formatFinalRes, 0);
1022: }
1023: if (ksp->viewReason) {
1024: PetscViewerPushFormat(ksp->viewerReason,PETSC_VIEWER_DEFAULT);
1025: KSPConvergedReasonView(ksp, ksp->viewerReason);
1026: PetscViewerPopFormat(ksp->viewerReason);
1027: }
1028: } else {
1029: for (n2 = 0; n2 < N2; n2 += Bbn) {
1030: MatDenseGetSubMatrix(B, n2, PetscMin(n2+Bbn, N2), &vB);
1031: MatDenseGetSubMatrix(X, n2, PetscMin(n2+Bbn, N2), &vX);
1032: (*ksp->ops->matsolve)(ksp, vB, vX);
1033: if (ksp->viewFinalRes) {
1034: KSPViewFinalMatResidual_Internal(ksp, vB, vX, ksp->viewerFinalRes, ksp->formatFinalRes, n2);
1035: }
1036: if (ksp->viewReason) {
1037: PetscViewerPushFormat(ksp->viewerReason,PETSC_VIEWER_DEFAULT);
1038: KSPConvergedReasonView(ksp, ksp->viewerReason);
1039: PetscViewerPopFormat(ksp->viewerReason);
1040: }
1041: MatDenseRestoreSubMatrix(B, &vB);
1042: MatDenseRestoreSubMatrix(X, &vX);
1043: }
1044: }
1045: if (ksp->view) {
1046: KSPView(ksp, ksp->viewer);
1047: }
1048: PetscLogEventEnd(KSP_MatSolve, ksp, B, X, 0);
1049: } else {
1050: PetscInfo1(ksp, "KSP type %s solving column by column\n", ((PetscObject)ksp)->type_name);
1051: for (n2 = 0; n2 < N2; ++n2) {
1052: MatDenseGetColumnVecRead(B, n2, &cb);
1053: MatDenseGetColumnVecWrite(X, n2, &cx);
1054: KSPSolve(ksp, cb, cx);
1055: MatDenseRestoreColumnVecWrite(X, n2, &cx);
1056: MatDenseRestoreColumnVecRead(B, n2, &cb);
1057: }
1058: }
1059: return(0);
1060: }
1062: /*@
1063: KSPSetMatSolveBlockSize - Sets the maximum number of columns treated simultaneously in KSPMatSolve().
1065: Logically collective
1067: Input Parameters:
1068: + ksp - iterative context
1069: - bs - block size
1071: Level: advanced
1073: .seealso: KSPMatSolve(), KSPGetMatSolveBlockSize(), -mat_mumps_icntl_27, -matmatmult_Bbn
1074: @*/
1075: PetscErrorCode KSPSetMatSolveBlockSize(KSP ksp, PetscInt bs)
1076: {
1082: PetscTryMethod(ksp, "KSPSetMatSolveBlockSize_C", (KSP, PetscInt), (ksp, bs));
1083: return(0);
1084: }
1086: /*@
1087: KSPGetMatSolveBlockSize - Gets the maximum number of columns treated simultaneously in KSPMatSolve().
1089: Input Parameter:
1090: . ksp - iterative context
1092: Output Parameter:
1093: . bs - block size
1095: Level: advanced
1097: .seealso: KSPMatSolve(), KSPSetMatSolveBlockSize(), -mat_mumps_icntl_27, -matmatmult_Bbn
1098: @*/
1099: PetscErrorCode KSPGetMatSolveBlockSize(KSP ksp, PetscInt *bs)
1100: {
1105: *bs = PETSC_DECIDE;
1106: PetscTryMethod(ksp, "KSPGetMatSolveBlockSize_C", (KSP, PetscInt*), (ksp, bs));
1107: return(0);
1108: }
1110: /*@
1111: KSPResetViewers - Resets all the viewers set from the options database during KSPSetFromOptions()
1113: Collective on ksp
1115: Input Parameter:
1116: . ksp - iterative context obtained from KSPCreate()
1118: Level: beginner
1120: .seealso: KSPCreate(), KSPSetUp(), KSPSolve(), KSPSetFromOptions(), KSP
1121: @*/
1122: PetscErrorCode KSPResetViewers(KSP ksp)
1123: {
1128: if (!ksp) return(0);
1129: PetscViewerDestroy(&ksp->viewer);
1130: PetscViewerDestroy(&ksp->viewerPre);
1131: PetscViewerDestroy(&ksp->viewerReason);
1132: PetscViewerDestroy(&ksp->viewerMat);
1133: PetscViewerDestroy(&ksp->viewerPMat);
1134: PetscViewerDestroy(&ksp->viewerRhs);
1135: PetscViewerDestroy(&ksp->viewerSol);
1136: PetscViewerDestroy(&ksp->viewerMatExp);
1137: PetscViewerDestroy(&ksp->viewerEV);
1138: PetscViewerDestroy(&ksp->viewerSV);
1139: PetscViewerDestroy(&ksp->viewerEVExp);
1140: PetscViewerDestroy(&ksp->viewerFinalRes);
1141: PetscViewerDestroy(&ksp->viewerPOpExp);
1142: PetscViewerDestroy(&ksp->viewerDScale);
1143: ksp->view = PETSC_FALSE;
1144: ksp->viewPre = PETSC_FALSE;
1145: ksp->viewReason = PETSC_FALSE;
1146: ksp->viewMat = PETSC_FALSE;
1147: ksp->viewPMat = PETSC_FALSE;
1148: ksp->viewRhs = PETSC_FALSE;
1149: ksp->viewSol = PETSC_FALSE;
1150: ksp->viewMatExp = PETSC_FALSE;
1151: ksp->viewEV = PETSC_FALSE;
1152: ksp->viewSV = PETSC_FALSE;
1153: ksp->viewEVExp = PETSC_FALSE;
1154: ksp->viewFinalRes = PETSC_FALSE;
1155: ksp->viewPOpExp = PETSC_FALSE;
1156: ksp->viewDScale = PETSC_FALSE;
1157: return(0);
1158: }
1160: /*@
1161: KSPReset - Resets a KSP context to the kspsetupcalled = 0 state and removes any allocated Vecs and Mats
1163: Collective on ksp
1165: Input Parameter:
1166: . ksp - iterative context obtained from KSPCreate()
1168: Level: beginner
1170: .seealso: KSPCreate(), KSPSetUp(), KSPSolve(), KSP
1171: @*/
1172: PetscErrorCode KSPReset(KSP ksp)
1173: {
1178: if (!ksp) return(0);
1179: if (ksp->ops->reset) {
1180: (*ksp->ops->reset)(ksp);
1181: }
1182: if (ksp->pc) {PCReset(ksp->pc);}
1183: if (ksp->guess) {
1184: KSPGuess guess = ksp->guess;
1185: if (guess->ops->reset) { (*guess->ops->reset)(guess); }
1186: }
1187: VecDestroyVecs(ksp->nwork,&ksp->work);
1188: VecDestroy(&ksp->vec_rhs);
1189: VecDestroy(&ksp->vec_sol);
1190: VecDestroy(&ksp->diagonal);
1191: VecDestroy(&ksp->truediagonal);
1193: KSPResetViewers(ksp);
1195: ksp->setupstage = KSP_SETUP_NEW;
1196: return(0);
1197: }
1199: /*@
1200: KSPDestroy - Destroys KSP context.
1202: Collective on ksp
1204: Input Parameter:
1205: . ksp - iterative context obtained from KSPCreate()
1207: Level: beginner
1209: .seealso: KSPCreate(), KSPSetUp(), KSPSolve(), KSP
1210: @*/
1211: PetscErrorCode KSPDestroy(KSP *ksp)
1212: {
1214: PC pc;
1217: if (!*ksp) return(0);
1219: if (--((PetscObject)(*ksp))->refct > 0) {*ksp = NULL; return(0);}
1221: PetscObjectSAWsViewOff((PetscObject)*ksp);
1223: /*
1224: Avoid a cascading call to PCReset(ksp->pc) from the following call:
1225: PCReset() shouldn't be called from KSPDestroy() as it is unprotected by pc's
1226: refcount (and may be shared, e.g., by other ksps).
1227: */
1228: pc = (*ksp)->pc;
1229: (*ksp)->pc = NULL;
1230: KSPReset((*ksp));
1231: (*ksp)->pc = pc;
1232: if ((*ksp)->ops->destroy) {(*(*ksp)->ops->destroy)(*ksp);}
1234: KSPGuessDestroy(&(*ksp)->guess);
1235: DMDestroy(&(*ksp)->dm);
1236: PCDestroy(&(*ksp)->pc);
1237: PetscFree((*ksp)->res_hist_alloc);
1238: if ((*ksp)->convergeddestroy) {
1239: (*(*ksp)->convergeddestroy)((*ksp)->cnvP);
1240: }
1241: KSPMonitorCancel((*ksp));
1242: PetscViewerDestroy(&(*ksp)->eigviewer);
1243: PetscHeaderDestroy(ksp);
1244: return(0);
1245: }
1247: /*@
1248: KSPSetPCSide - Sets the preconditioning side.
1250: Logically Collective on ksp
1252: Input Parameter:
1253: . ksp - iterative context obtained from KSPCreate()
1255: Output Parameter:
1256: . side - the preconditioning side, where side is one of
1257: .vb
1258: PC_LEFT - left preconditioning (default)
1259: PC_RIGHT - right preconditioning
1260: PC_SYMMETRIC - symmetric preconditioning
1261: .ve
1263: Options Database Keys:
1264: . -ksp_pc_side <right,left,symmetric>
1266: Notes:
1267: Left preconditioning is used by default for most Krylov methods except KSPFGMRES which only supports right preconditioning.
1269: For methods changing the side of the preconditioner changes the norm type that is used, see KSPSetNormType().
1271: Symmetric preconditioning is currently available only for the KSPQCG method. Note, however, that
1272: symmetric preconditioning can be emulated by using either right or left
1273: preconditioning and a pre or post processing step.
1275: Setting the PC side often affects the default norm type. See KSPSetNormType() for details.
1277: Level: intermediate
1279: .seealso: KSPGetPCSide(), KSPSetNormType(), KSPGetNormType(), KSP
1280: @*/
1281: PetscErrorCode KSPSetPCSide(KSP ksp,PCSide side)
1282: {
1286: ksp->pc_side = ksp->pc_side_set = side;
1287: return(0);
1288: }
1290: /*@
1291: KSPGetPCSide - Gets the preconditioning side.
1293: Not Collective
1295: Input Parameter:
1296: . ksp - iterative context obtained from KSPCreate()
1298: Output Parameter:
1299: . side - the preconditioning side, where side is one of
1300: .vb
1301: PC_LEFT - left preconditioning (default)
1302: PC_RIGHT - right preconditioning
1303: PC_SYMMETRIC - symmetric preconditioning
1304: .ve
1306: Level: intermediate
1308: .seealso: KSPSetPCSide(), KSP
1309: @*/
1310: PetscErrorCode KSPGetPCSide(KSP ksp,PCSide *side)
1311: {
1317: KSPSetUpNorms_Private(ksp,PETSC_TRUE,&ksp->normtype,&ksp->pc_side);
1318: *side = ksp->pc_side;
1319: return(0);
1320: }
1322: /*@
1323: KSPGetTolerances - Gets the relative, absolute, divergence, and maximum
1324: iteration tolerances used by the default KSP convergence tests.
1326: Not Collective
1328: Input Parameter:
1329: . ksp - the Krylov subspace context
1331: Output Parameters:
1332: + rtol - the relative convergence tolerance
1333: . abstol - the absolute convergence tolerance
1334: . dtol - the divergence tolerance
1335: - maxits - maximum number of iterations
1337: Notes:
1338: The user can specify NULL for any parameter that is not needed.
1340: Level: intermediate
1342: maximum, iterations
1344: .seealso: KSPSetTolerances(), KSP
1345: @*/
1346: PetscErrorCode KSPGetTolerances(KSP ksp,PetscReal *rtol,PetscReal *abstol,PetscReal *dtol,PetscInt *maxits)
1347: {
1350: if (abstol) *abstol = ksp->abstol;
1351: if (rtol) *rtol = ksp->rtol;
1352: if (dtol) *dtol = ksp->divtol;
1353: if (maxits) *maxits = ksp->max_it;
1354: return(0);
1355: }
1357: /*@
1358: KSPSetTolerances - Sets the relative, absolute, divergence, and maximum
1359: iteration tolerances used by the default KSP convergence testers.
1361: Logically Collective on ksp
1363: Input Parameters:
1364: + ksp - the Krylov subspace context
1365: . rtol - the relative convergence tolerance, relative decrease in the (possibly preconditioned) residual norm
1366: . abstol - the absolute convergence tolerance absolute size of the (possibly preconditioned) residual norm
1367: . dtol - the divergence tolerance, amount (possibly preconditioned) residual norm can increase before KSPConvergedDefault() concludes that the method is diverging
1368: - maxits - maximum number of iterations to use
1370: Options Database Keys:
1371: + -ksp_atol <abstol> - Sets abstol
1372: . -ksp_rtol <rtol> - Sets rtol
1373: . -ksp_divtol <dtol> - Sets dtol
1374: - -ksp_max_it <maxits> - Sets maxits
1376: Notes:
1377: Use PETSC_DEFAULT to retain the default value of any of the tolerances.
1379: See KSPConvergedDefault() for details how these parameters are used in the default convergence test. See also KSPSetConvergenceTest()
1380: for setting user-defined stopping criteria.
1382: Level: intermediate
1384: convergence, maximum, iterations
1386: .seealso: KSPGetTolerances(), KSPConvergedDefault(), KSPSetConvergenceTest(), KSP
1387: @*/
1388: PetscErrorCode KSPSetTolerances(KSP ksp,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt maxits)
1389: {
1397: if (rtol != PETSC_DEFAULT) {
1398: if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %g must be non-negative and less than 1.0",(double)rtol);
1399: ksp->rtol = rtol;
1400: }
1401: if (abstol != PETSC_DEFAULT) {
1402: if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol);
1403: ksp->abstol = abstol;
1404: }
1405: if (dtol != PETSC_DEFAULT) {
1406: if (dtol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Divergence tolerance %g must be larger than 1.0",(double)dtol);
1407: ksp->divtol = dtol;
1408: }
1409: if (maxits != PETSC_DEFAULT) {
1410: if (maxits < 0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxits);
1411: ksp->max_it = maxits;
1412: }
1413: return(0);
1414: }
1416: /*@
1417: KSPSetInitialGuessNonzero - Tells the iterative solver that the
1418: initial guess is nonzero; otherwise KSP assumes the initial guess
1419: is to be zero (and thus zeros it out before solving).
1421: Logically Collective on ksp
1423: Input Parameters:
1424: + ksp - iterative context obtained from KSPCreate()
1425: - flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero
1427: Options database keys:
1428: . -ksp_initial_guess_nonzero : use nonzero initial guess; this takes an optional truth value (0/1/no/yes/true/false)
1430: Level: beginner
1432: Notes:
1433: If this is not called the X vector is zeroed in the call to KSPSolve().
1435: .seealso: KSPGetInitialGuessNonzero(), KSPSetGuessType(), KSPGuessType, KSP
1436: @*/
1437: PetscErrorCode KSPSetInitialGuessNonzero(KSP ksp,PetscBool flg)
1438: {
1442: ksp->guess_zero = (PetscBool) !(int)flg;
1443: return(0);
1444: }
1446: /*@
1447: KSPGetInitialGuessNonzero - Determines whether the KSP solver is using
1448: a zero initial guess.
1450: Not Collective
1452: Input Parameter:
1453: . ksp - iterative context obtained from KSPCreate()
1455: Output Parameter:
1456: . flag - PETSC_TRUE if guess is nonzero, else PETSC_FALSE
1458: Level: intermediate
1460: .seealso: KSPSetInitialGuessNonzero(), KSP
1461: @*/
1462: PetscErrorCode KSPGetInitialGuessNonzero(KSP ksp,PetscBool *flag)
1463: {
1467: if (ksp->guess_zero) *flag = PETSC_FALSE;
1468: else *flag = PETSC_TRUE;
1469: return(0);
1470: }
1472: /*@
1473: KSPSetErrorIfNotConverged - Causes KSPSolve() to generate an error if the solver has not converged.
1475: Logically Collective on ksp
1477: Input Parameters:
1478: + ksp - iterative context obtained from KSPCreate()
1479: - flg - PETSC_TRUE indicates you want the error generated
1481: Options database keys:
1482: . -ksp_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false)
1484: Level: intermediate
1486: Notes:
1487: Normally PETSc continues if a linear solver fails to converge, you can call KSPGetConvergedReason() after a KSPSolve()
1488: to determine if it has converged.
1491: .seealso: KSPGetErrorIfNotConverged(), KSP
1492: @*/
1493: PetscErrorCode KSPSetErrorIfNotConverged(KSP ksp,PetscBool flg)
1494: {
1498: ksp->errorifnotconverged = flg;
1499: return(0);
1500: }
1502: /*@
1503: KSPGetErrorIfNotConverged - Will KSPSolve() generate an error if the solver does not converge?
1505: Not Collective
1507: Input Parameter:
1508: . ksp - iterative context obtained from KSPCreate()
1510: Output Parameter:
1511: . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE
1513: Level: intermediate
1515: .seealso: KSPSetErrorIfNotConverged(), KSP
1516: @*/
1517: PetscErrorCode KSPGetErrorIfNotConverged(KSP ksp,PetscBool *flag)
1518: {
1522: *flag = ksp->errorifnotconverged;
1523: return(0);
1524: }
1526: /*@
1527: KSPSetInitialGuessKnoll - Tells the iterative solver to use PCApply(pc,b,..) to compute the initial guess (The Knoll trick)
1529: Logically Collective on ksp
1531: Input Parameters:
1532: + ksp - iterative context obtained from KSPCreate()
1533: - flg - PETSC_TRUE or PETSC_FALSE
1535: Level: advanced
1537: Developer Note: the Knoll trick is not currently implemented using the KSPGuess class
1539: .seealso: KSPGetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero(), KSP
1540: @*/
1541: PetscErrorCode KSPSetInitialGuessKnoll(KSP ksp,PetscBool flg)
1542: {
1546: ksp->guess_knoll = flg;
1547: return(0);
1548: }
1550: /*@
1551: KSPGetInitialGuessKnoll - Determines whether the KSP solver is using the Knoll trick (using PCApply(pc,b,...) to compute
1552: the initial guess
1554: Not Collective
1556: Input Parameter:
1557: . ksp - iterative context obtained from KSPCreate()
1559: Output Parameter:
1560: . flag - PETSC_TRUE if using Knoll trick, else PETSC_FALSE
1562: Level: advanced
1564: .seealso: KSPSetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero(), KSP
1565: @*/
1566: PetscErrorCode KSPGetInitialGuessKnoll(KSP ksp,PetscBool *flag)
1567: {
1571: *flag = ksp->guess_knoll;
1572: return(0);
1573: }
1575: /*@
1576: KSPGetComputeSingularValues - Gets the flag indicating whether the extreme singular
1577: values will be calculated via a Lanczos or Arnoldi process as the linear
1578: system is solved.
1580: Not Collective
1582: Input Parameter:
1583: . ksp - iterative context obtained from KSPCreate()
1585: Output Parameter:
1586: . flg - PETSC_TRUE or PETSC_FALSE
1588: Options Database Key:
1589: . -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()
1591: Notes:
1592: Currently this option is not valid for all iterative methods.
1594: Many users may just want to use the monitoring routine
1595: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1596: to print the singular values at each iteration of the linear solve.
1598: Level: advanced
1600: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue(), KSP
1601: @*/
1602: PetscErrorCode KSPGetComputeSingularValues(KSP ksp,PetscBool *flg)
1603: {
1607: *flg = ksp->calc_sings;
1608: return(0);
1609: }
1611: /*@
1612: KSPSetComputeSingularValues - Sets a flag so that the extreme singular
1613: values will be calculated via a Lanczos or Arnoldi process as the linear
1614: system is solved.
1616: Logically Collective on ksp
1618: Input Parameters:
1619: + ksp - iterative context obtained from KSPCreate()
1620: - flg - PETSC_TRUE or PETSC_FALSE
1622: Options Database Key:
1623: . -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()
1625: Notes:
1626: Currently this option is not valid for all iterative methods.
1628: Many users may just want to use the monitoring routine
1629: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1630: to print the singular values at each iteration of the linear solve.
1632: Level: advanced
1634: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue(), KSP
1635: @*/
1636: PetscErrorCode KSPSetComputeSingularValues(KSP ksp,PetscBool flg)
1637: {
1641: ksp->calc_sings = flg;
1642: return(0);
1643: }
1645: /*@
1646: KSPGetComputeEigenvalues - Gets the flag indicating that the extreme eigenvalues
1647: values will be calculated via a Lanczos or Arnoldi process as the linear
1648: system is solved.
1650: Not Collective
1652: Input Parameter:
1653: . ksp - iterative context obtained from KSPCreate()
1655: Output Parameter:
1656: . flg - PETSC_TRUE or PETSC_FALSE
1658: Notes:
1659: Currently this option is not valid for all iterative methods.
1661: Level: advanced
1663: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly(), KSP
1664: @*/
1665: PetscErrorCode KSPGetComputeEigenvalues(KSP ksp,PetscBool *flg)
1666: {
1670: *flg = ksp->calc_sings;
1671: return(0);
1672: }
1674: /*@
1675: KSPSetComputeEigenvalues - Sets a flag so that the extreme eigenvalues
1676: values will be calculated via a Lanczos or Arnoldi process as the linear
1677: system is solved.
1679: Logically Collective on ksp
1681: Input Parameters:
1682: + ksp - iterative context obtained from KSPCreate()
1683: - flg - PETSC_TRUE or PETSC_FALSE
1685: Notes:
1686: Currently this option is not valid for all iterative methods.
1688: Level: advanced
1690: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly(), KSP
1691: @*/
1692: PetscErrorCode KSPSetComputeEigenvalues(KSP ksp,PetscBool flg)
1693: {
1697: ksp->calc_sings = flg;
1698: return(0);
1699: }
1701: /*@
1702: KSPSetComputeRitz - Sets a flag so that the Ritz or harmonic Ritz pairs
1703: will be calculated via a Lanczos or Arnoldi process as the linear
1704: system is solved.
1706: Logically Collective on ksp
1708: Input Parameters:
1709: + ksp - iterative context obtained from KSPCreate()
1710: - flg - PETSC_TRUE or PETSC_FALSE
1712: Notes:
1713: Currently this option is only valid for the GMRES method.
1715: Level: advanced
1717: .seealso: KSPComputeRitz(), KSP
1718: @*/
1719: PetscErrorCode KSPSetComputeRitz(KSP ksp, PetscBool flg)
1720: {
1724: ksp->calc_ritz = flg;
1725: return(0);
1726: }
1728: /*@
1729: KSPGetRhs - Gets the right-hand-side vector for the linear system to
1730: be solved.
1732: Not Collective
1734: Input Parameter:
1735: . ksp - iterative context obtained from KSPCreate()
1737: Output Parameter:
1738: . r - right-hand-side vector
1740: Level: developer
1742: .seealso: KSPGetSolution(), KSPSolve(), KSP
1743: @*/
1744: PetscErrorCode KSPGetRhs(KSP ksp,Vec *r)
1745: {
1749: *r = ksp->vec_rhs;
1750: return(0);
1751: }
1753: /*@
1754: KSPGetSolution - Gets the location of the solution for the
1755: linear system to be solved. Note that this may not be where the solution
1756: is stored during the iterative process; see KSPBuildSolution().
1758: Not Collective
1760: Input Parameters:
1761: . ksp - iterative context obtained from KSPCreate()
1763: Output Parameters:
1764: . v - solution vector
1766: Level: developer
1768: .seealso: KSPGetRhs(), KSPBuildSolution(), KSPSolve(), KSP
1769: @*/
1770: PetscErrorCode KSPGetSolution(KSP ksp,Vec *v)
1771: {
1775: *v = ksp->vec_sol;
1776: return(0);
1777: }
1779: /*@
1780: KSPSetPC - Sets the preconditioner to be used to calculate the
1781: application of the preconditioner on a vector.
1783: Collective on ksp
1785: Input Parameters:
1786: + ksp - iterative context obtained from KSPCreate()
1787: - pc - the preconditioner object (can be NULL)
1789: Notes:
1790: Use KSPGetPC() to retrieve the preconditioner context.
1792: Level: developer
1794: .seealso: KSPGetPC(), KSP
1795: @*/
1796: PetscErrorCode KSPSetPC(KSP ksp,PC pc)
1797: {
1802: if (pc) {
1805: }
1806: PetscObjectReference((PetscObject)pc);
1807: PCDestroy(&ksp->pc);
1808: ksp->pc = pc;
1809: PetscLogObjectParent((PetscObject)ksp,(PetscObject)ksp->pc);
1810: return(0);
1811: }
1813: /*@
1814: KSPGetPC - Returns a pointer to the preconditioner context
1815: set with KSPSetPC().
1817: Not Collective
1819: Input Parameters:
1820: . ksp - iterative context obtained from KSPCreate()
1822: Output Parameter:
1823: . pc - preconditioner context
1825: Level: developer
1827: .seealso: KSPSetPC(), KSP
1828: @*/
1829: PetscErrorCode KSPGetPC(KSP ksp,PC *pc)
1830: {
1836: if (!ksp->pc) {
1837: PCCreate(PetscObjectComm((PetscObject)ksp),&ksp->pc);
1838: PetscObjectIncrementTabLevel((PetscObject)ksp->pc,(PetscObject)ksp,0);
1839: PetscLogObjectParent((PetscObject)ksp,(PetscObject)ksp->pc);
1840: PetscObjectSetOptions((PetscObject)ksp->pc,((PetscObject)ksp)->options);
1841: }
1842: *pc = ksp->pc;
1843: return(0);
1844: }
1846: /*@
1847: KSPMonitor - runs the user provided monitor routines, if they exist
1849: Collective on ksp
1851: Input Parameters:
1852: + ksp - iterative context obtained from KSPCreate()
1853: . it - iteration number
1854: - rnorm - relative norm of the residual
1856: Notes:
1857: This routine is called by the KSP implementations.
1858: It does not typically need to be called by the user.
1860: Level: developer
1862: .seealso: KSPMonitorSet()
1863: @*/
1864: PetscErrorCode KSPMonitor(KSP ksp,PetscInt it,PetscReal rnorm)
1865: {
1866: PetscInt i, n = ksp->numbermonitors;
1870: for (i=0; i<n; i++) {
1871: (*ksp->monitor[i])(ksp,it,rnorm,ksp->monitorcontext[i]);
1872: }
1873: return(0);
1874: }
1876: /*@C
1877: KSPMonitorSet - Sets an ADDITIONAL function to be called at every iteration to monitor
1878: the residual/error etc.
1880: Logically Collective on ksp
1882: Input Parameters:
1883: + ksp - iterative context obtained from KSPCreate()
1884: . monitor - pointer to function (if this is NULL, it turns off monitoring
1885: . mctx - [optional] context for private data for the
1886: monitor routine (use NULL if no context is desired)
1887: - monitordestroy - [optional] routine that frees monitor context
1888: (may be NULL)
1890: Calling Sequence of monitor:
1891: $ monitor (KSP ksp, PetscInt it, PetscReal rnorm, void *mctx)
1893: + ksp - iterative context obtained from KSPCreate()
1894: . it - iteration number
1895: . rnorm - (estimated) 2-norm of (preconditioned) residual
1896: - mctx - optional monitoring context, as set by KSPMonitorSet()
1898: Options Database Keys:
1899: + -ksp_monitor - sets KSPMonitorDefault()
1900: . -ksp_monitor_true_residual - sets KSPMonitorTrueResidualNorm()
1901: . -ksp_monitor_max - sets KSPMonitorTrueResidualMaxNorm()
1902: . -ksp_monitor_lg_residualnorm - sets line graph monitor,
1903: uses KSPMonitorLGResidualNormCreate()
1904: . -ksp_monitor_lg_true_residualnorm - sets line graph monitor,
1905: uses KSPMonitorLGResidualNormCreate()
1906: . -ksp_monitor_singular_value - sets KSPMonitorSingularValue()
1907: - -ksp_monitor_cancel - cancels all monitors that have
1908: been hardwired into a code by
1909: calls to KSPMonitorSet(), but
1910: does not cancel those set via
1911: the options database.
1913: Notes:
1914: The default is to do nothing. To print the residual, or preconditioned
1915: residual if KSPSetNormType(ksp,KSP_NORM_PRECONDITIONED) was called, use
1916: KSPMonitorDefault() as the monitoring routine, with a ASCII viewer as the
1917: context.
1919: Several different monitoring routines may be set by calling
1920: KSPMonitorSet() multiple times; all will be called in the
1921: order in which they were set.
1923: Fortran Notes:
1924: Only a single monitor function can be set for each KSP object
1926: Level: beginner
1928: .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSPMonitorCancel(), KSP
1929: @*/
1930: PetscErrorCode KSPMonitorSet(KSP ksp,PetscErrorCode (*monitor)(KSP,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
1931: {
1932: PetscInt i;
1934: PetscBool identical;
1938: for (i=0; i<ksp->numbermonitors;i++) {
1939: PetscMonitorCompare((PetscErrorCode (*)(void))monitor,mctx,monitordestroy,(PetscErrorCode (*)(void))ksp->monitor[i],ksp->monitorcontext[i],ksp->monitordestroy[i],&identical);
1940: if (identical) return(0);
1941: }
1942: if (ksp->numbermonitors >= MAXKSPMONITORS) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Too many KSP monitors set");
1943: ksp->monitor[ksp->numbermonitors] = monitor;
1944: ksp->monitordestroy[ksp->numbermonitors] = monitordestroy;
1945: ksp->monitorcontext[ksp->numbermonitors++] = (void*)mctx;
1946: return(0);
1947: }
1949: /*@
1950: KSPMonitorCancel - Clears all monitors for a KSP object.
1952: Logically Collective on ksp
1954: Input Parameters:
1955: . ksp - iterative context obtained from KSPCreate()
1957: Options Database Key:
1958: . -ksp_monitor_cancel - Cancels all monitors that have
1959: been hardwired into a code by calls to KSPMonitorSet(),
1960: but does not cancel those set via the options database.
1962: Level: intermediate
1964: .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSPMonitorSet(), KSP
1965: @*/
1966: PetscErrorCode KSPMonitorCancel(KSP ksp)
1967: {
1969: PetscInt i;
1973: for (i=0; i<ksp->numbermonitors; i++) {
1974: if (ksp->monitordestroy[i]) {
1975: (*ksp->monitordestroy[i])(&ksp->monitorcontext[i]);
1976: }
1977: }
1978: ksp->numbermonitors = 0;
1979: return(0);
1980: }
1982: /*@C
1983: KSPGetMonitorContext - Gets the monitoring context, as set by
1984: KSPMonitorSet() for the FIRST monitor only.
1986: Not Collective
1988: Input Parameter:
1989: . ksp - iterative context obtained from KSPCreate()
1991: Output Parameter:
1992: . ctx - monitoring context
1994: Level: intermediate
1996: .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSP
1997: @*/
1998: PetscErrorCode KSPGetMonitorContext(KSP ksp,void **ctx)
1999: {
2002: *ctx = (ksp->monitorcontext[0]);
2003: return(0);
2004: }
2006: /*@
2007: KSPSetResidualHistory - Sets the array used to hold the residual history.
2008: If set, this array will contain the residual norms computed at each
2009: iteration of the solver.
2011: Not Collective
2013: Input Parameters:
2014: + ksp - iterative context obtained from KSPCreate()
2015: . a - array to hold history
2016: . na - size of a
2017: - reset - PETSC_TRUE indicates the history counter is reset to zero
2018: for each new linear solve
2020: Level: advanced
2022: Notes:
2023: The array is NOT freed by PETSc so the user needs to keep track of
2024: it and destroy once the KSP object is destroyed.
2026: If 'a' is NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
2027: default array of length 10000 is allocated.
2029: .seealso: KSPGetResidualHistory(), KSP
2031: @*/
2032: PetscErrorCode KSPSetResidualHistory(KSP ksp,PetscReal a[],PetscInt na,PetscBool reset)
2033: {
2039: PetscFree(ksp->res_hist_alloc);
2040: if (na != PETSC_DECIDE && na != PETSC_DEFAULT && a) {
2041: ksp->res_hist = a;
2042: ksp->res_hist_max = na;
2043: } else {
2044: if (na != PETSC_DECIDE && na != PETSC_DEFAULT) ksp->res_hist_max = na;
2045: else ksp->res_hist_max = 10000; /* like default ksp->max_it */
2046: PetscCalloc1(ksp->res_hist_max,&ksp->res_hist_alloc);
2048: ksp->res_hist = ksp->res_hist_alloc;
2049: }
2050: ksp->res_hist_len = 0;
2051: ksp->res_hist_reset = reset;
2052: return(0);
2053: }
2055: /*@C
2056: KSPGetResidualHistory - Gets the array used to hold the residual history
2057: and the number of residuals it contains.
2059: Not Collective
2061: Input Parameter:
2062: . ksp - iterative context obtained from KSPCreate()
2064: Output Parameters:
2065: + a - pointer to array to hold history (or NULL)
2066: - na - number of used entries in a (or NULL)
2068: Level: advanced
2070: Notes:
2071: Can only be called after a KSPSetResidualHistory() otherwise a and na are set to zero
2073: The Fortran version of this routine has a calling sequence
2074: $ call KSPGetResidualHistory(KSP ksp, integer na, integer ierr)
2075: note that you have passed a Fortran array into KSPSetResidualHistory() and you need
2076: to access the residual values from this Fortran array you provided. Only the na (number of
2077: residual norms currently held) is set.
2079: .seealso: KSPGetResidualHistory(), KSP
2081: @*/
2082: PetscErrorCode KSPGetResidualHistory(KSP ksp,PetscReal *a[],PetscInt *na)
2083: {
2086: if (a) *a = ksp->res_hist;
2087: if (na) *na = ksp->res_hist_len;
2088: return(0);
2089: }
2091: /*@C
2092: KSPSetConvergenceTest - Sets the function to be used to determine
2093: convergence.
2095: Logically Collective on ksp
2097: Input Parameters:
2098: + ksp - iterative context obtained from KSPCreate()
2099: . converge - pointer to the function
2100: . cctx - context for private data for the convergence routine (may be null)
2101: - destroy - a routine for destroying the context (may be null)
2103: Calling sequence of converge:
2104: $ converge (KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)
2106: + ksp - iterative context obtained from KSPCreate()
2107: . it - iteration number
2108: . rnorm - (estimated) 2-norm of (preconditioned) residual
2109: . reason - the reason why it has converged or diverged
2110: - cctx - optional convergence context, as set by KSPSetConvergenceTest()
2113: Notes:
2114: Must be called after the KSP type has been set so put this after
2115: a call to KSPSetType(), or KSPSetFromOptions().
2117: The default convergence test, KSPConvergedDefault(), aborts if the
2118: residual grows to more than 10000 times the initial residual.
2120: The default is a combination of relative and absolute tolerances.
2121: The residual value that is tested may be an approximation; routines
2122: that need exact values should compute them.
2124: In the default PETSc convergence test, the precise values of reason
2125: are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
2127: Level: advanced
2129: .seealso: KSPConvergedDefault(), KSPGetConvergenceContext(), KSPSetTolerances(), KSP, KSPGetConvergenceTest(), KSPGetAndClearConvergenceTest()
2130: @*/
2131: PetscErrorCode KSPSetConvergenceTest(KSP ksp,PetscErrorCode (*converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
2132: {
2137: if (ksp->convergeddestroy) {
2138: (*ksp->convergeddestroy)(ksp->cnvP);
2139: }
2140: ksp->converged = converge;
2141: ksp->convergeddestroy = destroy;
2142: ksp->cnvP = (void*)cctx;
2143: return(0);
2144: }
2146: /*@C
2147: KSPGetConvergenceTest - Gets the function to be used to determine
2148: convergence.
2150: Logically Collective on ksp
2152: Input Parameter:
2153: . ksp - iterative context obtained from KSPCreate()
2155: Output Parameter:
2156: + converge - pointer to convergence test function
2157: . cctx - context for private data for the convergence routine (may be null)
2158: - destroy - a routine for destroying the context (may be null)
2160: Calling sequence of converge:
2161: $ converge (KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)
2163: + ksp - iterative context obtained from KSPCreate()
2164: . it - iteration number
2165: . rnorm - (estimated) 2-norm of (preconditioned) residual
2166: . reason - the reason why it has converged or diverged
2167: - cctx - optional convergence context, as set by KSPSetConvergenceTest()
2169: Level: advanced
2171: .seealso: KSPConvergedDefault(), KSPGetConvergenceContext(), KSPSetTolerances(), KSP, KSPSetConvergenceTest(), KSPGetAndClearConvergenceTest()
2172: @*/
2173: PetscErrorCode KSPGetConvergenceTest(KSP ksp,PetscErrorCode (**converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void **cctx,PetscErrorCode (**destroy)(void*))
2174: {
2177: if (converge) *converge = ksp->converged;
2178: if (destroy) *destroy = ksp->convergeddestroy;
2179: if (cctx) *cctx = ksp->cnvP;
2180: return(0);
2181: }
2183: /*@C
2184: KSPGetAndClearConvergenceTest - Gets the function to be used to determine convergence. Removes the current test without calling destroy on the test context
2186: Logically Collective on ksp
2188: Input Parameter:
2189: . ksp - iterative context obtained from KSPCreate()
2191: Output Parameter:
2192: + converge - pointer to convergence test function
2193: . cctx - context for private data for the convergence routine
2194: - destroy - a routine for destroying the context
2196: Calling sequence of converge:
2197: $ converge (KSP ksp, PetscInt it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)
2199: + ksp - iterative context obtained from KSPCreate()
2200: . it - iteration number
2201: . rnorm - (estimated) 2-norm of (preconditioned) residual
2202: . reason - the reason why it has converged or diverged
2203: - cctx - optional convergence context, as set by KSPSetConvergenceTest()
2205: Level: advanced
2207: Notes: This is intended to be used to allow transferring the convergence test (and its context) to another testing object (for example another KSP) and then calling
2208: KSPSetConvergenceTest() on this original KSP. If you just called KSPGetConvergenceTest() followed by KSPSetConvergenceTest() the original context information
2209: would be destroyed and hence the transferred context would be invalid and trigger a crash on use
2211: .seealso: KSPConvergedDefault(), KSPGetConvergenceContext(), KSPSetTolerances(), KSP, KSPSetConvergenceTest(), KSPGetConvergenceTest()
2212: @*/
2213: PetscErrorCode KSPGetAndClearConvergenceTest(KSP ksp,PetscErrorCode (**converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void **cctx,PetscErrorCode (**destroy)(void*))
2214: {
2217: *converge = ksp->converged;
2218: *destroy = ksp->convergeddestroy;
2219: *cctx = ksp->cnvP;
2220: ksp->converged = NULL;
2221: ksp->cnvP = NULL;
2222: ksp->convergeddestroy = NULL;
2223: return(0);
2224: }
2226: /*@C
2227: KSPGetConvergenceContext - Gets the convergence context set with
2228: KSPSetConvergenceTest().
2230: Not Collective
2232: Input Parameter:
2233: . ksp - iterative context obtained from KSPCreate()
2235: Output Parameter:
2236: . ctx - monitoring context
2238: Level: advanced
2240: .seealso: KSPConvergedDefault(), KSPSetConvergenceTest(), KSP
2241: @*/
2242: PetscErrorCode KSPGetConvergenceContext(KSP ksp,void **ctx)
2243: {
2246: *ctx = ksp->cnvP;
2247: return(0);
2248: }
2250: /*@C
2251: KSPBuildSolution - Builds the approximate solution in a vector provided.
2252: This routine is NOT commonly needed (see KSPSolve()).
2254: Collective on ksp
2256: Input Parameter:
2257: . ctx - iterative context obtained from KSPCreate()
2259: Output Parameter:
2260: Provide exactly one of
2261: + v - location to stash solution.
2262: - V - the solution is returned in this location. This vector is created
2263: internally. This vector should NOT be destroyed by the user with
2264: VecDestroy().
2266: Notes:
2267: This routine can be used in one of two ways
2268: .vb
2269: KSPBuildSolution(ksp,NULL,&V);
2270: or
2271: KSPBuildSolution(ksp,v,NULL); or KSPBuildSolution(ksp,v,&v);
2272: .ve
2273: In the first case an internal vector is allocated to store the solution
2274: (the user cannot destroy this vector). In the second case the solution
2275: is generated in the vector that the user provides. Note that for certain
2276: methods, such as KSPCG, the second case requires a copy of the solution,
2277: while in the first case the call is essentially free since it simply
2278: returns the vector where the solution already is stored. For some methods
2279: like GMRES this is a reasonably expensive operation and should only be
2280: used in truly needed.
2282: Level: advanced
2284: .seealso: KSPGetSolution(), KSPBuildResidual(), KSP
2285: @*/
2286: PetscErrorCode KSPBuildSolution(KSP ksp,Vec v,Vec *V)
2287: {
2292: if (!V && !v) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONG,"Must provide either v or V");
2293: if (!V) V = &v;
2294: (*ksp->ops->buildsolution)(ksp,v,V);
2295: return(0);
2296: }
2298: /*@C
2299: KSPBuildResidual - Builds the residual in a vector provided.
2301: Collective on ksp
2303: Input Parameter:
2304: . ksp - iterative context obtained from KSPCreate()
2306: Output Parameters:
2307: + v - optional location to stash residual. If v is not provided,
2308: then a location is generated.
2309: . t - work vector. If not provided then one is generated.
2310: - V - the residual
2312: Notes:
2313: Regardless of whether or not v is provided, the residual is
2314: returned in V.
2316: Level: advanced
2318: .seealso: KSPBuildSolution()
2319: @*/
2320: PetscErrorCode KSPBuildResidual(KSP ksp,Vec t,Vec v,Vec *V)
2321: {
2323: PetscBool flag = PETSC_FALSE;
2324: Vec w = v,tt = t;
2328: if (!w) {
2329: VecDuplicate(ksp->vec_rhs,&w);
2330: PetscLogObjectParent((PetscObject)ksp,(PetscObject)w);
2331: }
2332: if (!tt) {
2333: VecDuplicate(ksp->vec_sol,&tt); flag = PETSC_TRUE;
2334: PetscLogObjectParent((PetscObject)ksp,(PetscObject)tt);
2335: }
2336: (*ksp->ops->buildresidual)(ksp,tt,w,V);
2337: if (flag) {VecDestroy(&tt);}
2338: return(0);
2339: }
2341: /*@
2342: KSPSetDiagonalScale - Tells KSP to symmetrically diagonally scale the system
2343: before solving. This actually CHANGES the matrix (and right hand side).
2345: Logically Collective on ksp
2347: Input Parameter:
2348: + ksp - the KSP context
2349: - scale - PETSC_TRUE or PETSC_FALSE
2351: Options Database Key:
2352: + -ksp_diagonal_scale -
2353: - -ksp_diagonal_scale_fix - scale the matrix back AFTER the solve
2356: Notes:
2357: Scales the matrix by D^(-1/2) A D^(-1/2) [D^(1/2) x ] = D^(-1/2) b
2358: where D_{ii} is 1/abs(A_{ii}) unless A_{ii} is zero and then it is 1.
2360: BE CAREFUL with this routine: it actually scales the matrix and right
2361: hand side that define the system. After the system is solved the matrix
2362: and right hand side remain scaled unless you use KSPSetDiagonalScaleFix()
2364: This should NOT be used within the SNES solves if you are using a line
2365: search.
2367: If you use this with the PCType Eisenstat preconditioner than you can
2368: use the PCEisenstatSetNoDiagonalScaling() option, or -pc_eisenstat_no_diagonal_scaling
2369: to save some unneeded, redundant flops.
2371: Level: intermediate
2373: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScaleFix(), KSP
2374: @*/
2375: PetscErrorCode KSPSetDiagonalScale(KSP ksp,PetscBool scale)
2376: {
2380: ksp->dscale = scale;
2381: return(0);
2382: }
2384: /*@
2385: KSPGetDiagonalScale - Checks if KSP solver scales the matrix and
2386: right hand side
2388: Not Collective
2390: Input Parameter:
2391: . ksp - the KSP context
2393: Output Parameter:
2394: . scale - PETSC_TRUE or PETSC_FALSE
2396: Notes:
2397: BE CAREFUL with this routine: it actually scales the matrix and right
2398: hand side that define the system. After the system is solved the matrix
2399: and right hand side remain scaled unless you use KSPSetDiagonalScaleFix()
2401: Level: intermediate
2403: .seealso: KSPSetDiagonalScale(), KSPSetDiagonalScaleFix(), KSP
2404: @*/
2405: PetscErrorCode KSPGetDiagonalScale(KSP ksp,PetscBool *scale)
2406: {
2410: *scale = ksp->dscale;
2411: return(0);
2412: }
2414: /*@
2415: KSPSetDiagonalScaleFix - Tells KSP to diagonally scale the system
2416: back after solving.
2418: Logically Collective on ksp
2420: Input Parameter:
2421: + ksp - the KSP context
2422: - fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
2423: rescale (default)
2425: Notes:
2426: Must be called after KSPSetDiagonalScale()
2428: Using this will slow things down, because it rescales the matrix before and
2429: after each linear solve. This is intended mainly for testing to allow one
2430: to easily get back the original system to make sure the solution computed is
2431: accurate enough.
2433: Level: intermediate
2435: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPGetDiagonalScaleFix(), KSP
2436: @*/
2437: PetscErrorCode KSPSetDiagonalScaleFix(KSP ksp,PetscBool fix)
2438: {
2442: ksp->dscalefix = fix;
2443: return(0);
2444: }
2446: /*@
2447: KSPGetDiagonalScaleFix - Determines if KSP diagonally scales the system
2448: back after solving.
2450: Not Collective
2452: Input Parameter:
2453: . ksp - the KSP context
2455: Output Parameter:
2456: . fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
2457: rescale (default)
2459: Notes:
2460: Must be called after KSPSetDiagonalScale()
2462: If PETSC_TRUE will slow things down, because it rescales the matrix before and
2463: after each linear solve. This is intended mainly for testing to allow one
2464: to easily get back the original system to make sure the solution computed is
2465: accurate enough.
2467: Level: intermediate
2469: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPSetDiagonalScaleFix(), KSP
2470: @*/
2471: PetscErrorCode KSPGetDiagonalScaleFix(KSP ksp,PetscBool *fix)
2472: {
2476: *fix = ksp->dscalefix;
2477: return(0);
2478: }
2480: /*@C
2481: KSPSetComputeOperators - set routine to compute the linear operators
2483: Logically Collective
2485: Input Arguments:
2486: + ksp - the KSP context
2487: . func - function to compute the operators
2488: - ctx - optional context
2490: Calling sequence of func:
2491: $ func(KSP ksp,Mat A,Mat B,void *ctx)
2493: + ksp - the KSP context
2494: . A - the linear operator
2495: . B - preconditioning matrix
2496: - ctx - optional user-provided context
2498: Notes:
2499: The user provided func() will be called automatically at the very next call to KSPSolve(). It will not be called at future KSPSolve() calls
2500: unless either KSPSetComputeOperators() or KSPSetOperators() is called before that KSPSolve() is called.
2502: To reuse the same preconditioner for the next KSPSolve() and not compute a new one based on the most recently computed matrix call KSPSetReusePreconditioner()
2504: Level: beginner
2506: .seealso: KSPSetOperators(), KSPSetComputeRHS(), DMKSPSetComputeOperators(), KSPSetComputeInitialGuess()
2507: @*/
2508: PetscErrorCode KSPSetComputeOperators(KSP ksp,PetscErrorCode (*func)(KSP,Mat,Mat,void*),void *ctx)
2509: {
2511: DM dm;
2515: KSPGetDM(ksp,&dm);
2516: DMKSPSetComputeOperators(dm,func,ctx);
2517: if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX;
2518: return(0);
2519: }
2521: /*@C
2522: KSPSetComputeRHS - set routine to compute the right hand side of the linear system
2524: Logically Collective
2526: Input Arguments:
2527: + ksp - the KSP context
2528: . func - function to compute the right hand side
2529: - ctx - optional context
2531: Calling sequence of func:
2532: $ func(KSP ksp,Vec b,void *ctx)
2534: + ksp - the KSP context
2535: . b - right hand side of linear system
2536: - ctx - optional user-provided context
2538: Notes:
2539: The routine you provide will be called EACH you call KSPSolve() to prepare the new right hand side for that solve
2541: Level: beginner
2543: .seealso: KSPSolve(), DMKSPSetComputeRHS(), KSPSetComputeOperators()
2544: @*/
2545: PetscErrorCode KSPSetComputeRHS(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2546: {
2548: DM dm;
2552: KSPGetDM(ksp,&dm);
2553: DMKSPSetComputeRHS(dm,func,ctx);
2554: return(0);
2555: }
2557: /*@C
2558: KSPSetComputeInitialGuess - set routine to compute the initial guess of the linear system
2560: Logically Collective
2562: Input Arguments:
2563: + ksp - the KSP context
2564: . func - function to compute the initial guess
2565: - ctx - optional context
2567: Calling sequence of func:
2568: $ func(KSP ksp,Vec x,void *ctx)
2570: + ksp - the KSP context
2571: . x - solution vector
2572: - ctx - optional user-provided context
2574: Notes: This should only be used in conjunction with KSPSetComputeRHS(), KSPSetComputeOperators(), otherwise
2575: call KSPSetInitialGuessNonzero() and set the initial guess values in the solution vector passed to KSPSolve().
2577: Level: beginner
2579: .seealso: KSPSolve(), KSPSetComputeRHS(), KSPSetComputeOperators(), DMKSPSetComputeInitialGuess()
2580: @*/
2581: PetscErrorCode KSPSetComputeInitialGuess(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2582: {
2584: DM dm;
2588: KSPGetDM(ksp,&dm);
2589: DMKSPSetComputeInitialGuess(dm,func,ctx);
2590: return(0);
2591: }