Actual source code: itfunc.c
2: /*
3: Interface KSP routines that the user calls.
4: */
6: #include <private/kspimpl.h> /*I "petscksp.h" I*/
10: /*@
11: KSPComputeExtremeSingularValues - Computes the extreme singular values
12: for the preconditioned operator. Called after or during KSPSolve().
14: Not Collective
16: Input Parameter:
17: . ksp - iterative context obtained from KSPCreate()
19: Output Parameters:
20: . emin, emax - extreme singular values
22: Notes:
23: One must call KSPSetComputeSingularValues() before calling KSPSetUp()
24: (or use the option -ksp_compute_eigenvalues) in order for this routine to work correctly.
26: Many users may just want to use the monitoring routine
27: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
28: to print the extreme singular values at each iteration of the linear solve.
30: Level: advanced
32: .keywords: KSP, compute, extreme, singular, values
34: .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeEigenvalues()
35: @*/
36: PetscErrorCode KSPComputeExtremeSingularValues(KSP ksp,PetscReal *emax,PetscReal *emin)
37: {
44: if (!ksp->calc_sings) SETERRQ(((PetscObject)ksp)->comm,4,"Singular values not requested before KSPSetUp()");
46: if (ksp->ops->computeextremesingularvalues) {
47: (*ksp->ops->computeextremesingularvalues)(ksp,emax,emin);
48: } else {
49: *emin = -1.0;
50: *emax = -1.0;
51: }
52: return(0);
53: }
57: /*@
58: KSPComputeEigenvalues - Computes the extreme eigenvalues for the
59: preconditioned operator. Called after or during KSPSolve().
61: Not Collective
63: Input Parameter:
64: + ksp - iterative context obtained from KSPCreate()
65: - n - size of arrays r and c. The number of eigenvalues computed (neig) will, in
66: general, be less than this.
68: Output Parameters:
69: + r - real part of computed eigenvalues
70: . c - complex part of computed eigenvalues
71: - neig - number of eigenvalues computed (will be less than or equal to n)
73: Options Database Keys:
74: + -ksp_compute_eigenvalues - Prints eigenvalues to stdout
75: - -ksp_plot_eigenvalues - Plots eigenvalues in an x-window display
77: Notes:
78: The number of eigenvalues estimated depends on the size of the Krylov space
79: generated during the KSPSolve() ; for example, with
80: CG it corresponds to the number of CG iterations, for GMRES it is the number
81: of GMRES iterations SINCE the last restart. Any extra space in r[] and c[]
82: will be ignored.
84: KSPComputeEigenvalues() does not usually provide accurate estimates; it is
85: intended only for assistance in understanding the convergence of iterative
86: methods, not for eigenanalysis.
88: One must call KSPSetComputeEigenvalues() before calling KSPSetUp()
89: in order for this routine to work correctly.
91: Many users may just want to use the monitoring routine
92: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
93: to print the singular values at each iteration of the linear solve.
95: Level: advanced
97: .keywords: KSP, compute, extreme, singular, values
99: .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues()
100: @*/
101: PetscErrorCode KSPComputeEigenvalues(KSP ksp,PetscInt n,PetscReal *r,PetscReal *c,PetscInt *neig)
102: {
110: if (!ksp->calc_sings) SETERRQ(((PetscObject)ksp)->comm,4,"Eigenvalues not requested before KSPSetUp()");
112: if (ksp->ops->computeeigenvalues) {
113: (*ksp->ops->computeeigenvalues)(ksp,n,r,c,neig);
114: } else {
115: *neig = 0;
116: }
117: return(0);
118: }
122: /*@
123: KSPSetUpOnBlocks - Sets up the preconditioner for each block in
124: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
125: methods.
127: Collective on KSP
129: Input Parameter:
130: . ksp - the KSP context
132: Notes:
133: KSPSetUpOnBlocks() is a routine that the user can optinally call for
134: more precise profiling (via -log_summary) of the setup phase for these
135: block preconditioners. If the user does not call KSPSetUpOnBlocks(),
136: it will automatically be called from within KSPSolve().
137:
138: Calling KSPSetUpOnBlocks() is the same as calling PCSetUpOnBlocks()
139: on the PC context within the KSP context.
141: Level: advanced
143: .keywords: KSP, setup, blocks
145: .seealso: PCSetUpOnBlocks(), KSPSetUp(), PCSetUp()
146: @*/
147: PetscErrorCode KSPSetUpOnBlocks(KSP ksp)
148: {
153: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
154: PCSetUpOnBlocks(ksp->pc);
155: return(0);
156: }
160: /*@
161: KSPSetUp - Sets up the internal data structures for the
162: later use of an iterative solver.
164: Collective on KSP
166: Input Parameter:
167: . ksp - iterative context obtained from KSPCreate()
169: Level: developer
171: .keywords: KSP, setup
173: .seealso: KSPCreate(), KSPSolve(), KSPDestroy()
174: @*/
175: PetscErrorCode KSPSetUp(KSP ksp)
176: {
178: PetscBool ir = PETSC_FALSE,ig = PETSC_FALSE;
179: Mat A;
180: MatStructure stflg = SAME_NONZERO_PATTERN;
182: /* PetscBool im = PETSC_FALSE; */
187: /* reset the convergence flag from the previous solves */
188: ksp->reason = KSP_CONVERGED_ITERATING;
190: if (!((PetscObject)ksp)->type_name){
191: KSPSetType(ksp,KSPGMRES);
192: }
193: KSPSetUpNorms_Private(ksp);
195: if (ksp->dmActive && !ksp->setupstage) {
196: /* first time in so build matrix and vector data structures using DM */
197: if (!ksp->vec_rhs) {DMCreateGlobalVector(ksp->dm,&ksp->vec_rhs);}
198: if (!ksp->vec_sol) {DMCreateGlobalVector(ksp->dm,&ksp->vec_sol);}
199: DMCreateMatrix(ksp->dm,MATAIJ,&A);
200: KSPSetOperators(ksp,A,A,stflg);
201: PetscObjectDereference((PetscObject)A);
202: }
204: if (ksp->dmActive) {
205: DMHasInitialGuess(ksp->dm,&ig);
206: if (ig && ksp->setupstage != KSP_SETUP_NEWRHS) {
207: /* only computes initial guess the first time through */
208: DMComputeInitialGuess(ksp->dm,ksp->vec_sol);
209: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
210: }
211: DMHasFunction(ksp->dm,&ir);
212: if (ir) {
213: DMComputeFunction(ksp->dm,PETSC_NULL,ksp->vec_rhs);
214: }
216: if (ksp->setupstage != KSP_SETUP_NEWRHS) {
217: KSPGetOperators(ksp,&A,&A,PETSC_NULL);
218: DMComputeJacobian(ksp->dm,PETSC_NULL,A,A,&stflg);
219: KSPSetOperators(ksp,A,A,stflg);
220: }
221: }
223: if (ksp->setupstage == KSP_SETUP_NEWRHS) return(0);
224: PetscLogEventBegin(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);
226: if (!ksp->setupstage) {
227: (*ksp->ops->setup)(ksp);
228: }
230: /* scale the matrix if requested */
231: if (ksp->dscale) {
232: Mat mat,pmat;
233: PetscScalar *xx;
234: PetscInt i,n;
235: PetscBool zeroflag = PETSC_FALSE;
236: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
237: PCGetOperators(ksp->pc,&mat,&pmat,PETSC_NULL);
238: if (!ksp->diagonal) { /* allocate vector to hold diagonal */
239: MatGetVecs(pmat,&ksp->diagonal,0);
240: }
241: MatGetDiagonal(pmat,ksp->diagonal);
242: VecGetLocalSize(ksp->diagonal,&n);
243: VecGetArray(ksp->diagonal,&xx);
244: for (i=0; i<n; i++) {
245: if (xx[i] != 0.0) xx[i] = 1.0/PetscSqrtReal(PetscAbsScalar(xx[i]));
246: else {
247: xx[i] = 1.0;
248: zeroflag = PETSC_TRUE;
249: }
250: }
251: VecRestoreArray(ksp->diagonal,&xx);
252: if (zeroflag) {
253: PetscInfo(ksp,"Zero detected in diagonal of matrix, using 1 at those locations\n");
254: }
255: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
256: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
257: ksp->dscalefix2 = PETSC_FALSE;
258: }
259: PetscLogEventEnd(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);
260: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
261: PCSetUp(ksp->pc);
262: if (ksp->nullsp) {
263: PetscBool test = PETSC_FALSE;
264: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_test_null_space",&test,PETSC_NULL);
265: if (test) {
266: Mat mat;
267: PCGetOperators(ksp->pc,&mat,PETSC_NULL,PETSC_NULL);
268: MatNullSpaceTest(ksp->nullsp,mat,PETSC_NULL);
269: }
270: }
271: ksp->setupstage = KSP_SETUP_NEWRHS;
272: return(0);
273: }
277: /*@
278: KSPSolve - Solves linear system.
280: Collective on KSP
282: Parameter:
283: + ksp - iterative context obtained from KSPCreate()
284: . b - the right hand side vector
285: - x - the solution (this may be the same vector as b, then b will be overwritten with answer)
287: Options Database Keys:
288: + -ksp_compute_eigenvalues - compute preconditioned operators eigenvalues
289: . -ksp_plot_eigenvalues - plot the computed eigenvalues in an X-window
290: . -ksp_compute_eigenvalues_explicitly - compute the eigenvalues by forming the dense operator and useing LAPACK
291: . -ksp_plot_eigenvalues_explicitly - plot the explicitly computing eigenvalues
292: . -ksp_view_binary - save matrix and right hand side that define linear system to the default binary viewer (can be
293: read later with src/ksp/examples/tutorials/ex10.c for testing solvers)
294: . -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations
295: . -ksp_final_residual - print 2-norm of true linear system residual at the end of the solution process
296: - -ksp_view - print the ksp data structure at the end of the system solution
298: Notes:
300: If one uses KSPSetDM() then x or b need not be passed. Use KSPGetSolution() to access the solution in this case.
302: The operator is specified with KSPSetOperators().
304: Call KSPGetConvergedReason() to determine if the solver converged or failed and
305: why. The number of iterations can be obtained from KSPGetIterationNumber().
306:
307: If using a direct method (e.g., via the KSP solver
308: KSPPREONLY and a preconditioner such as PCLU/PCILU),
309: then its=1. See KSPSetTolerances() and KSPDefaultConverged()
310: for more details.
312: Understanding Convergence:
313: The routines KSPMonitorSet(), KSPComputeEigenvalues(), and
314: KSPComputeEigenvaluesExplicitly() provide information on additional
315: options to monitor convergence and print eigenvalue information.
317: Level: beginner
319: .keywords: KSP, solve, linear system
321: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPDefaultConverged(),
322: KSPSolveTranspose(), KSPGetIterationNumber()
323: @*/
324: PetscErrorCode KSPSolve(KSP ksp,Vec b,Vec x)
325: {
327: PetscMPIInt rank;
328: PetscBool flag1,flag2,flg = PETSC_FALSE,inXisinB=PETSC_FALSE,guess_zero;
329: char view[10];
330: char filename[PETSC_MAX_PATH_LEN];
331: PetscViewer viewer;
332:
339: if (x && x == b) {
340: if (!ksp->guess_zero) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_INCOMP,"Cannot use x == b with nonzero initial guess");
341: VecDuplicate(b,&x);
342: inXisinB = PETSC_TRUE;
343: }
344: if (b) {
345: PetscObjectReference((PetscObject)b);
346: VecDestroy(&ksp->vec_rhs);
347: ksp->vec_rhs = b;
348: }
349: if (x) {
350: PetscObjectReference((PetscObject)x);
351: VecDestroy(&ksp->vec_sol);
352: ksp->vec_sol = x;
353: }
355: flag1 = flag2 = PETSC_FALSE;
356: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_view_binary",&flag1,PETSC_NULL);
357: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_view_binary_pre",&flag2,PETSC_NULL);
358: if (flag1 || flag2) {
359: Mat mat,premat;
360: PetscViewer viewer = PETSC_VIEWER_BINARY_(((PetscObject)ksp)->comm);
361: PCGetOperators(ksp->pc,&mat,&premat,PETSC_NULL);
362: if (flag1) {MatView(mat,viewer);}
363: if (flag2) {MatView(premat,viewer);}
364: VecView(ksp->vec_rhs,viewer);
365: }
366: PetscLogEventBegin(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
368: /* reset the residual history list if requested */
369: if (ksp->res_hist_reset) ksp->res_hist_len = 0;
371: PetscOptionsGetString(((PetscObject)ksp)->prefix,"-ksp_view_before",view,10,&flg);
372: if (flg) {
373: PetscViewer viewer;
374: PetscViewerASCIIGetStdout(((PetscObject)ksp)->comm,&viewer);
375: KSPView(ksp,viewer);
376: }
378: ksp->transpose_solve = PETSC_FALSE;
380: if (ksp->guess) {
381: KSPFischerGuessFormGuess(ksp->guess,ksp->vec_rhs,ksp->vec_sol);
382: ksp->guess_zero = PETSC_FALSE;
383: }
384: /* KSPSetUp() scales the matrix if needed */
385: KSPSetUp(ksp);
386: KSPSetUpOnBlocks(ksp);
388: /* diagonal scale RHS if called for */
389: if (ksp->dscale) {
390: VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
391: /* second time in, but matrix was scaled back to original */
392: if (ksp->dscalefix && ksp->dscalefix2) {
393: Mat mat,pmat;
395: PCGetOperators(ksp->pc,&mat,&pmat,PETSC_NULL);
396: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
397: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
398: }
400: /* scale initial guess */
401: if (!ksp->guess_zero) {
402: if (!ksp->truediagonal) {
403: VecDuplicate(ksp->diagonal,&ksp->truediagonal);
404: VecCopy(ksp->diagonal,ksp->truediagonal);
405: VecReciprocal(ksp->truediagonal);
406: }
407: VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->truediagonal);
408: }
409: }
410: PCPreSolve(ksp->pc,ksp);
412: if (ksp->guess_zero) { VecSet(ksp->vec_sol,0.0);}
413: if (ksp->guess_knoll) {
414: PCApply(ksp->pc,ksp->vec_rhs,ksp->vec_sol);
415: KSP_RemoveNullSpace(ksp,ksp->vec_sol);
416: ksp->guess_zero = PETSC_FALSE;
417: }
419: /* can we mark the initial guess as zero for this solve? */
420: guess_zero = ksp->guess_zero;
421: if (!ksp->guess_zero) {
422: PetscReal norm;
424: VecNormAvailable(ksp->vec_sol,NORM_2,&flg,&norm);
425: if (flg && !norm) {
426: ksp->guess_zero = PETSC_TRUE;
427: }
428: }
429: (*ksp->ops->solve)(ksp);
430: ksp->guess_zero = guess_zero;
432: if (!ksp->reason) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
433: if (ksp->printreason) {
434: PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm),((PetscObject)ksp)->tablevel);
435: if (ksp->reason > 0) {
436: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm),"Linear solve converged due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
437: } else {
438: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm),"Linear solve did not converge due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
439: }
440: PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm),((PetscObject)ksp)->tablevel);
441: }
442: PCPostSolve(ksp->pc,ksp);
444: /* diagonal scale solution if called for */
445: if (ksp->dscale) {
446: VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->diagonal);
447: /* unscale right hand side and matrix */
448: if (ksp->dscalefix) {
449: Mat mat,pmat;
451: VecReciprocal(ksp->diagonal);
452: VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
453: PCGetOperators(ksp->pc,&mat,&pmat,PETSC_NULL);
454: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
455: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
456: VecReciprocal(ksp->diagonal);
457: ksp->dscalefix2 = PETSC_TRUE;
458: }
459: }
460: PetscLogEventEnd(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
462: if (ksp->guess) {
463: KSPFischerGuessUpdate(ksp->guess,ksp->vec_sol);
464: }
466: MPI_Comm_rank(((PetscObject)ksp)->comm,&rank);
468: flag1 = PETSC_FALSE;
469: flag2 = PETSC_FALSE;
470: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_compute_eigenvalues",&flag1,PETSC_NULL);
471: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_plot_eigenvalues",&flag2,PETSC_NULL);
472: if (flag1 || flag2) {
473: PetscInt nits,n,i,neig;
474: PetscReal *r,*c;
475:
476: KSPGetIterationNumber(ksp,&nits);
477: n = nits+2;
479: if (!nits) {
480: PetscPrintf(((PetscObject)ksp)->comm,"Zero iterations in solver, cannot approximate any eigenvalues\n");
481: } else {
482: PetscMalloc(2*n*sizeof(PetscReal),&r);
483: c = r + n;
484: KSPComputeEigenvalues(ksp,n,r,c,&neig);
485: if (flag1) {
486: PetscPrintf(((PetscObject)ksp)->comm,"Iteratively computed eigenvalues\n");
487: for (i=0; i<neig; i++) {
488: if (c[i] >= 0.0) {PetscPrintf(((PetscObject)ksp)->comm,"%G + %Gi\n",r[i],c[i]);}
489: else {PetscPrintf(((PetscObject)ksp)->comm,"%G - %Gi\n",r[i],-c[i]);}
490: }
491: }
492: if (flag2 && !rank) {
493: PetscDraw draw;
494: PetscDrawSP drawsp;
496: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Iteratively Computed Eigenvalues",PETSC_DECIDE,PETSC_DECIDE,300,300,&viewer);
497: PetscViewerDrawGetDraw(viewer,0,&draw);
498: PetscDrawSPCreate(draw,1,&drawsp);
499: for (i=0; i<neig; i++) {
500: PetscDrawSPAddPoint(drawsp,r+i,c+i);
501: }
502: PetscDrawSPDraw(drawsp);
503: PetscDrawSPDestroy(&drawsp);
504: PetscViewerDestroy(&viewer);
505: }
506: PetscFree(r);
507: }
508: }
510: flag1 = PETSC_FALSE;
511: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_compute_singularvalues",&flag1,PETSC_NULL);
512: if (flag1) {
513: PetscInt nits;
515: KSPGetIterationNumber(ksp,&nits);
516: if (!nits) {
517: PetscPrintf(((PetscObject)ksp)->comm,"Zero iterations in solver, cannot approximate any singular values\n");
518: } else {
519: PetscReal emax,emin;
521: KSPComputeExtremeSingularValues(ksp,&emax,&emin);
522: PetscPrintf(((PetscObject)ksp)->comm,"Iteratively computed extreme singular values: max %G min %G max/min %G\n",emax,emin,emax/emin);
523: }
524: }
527: flag1 = PETSC_FALSE;
528: flag2 = PETSC_FALSE;
529: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_compute_eigenvalues_explicitly",&flag1,PETSC_NULL);
530: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_plot_eigenvalues_explicitly",&flag2,PETSC_NULL);
531: if (flag1 || flag2) {
532: PetscInt n,i;
533: PetscReal *r,*c;
534: VecGetSize(ksp->vec_sol,&n);
535: PetscMalloc2(n,PetscReal,&r,n,PetscReal,&c);
536: KSPComputeEigenvaluesExplicitly(ksp,n,r,c);
537: if (flag1) {
538: PetscPrintf(((PetscObject)ksp)->comm,"Explicitly computed eigenvalues\n");
539: for (i=0; i<n; i++) {
540: if (c[i] >= 0.0) {PetscPrintf(((PetscObject)ksp)->comm,"%G + %Gi\n",r[i],c[i]);}
541: else {PetscPrintf(((PetscObject)ksp)->comm,"%G - %Gi\n",r[i],-c[i]);}
542: }
543: }
544: if (flag2 && !rank) {
545: PetscDraw draw;
546: PetscDrawSP drawsp;
548: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Explicitly Computed Eigenvalues",0,320,300,300,&viewer);
549: PetscViewerDrawGetDraw(viewer,0,&draw);
550: PetscDrawSPCreate(draw,1,&drawsp);
551: for (i=0; i<n; i++) {
552: PetscDrawSPAddPoint(drawsp,r+i,c+i);
553: }
554: PetscDrawSPDraw(drawsp);
555: PetscDrawSPDestroy(&drawsp);
556: PetscViewerDestroy(&viewer);
557: }
558: PetscFree2(r,c);
559: }
561: flag2 = PETSC_FALSE;
562: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_view_operator",&flag2,PETSC_NULL);
563: if (flag2) {
564: Mat A,B;
565: PetscViewer viewer;
567: PCGetOperators(ksp->pc,&A,PETSC_NULL,PETSC_NULL);
568: MatComputeExplicitOperator(A,&B);
569: PetscViewerASCIIGetStdout(((PetscObject)ksp)->comm,&viewer);
570: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
571: MatView(B,viewer);
572: PetscViewerPopFormat(viewer);
573: MatDestroy(&B);
574: }
575: flag2 = PETSC_FALSE;
576: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_view_operator_binary",&flag2,PETSC_NULL);
577: if (flag2) {
578: Mat A,B;
579: PCGetOperators(ksp->pc,&A,PETSC_NULL,PETSC_NULL);
580: MatComputeExplicitOperator(A,&B);
581: MatView(B,PETSC_VIEWER_BINARY_(((PetscObject)ksp)->comm));
582: MatDestroy(&B);
583: }
584: flag2 = PETSC_FALSE;
585: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_view_preconditioned_operator_binary",&flag2,PETSC_NULL);
586: if (flag2) {
587: Mat B;
588: KSPComputeExplicitOperator(ksp,&B);
589: MatView(B,PETSC_VIEWER_BINARY_(((PetscObject)ksp)->comm));
590: MatDestroy(&B);
591: }
592: flag2 = PETSC_FALSE;
593: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_view_preconditioner_binary",&flag2,PETSC_NULL);
594: if (flag2) {
595: Mat B;
596: PCComputeExplicitOperator(ksp->pc,&B);
597: MatView(B,PETSC_VIEWER_BINARY_(((PetscObject)ksp)->comm));
598: MatDestroy(&B);
599: }
600: PetscOptionsGetString(((PetscObject)ksp)->prefix,"-ksp_view",filename,PETSC_MAX_PATH_LEN,&flg);
601: if (flg && !PetscPreLoadingOn) {
602: PetscViewerASCIIOpen(((PetscObject)ksp)->comm,filename,&viewer);
603: KSPView(ksp,viewer);
604: PetscViewerDestroy(&viewer);
605: }
606: flg = PETSC_FALSE;
607: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_final_residual",&flg,PETSC_NULL);
608: if (flg) {
609: Mat A;
610: Vec t;
611: PetscReal norm;
612: if (ksp->dscale && !ksp->dscalefix) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot compute final scale with -ksp_diagonal_scale except also with -ksp_diagonal_scale_fix");
613: PCGetOperators(ksp->pc,&A,0,0);
614: VecDuplicate(ksp->vec_sol,&t);
615: KSP_MatMult(ksp,A,ksp->vec_sol,t);
616: VecAYPX(t, -1.0, ksp->vec_rhs);
617: VecNorm(t,NORM_2,&norm);
618: VecDestroy(&t);
619: PetscPrintf(((PetscObject)ksp)->comm,"KSP final norm of residual %G\n",norm);
620: }
621: if (inXisinB) {
622: VecCopy(x,b);
623: VecDestroy(&x);
624: }
625: if (ksp->errorifnotconverged && ksp->reason < 0) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged");
626: return(0);
627: }
631: /*@
632: KSPSolveTranspose - Solves the transpose of a linear system.
634: Collective on KSP
636: Input Parameter:
637: + ksp - iterative context obtained from KSPCreate()
638: . b - right hand side vector
639: - x - solution vector
641: Notes: For complex numbers this solve the non-Hermitian transpose system.
643: Developer Notes: We need to implement a KSPSolveHermitianTranspose()
645: Level: developer
647: .keywords: KSP, solve, linear system
649: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPDefaultConverged(),
650: KSPSolve()
651: @*/
653: PetscErrorCode KSPSolveTranspose(KSP ksp,Vec b,Vec x)
654: {
656: PetscBool inXisinB=PETSC_FALSE;
662: if (x == b) {
663: VecDuplicate(b,&x);
664: inXisinB = PETSC_TRUE;
665: }
666: PetscObjectReference((PetscObject)b);
667: PetscObjectReference((PetscObject)x);
668: VecDestroy(&ksp->vec_rhs);
669: VecDestroy(&ksp->vec_sol);
670: ksp->vec_rhs = b;
671: ksp->vec_sol = x;
672: ksp->transpose_solve = PETSC_TRUE;
673: KSPSetUp(ksp);
674: if (ksp->guess_zero) { VecSet(ksp->vec_sol,0.0);}
675: (*ksp->ops->solve)(ksp);
676: if (!ksp->reason) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
677: if (ksp->printreason) {
678: if (ksp->reason > 0) {
679: PetscPrintf(((PetscObject)ksp)->comm,"Linear solve converged due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
680: } else {
681: PetscPrintf(((PetscObject)ksp)->comm,"Linear solve did not converge due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
682: }
683: }
684: if (inXisinB) {
685: VecCopy(x,b);
686: VecDestroy(&x);
687: }
688: if (ksp->errorifnotconverged && ksp->reason < 0) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged");
689: return(0);
690: }
694: /*@
695: KSPReset - Resets a KSP context to the kspsetupcalled = 0 state and removes any allocated Vecs and Mats
697: Collective on KSP
699: Input Parameter:
700: . ksp - iterative context obtained from KSPCreate()
702: Level: beginner
704: .keywords: KSP, destroy
706: .seealso: KSPCreate(), KSPSetUp(), KSPSolve()
707: @*/
708: PetscErrorCode KSPReset(KSP ksp)
709: {
714: if (!ksp) return(0);
715: if (ksp->ops->reset) {
716: (*ksp->ops->reset)(ksp);
717: }
718: if (ksp->pc) {PCReset(ksp->pc);}
719: KSPFischerGuessDestroy(&ksp->guess);
720: VecDestroyVecs(ksp->nwork,&ksp->work);
721: VecDestroy(&ksp->vec_rhs);
722: VecDestroy(&ksp->vec_sol);
723: VecDestroy(&ksp->diagonal);
724: VecDestroy(&ksp->truediagonal);
725: MatNullSpaceDestroy(&ksp->nullsp);
726: ksp->setupstage = KSP_SETUP_NEW;
727: return(0);
728: }
732: /*@
733: KSPDestroy - Destroys KSP context.
735: Collective on KSP
737: Input Parameter:
738: . ksp - iterative context obtained from KSPCreate()
740: Level: beginner
742: .keywords: KSP, destroy
744: .seealso: KSPCreate(), KSPSetUp(), KSPSolve()
745: @*/
746: PetscErrorCode KSPDestroy(KSP *ksp)
747: {
751: if (!*ksp) return(0);
753: if (--((PetscObject)(*ksp))->refct > 0) {*ksp = 0; return(0);}
755: KSPReset((*ksp));
757: PetscObjectDepublish((*ksp));
758: if ((*ksp)->ops->destroy) {(*(*ksp)->ops->destroy)(*ksp);}
760: DMDestroy(&(*ksp)->dm);
761: PCDestroy(&(*ksp)->pc);
762: PetscFree((*ksp)->res_hist_alloc);
763: if ((*ksp)->convergeddestroy) {
764: (*(*ksp)->convergeddestroy)((*ksp)->cnvP);
765: }
766: KSPMonitorCancel((*ksp));
767: PetscHeaderDestroy(ksp);
768: return(0);
769: }
773: /*@
774: KSPSetPCSide - Sets the preconditioning side.
776: Logically Collective on KSP
778: Input Parameter:
779: . ksp - iterative context obtained from KSPCreate()
781: Output Parameter:
782: . side - the preconditioning side, where side is one of
783: .vb
784: PC_LEFT - left preconditioning (default)
785: PC_RIGHT - right preconditioning
786: PC_SYMMETRIC - symmetric preconditioning
787: .ve
789: Options Database Keys:
790: . -ksp_pc_side <right,left,symmetric>
792: Notes:
793: Left preconditioning is used by default for most Krylov methods except KSPFGMRES which only supports right preconditioning.
794: Symmetric preconditioning is currently available only for the KSPQCG method. Note, however, that
795: symmetric preconditioning can be emulated by using either right or left
796: preconditioning and a pre or post processing step.
798: Level: intermediate
800: .keywords: KSP, set, right, left, symmetric, side, preconditioner, flag
802: .seealso: KSPGetPCSide()
803: @*/
804: PetscErrorCode KSPSetPCSide(KSP ksp,PCSide side)
805: {
809: ksp->pc_side = side;
810: return(0);
811: }
815: /*@
816: KSPGetPCSide - Gets the preconditioning side.
818: Not Collective
820: Input Parameter:
821: . ksp - iterative context obtained from KSPCreate()
823: Output Parameter:
824: . side - the preconditioning side, where side is one of
825: .vb
826: PC_LEFT - left preconditioning (default)
827: PC_RIGHT - right preconditioning
828: PC_SYMMETRIC - symmetric preconditioning
829: .ve
831: Level: intermediate
833: .keywords: KSP, get, right, left, symmetric, side, preconditioner, flag
835: .seealso: KSPSetPCSide()
836: @*/
837: PetscErrorCode KSPGetPCSide(KSP ksp,PCSide *side)
838: {
844: KSPSetUpNorms_Private(ksp);
845: *side = ksp->pc_side;
846: return(0);
847: }
851: /*@
852: KSPGetTolerances - Gets the relative, absolute, divergence, and maximum
853: iteration tolerances used by the default KSP convergence tests.
855: Not Collective
857: Input Parameter:
858: . ksp - the Krylov subspace context
859:
860: Output Parameters:
861: + rtol - the relative convergence tolerance
862: . abstol - the absolute convergence tolerance
863: . dtol - the divergence tolerance
864: - maxits - maximum number of iterations
866: Notes:
867: The user can specify PETSC_NULL for any parameter that is not needed.
869: Level: intermediate
871: .keywords: KSP, get, tolerance, absolute, relative, divergence, convergence,
872: maximum, iterations
874: .seealso: KSPSetTolerances()
875: @*/
876: PetscErrorCode KSPGetTolerances(KSP ksp,PetscReal *rtol,PetscReal *abstol,PetscReal *dtol,PetscInt *maxits)
877: {
880: if (abstol) *abstol = ksp->abstol;
881: if (rtol) *rtol = ksp->rtol;
882: if (dtol) *dtol = ksp->divtol;
883: if (maxits) *maxits = ksp->max_it;
884: return(0);
885: }
889: /*@
890: KSPSetTolerances - Sets the relative, absolute, divergence, and maximum
891: iteration tolerances used by the default KSP convergence testers.
893: Logically Collective on KSP
895: Input Parameters:
896: + ksp - the Krylov subspace context
897: . rtol - the relative convergence tolerance
898: (relative decrease in the residual norm)
899: . abstol - the absolute convergence tolerance
900: (absolute size of the residual norm)
901: . dtol - the divergence tolerance
902: (amount residual can increase before KSPDefaultConverged()
903: concludes that the method is diverging)
904: - maxits - maximum number of iterations to use
906: Options Database Keys:
907: + -ksp_atol <abstol> - Sets abstol
908: . -ksp_rtol <rtol> - Sets rtol
909: . -ksp_divtol <dtol> - Sets dtol
910: - -ksp_max_it <maxits> - Sets maxits
912: Notes:
913: Use PETSC_DEFAULT to retain the default value of any of the tolerances.
915: See KSPDefaultConverged() for details on the use of these parameters
916: in the default convergence test. See also KSPSetConvergenceTest()
917: for setting user-defined stopping criteria.
919: Level: intermediate
921: .keywords: KSP, set, tolerance, absolute, relative, divergence,
922: convergence, maximum, iterations
924: .seealso: KSPGetTolerances(), KSPDefaultConverged(), KSPSetConvergenceTest()
925: @*/
926: PetscErrorCode KSPSetTolerances(KSP ksp,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt maxits)
927: {
935: if (rtol != PETSC_DEFAULT) {
936: if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %G must be non-negative and less than 1.0",rtol); ksp->rtol = rtol;
937: }
938: if (abstol != PETSC_DEFAULT) {
939: if (abstol < 0.0) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol);
940: ksp->abstol = abstol;
941: }
942: if (dtol != PETSC_DEFAULT) {
943: if (dtol < 0.0) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Divergence tolerance %G must be larger than 1.0",dtol);
944: ksp->divtol = dtol;
945: }
946: if (maxits != PETSC_DEFAULT) {
947: if (maxits < 0) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxits);
948: ksp->max_it = maxits;
949: }
950: return(0);
951: }
955: /*@
956: KSPSetInitialGuessNonzero - Tells the iterative solver that the
957: initial guess is nonzero; otherwise KSP assumes the initial guess
958: is to be zero (and thus zeros it out before solving).
960: Logically Collective on KSP
962: Input Parameters:
963: + ksp - iterative context obtained from KSPCreate()
964: - flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero
966: Options database keys:
967: . -ksp_initial_guess_nonzero : use nonzero initial guess; this takes an optional truth value (0/1/no/yes/true/false)
969: Level: beginner
971: Notes:
972: If this is not called the X vector is zeroed in the call to KSPSolve().
974: .keywords: KSP, set, initial guess, nonzero
976: .seealso: KSPGetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll()
977: @*/
978: PetscErrorCode KSPSetInitialGuessNonzero(KSP ksp,PetscBool flg)
979: {
983: ksp->guess_zero = (PetscBool)!(int)flg;
984: return(0);
985: }
989: /*@
990: KSPGetInitialGuessNonzero - Determines whether the KSP solver is using
991: a zero initial guess.
993: Not Collective
995: Input Parameter:
996: . ksp - iterative context obtained from KSPCreate()
998: Output Parameter:
999: . flag - PETSC_TRUE if guess is nonzero, else PETSC_FALSE
1001: Level: intermediate
1003: .keywords: KSP, set, initial guess, nonzero
1005: .seealso: KSPSetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll()
1006: @*/
1007: PetscErrorCode KSPGetInitialGuessNonzero(KSP ksp,PetscBool *flag)
1008: {
1012: if (ksp->guess_zero) *flag = PETSC_FALSE;
1013: else *flag = PETSC_TRUE;
1014: return(0);
1015: }
1019: /*@
1020: KSPSetErrorIfNotConverged - Causes KSPSolve() to generate an error if the solver has not converged.
1022: Logically Collective on KSP
1024: Input Parameters:
1025: + ksp - iterative context obtained from KSPCreate()
1026: - flg - PETSC_TRUE indicates you want the error generated
1028: Options database keys:
1029: . -ksp_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false)
1031: Level: intermediate
1033: Notes:
1034: Normally PETSc continues if a linear solver fails to converge, you can call KSPGetConvergedReason() after a KSPSolve()
1035: to determine if it has converged.
1037: .keywords: KSP, set, initial guess, nonzero
1039: .seealso: KSPGetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll(), KSPGetErrorIfNotConverged()
1040: @*/
1041: PetscErrorCode KSPSetErrorIfNotConverged(KSP ksp,PetscBool flg)
1042: {
1046: ksp->errorifnotconverged = flg;
1047: return(0);
1048: }
1052: /*@
1053: KSPGetErrorIfNotConverged - Will KSPSolve() generate an error if the solver does not converge?
1055: Not Collective
1057: Input Parameter:
1058: . ksp - iterative context obtained from KSPCreate()
1060: Output Parameter:
1061: . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE
1063: Level: intermediate
1065: .keywords: KSP, set, initial guess, nonzero
1067: .seealso: KSPSetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll(), KSPSetErrorIfNotConverged()
1068: @*/
1069: PetscErrorCode KSPGetErrorIfNotConverged(KSP ksp,PetscBool *flag)
1070: {
1074: *flag = ksp->errorifnotconverged;
1075: return(0);
1076: }
1080: /*@
1081: KSPSetInitialGuessKnoll - Tells the iterative solver to use PCApply(pc,b,..) to compute the initial guess (The Knoll trick)
1083: Logically Collective on KSP
1085: Input Parameters:
1086: + ksp - iterative context obtained from KSPCreate()
1087: - flg - PETSC_TRUE or PETSC_FALSE
1089: Level: advanced
1092: .keywords: KSP, set, initial guess, nonzero
1094: .seealso: KSPGetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero()
1095: @*/
1096: PetscErrorCode KSPSetInitialGuessKnoll(KSP ksp,PetscBool flg)
1097: {
1101: ksp->guess_knoll = flg;
1102: return(0);
1103: }
1107: /*@
1108: KSPGetInitialGuessKnoll - Determines whether the KSP solver is using the Knoll trick (using PCApply(pc,b,...) to compute
1109: the initial guess
1111: Not Collective
1113: Input Parameter:
1114: . ksp - iterative context obtained from KSPCreate()
1116: Output Parameter:
1117: . flag - PETSC_TRUE if using Knoll trick, else PETSC_FALSE
1119: Level: advanced
1121: .keywords: KSP, set, initial guess, nonzero
1123: .seealso: KSPSetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero()
1124: @*/
1125: PetscErrorCode KSPGetInitialGuessKnoll(KSP ksp,PetscBool *flag)
1126: {
1130: *flag = ksp->guess_knoll;
1131: return(0);
1132: }
1136: /*@
1137: KSPGetComputeSingularValues - Gets the flag indicating whether the extreme singular
1138: values will be calculated via a Lanczos or Arnoldi process as the linear
1139: system is solved.
1141: Not Collective
1143: Input Parameter:
1144: . ksp - iterative context obtained from KSPCreate()
1146: Output Parameter:
1147: . flg - PETSC_TRUE or PETSC_FALSE
1149: Options Database Key:
1150: . -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()
1152: Notes:
1153: Currently this option is not valid for all iterative methods.
1155: Many users may just want to use the monitoring routine
1156: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1157: to print the singular values at each iteration of the linear solve.
1159: Level: advanced
1161: .keywords: KSP, set, compute, singular values
1163: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue()
1164: @*/
1165: PetscErrorCode KSPGetComputeSingularValues(KSP ksp,PetscBool *flg)
1166: {
1170: *flg = ksp->calc_sings;
1171: return(0);
1172: }
1176: /*@
1177: KSPSetComputeSingularValues - Sets a flag so that the extreme singular
1178: values will be calculated via a Lanczos or Arnoldi process as the linear
1179: system is solved.
1181: Logically Collective on KSP
1183: Input Parameters:
1184: + ksp - iterative context obtained from KSPCreate()
1185: - flg - PETSC_TRUE or PETSC_FALSE
1187: Options Database Key:
1188: . -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()
1190: Notes:
1191: Currently this option is not valid for all iterative methods.
1193: Many users may just want to use the monitoring routine
1194: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1195: to print the singular values at each iteration of the linear solve.
1197: Level: advanced
1199: .keywords: KSP, set, compute, singular values
1201: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue()
1202: @*/
1203: PetscErrorCode KSPSetComputeSingularValues(KSP ksp,PetscBool flg)
1204: {
1208: ksp->calc_sings = flg;
1209: return(0);
1210: }
1214: /*@
1215: KSPGetComputeEigenvalues - Gets the flag indicating that the extreme eigenvalues
1216: values will be calculated via a Lanczos or Arnoldi process as the linear
1217: system is solved.
1219: Not Collective
1221: Input Parameter:
1222: . ksp - iterative context obtained from KSPCreate()
1224: Output Parameter:
1225: . flg - PETSC_TRUE or PETSC_FALSE
1227: Notes:
1228: Currently this option is not valid for all iterative methods.
1230: Level: advanced
1232: .keywords: KSP, set, compute, eigenvalues
1234: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly()
1235: @*/
1236: PetscErrorCode KSPGetComputeEigenvalues(KSP ksp,PetscBool *flg)
1237: {
1241: *flg = ksp->calc_sings;
1242: return(0);
1243: }
1247: /*@
1248: KSPSetComputeEigenvalues - Sets a flag so that the extreme eigenvalues
1249: values will be calculated via a Lanczos or Arnoldi process as the linear
1250: system is solved.
1252: Logically Collective on KSP
1254: Input Parameters:
1255: + ksp - iterative context obtained from KSPCreate()
1256: - flg - PETSC_TRUE or PETSC_FALSE
1258: Notes:
1259: Currently this option is not valid for all iterative methods.
1261: Level: advanced
1263: .keywords: KSP, set, compute, eigenvalues
1265: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly()
1266: @*/
1267: PetscErrorCode KSPSetComputeEigenvalues(KSP ksp,PetscBool flg)
1268: {
1272: ksp->calc_sings = flg;
1273: return(0);
1274: }
1278: /*@
1279: KSPGetRhs - Gets the right-hand-side vector for the linear system to
1280: be solved.
1282: Not Collective
1284: Input Parameter:
1285: . ksp - iterative context obtained from KSPCreate()
1287: Output Parameter:
1288: . r - right-hand-side vector
1290: Level: developer
1292: .keywords: KSP, get, right-hand-side, rhs
1294: .seealso: KSPGetSolution(), KSPSolve()
1295: @*/
1296: PetscErrorCode KSPGetRhs(KSP ksp,Vec *r)
1297: {
1301: *r = ksp->vec_rhs;
1302: return(0);
1303: }
1307: /*@
1308: KSPGetSolution - Gets the location of the solution for the
1309: linear system to be solved. Note that this may not be where the solution
1310: is stored during the iterative process; see KSPBuildSolution().
1312: Not Collective
1314: Input Parameters:
1315: . ksp - iterative context obtained from KSPCreate()
1317: Output Parameters:
1318: . v - solution vector
1320: Level: developer
1322: .keywords: KSP, get, solution
1324: .seealso: KSPGetRhs(), KSPBuildSolution(), KSPSolve()
1325: @*/
1326: PetscErrorCode KSPGetSolution(KSP ksp,Vec *v)
1327: {
1331: *v = ksp->vec_sol;
1332: return(0);
1333: }
1337: /*@
1338: KSPSetPC - Sets the preconditioner to be used to calculate the
1339: application of the preconditioner on a vector.
1341: Collective on KSP
1343: Input Parameters:
1344: + ksp - iterative context obtained from KSPCreate()
1345: - pc - the preconditioner object
1347: Notes:
1348: Use KSPGetPC() to retrieve the preconditioner context (for example,
1349: to free it at the end of the computations).
1351: Level: developer
1353: .keywords: KSP, set, precondition, Binv
1355: .seealso: KSPGetPC()
1356: @*/
1357: PetscErrorCode KSPSetPC(KSP ksp,PC pc)
1358: {
1365: PetscObjectReference((PetscObject)pc);
1366: PCDestroy(&ksp->pc);
1367: ksp->pc = pc;
1368: PetscLogObjectParent(ksp,ksp->pc);
1369: return(0);
1370: }
1374: /*@
1375: KSPGetPC - Returns a pointer to the preconditioner context
1376: set with KSPSetPC().
1378: Not Collective
1380: Input Parameters:
1381: . ksp - iterative context obtained from KSPCreate()
1383: Output Parameter:
1384: . pc - preconditioner context
1386: Level: developer
1388: .keywords: KSP, get, preconditioner, Binv
1390: .seealso: KSPSetPC()
1391: @*/
1392: PetscErrorCode KSPGetPC(KSP ksp,PC *pc)
1393: {
1399: if (!ksp->pc) {
1400: PCCreate(((PetscObject)ksp)->comm,&ksp->pc);
1401: PetscObjectIncrementTabLevel((PetscObject)ksp->pc,(PetscObject)ksp,0);
1402: PetscLogObjectParent(ksp,ksp->pc);
1403: }
1404: *pc = ksp->pc;
1405: return(0);
1406: }
1410: /*@
1411: KSPMonitor - runs the user provided monitor routines, if they exist
1413: Collective on KSP
1415: Input Parameters:
1416: + ksp - iterative context obtained from KSPCreate()
1417: . it - iteration number
1418: - rnorm - relative norm of the residual
1420: Notes:
1421: This routine is called by the KSP implementations.
1422: It does not typically need to be called by the user.
1424: Level: developer
1426: .seealso: KSPMonitorSet()
1427: @*/
1428: PetscErrorCode KSPMonitor(KSP ksp,PetscInt it,PetscReal rnorm)
1429: {
1430: PetscInt i, n = ksp->numbermonitors;
1434: for (i=0; i<n; i++) {
1435: (*ksp->monitor[i])(ksp,it,rnorm,ksp->monitorcontext[i]);
1436: }
1437: return(0);
1438: }
1442: /*@C
1443: KSPMonitorSet - Sets an ADDITIONAL function to be called at every iteration to monitor
1444: the residual/error etc.
1445:
1446: Logically Collective on KSP
1448: Input Parameters:
1449: + ksp - iterative context obtained from KSPCreate()
1450: . monitor - pointer to function (if this is PETSC_NULL, it turns off monitoring
1451: . mctx - [optional] context for private data for the
1452: monitor routine (use PETSC_NULL if no context is desired)
1453: - monitordestroy - [optional] routine that frees monitor context
1454: (may be PETSC_NULL)
1456: Calling Sequence of monitor:
1457: $ monitor (KSP ksp, int it, PetscReal rnorm, void *mctx)
1459: + ksp - iterative context obtained from KSPCreate()
1460: . it - iteration number
1461: . rnorm - (estimated) 2-norm of (preconditioned) residual
1462: - mctx - optional monitoring context, as set by KSPMonitorSet()
1464: Options Database Keys:
1465: + -ksp_monitor - sets KSPMonitorDefault()
1466: . -ksp_monitor_true_residual - sets KSPMonitorTrueResidualNorm()
1467: . -ksp_monitor_draw - sets line graph monitor,
1468: uses KSPMonitorLGCreate()
1469: . -ksp_monitor_draw_true_residual - sets line graph monitor,
1470: uses KSPMonitorLGCreate()
1471: . -ksp_monitor_singular_value - sets KSPMonitorSingularValue()
1472: - -ksp_monitor_cancel - cancels all monitors that have
1473: been hardwired into a code by
1474: calls to KSPMonitorSet(), but
1475: does not cancel those set via
1476: the options database.
1478: Notes:
1479: The default is to do nothing. To print the residual, or preconditioned
1480: residual if KSPSetNormType(ksp,KSP_NORM_PRECONDITIONED) was called, use
1481: KSPMonitorDefault() as the monitoring routine, with a null monitoring
1482: context.
1484: Several different monitoring routines may be set by calling
1485: KSPMonitorSet() multiple times; all will be called in the
1486: order in which they were set.
1488: Fortran notes: Only a single monitor function can be set for each KSP object
1490: Level: beginner
1492: .keywords: KSP, set, monitor
1494: .seealso: KSPMonitorDefault(), KSPMonitorLGCreate(), KSPMonitorCancel()
1495: @*/
1496: PetscErrorCode KSPMonitorSet(KSP ksp,PetscErrorCode (*monitor)(KSP,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
1497: {
1498: PetscInt i;
1503: if (ksp->numbermonitors >= MAXKSPMONITORS) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many KSP monitors set");
1504: for (i=0; i<ksp->numbermonitors;i++) {
1505: if (monitor == ksp->monitor[i] && monitordestroy == ksp->monitordestroy[i] && mctx == ksp->monitorcontext[i]) {
1506: if (monitordestroy) {
1507: (*monitordestroy)(&mctx);
1508: }
1509: return(0);
1510: }
1511: }
1512: ksp->monitor[ksp->numbermonitors] = monitor;
1513: ksp->monitordestroy[ksp->numbermonitors] = monitordestroy;
1514: ksp->monitorcontext[ksp->numbermonitors++] = (void*)mctx;
1515: return(0);
1516: }
1520: /*@
1521: KSPMonitorCancel - Clears all monitors for a KSP object.
1523: Logically Collective on KSP
1525: Input Parameters:
1526: . ksp - iterative context obtained from KSPCreate()
1528: Options Database Key:
1529: . -ksp_monitor_cancel - Cancels all monitors that have
1530: been hardwired into a code by calls to KSPMonitorSet(),
1531: but does not cancel those set via the options database.
1533: Level: intermediate
1535: .keywords: KSP, set, monitor
1537: .seealso: KSPMonitorDefault(), KSPMonitorLGCreate(), KSPMonitorSet()
1538: @*/
1539: PetscErrorCode KSPMonitorCancel(KSP ksp)
1540: {
1542: PetscInt i;
1546: for (i=0; i<ksp->numbermonitors; i++) {
1547: if (ksp->monitordestroy[i]) {
1548: (*ksp->monitordestroy[i])(&ksp->monitorcontext[i]);
1549: }
1550: }
1551: ksp->numbermonitors = 0;
1552: return(0);
1553: }
1557: /*@C
1558: KSPGetMonitorContext - Gets the monitoring context, as set by
1559: KSPMonitorSet() for the FIRST monitor only.
1561: Not Collective
1563: Input Parameter:
1564: . ksp - iterative context obtained from KSPCreate()
1566: Output Parameter:
1567: . ctx - monitoring context
1569: Level: intermediate
1571: .keywords: KSP, get, monitor, context
1573: .seealso: KSPMonitorDefault(), KSPMonitorLGCreate()
1574: @*/
1575: PetscErrorCode KSPGetMonitorContext(KSP ksp,void **ctx)
1576: {
1579: *ctx = (ksp->monitorcontext[0]);
1580: return(0);
1581: }
1585: /*@
1586: KSPSetResidualHistory - Sets the array used to hold the residual history.
1587: If set, this array will contain the residual norms computed at each
1588: iteration of the solver.
1590: Not Collective
1592: Input Parameters:
1593: + ksp - iterative context obtained from KSPCreate()
1594: . a - array to hold history
1595: . na - size of a
1596: - reset - PETSC_TRUE indicates the history counter is reset to zero
1597: for each new linear solve
1599: Level: advanced
1601: Notes: The array is NOT freed by PETSc so the user needs to keep track of
1602: it and destroy once the KSP object is destroyed.
1604: If 'a' is PETSC_NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
1605: default array of length 10000 is allocated.
1607: .keywords: KSP, set, residual, history, norm
1609: .seealso: KSPGetResidualHistory()
1611: @*/
1612: PetscErrorCode KSPSetResidualHistory(KSP ksp,PetscReal a[],PetscInt na,PetscBool reset)
1613: {
1619: PetscFree(ksp->res_hist_alloc);
1620: if (na != PETSC_DECIDE && na != PETSC_DEFAULT && a) {
1621: ksp->res_hist = a;
1622: ksp->res_hist_max = na;
1623: } else {
1624: if (na != PETSC_DECIDE && na != PETSC_DEFAULT)
1625: ksp->res_hist_max = na;
1626: else
1627: ksp->res_hist_max = 10000; /* like default ksp->max_it */
1628: PetscMalloc(ksp->res_hist_max*sizeof(PetscReal),&ksp->res_hist_alloc);
1629: ksp->res_hist = ksp->res_hist_alloc;
1630: }
1631: ksp->res_hist_len = 0;
1632: ksp->res_hist_reset = reset;
1634: return(0);
1635: }
1639: /*@C
1640: KSPGetResidualHistory - Gets the array used to hold the residual history
1641: and the number of residuals it contains.
1643: Not Collective
1645: Input Parameter:
1646: . ksp - iterative context obtained from KSPCreate()
1648: Output Parameters:
1649: + a - pointer to array to hold history (or PETSC_NULL)
1650: - na - number of used entries in a (or PETSC_NULL)
1652: Level: advanced
1654: Notes:
1655: Can only be called after a KSPSetResidualHistory() otherwise a and na are set to zero
1657: The Fortran version of this routine has a calling sequence
1658: $ call KSPGetResidualHistory(KSP ksp, integer na, integer ierr)
1659: note that you have passed a Fortran array into KSPSetResidualHistory() and you need
1660: to access the residual values from this Fortran array you provided. Only the na (number of
1661: residual norms currently held) is set.
1663: .keywords: KSP, get, residual, history, norm
1665: .seealso: KSPGetResidualHistory()
1667: @*/
1668: PetscErrorCode KSPGetResidualHistory(KSP ksp,PetscReal *a[],PetscInt *na)
1669: {
1672: if (a) *a = ksp->res_hist;
1673: if (na) *na = ksp->res_hist_len;
1674: return(0);
1675: }
1679: /*@C
1680: KSPSetConvergenceTest - Sets the function to be used to determine
1681: convergence.
1683: Logically Collective on KSP
1685: Input Parameters:
1686: + ksp - iterative context obtained from KSPCreate()
1687: . converge - pointer to int function
1688: . cctx - context for private data for the convergence routine (may be null)
1689: - destroy - a routine for destroying the context (may be null)
1691: Calling sequence of converge:
1692: $ converge (KSP ksp, int it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)
1694: + ksp - iterative context obtained from KSPCreate()
1695: . it - iteration number
1696: . rnorm - (estimated) 2-norm of (preconditioned) residual
1697: . reason - the reason why it has converged or diverged
1698: - cctx - optional convergence context, as set by KSPSetConvergenceTest()
1701: Notes:
1702: Must be called after the KSP type has been set so put this after
1703: a call to KSPSetType(), or KSPSetFromOptions().
1705: The default convergence test, KSPDefaultConverged(), aborts if the
1706: residual grows to more than 10000 times the initial residual.
1708: The default is a combination of relative and absolute tolerances.
1709: The residual value that is tested may be an approximation; routines
1710: that need exact values should compute them.
1712: In the default PETSc convergence test, the precise values of reason
1713: are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
1715: Level: advanced
1717: .keywords: KSP, set, convergence, test, context
1719: .seealso: KSPDefaultConverged(), KSPGetConvergenceContext(), KSPSetTolerances()
1720: @*/
1721: PetscErrorCode KSPSetConvergenceTest(KSP ksp,PetscErrorCode (*converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
1722: {
1727: if (ksp->convergeddestroy) {
1728: (*ksp->convergeddestroy)(ksp->cnvP);
1729: }
1730: ksp->converged = converge;
1731: ksp->convergeddestroy = destroy;
1732: ksp->cnvP = (void*)cctx;
1733: return(0);
1734: }
1738: /*@C
1739: KSPGetConvergenceContext - Gets the convergence context set with
1740: KSPSetConvergenceTest().
1742: Not Collective
1744: Input Parameter:
1745: . ksp - iterative context obtained from KSPCreate()
1747: Output Parameter:
1748: . ctx - monitoring context
1750: Level: advanced
1752: .keywords: KSP, get, convergence, test, context
1754: .seealso: KSPDefaultConverged(), KSPSetConvergenceTest()
1755: @*/
1756: PetscErrorCode KSPGetConvergenceContext(KSP ksp,void **ctx)
1757: {
1760: *ctx = ksp->cnvP;
1761: return(0);
1762: }
1766: /*@C
1767: KSPBuildSolution - Builds the approximate solution in a vector provided.
1768: This routine is NOT commonly needed (see KSPSolve()).
1770: Collective on KSP
1772: Input Parameter:
1773: . ctx - iterative context obtained from KSPCreate()
1775: Output Parameter:
1776: Provide exactly one of
1777: + v - location to stash solution.
1778: - V - the solution is returned in this location. This vector is created
1779: internally. This vector should NOT be destroyed by the user with
1780: VecDestroy().
1782: Notes:
1783: This routine can be used in one of two ways
1784: .vb
1785: KSPBuildSolution(ksp,PETSC_NULL,&V);
1786: or
1787: KSPBuildSolution(ksp,v,PETSC_NULL); or KSPBuildSolution(ksp,v,&v);
1788: .ve
1789: In the first case an internal vector is allocated to store the solution
1790: (the user cannot destroy this vector). In the second case the solution
1791: is generated in the vector that the user provides. Note that for certain
1792: methods, such as KSPCG, the second case requires a copy of the solution,
1793: while in the first case the call is essentially free since it simply
1794: returns the vector where the solution already is stored. For some methods
1795: like GMRES this is a reasonably expensive operation and should only be
1796: used in truly needed.
1798: Level: advanced
1800: .keywords: KSP, build, solution
1802: .seealso: KSPGetSolution(), KSPBuildResidual()
1803: @*/
1804: PetscErrorCode KSPBuildSolution(KSP ksp,Vec v,Vec *V)
1805: {
1810: if (!V && !v) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_WRONG,"Must provide either v or V");
1811: if (!V) V = &v;
1812: (*ksp->ops->buildsolution)(ksp,v,V);
1813: return(0);
1814: }
1818: /*@C
1819: KSPBuildResidual - Builds the residual in a vector provided.
1821: Collective on KSP
1823: Input Parameter:
1824: . ksp - iterative context obtained from KSPCreate()
1826: Output Parameters:
1827: + v - optional location to stash residual. If v is not provided,
1828: then a location is generated.
1829: . t - work vector. If not provided then one is generated.
1830: - V - the residual
1832: Notes:
1833: Regardless of whether or not v is provided, the residual is
1834: returned in V.
1836: Level: advanced
1838: .keywords: KSP, build, residual
1840: .seealso: KSPBuildSolution()
1841: @*/
1842: PetscErrorCode KSPBuildResidual(KSP ksp,Vec t,Vec v,Vec *V)
1843: {
1845: PetscBool flag = PETSC_FALSE;
1846: Vec w = v,tt = t;
1850: if (!w) {
1851: VecDuplicate(ksp->vec_rhs,&w);
1852: PetscLogObjectParent((PetscObject)ksp,w);
1853: }
1854: if (!tt) {
1855: VecDuplicate(ksp->vec_sol,&tt); flag = PETSC_TRUE;
1856: PetscLogObjectParent((PetscObject)ksp,tt);
1857: }
1858: (*ksp->ops->buildresidual)(ksp,tt,w,V);
1859: if (flag) {VecDestroy(&tt);}
1860: return(0);
1861: }
1865: /*@
1866: KSPSetDiagonalScale - Tells KSP to symmetrically diagonally scale the system
1867: before solving. This actually CHANGES the matrix (and right hand side).
1869: Logically Collective on KSP
1871: Input Parameter:
1872: + ksp - the KSP context
1873: - scale - PETSC_TRUE or PETSC_FALSE
1875: Options Database Key:
1876: + -ksp_diagonal_scale -
1877: - -ksp_diagonal_scale_fix - scale the matrix back AFTER the solve
1880: BE CAREFUL with this routine: it actually scales the matrix and right
1881: hand side that define the system. After the system is solved the matrix
1882: and right hand side remain scaled.
1884: This routine is only used if the matrix and preconditioner matrix are
1885: the same thing.
1887: This should NOT be used within the SNES solves if you are using a line
1888: search.
1889:
1890: If you use this with the PCType Eisenstat preconditioner than you can
1891: use the PCEisenstatNoDiagonalScaling() option, or -pc_eisenstat_no_diagonal_scaling
1892: to save some unneeded, redundant flops.
1894: Level: intermediate
1896: .keywords: KSP, set, options, prefix, database
1898: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScaleFix()
1899: @*/
1900: PetscErrorCode KSPSetDiagonalScale(KSP ksp,PetscBool scale)
1901: {
1905: ksp->dscale = scale;
1906: return(0);
1907: }
1911: /*@
1912: KSPGetDiagonalScale - Checks if KSP solver scales the matrix and
1913: right hand side
1915: Not Collective
1917: Input Parameter:
1918: . ksp - the KSP context
1920: Output Parameter:
1921: . scale - PETSC_TRUE or PETSC_FALSE
1923: Notes:
1924: BE CAREFUL with this routine: it actually scales the matrix and right
1925: hand side that define the system. After the system is solved the matrix
1926: and right hand side remain scaled.
1928: This routine is only used if the matrix and preconditioner matrix are
1929: the same thing.
1931: Level: intermediate
1933: .keywords: KSP, set, options, prefix, database
1935: .seealso: KSPSetDiagonalScale(), KSPSetDiagonalScaleFix()
1936: @*/
1937: PetscErrorCode KSPGetDiagonalScale(KSP ksp,PetscBool *scale)
1938: {
1942: *scale = ksp->dscale;
1943: return(0);
1944: }
1948: /*@
1949: KSPSetDiagonalScaleFix - Tells KSP to diagonally scale the system
1950: back after solving.
1952: Logically Collective on KSP
1954: Input Parameter:
1955: + ksp - the KSP context
1956: - fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
1957: rescale (default)
1959: Notes:
1960: Must be called after KSPSetDiagonalScale()
1962: Using this will slow things down, because it rescales the matrix before and
1963: after each linear solve. This is intended mainly for testing to allow one
1964: to easily get back the original system to make sure the solution computed is
1965: accurate enough.
1967: This routine is only used if the matrix and preconditioner matrix are
1968: the same thing.
1970: Level: intermediate
1972: .keywords: KSP, set, options, prefix, database
1974: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPGetDiagonalScaleFix()
1975: @*/
1976: PetscErrorCode KSPSetDiagonalScaleFix(KSP ksp,PetscBool fix)
1977: {
1981: ksp->dscalefix = fix;
1982: return(0);
1983: }
1987: /*@
1988: KSPGetDiagonalScaleFix - Determines if KSP diagonally scales the system
1989: back after solving.
1991: Not Collective
1993: Input Parameter:
1994: . ksp - the KSP context
1996: Output Parameter:
1997: . fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
1998: rescale (default)
2000: Notes:
2001: Must be called after KSPSetDiagonalScale()
2003: If PETSC_TRUE will slow things down, because it rescales the matrix before and
2004: after each linear solve. This is intended mainly for testing to allow one
2005: to easily get back the original system to make sure the solution computed is
2006: accurate enough.
2008: This routine is only used if the matrix and preconditioner matrix are
2009: the same thing.
2011: Level: intermediate
2013: .keywords: KSP, set, options, prefix, database
2015: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPSetDiagonalScaleFix()
2016: @*/
2017: PetscErrorCode KSPGetDiagonalScaleFix(KSP ksp,PetscBool *fix)
2018: {
2022: *fix = ksp->dscalefix;
2023: return(0);
2024: }