Actual source code: itfunc.c
petsc-3.4.5 2014-06-29
2: /*
3: Interface KSP routines that the user calls.
4: */
6: #include <petsc-private/kspimpl.h> /*I "petscksp.h" I*/
7: #include <petscdm.h>
11: /*@
12: KSPComputeExtremeSingularValues - Computes the extreme singular values
13: for the preconditioned operator. Called after or during KSPSolve().
15: Not Collective
17: Input Parameter:
18: . ksp - iterative context obtained from KSPCreate()
20: Output Parameters:
21: . emin, emax - extreme singular values
23: Options Database Keys:
24: . -ksp_compute_singularvalues - compute extreme singular values and print when KSPSolve completes.
26: Notes:
27: One must call KSPSetComputeSingularValues() before calling KSPSetUp()
28: (or use the option -ksp_compute_eigenvalues) in order for this routine to work correctly.
30: Many users may just want to use the monitoring routine
31: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
32: to print the extreme singular values at each iteration of the linear solve.
34: Estimates of the smallest singular value may be very inaccurate, especially if the Krylov method has not converged.
35: The largest singular value is usually accurate to within a few percent if the method has converged, but is still not
36: intended for eigenanalysis.
38: Disable restarts if using KSPGMRES, otherwise this estimate will only be using those iterations after the last
39: restart. See KSPGMRESSetRestart() for more details.
41: Level: advanced
43: .keywords: KSP, compute, extreme, singular, values
45: .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeEigenvalues()
46: @*/
47: PetscErrorCode KSPComputeExtremeSingularValues(KSP ksp,PetscReal *emax,PetscReal *emin)
48: {
55: if (!ksp->calc_sings) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Singular values not requested before KSPSetUp()");
57: if (ksp->ops->computeextremesingularvalues) {
58: (*ksp->ops->computeextremesingularvalues)(ksp,emax,emin);
59: } else {
60: *emin = -1.0;
61: *emax = -1.0;
62: }
63: return(0);
64: }
68: /*@
69: KSPComputeEigenvalues - Computes the extreme eigenvalues for the
70: preconditioned operator. Called after or during KSPSolve().
72: Not Collective
74: Input Parameter:
75: + ksp - iterative context obtained from KSPCreate()
76: - n - size of arrays r and c. The number of eigenvalues computed (neig) will, in
77: general, be less than this.
79: Output Parameters:
80: + r - real part of computed eigenvalues
81: . c - complex part of computed eigenvalues
82: - neig - number of eigenvalues computed (will be less than or equal to n)
84: Options Database Keys:
85: + -ksp_compute_eigenvalues - Prints eigenvalues to stdout
86: - -ksp_plot_eigenvalues - Plots eigenvalues in an x-window display
88: Notes:
89: The number of eigenvalues estimated depends on the size of the Krylov space
90: generated during the KSPSolve() ; for example, with
91: CG it corresponds to the number of CG iterations, for GMRES it is the number
92: of GMRES iterations SINCE the last restart. Any extra space in r[] and c[]
93: will be ignored.
95: KSPComputeEigenvalues() does not usually provide accurate estimates; it is
96: intended only for assistance in understanding the convergence of iterative
97: methods, not for eigenanalysis.
99: One must call KSPSetComputeEigenvalues() before calling KSPSetUp()
100: in order for this routine to work correctly.
102: Many users may just want to use the monitoring routine
103: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
104: to print the singular values at each iteration of the linear solve.
106: Level: advanced
108: .keywords: KSP, compute, extreme, singular, values
110: .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues()
111: @*/
112: PetscErrorCode KSPComputeEigenvalues(KSP ksp,PetscInt n,PetscReal *r,PetscReal *c,PetscInt *neig)
113: {
121: if (!ksp->calc_sings) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Eigenvalues not requested before KSPSetUp()");
123: if (ksp->ops->computeeigenvalues) {
124: (*ksp->ops->computeeigenvalues)(ksp,n,r,c,neig);
125: } else {
126: *neig = 0;
127: }
128: return(0);
129: }
133: /*@
134: KSPSetUpOnBlocks - Sets up the preconditioner for each block in
135: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
136: methods.
138: Collective on KSP
140: Input Parameter:
141: . ksp - the KSP context
143: Notes:
144: KSPSetUpOnBlocks() is a routine that the user can optinally call for
145: more precise profiling (via -log_summary) of the setup phase for these
146: block preconditioners. If the user does not call KSPSetUpOnBlocks(),
147: it will automatically be called from within KSPSolve().
149: Calling KSPSetUpOnBlocks() is the same as calling PCSetUpOnBlocks()
150: on the PC context within the KSP context.
152: Level: advanced
154: .keywords: KSP, setup, blocks
156: .seealso: PCSetUpOnBlocks(), KSPSetUp(), PCSetUp()
157: @*/
158: PetscErrorCode KSPSetUpOnBlocks(KSP ksp)
159: {
164: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
165: PCSetUpOnBlocks(ksp->pc);
166: return(0);
167: }
171: /*@
172: KSPSetUp - Sets up the internal data structures for the
173: later use of an iterative solver.
175: Collective on KSP
177: Input Parameter:
178: . ksp - iterative context obtained from KSPCreate()
180: Level: developer
182: .keywords: KSP, setup
184: .seealso: KSPCreate(), KSPSolve(), KSPDestroy()
185: @*/
186: PetscErrorCode KSPSetUp(KSP ksp)
187: {
189: Mat A,B;
190: MatStructure stflg = SAME_NONZERO_PATTERN;
195: /* reset the convergence flag from the previous solves */
196: ksp->reason = KSP_CONVERGED_ITERATING;
198: if (!((PetscObject)ksp)->type_name) {
199: KSPSetType(ksp,KSPGMRES);
200: }
201: KSPSetUpNorms_Private(ksp,&ksp->normtype,&ksp->pc_side);
203: if (ksp->dmActive && !ksp->setupstage) {
204: /* first time in so build matrix and vector data structures using DM */
205: if (!ksp->vec_rhs) {DMCreateGlobalVector(ksp->dm,&ksp->vec_rhs);}
206: if (!ksp->vec_sol) {DMCreateGlobalVector(ksp->dm,&ksp->vec_sol);}
207: DMCreateMatrix(ksp->dm,MATAIJ,&A);
208: KSPSetOperators(ksp,A,A,stflg);
209: PetscObjectDereference((PetscObject)A);
210: }
212: if (ksp->dmActive) {
213: DMKSP kdm;
214: DMGetDMKSP(ksp->dm,&kdm);
216: if (kdm->ops->computeinitialguess && ksp->setupstage != KSP_SETUP_NEWRHS) {
217: /* only computes initial guess the first time through */
218: (*kdm->ops->computeinitialguess)(ksp,ksp->vec_sol,kdm->initialguessctx);
219: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
220: }
221: if (kdm->ops->computerhs) {
222: (*kdm->ops->computerhs)(ksp,ksp->vec_rhs,kdm->rhsctx);
223: }
225: if (ksp->setupstage != KSP_SETUP_NEWRHS) {
226: if (kdm->ops->computeoperators) {
227: KSPGetOperators(ksp,&A,&B,NULL);
228: (*kdm->ops->computeoperators)(ksp,A,B,&stflg,kdm->operatorsctx);
229: KSPSetOperators(ksp,A,B,stflg);
230: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You called KSPSetDM() but did not use DMKSPSetComputeOperators() or KSPSetDMActive(dm,PETSC_FALSE);");
231: }
232: }
234: if (ksp->setupstage == KSP_SETUP_NEWRHS) return(0);
235: PetscLogEventBegin(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);
237: switch (ksp->setupstage) {
238: case KSP_SETUP_NEW:
239: (*ksp->ops->setup)(ksp);
240: break;
241: case KSP_SETUP_NEWMATRIX: { /* This should be replaced with a more general mechanism */
242: KSPChebyshevSetNewMatrix(ksp);
243: } break;
244: default: break;
245: }
247: /* scale the matrix if requested */
248: if (ksp->dscale) {
249: Mat mat,pmat;
250: PetscScalar *xx;
251: PetscInt i,n;
252: PetscBool zeroflag = PETSC_FALSE;
253: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
254: PCGetOperators(ksp->pc,&mat,&pmat,NULL);
255: if (!ksp->diagonal) { /* allocate vector to hold diagonal */
256: MatGetVecs(pmat,&ksp->diagonal,0);
257: }
258: MatGetDiagonal(pmat,ksp->diagonal);
259: VecGetLocalSize(ksp->diagonal,&n);
260: VecGetArray(ksp->diagonal,&xx);
261: for (i=0; i<n; i++) {
262: if (xx[i] != 0.0) xx[i] = 1.0/PetscSqrtReal(PetscAbsScalar(xx[i]));
263: else {
264: xx[i] = 1.0;
265: zeroflag = PETSC_TRUE;
266: }
267: }
268: VecRestoreArray(ksp->diagonal,&xx);
269: if (zeroflag) {
270: PetscInfo(ksp,"Zero detected in diagonal of matrix, using 1 at those locations\n");
271: }
272: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
273: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
274: ksp->dscalefix2 = PETSC_FALSE;
275: }
276: PetscLogEventEnd(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);
277: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
278: PCSetUp(ksp->pc);
279: if (ksp->nullsp) {
280: PetscBool test = PETSC_FALSE;
281: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_test_null_space",&test,NULL);
282: if (test) {
283: Mat mat;
284: PCGetOperators(ksp->pc,&mat,NULL,NULL);
285: MatNullSpaceTest(ksp->nullsp,mat,NULL);
286: }
287: }
288: ksp->setupstage = KSP_SETUP_NEWRHS;
289: return(0);
290: }
292: #include <petscdraw.h>
295: /*@
296: KSPSolve - Solves linear system.
298: Collective on KSP
300: Parameter:
301: + ksp - iterative context obtained from KSPCreate()
302: . b - the right hand side vector
303: - x - the solution (this may be the same vector as b, then b will be overwritten with answer)
305: Options Database Keys:
306: + -ksp_compute_eigenvalues - compute preconditioned operators eigenvalues
307: . -ksp_plot_eigenvalues - plot the computed eigenvalues in an X-window
308: . -ksp_compute_eigenvalues_explicitly - compute the eigenvalues by forming the dense operator and useing LAPACK
309: . -ksp_plot_eigenvalues_explicitly - plot the explicitly computing eigenvalues
310: . -ksp_view_mat binary - save matrix to the default binary viewer
311: . -ksp_view_pmat binary - save matrix used to build preconditioner to the default binary viewer
312: . -ksp_view_rhs binary - save right hand side vector to the default binary viewer
313: . -ksp_view_solution binary - save computed solution vector to the default binary viewer
314: (can be read later with src/ksp/examples/tutorials/ex10.c for testing solvers)
315: . -ksp_view_mat_explicit - for matrix-free operators, computes the matrix entries and views them
316: . -ksp_view_preconditioned_operator_explicit - computes the product of the preconditioner and matrix as an explicit matrix and views it
317: . -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations
318: . -ksp_final_residual - print 2-norm of true linear system residual at the end of the solution process
319: - -ksp_view - print the ksp data structure at the end of the system solution
321: Notes:
323: If one uses KSPSetDM() then x or b need not be passed. Use KSPGetSolution() to access the solution in this case.
325: The operator is specified with KSPSetOperators().
327: Call KSPGetConvergedReason() to determine if the solver converged or failed and
328: why. The number of iterations can be obtained from KSPGetIterationNumber().
330: If using a direct method (e.g., via the KSP solver
331: KSPPREONLY and a preconditioner such as PCLU/PCILU),
332: then its=1. See KSPSetTolerances() and KSPDefaultConverged()
333: for more details.
335: Understanding Convergence:
336: The routines KSPMonitorSet(), KSPComputeEigenvalues(), and
337: KSPComputeEigenvaluesExplicitly() provide information on additional
338: options to monitor convergence and print eigenvalue information.
340: Level: beginner
342: .keywords: KSP, solve, linear system
344: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPDefaultConverged(),
345: KSPSolveTranspose(), KSPGetIterationNumber()
346: @*/
347: PetscErrorCode KSPSolve(KSP ksp,Vec b,Vec x)
348: {
349: PetscErrorCode ierr;
350: PetscMPIInt rank;
351: PetscBool flag1,flag2,flag3,flg = PETSC_FALSE,inXisinB=PETSC_FALSE,guess_zero;
352: PetscViewer viewer;
353: Mat mat,premat;
354: PetscViewerFormat format;
361: if (x && x == b) {
362: if (!ksp->guess_zero) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Cannot use x == b with nonzero initial guess");
363: VecDuplicate(b,&x);
364: inXisinB = PETSC_TRUE;
365: }
366: if (b) {
367: PetscObjectReference((PetscObject)b);
368: VecDestroy(&ksp->vec_rhs);
369: ksp->vec_rhs = b;
370: }
371: if (x) {
372: PetscObjectReference((PetscObject)x);
373: VecDestroy(&ksp->vec_sol);
374: ksp->vec_sol = x;
375: }
377: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_view_pre",&viewer,&format,&flg);
378: if (flg && !PetscPreLoadingOn) {
379: PetscViewerPushFormat(viewer,format);
380: KSPView(ksp,viewer);
381: PetscViewerPopFormat(viewer);
382: PetscViewerDestroy(&viewer);
383: }
385: if (ksp->presolve) {
386: (*ksp->presolve)(ksp,ksp->vec_rhs,ksp->vec_sol,ksp->prectx);
387: }
388: PetscLogEventBegin(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
390: /* reset the residual history list if requested */
391: if (ksp->res_hist_reset) ksp->res_hist_len = 0;
392: ksp->transpose_solve = PETSC_FALSE;
394: if (ksp->guess) {
395: KSPFischerGuessFormGuess(ksp->guess,ksp->vec_rhs,ksp->vec_sol);
396: ksp->guess_zero = PETSC_FALSE;
397: }
398: /* KSPSetUp() scales the matrix if needed */
399: KSPSetUp(ksp);
400: KSPSetUpOnBlocks(ksp);
402: /* diagonal scale RHS if called for */
403: if (ksp->dscale) {
404: VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
405: /* second time in, but matrix was scaled back to original */
406: if (ksp->dscalefix && ksp->dscalefix2) {
407: Mat mat,pmat;
409: PCGetOperators(ksp->pc,&mat,&pmat,NULL);
410: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
411: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
412: }
414: /* scale initial guess */
415: if (!ksp->guess_zero) {
416: if (!ksp->truediagonal) {
417: VecDuplicate(ksp->diagonal,&ksp->truediagonal);
418: VecCopy(ksp->diagonal,ksp->truediagonal);
419: VecReciprocal(ksp->truediagonal);
420: }
421: VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->truediagonal);
422: }
423: }
424: PCPreSolve(ksp->pc,ksp);
426: if (ksp->guess_zero) { VecSet(ksp->vec_sol,0.0);}
427: if (ksp->guess_knoll) {
428: PCApply(ksp->pc,ksp->vec_rhs,ksp->vec_sol);
429: KSP_RemoveNullSpace(ksp,ksp->vec_sol);
430: ksp->guess_zero = PETSC_FALSE;
431: }
433: /* can we mark the initial guess as zero for this solve? */
434: guess_zero = ksp->guess_zero;
435: if (!ksp->guess_zero) {
436: PetscReal norm;
438: VecNormAvailable(ksp->vec_sol,NORM_2,&flg,&norm);
439: if (flg && !norm) ksp->guess_zero = PETSC_TRUE;
440: }
441: (*ksp->ops->solve)(ksp);
442: ksp->guess_zero = guess_zero;
444: if (!ksp->reason) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
445: if (ksp->printreason) {
446: PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)),((PetscObject)ksp)->tablevel);
447: if (ksp->reason > 0) {
448: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)),"Linear solve converged due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
449: } else {
450: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)),"Linear solve did not converge due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
451: }
452: PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp)),((PetscObject)ksp)->tablevel);
453: }
454: PCPostSolve(ksp->pc,ksp);
456: /* diagonal scale solution if called for */
457: if (ksp->dscale) {
458: VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->diagonal);
459: /* unscale right hand side and matrix */
460: if (ksp->dscalefix) {
461: Mat mat,pmat;
463: VecReciprocal(ksp->diagonal);
464: VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
465: PCGetOperators(ksp->pc,&mat,&pmat,NULL);
466: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
467: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
468: VecReciprocal(ksp->diagonal);
469: ksp->dscalefix2 = PETSC_TRUE;
470: }
471: }
472: PetscLogEventEnd(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
473: if (ksp->postsolve) {
474: (*ksp->postsolve)(ksp,ksp->vec_rhs,ksp->vec_sol,ksp->postctx);
475: }
477: if (ksp->guess) {
478: KSPFischerGuessUpdate(ksp->guess,ksp->vec_sol);
479: }
481: MPI_Comm_rank(PetscObjectComm((PetscObject)ksp),&rank);
483: PCGetOperators(ksp->pc,&mat,&premat,NULL);
484: MatViewFromOptions(mat,((PetscObject)ksp)->prefix,"-ksp_view_mat");
485: MatViewFromOptions(premat,((PetscObject)ksp)->prefix,"-ksp_view_pmat");
486: VecViewFromOptions(ksp->vec_rhs,((PetscObject)ksp)->prefix,"-ksp_view_rhs");
488: flag1 = PETSC_FALSE;
489: flag2 = PETSC_FALSE;
490: flag3 = PETSC_FALSE;
491: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_compute_eigenvalues",&flag1,NULL);
492: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_plot_eigenvalues",&flag2,NULL);
493: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_plot_eigencontours",&flag3,NULL);
494: if (flag1 || flag2 || flag3) {
495: PetscInt nits,n,i,neig;
496: PetscReal *r,*c;
498: KSPGetIterationNumber(ksp,&nits);
499: n = nits+2;
501: if (!nits) {
502: PetscPrintf(PetscObjectComm((PetscObject)ksp),"Zero iterations in solver, cannot approximate any eigenvalues\n");
503: } else {
504: PetscMalloc(2*n*sizeof(PetscReal),&r);
505: c = r + n;
506: KSPComputeEigenvalues(ksp,n,r,c,&neig);
507: if (flag1) {
508: PetscPrintf(PetscObjectComm((PetscObject)ksp),"Iteratively computed eigenvalues\n");
509: for (i=0; i<neig; i++) {
510: if (c[i] >= 0.0) {
511: PetscPrintf(PetscObjectComm((PetscObject)ksp),"%G + %Gi\n",r[i],c[i]);
512: } else {
513: PetscPrintf(PetscObjectComm((PetscObject)ksp),"%G - %Gi\n",r[i],-c[i]);
514: }
515: }
516: }
517: if (flag2 && !rank) {
518: PetscDraw draw;
519: PetscDrawSP drawsp;
521: if (!ksp->eigviewer) {
522: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Iteratively Computed Eigenvalues",PETSC_DECIDE,PETSC_DECIDE,400,400,&ksp->eigviewer);
523: }
524: PetscViewerDrawGetDraw(ksp->eigviewer,0,&draw);
525: PetscDrawSPCreate(draw,1,&drawsp);
526: PetscDrawSPReset(drawsp);
527: for (i=0; i<neig; i++) {
528: PetscDrawSPAddPoint(drawsp,r+i,c+i);
529: }
530: PetscDrawSPDraw(drawsp,PETSC_TRUE);
531: PetscDrawSPDestroy(&drawsp);
532: }
533: if (flag3 && !rank) {
534: KSPPlotEigenContours_Private(ksp,neig,r,c);
535: }
536: PetscFree(r);
537: }
538: }
540: flag1 = PETSC_FALSE;
541: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_compute_singularvalues",&flag1,NULL);
542: if (flag1) {
543: PetscInt nits;
545: KSPGetIterationNumber(ksp,&nits);
546: if (!nits) {
547: PetscPrintf(PetscObjectComm((PetscObject)ksp),"Zero iterations in solver, cannot approximate any singular values\n");
548: } else {
549: PetscReal emax,emin;
551: KSPComputeExtremeSingularValues(ksp,&emax,&emin);
552: PetscPrintf(PetscObjectComm((PetscObject)ksp),"Iteratively computed extreme singular values: max %G min %G max/min %G\n",emax,emin,emax/emin);
553: }
554: }
557: flag1 = PETSC_FALSE;
558: flag2 = PETSC_FALSE;
559: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_compute_eigenvalues_explicitly",&flag1,NULL);
560: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_plot_eigenvalues_explicitly",&flag2,NULL);
561: if (flag1 || flag2) {
562: PetscInt n,i;
563: PetscReal *r,*c;
564: VecGetSize(ksp->vec_sol,&n);
565: PetscMalloc2(n,PetscReal,&r,n,PetscReal,&c);
566: KSPComputeEigenvaluesExplicitly(ksp,n,r,c);
567: if (flag1) {
568: PetscPrintf(PetscObjectComm((PetscObject)ksp),"Explicitly computed eigenvalues\n");
569: for (i=0; i<n; i++) {
570: if (c[i] >= 0.0) {
571: PetscPrintf(PetscObjectComm((PetscObject)ksp),"%G + %Gi\n",r[i],c[i]);
572: } else {
573: PetscPrintf(PetscObjectComm((PetscObject)ksp),"%G - %Gi\n",r[i],-c[i]);
574: }
575: }
576: }
577: if (flag2 && !rank) {
578: PetscDraw draw;
579: PetscDrawSP drawsp;
581: if (!ksp->eigviewer) {
582: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Explicitly Computed Eigenvalues",0,320,400,400,&ksp->eigviewer);
583: }
584: PetscViewerDrawGetDraw(ksp->eigviewer,0,&draw);
585: PetscDrawSPCreate(draw,1,&drawsp);
586: PetscDrawSPReset(drawsp);
587: for (i=0; i<n; i++) {
588: PetscDrawSPAddPoint(drawsp,r+i,c+i);
589: }
590: PetscDrawSPDraw(drawsp,PETSC_TRUE);
591: PetscDrawSPDestroy(&drawsp);
592: }
593: PetscFree2(r,c);
594: }
596: PetscOptionsHasName(((PetscObject)ksp)->prefix,"-ksp_view_mat_explicit",&flag2);
597: if (flag2) {
598: Mat A,B;
599: PCGetOperators(ksp->pc,&A,NULL,NULL);
600: MatComputeExplicitOperator(A,&B);
601: MatViewFromOptions(B,((PetscObject)ksp)->prefix,"-ksp_view_mat_explicit");
602: MatDestroy(&B);
603: }
604: PetscOptionsHasName(((PetscObject)ksp)->prefix,"-ksp_view_preconditioned_operator_explicit",&flag2);
605: if (flag2) {
606: Mat B;
607: KSPComputeExplicitOperator(ksp,&B);
608: MatViewFromOptions(B,((PetscObject)ksp)->prefix,"-ksp_view_preconditioned_operator_explicit");
609: MatDestroy(&B);
610: }
611: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_view",&viewer,&format,&flg);
612: if (flg && !PetscPreLoadingOn) {
613: PetscViewerPushFormat(viewer,format);
614: KSPView(ksp,viewer);
615: PetscViewerPopFormat(viewer);
616: PetscViewerDestroy(&viewer);
617: }
619: flg = PETSC_FALSE;
620: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_final_residual",&flg,NULL);
621: if (flg) {
622: Mat A;
623: Vec t;
624: PetscReal norm;
625: 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");
626: PCGetOperators(ksp->pc,&A,0,0);
627: VecDuplicate(ksp->vec_rhs,&t);
628: KSP_MatMult(ksp,A,ksp->vec_sol,t);
629: VecAYPX(t, -1.0, ksp->vec_rhs);
630: VecNorm(t,NORM_2,&norm);
631: VecDestroy(&t);
632: PetscPrintf(PetscObjectComm((PetscObject)ksp),"KSP final norm of residual %G\n",norm);
633: }
634: VecViewFromOptions(ksp->vec_sol,((PetscObject)ksp)->prefix,"-ksp_view_solution");
636: if (inXisinB) {
637: VecCopy(x,b);
638: VecDestroy(&x);
639: }
640: PetscObjectAMSBlock((PetscObject)ksp);
641: if (ksp->errorifnotconverged && ksp->reason < 0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged");
642: return(0);
643: }
647: /*@
648: KSPSolveTranspose - Solves the transpose of a linear system.
650: Collective on KSP
652: Input Parameter:
653: + ksp - iterative context obtained from KSPCreate()
654: . b - right hand side vector
655: - x - solution vector
657: Notes: For complex numbers this solve the non-Hermitian transpose system.
659: Developer Notes: We need to implement a KSPSolveHermitianTranspose()
661: Level: developer
663: .keywords: KSP, solve, linear system
665: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPDefaultConverged(),
666: KSPSolve()
667: @*/
669: PetscErrorCode KSPSolveTranspose(KSP ksp,Vec b,Vec x)
670: {
672: PetscBool inXisinB=PETSC_FALSE;
678: if (x == b) {
679: VecDuplicate(b,&x);
680: inXisinB = PETSC_TRUE;
681: }
682: PetscObjectReference((PetscObject)b);
683: PetscObjectReference((PetscObject)x);
684: VecDestroy(&ksp->vec_rhs);
685: VecDestroy(&ksp->vec_sol);
687: ksp->vec_rhs = b;
688: ksp->vec_sol = x;
689: ksp->transpose_solve = PETSC_TRUE;
691: KSPSetUp(ksp);
692: if (ksp->guess_zero) { VecSet(ksp->vec_sol,0.0);}
693: (*ksp->ops->solve)(ksp);
694: if (!ksp->reason) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
695: if (ksp->printreason) {
696: if (ksp->reason > 0) {
697: PetscPrintf(PetscObjectComm((PetscObject)ksp),"Linear solve converged due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
698: } else {
699: PetscPrintf(PetscObjectComm((PetscObject)ksp),"Linear solve did not converge due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
700: }
701: }
702: if (inXisinB) {
703: VecCopy(x,b);
704: VecDestroy(&x);
705: }
706: if (ksp->errorifnotconverged && ksp->reason < 0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged");
707: return(0);
708: }
712: /*@
713: KSPReset - Resets a KSP context to the kspsetupcalled = 0 state and removes any allocated Vecs and Mats
715: Collective on KSP
717: Input Parameter:
718: . ksp - iterative context obtained from KSPCreate()
720: Level: beginner
722: .keywords: KSP, destroy
724: .seealso: KSPCreate(), KSPSetUp(), KSPSolve()
725: @*/
726: PetscErrorCode KSPReset(KSP ksp)
727: {
732: if (!ksp) return(0);
733: if (ksp->ops->reset) {
734: (*ksp->ops->reset)(ksp);
735: }
736: if (ksp->pc) {PCReset(ksp->pc);}
737: KSPFischerGuessDestroy(&ksp->guess);
738: VecDestroyVecs(ksp->nwork,&ksp->work);
739: VecDestroy(&ksp->vec_rhs);
740: VecDestroy(&ksp->vec_sol);
741: VecDestroy(&ksp->diagonal);
742: VecDestroy(&ksp->truediagonal);
743: MatNullSpaceDestroy(&ksp->nullsp);
745: ksp->setupstage = KSP_SETUP_NEW;
746: return(0);
747: }
751: /*@
752: KSPDestroy - Destroys KSP context.
754: Collective on KSP
756: Input Parameter:
757: . ksp - iterative context obtained from KSPCreate()
759: Level: beginner
761: .keywords: KSP, destroy
763: .seealso: KSPCreate(), KSPSetUp(), KSPSolve()
764: @*/
765: PetscErrorCode KSPDestroy(KSP *ksp)
766: {
768: PC pc;
771: if (!*ksp) return(0);
773: if (--((PetscObject)(*ksp))->refct > 0) {*ksp = 0; return(0);}
775: PetscObjectAMSViewOff((PetscObject)*ksp);
776: /*
777: Avoid a cascading call to PCReset(ksp->pc) from the following call:
778: PCReset() shouldn't be called from KSPDestroy() as it is unprotected by pc's
779: refcount (and may be shared, e.g., by other ksps).
780: */
781: pc = (*ksp)->pc;
782: (*ksp)->pc = NULL;
783: KSPReset((*ksp));
784: (*ksp)->pc = pc;
785: if ((*ksp)->ops->destroy) {(*(*ksp)->ops->destroy)(*ksp);}
787: DMDestroy(&(*ksp)->dm);
788: PCDestroy(&(*ksp)->pc);
789: PetscFree((*ksp)->res_hist_alloc);
790: if ((*ksp)->convergeddestroy) {
791: (*(*ksp)->convergeddestroy)((*ksp)->cnvP);
792: }
793: KSPMonitorCancel((*ksp));
794: PetscViewerDestroy(&(*ksp)->eigviewer);
795: PetscHeaderDestroy(ksp);
796: return(0);
797: }
801: /*@
802: KSPSetPCSide - Sets the preconditioning side.
804: Logically Collective on KSP
806: Input Parameter:
807: . ksp - iterative context obtained from KSPCreate()
809: Output Parameter:
810: . side - the preconditioning side, where side is one of
811: .vb
812: PC_LEFT - left preconditioning (default)
813: PC_RIGHT - right preconditioning
814: PC_SYMMETRIC - symmetric preconditioning
815: .ve
817: Options Database Keys:
818: . -ksp_pc_side <right,left,symmetric>
820: Notes:
821: Left preconditioning is used by default for most Krylov methods except KSPFGMRES which only supports right preconditioning.
822: Symmetric preconditioning is currently available only for the KSPQCG method. Note, however, that
823: symmetric preconditioning can be emulated by using either right or left
824: preconditioning and a pre or post processing step.
826: Level: intermediate
828: .keywords: KSP, set, right, left, symmetric, side, preconditioner, flag
830: .seealso: KSPGetPCSide()
831: @*/
832: PetscErrorCode KSPSetPCSide(KSP ksp,PCSide side)
833: {
837: ksp->pc_side = side;
838: return(0);
839: }
843: /*@
844: KSPGetPCSide - Gets the preconditioning side.
846: Not Collective
848: Input Parameter:
849: . ksp - iterative context obtained from KSPCreate()
851: Output Parameter:
852: . side - the preconditioning side, where side is one of
853: .vb
854: PC_LEFT - left preconditioning (default)
855: PC_RIGHT - right preconditioning
856: PC_SYMMETRIC - symmetric preconditioning
857: .ve
859: Level: intermediate
861: .keywords: KSP, get, right, left, symmetric, side, preconditioner, flag
863: .seealso: KSPSetPCSide()
864: @*/
865: PetscErrorCode KSPGetPCSide(KSP ksp,PCSide *side)
866: {
872: KSPSetUpNorms_Private(ksp,&ksp->normtype,&ksp->pc_side);
873: *side = ksp->pc_side;
874: return(0);
875: }
879: /*@
880: KSPGetTolerances - Gets the relative, absolute, divergence, and maximum
881: iteration tolerances used by the default KSP convergence tests.
883: Not Collective
885: Input Parameter:
886: . ksp - the Krylov subspace context
888: Output Parameters:
889: + rtol - the relative convergence tolerance
890: . abstol - the absolute convergence tolerance
891: . dtol - the divergence tolerance
892: - maxits - maximum number of iterations
894: Notes:
895: The user can specify NULL for any parameter that is not needed.
897: Level: intermediate
899: .keywords: KSP, get, tolerance, absolute, relative, divergence, convergence,
900: maximum, iterations
902: .seealso: KSPSetTolerances()
903: @*/
904: PetscErrorCode KSPGetTolerances(KSP ksp,PetscReal *rtol,PetscReal *abstol,PetscReal *dtol,PetscInt *maxits)
905: {
908: if (abstol) *abstol = ksp->abstol;
909: if (rtol) *rtol = ksp->rtol;
910: if (dtol) *dtol = ksp->divtol;
911: if (maxits) *maxits = ksp->max_it;
912: return(0);
913: }
917: /*@
918: KSPSetTolerances - Sets the relative, absolute, divergence, and maximum
919: iteration tolerances used by the default KSP convergence testers.
921: Logically Collective on KSP
923: Input Parameters:
924: + ksp - the Krylov subspace context
925: . rtol - the relative convergence tolerance
926: (relative decrease in the residual norm)
927: . abstol - the absolute convergence tolerance
928: (absolute size of the residual norm)
929: . dtol - the divergence tolerance
930: (amount residual can increase before KSPDefaultConverged()
931: concludes that the method is diverging)
932: - maxits - maximum number of iterations to use
934: Options Database Keys:
935: + -ksp_atol <abstol> - Sets abstol
936: . -ksp_rtol <rtol> - Sets rtol
937: . -ksp_divtol <dtol> - Sets dtol
938: - -ksp_max_it <maxits> - Sets maxits
940: Notes:
941: Use PETSC_DEFAULT to retain the default value of any of the tolerances.
943: See KSPDefaultConverged() for details on the use of these parameters
944: in the default convergence test. See also KSPSetConvergenceTest()
945: for setting user-defined stopping criteria.
947: Level: intermediate
949: .keywords: KSP, set, tolerance, absolute, relative, divergence,
950: convergence, maximum, iterations
952: .seealso: KSPGetTolerances(), KSPDefaultConverged(), KSPSetConvergenceTest()
953: @*/
954: PetscErrorCode KSPSetTolerances(KSP ksp,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt maxits)
955: {
963: if (rtol != PETSC_DEFAULT) {
964: 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",rtol);
965: ksp->rtol = rtol;
966: }
967: if (abstol != PETSC_DEFAULT) {
968: if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol);
969: ksp->abstol = abstol;
970: }
971: if (dtol != PETSC_DEFAULT) {
972: if (dtol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Divergence tolerance %G must be larger than 1.0",dtol);
973: ksp->divtol = dtol;
974: }
975: if (maxits != PETSC_DEFAULT) {
976: if (maxits < 0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxits);
977: ksp->max_it = maxits;
978: }
979: return(0);
980: }
984: /*@
985: KSPSetInitialGuessNonzero - Tells the iterative solver that the
986: initial guess is nonzero; otherwise KSP assumes the initial guess
987: is to be zero (and thus zeros it out before solving).
989: Logically Collective on KSP
991: Input Parameters:
992: + ksp - iterative context obtained from KSPCreate()
993: - flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero
995: Options database keys:
996: . -ksp_initial_guess_nonzero : use nonzero initial guess; this takes an optional truth value (0/1/no/yes/true/false)
998: Level: beginner
1000: Notes:
1001: If this is not called the X vector is zeroed in the call to KSPSolve().
1003: .keywords: KSP, set, initial guess, nonzero
1005: .seealso: KSPGetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll()
1006: @*/
1007: PetscErrorCode KSPSetInitialGuessNonzero(KSP ksp,PetscBool flg)
1008: {
1012: ksp->guess_zero = (PetscBool) !(int)flg;
1013: return(0);
1014: }
1018: /*@
1019: KSPGetInitialGuessNonzero - Determines whether the KSP solver is using
1020: a zero initial guess.
1022: Not Collective
1024: Input Parameter:
1025: . ksp - iterative context obtained from KSPCreate()
1027: Output Parameter:
1028: . flag - PETSC_TRUE if guess is nonzero, else PETSC_FALSE
1030: Level: intermediate
1032: .keywords: KSP, set, initial guess, nonzero
1034: .seealso: KSPSetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll()
1035: @*/
1036: PetscErrorCode KSPGetInitialGuessNonzero(KSP ksp,PetscBool *flag)
1037: {
1041: if (ksp->guess_zero) *flag = PETSC_FALSE;
1042: else *flag = PETSC_TRUE;
1043: return(0);
1044: }
1048: /*@
1049: KSPSetErrorIfNotConverged - Causes KSPSolve() to generate an error if the solver has not converged.
1051: Logically Collective on KSP
1053: Input Parameters:
1054: + ksp - iterative context obtained from KSPCreate()
1055: - flg - PETSC_TRUE indicates you want the error generated
1057: Options database keys:
1058: . -ksp_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false)
1060: Level: intermediate
1062: Notes:
1063: Normally PETSc continues if a linear solver fails to converge, you can call KSPGetConvergedReason() after a KSPSolve()
1064: to determine if it has converged.
1066: .keywords: KSP, set, initial guess, nonzero
1068: .seealso: KSPGetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll(), KSPGetErrorIfNotConverged()
1069: @*/
1070: PetscErrorCode KSPSetErrorIfNotConverged(KSP ksp,PetscBool flg)
1071: {
1075: ksp->errorifnotconverged = flg;
1076: return(0);
1077: }
1081: /*@
1082: KSPGetErrorIfNotConverged - Will KSPSolve() generate an error if the solver does not converge?
1084: Not Collective
1086: Input Parameter:
1087: . ksp - iterative context obtained from KSPCreate()
1089: Output Parameter:
1090: . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE
1092: Level: intermediate
1094: .keywords: KSP, set, initial guess, nonzero
1096: .seealso: KSPSetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll(), KSPSetErrorIfNotConverged()
1097: @*/
1098: PetscErrorCode KSPGetErrorIfNotConverged(KSP ksp,PetscBool *flag)
1099: {
1103: *flag = ksp->errorifnotconverged;
1104: return(0);
1105: }
1109: /*@
1110: KSPSetInitialGuessKnoll - Tells the iterative solver to use PCApply(pc,b,..) to compute the initial guess (The Knoll trick)
1112: Logically Collective on KSP
1114: Input Parameters:
1115: + ksp - iterative context obtained from KSPCreate()
1116: - flg - PETSC_TRUE or PETSC_FALSE
1118: Level: advanced
1121: .keywords: KSP, set, initial guess, nonzero
1123: .seealso: KSPGetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero()
1124: @*/
1125: PetscErrorCode KSPSetInitialGuessKnoll(KSP ksp,PetscBool flg)
1126: {
1130: ksp->guess_knoll = flg;
1131: return(0);
1132: }
1136: /*@
1137: KSPGetInitialGuessKnoll - Determines whether the KSP solver is using the Knoll trick (using PCApply(pc,b,...) to compute
1138: the initial guess
1140: Not Collective
1142: Input Parameter:
1143: . ksp - iterative context obtained from KSPCreate()
1145: Output Parameter:
1146: . flag - PETSC_TRUE if using Knoll trick, else PETSC_FALSE
1148: Level: advanced
1150: .keywords: KSP, set, initial guess, nonzero
1152: .seealso: KSPSetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero()
1153: @*/
1154: PetscErrorCode KSPGetInitialGuessKnoll(KSP ksp,PetscBool *flag)
1155: {
1159: *flag = ksp->guess_knoll;
1160: return(0);
1161: }
1165: /*@
1166: KSPGetComputeSingularValues - Gets the flag indicating whether the extreme singular
1167: values will be calculated via a Lanczos or Arnoldi process as the linear
1168: system is solved.
1170: Not Collective
1172: Input Parameter:
1173: . ksp - iterative context obtained from KSPCreate()
1175: Output Parameter:
1176: . flg - PETSC_TRUE or PETSC_FALSE
1178: Options Database Key:
1179: . -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()
1181: Notes:
1182: Currently this option is not valid for all iterative methods.
1184: Many users may just want to use the monitoring routine
1185: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1186: to print the singular values at each iteration of the linear solve.
1188: Level: advanced
1190: .keywords: KSP, set, compute, singular values
1192: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue()
1193: @*/
1194: PetscErrorCode KSPGetComputeSingularValues(KSP ksp,PetscBool *flg)
1195: {
1199: *flg = ksp->calc_sings;
1200: return(0);
1201: }
1205: /*@
1206: KSPSetComputeSingularValues - Sets a flag so that the extreme singular
1207: values will be calculated via a Lanczos or Arnoldi process as the linear
1208: system is solved.
1210: Logically Collective on KSP
1212: Input Parameters:
1213: + ksp - iterative context obtained from KSPCreate()
1214: - flg - PETSC_TRUE or PETSC_FALSE
1216: Options Database Key:
1217: . -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()
1219: Notes:
1220: Currently this option is not valid for all iterative methods.
1222: Many users may just want to use the monitoring routine
1223: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1224: to print the singular values at each iteration of the linear solve.
1226: Level: advanced
1228: .keywords: KSP, set, compute, singular values
1230: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue()
1231: @*/
1232: PetscErrorCode KSPSetComputeSingularValues(KSP ksp,PetscBool flg)
1233: {
1237: ksp->calc_sings = flg;
1238: return(0);
1239: }
1243: /*@
1244: KSPGetComputeEigenvalues - Gets the flag indicating that the extreme eigenvalues
1245: values will be calculated via a Lanczos or Arnoldi process as the linear
1246: system is solved.
1248: Not Collective
1250: Input Parameter:
1251: . ksp - iterative context obtained from KSPCreate()
1253: Output Parameter:
1254: . flg - PETSC_TRUE or PETSC_FALSE
1256: Notes:
1257: Currently this option is not valid for all iterative methods.
1259: Level: advanced
1261: .keywords: KSP, set, compute, eigenvalues
1263: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly()
1264: @*/
1265: PetscErrorCode KSPGetComputeEigenvalues(KSP ksp,PetscBool *flg)
1266: {
1270: *flg = ksp->calc_sings;
1271: return(0);
1272: }
1276: /*@
1277: KSPSetComputeEigenvalues - Sets a flag so that the extreme eigenvalues
1278: values will be calculated via a Lanczos or Arnoldi process as the linear
1279: system is solved.
1281: Logically Collective on KSP
1283: Input Parameters:
1284: + ksp - iterative context obtained from KSPCreate()
1285: - flg - PETSC_TRUE or PETSC_FALSE
1287: Notes:
1288: Currently this option is not valid for all iterative methods.
1290: Level: advanced
1292: .keywords: KSP, set, compute, eigenvalues
1294: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly()
1295: @*/
1296: PetscErrorCode KSPSetComputeEigenvalues(KSP ksp,PetscBool flg)
1297: {
1301: ksp->calc_sings = flg;
1302: return(0);
1303: }
1307: /*@
1308: KSPGetRhs - Gets the right-hand-side vector for the linear system to
1309: be solved.
1311: Not Collective
1313: Input Parameter:
1314: . ksp - iterative context obtained from KSPCreate()
1316: Output Parameter:
1317: . r - right-hand-side vector
1319: Level: developer
1321: .keywords: KSP, get, right-hand-side, rhs
1323: .seealso: KSPGetSolution(), KSPSolve()
1324: @*/
1325: PetscErrorCode KSPGetRhs(KSP ksp,Vec *r)
1326: {
1330: *r = ksp->vec_rhs;
1331: return(0);
1332: }
1336: /*@
1337: KSPGetSolution - Gets the location of the solution for the
1338: linear system to be solved. Note that this may not be where the solution
1339: is stored during the iterative process; see KSPBuildSolution().
1341: Not Collective
1343: Input Parameters:
1344: . ksp - iterative context obtained from KSPCreate()
1346: Output Parameters:
1347: . v - solution vector
1349: Level: developer
1351: .keywords: KSP, get, solution
1353: .seealso: KSPGetRhs(), KSPBuildSolution(), KSPSolve()
1354: @*/
1355: PetscErrorCode KSPGetSolution(KSP ksp,Vec *v)
1356: {
1360: *v = ksp->vec_sol;
1361: return(0);
1362: }
1366: /*@
1367: KSPSetPC - Sets the preconditioner to be used to calculate the
1368: application of the preconditioner on a vector.
1370: Collective on KSP
1372: Input Parameters:
1373: + ksp - iterative context obtained from KSPCreate()
1374: - pc - the preconditioner object
1376: Notes:
1377: Use KSPGetPC() to retrieve the preconditioner context (for example,
1378: to free it at the end of the computations).
1380: Level: developer
1382: .keywords: KSP, set, precondition, Binv
1384: .seealso: KSPGetPC()
1385: @*/
1386: PetscErrorCode KSPSetPC(KSP ksp,PC pc)
1387: {
1394: PetscObjectReference((PetscObject)pc);
1395: PCDestroy(&ksp->pc);
1396: ksp->pc = pc;
1397: PetscLogObjectParent(ksp,ksp->pc);
1398: return(0);
1399: }
1403: /*@
1404: KSPGetPC - Returns a pointer to the preconditioner context
1405: set with KSPSetPC().
1407: Not Collective
1409: Input Parameters:
1410: . ksp - iterative context obtained from KSPCreate()
1412: Output Parameter:
1413: . pc - preconditioner context
1415: Level: developer
1417: .keywords: KSP, get, preconditioner, Binv
1419: .seealso: KSPSetPC()
1420: @*/
1421: PetscErrorCode KSPGetPC(KSP ksp,PC *pc)
1422: {
1428: if (!ksp->pc) {
1429: PCCreate(PetscObjectComm((PetscObject)ksp),&ksp->pc);
1430: PetscObjectIncrementTabLevel((PetscObject)ksp->pc,(PetscObject)ksp,0);
1431: PetscLogObjectParent(ksp,ksp->pc);
1432: }
1433: *pc = ksp->pc;
1434: return(0);
1435: }
1439: /*@
1440: KSPMonitor - runs the user provided monitor routines, if they exist
1442: Collective on KSP
1444: Input Parameters:
1445: + ksp - iterative context obtained from KSPCreate()
1446: . it - iteration number
1447: - rnorm - relative norm of the residual
1449: Notes:
1450: This routine is called by the KSP implementations.
1451: It does not typically need to be called by the user.
1453: Level: developer
1455: .seealso: KSPMonitorSet()
1456: @*/
1457: PetscErrorCode KSPMonitor(KSP ksp,PetscInt it,PetscReal rnorm)
1458: {
1459: PetscInt i, n = ksp->numbermonitors;
1463: for (i=0; i<n; i++) {
1464: (*ksp->monitor[i])(ksp,it,rnorm,ksp->monitorcontext[i]);
1465: }
1466: return(0);
1467: }
1471: /*@C
1472: KSPMonitorSet - Sets an ADDITIONAL function to be called at every iteration to monitor
1473: the residual/error etc.
1475: Logically Collective on KSP
1477: Input Parameters:
1478: + ksp - iterative context obtained from KSPCreate()
1479: . monitor - pointer to function (if this is NULL, it turns off monitoring
1480: . mctx - [optional] context for private data for the
1481: monitor routine (use NULL if no context is desired)
1482: - monitordestroy - [optional] routine that frees monitor context
1483: (may be NULL)
1485: Calling Sequence of monitor:
1486: $ monitor (KSP ksp, int it, PetscReal rnorm, void *mctx)
1488: + ksp - iterative context obtained from KSPCreate()
1489: . it - iteration number
1490: . rnorm - (estimated) 2-norm of (preconditioned) residual
1491: - mctx - optional monitoring context, as set by KSPMonitorSet()
1493: Options Database Keys:
1494: + -ksp_monitor - sets KSPMonitorDefault()
1495: . -ksp_monitor_true_residual - sets KSPMonitorTrueResidualNorm()
1496: . -ksp_monitor_max - sets KSPMonitorTrueResidualMaxNorm()
1497: . -ksp_monitor_lg_residualnorm - sets line graph monitor,
1498: uses KSPMonitorLGResidualNormCreate()
1499: . -ksp_monitor_lg_true_residualnorm - sets line graph monitor,
1500: uses KSPMonitorLGResidualNormCreate()
1501: . -ksp_monitor_singular_value - sets KSPMonitorSingularValue()
1502: - -ksp_monitor_cancel - cancels all monitors that have
1503: been hardwired into a code by
1504: calls to KSPMonitorSet(), but
1505: does not cancel those set via
1506: the options database.
1508: Notes:
1509: The default is to do nothing. To print the residual, or preconditioned
1510: residual if KSPSetNormType(ksp,KSP_NORM_PRECONDITIONED) was called, use
1511: KSPMonitorDefault() as the monitoring routine, with a null monitoring
1512: context.
1514: Several different monitoring routines may be set by calling
1515: KSPMonitorSet() multiple times; all will be called in the
1516: order in which they were set.
1518: Fortran notes: Only a single monitor function can be set for each KSP object
1520: Level: beginner
1522: .keywords: KSP, set, monitor
1524: .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSPMonitorCancel()
1525: @*/
1526: PetscErrorCode KSPMonitorSet(KSP ksp,PetscErrorCode (*monitor)(KSP,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
1527: {
1528: PetscInt i;
1533: if (ksp->numbermonitors >= MAXKSPMONITORS) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Too many KSP monitors set");
1534: for (i=0; i<ksp->numbermonitors;i++) {
1535: if (monitor == ksp->monitor[i] && monitordestroy == ksp->monitordestroy[i] && mctx == ksp->monitorcontext[i]) {
1536: if (monitordestroy) {
1537: (*monitordestroy)(&mctx);
1538: }
1539: return(0);
1540: }
1541: }
1542: ksp->monitor[ksp->numbermonitors] = monitor;
1543: ksp->monitordestroy[ksp->numbermonitors] = monitordestroy;
1544: ksp->monitorcontext[ksp->numbermonitors++] = (void*)mctx;
1545: return(0);
1546: }
1550: /*@
1551: KSPMonitorCancel - Clears all monitors for a KSP object.
1553: Logically Collective on KSP
1555: Input Parameters:
1556: . ksp - iterative context obtained from KSPCreate()
1558: Options Database Key:
1559: . -ksp_monitor_cancel - Cancels all monitors that have
1560: been hardwired into a code by calls to KSPMonitorSet(),
1561: but does not cancel those set via the options database.
1563: Level: intermediate
1565: .keywords: KSP, set, monitor
1567: .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSPMonitorSet()
1568: @*/
1569: PetscErrorCode KSPMonitorCancel(KSP ksp)
1570: {
1572: PetscInt i;
1576: for (i=0; i<ksp->numbermonitors; i++) {
1577: if (ksp->monitordestroy[i]) {
1578: (*ksp->monitordestroy[i])(&ksp->monitorcontext[i]);
1579: }
1580: }
1581: ksp->numbermonitors = 0;
1582: return(0);
1583: }
1587: /*@C
1588: KSPGetMonitorContext - Gets the monitoring context, as set by
1589: KSPMonitorSet() for the FIRST monitor only.
1591: Not Collective
1593: Input Parameter:
1594: . ksp - iterative context obtained from KSPCreate()
1596: Output Parameter:
1597: . ctx - monitoring context
1599: Level: intermediate
1601: .keywords: KSP, get, monitor, context
1603: .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate()
1604: @*/
1605: PetscErrorCode KSPGetMonitorContext(KSP ksp,void **ctx)
1606: {
1609: *ctx = (ksp->monitorcontext[0]);
1610: return(0);
1611: }
1615: /*@
1616: KSPSetResidualHistory - Sets the array used to hold the residual history.
1617: If set, this array will contain the residual norms computed at each
1618: iteration of the solver.
1620: Not Collective
1622: Input Parameters:
1623: + ksp - iterative context obtained from KSPCreate()
1624: . a - array to hold history
1625: . na - size of a
1626: - reset - PETSC_TRUE indicates the history counter is reset to zero
1627: for each new linear solve
1629: Level: advanced
1631: Notes: The array is NOT freed by PETSc so the user needs to keep track of
1632: it and destroy once the KSP object is destroyed.
1634: If 'a' is NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
1635: default array of length 10000 is allocated.
1637: .keywords: KSP, set, residual, history, norm
1639: .seealso: KSPGetResidualHistory()
1641: @*/
1642: PetscErrorCode KSPSetResidualHistory(KSP ksp,PetscReal a[],PetscInt na,PetscBool reset)
1643: {
1649: PetscFree(ksp->res_hist_alloc);
1650: if (na != PETSC_DECIDE && na != PETSC_DEFAULT && a) {
1651: ksp->res_hist = a;
1652: ksp->res_hist_max = na;
1653: } else {
1654: if (na != PETSC_DECIDE && na != PETSC_DEFAULT) ksp->res_hist_max = na;
1655: else ksp->res_hist_max = 10000; /* like default ksp->max_it */
1656: PetscMalloc(ksp->res_hist_max*sizeof(PetscReal),&ksp->res_hist_alloc);
1658: ksp->res_hist = ksp->res_hist_alloc;
1659: }
1660: ksp->res_hist_len = 0;
1661: ksp->res_hist_reset = reset;
1662: return(0);
1663: }
1667: /*@C
1668: KSPGetResidualHistory - Gets the array used to hold the residual history
1669: and the number of residuals it contains.
1671: Not Collective
1673: Input Parameter:
1674: . ksp - iterative context obtained from KSPCreate()
1676: Output Parameters:
1677: + a - pointer to array to hold history (or NULL)
1678: - na - number of used entries in a (or NULL)
1680: Level: advanced
1682: Notes:
1683: Can only be called after a KSPSetResidualHistory() otherwise a and na are set to zero
1685: The Fortran version of this routine has a calling sequence
1686: $ call KSPGetResidualHistory(KSP ksp, integer na, integer ierr)
1687: note that you have passed a Fortran array into KSPSetResidualHistory() and you need
1688: to access the residual values from this Fortran array you provided. Only the na (number of
1689: residual norms currently held) is set.
1691: .keywords: KSP, get, residual, history, norm
1693: .seealso: KSPGetResidualHistory()
1695: @*/
1696: PetscErrorCode KSPGetResidualHistory(KSP ksp,PetscReal *a[],PetscInt *na)
1697: {
1700: if (a) *a = ksp->res_hist;
1701: if (na) *na = ksp->res_hist_len;
1702: return(0);
1703: }
1707: /*@C
1708: KSPSetConvergenceTest - Sets the function to be used to determine
1709: convergence.
1711: Logically Collective on KSP
1713: Input Parameters:
1714: + ksp - iterative context obtained from KSPCreate()
1715: . converge - pointer to int function
1716: . cctx - context for private data for the convergence routine (may be null)
1717: - destroy - a routine for destroying the context (may be null)
1719: Calling sequence of converge:
1720: $ converge (KSP ksp, int it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)
1722: + ksp - iterative context obtained from KSPCreate()
1723: . it - iteration number
1724: . rnorm - (estimated) 2-norm of (preconditioned) residual
1725: . reason - the reason why it has converged or diverged
1726: - cctx - optional convergence context, as set by KSPSetConvergenceTest()
1729: Notes:
1730: Must be called after the KSP type has been set so put this after
1731: a call to KSPSetType(), or KSPSetFromOptions().
1733: The default convergence test, KSPDefaultConverged(), aborts if the
1734: residual grows to more than 10000 times the initial residual.
1736: The default is a combination of relative and absolute tolerances.
1737: The residual value that is tested may be an approximation; routines
1738: that need exact values should compute them.
1740: In the default PETSc convergence test, the precise values of reason
1741: are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
1743: Level: advanced
1745: .keywords: KSP, set, convergence, test, context
1747: .seealso: KSPDefaultConverged(), KSPGetConvergenceContext(), KSPSetTolerances()
1748: @*/
1749: PetscErrorCode KSPSetConvergenceTest(KSP ksp,PetscErrorCode (*converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
1750: {
1755: if (ksp->convergeddestroy) {
1756: (*ksp->convergeddestroy)(ksp->cnvP);
1757: }
1758: ksp->converged = converge;
1759: ksp->convergeddestroy = destroy;
1760: ksp->cnvP = (void*)cctx;
1761: return(0);
1762: }
1766: /*@C
1767: KSPGetConvergenceContext - Gets the convergence context set with
1768: KSPSetConvergenceTest().
1770: Not Collective
1772: Input Parameter:
1773: . ksp - iterative context obtained from KSPCreate()
1775: Output Parameter:
1776: . ctx - monitoring context
1778: Level: advanced
1780: .keywords: KSP, get, convergence, test, context
1782: .seealso: KSPDefaultConverged(), KSPSetConvergenceTest()
1783: @*/
1784: PetscErrorCode KSPGetConvergenceContext(KSP ksp,void **ctx)
1785: {
1788: *ctx = ksp->cnvP;
1789: return(0);
1790: }
1794: /*@C
1795: KSPBuildSolution - Builds the approximate solution in a vector provided.
1796: This routine is NOT commonly needed (see KSPSolve()).
1798: Collective on KSP
1800: Input Parameter:
1801: . ctx - iterative context obtained from KSPCreate()
1803: Output Parameter:
1804: Provide exactly one of
1805: + v - location to stash solution.
1806: - V - the solution is returned in this location. This vector is created
1807: internally. This vector should NOT be destroyed by the user with
1808: VecDestroy().
1810: Notes:
1811: This routine can be used in one of two ways
1812: .vb
1813: KSPBuildSolution(ksp,NULL,&V);
1814: or
1815: KSPBuildSolution(ksp,v,NULL); or KSPBuildSolution(ksp,v,&v);
1816: .ve
1817: In the first case an internal vector is allocated to store the solution
1818: (the user cannot destroy this vector). In the second case the solution
1819: is generated in the vector that the user provides. Note that for certain
1820: methods, such as KSPCG, the second case requires a copy of the solution,
1821: while in the first case the call is essentially free since it simply
1822: returns the vector where the solution already is stored. For some methods
1823: like GMRES this is a reasonably expensive operation and should only be
1824: used in truly needed.
1826: Level: advanced
1828: .keywords: KSP, build, solution
1830: .seealso: KSPGetSolution(), KSPBuildResidual()
1831: @*/
1832: PetscErrorCode KSPBuildSolution(KSP ksp,Vec v,Vec *V)
1833: {
1838: if (!V && !v) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONG,"Must provide either v or V");
1839: if (!V) V = &v;
1840: (*ksp->ops->buildsolution)(ksp,v,V);
1841: return(0);
1842: }
1846: /*@C
1847: KSPBuildResidual - Builds the residual in a vector provided.
1849: Collective on KSP
1851: Input Parameter:
1852: . ksp - iterative context obtained from KSPCreate()
1854: Output Parameters:
1855: + v - optional location to stash residual. If v is not provided,
1856: then a location is generated.
1857: . t - work vector. If not provided then one is generated.
1858: - V - the residual
1860: Notes:
1861: Regardless of whether or not v is provided, the residual is
1862: returned in V.
1864: Level: advanced
1866: .keywords: KSP, build, residual
1868: .seealso: KSPBuildSolution()
1869: @*/
1870: PetscErrorCode KSPBuildResidual(KSP ksp,Vec t,Vec v,Vec *V)
1871: {
1873: PetscBool flag = PETSC_FALSE;
1874: Vec w = v,tt = t;
1878: if (!w) {
1879: VecDuplicate(ksp->vec_rhs,&w);
1880: PetscLogObjectParent((PetscObject)ksp,w);
1881: }
1882: if (!tt) {
1883: VecDuplicate(ksp->vec_sol,&tt); flag = PETSC_TRUE;
1884: PetscLogObjectParent((PetscObject)ksp,tt);
1885: }
1886: (*ksp->ops->buildresidual)(ksp,tt,w,V);
1887: if (flag) {VecDestroy(&tt);}
1888: return(0);
1889: }
1893: /*@
1894: KSPSetDiagonalScale - Tells KSP to symmetrically diagonally scale the system
1895: before solving. This actually CHANGES the matrix (and right hand side).
1897: Logically Collective on KSP
1899: Input Parameter:
1900: + ksp - the KSP context
1901: - scale - PETSC_TRUE or PETSC_FALSE
1903: Options Database Key:
1904: + -ksp_diagonal_scale -
1905: - -ksp_diagonal_scale_fix - scale the matrix back AFTER the solve
1908: Notes: Scales the matrix by D^(-1/2) A D^(-1/2) [D^(1/2) x ] = D^(-1/2) b
1909: where D_{ii} is 1/abs(A_{ii}) unless A_{ii} is zero and then it is 1.
1911: BE CAREFUL with this routine: it actually scales the matrix and right
1912: hand side that define the system. After the system is solved the matrix
1913: and right hand side remain scaled unless you use KSPSetDiagonalScaleFix()
1915: This should NOT be used within the SNES solves if you are using a line
1916: search.
1918: If you use this with the PCType Eisenstat preconditioner than you can
1919: use the PCEisenstatNoDiagonalScaling() option, or -pc_eisenstat_no_diagonal_scaling
1920: to save some unneeded, redundant flops.
1922: Level: intermediate
1924: .keywords: KSP, set, options, prefix, database
1926: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScaleFix()
1927: @*/
1928: PetscErrorCode KSPSetDiagonalScale(KSP ksp,PetscBool scale)
1929: {
1933: ksp->dscale = scale;
1934: return(0);
1935: }
1939: /*@
1940: KSPGetDiagonalScale - Checks if KSP solver scales the matrix and
1941: right hand side
1943: Not Collective
1945: Input Parameter:
1946: . ksp - the KSP context
1948: Output Parameter:
1949: . scale - PETSC_TRUE or PETSC_FALSE
1951: Notes:
1952: BE CAREFUL with this routine: it actually scales the matrix and right
1953: hand side that define the system. After the system is solved the matrix
1954: and right hand side remain scaled unless you use KSPSetDiagonalScaleFix()
1956: Level: intermediate
1958: .keywords: KSP, set, options, prefix, database
1960: .seealso: KSPSetDiagonalScale(), KSPSetDiagonalScaleFix()
1961: @*/
1962: PetscErrorCode KSPGetDiagonalScale(KSP ksp,PetscBool *scale)
1963: {
1967: *scale = ksp->dscale;
1968: return(0);
1969: }
1973: /*@
1974: KSPSetDiagonalScaleFix - Tells KSP to diagonally scale the system
1975: back after solving.
1977: Logically Collective on KSP
1979: Input Parameter:
1980: + ksp - the KSP context
1981: - fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
1982: rescale (default)
1984: Notes:
1985: Must be called after KSPSetDiagonalScale()
1987: Using this will slow things down, because it rescales the matrix before and
1988: after each linear solve. This is intended mainly for testing to allow one
1989: to easily get back the original system to make sure the solution computed is
1990: accurate enough.
1992: Level: intermediate
1994: .keywords: KSP, set, options, prefix, database
1996: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPGetDiagonalScaleFix()
1997: @*/
1998: PetscErrorCode KSPSetDiagonalScaleFix(KSP ksp,PetscBool fix)
1999: {
2003: ksp->dscalefix = fix;
2004: return(0);
2005: }
2009: /*@
2010: KSPGetDiagonalScaleFix - Determines if KSP diagonally scales the system
2011: back after solving.
2013: Not Collective
2015: Input Parameter:
2016: . ksp - the KSP context
2018: Output Parameter:
2019: . fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
2020: rescale (default)
2022: Notes:
2023: Must be called after KSPSetDiagonalScale()
2025: If PETSC_TRUE will slow things down, because it rescales the matrix before and
2026: after each linear solve. This is intended mainly for testing to allow one
2027: to easily get back the original system to make sure the solution computed is
2028: accurate enough.
2030: Level: intermediate
2032: .keywords: KSP, set, options, prefix, database
2034: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPSetDiagonalScaleFix()
2035: @*/
2036: PetscErrorCode KSPGetDiagonalScaleFix(KSP ksp,PetscBool *fix)
2037: {
2041: *fix = ksp->dscalefix;
2042: return(0);
2043: }
2047: /*@C
2048: KSPSetComputeOperators - set routine to compute the linear operators
2050: Logically Collective
2052: Input Arguments:
2053: + ksp - the KSP context
2054: . func - function to compute the operators
2055: - ctx - optional context
2057: Calling sequence of func:
2058: $ func(KSP ksp,Mat *A,Mat *B,MatStructure *mstruct,void *ctx)
2060: + ksp - the KSP context
2061: . A - the linear operator
2062: . B - preconditioning matrix
2063: . mstruct - flag indicating structure, same as in KSPSetOperators(), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER
2064: - ctx - optional user-provided context
2066: Level: beginner
2068: .seealso: KSPSetOperators(), DMKSPSetComputeOperators()
2069: @*/
2070: PetscErrorCode KSPSetComputeOperators(KSP ksp,PetscErrorCode (*func)(KSP,Mat,Mat,MatStructure*,void*),void *ctx)
2071: {
2073: DM dm;
2077: KSPGetDM(ksp,&dm);
2078: DMKSPSetComputeOperators(dm,func,ctx);
2079: if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX;
2080: return(0);
2081: }
2085: /*@C
2086: KSPSetComputeRHS - set routine to compute the right hand side of the linear system
2088: Logically Collective
2090: Input Arguments:
2091: + ksp - the KSP context
2092: . func - function to compute the right hand side
2093: - ctx - optional context
2095: Calling sequence of func:
2096: $ func(KSP ksp,Vec b,void *ctx)
2098: + ksp - the KSP context
2099: . b - right hand side of linear system
2100: - ctx - optional user-provided context
2102: Level: beginner
2104: .seealso: KSPSolve(), DMKSPSetComputeRHS()
2105: @*/
2106: PetscErrorCode KSPSetComputeRHS(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2107: {
2109: DM dm;
2113: KSPGetDM(ksp,&dm);
2114: DMKSPSetComputeRHS(dm,func,ctx);
2115: return(0);
2116: }
2120: /*@C
2121: KSPSetComputeInitialGuess - set routine to compute the initial guess of the linear system
2123: Logically Collective
2125: Input Arguments:
2126: + ksp - the KSP context
2127: . func - function to compute the initial guess
2128: - ctx - optional context
2130: Calling sequence of func:
2131: $ func(KSP ksp,Vec x,void *ctx)
2133: + ksp - the KSP context
2134: . x - solution vector
2135: - ctx - optional user-provided context
2137: Level: beginner
2139: .seealso: KSPSolve(), DMKSPSetComputeInitialGuess()
2140: @*/
2141: PetscErrorCode KSPSetComputeInitialGuess(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2142: {
2144: DM dm;
2148: KSPGetDM(ksp,&dm);
2149: DMKSPSetComputeInitialGuess(dm,func,ctx);
2150: return(0);
2151: }