Actual source code: itfunc.c
petsc-3.9.4 2018-09-11
2: /*
3: Interface KSP routines that the user calls.
4: */
6: #include <petsc/private/kspimpl.h>
7: #include <petscdm.h>
9: /*@
10: KSPComputeExtremeSingularValues - Computes the extreme singular values
11: for the preconditioned operator. Called after or during KSPSolve().
13: Not Collective
15: Input Parameter:
16: . ksp - iterative context obtained from KSPCreate()
18: Output Parameters:
19: . emin, emax - extreme singular values
21: Options Database Keys:
22: . -ksp_compute_singularvalues - compute extreme singular values and print when KSPSolve completes.
24: Notes:
25: One must call KSPSetComputeSingularValues() before calling KSPSetUp()
26: (or use the option -ksp_compute_eigenvalues) in order for this routine to work correctly.
28: Many users may just want to use the monitoring routine
29: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
30: to print the extreme singular values at each iteration of the linear solve.
32: Estimates of the smallest singular value may be very inaccurate, especially if the Krylov method has not converged.
33: The largest singular value is usually accurate to within a few percent if the method has converged, but is still not
34: intended for eigenanalysis.
36: Disable restarts if using KSPGMRES, otherwise this estimate will only be using those iterations after the last
37: restart. See KSPGMRESSetRestart() for more details.
39: Level: advanced
41: .keywords: compute, extreme, singular, values
43: .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeEigenvalues(), KSP
44: @*/
45: PetscErrorCode KSPComputeExtremeSingularValues(KSP ksp,PetscReal *emax,PetscReal *emin)
46: {
53: if (!ksp->calc_sings) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Singular values not requested before KSPSetUp()");
55: if (ksp->ops->computeextremesingularvalues) {
56: (*ksp->ops->computeextremesingularvalues)(ksp,emax,emin);
57: } else {
58: *emin = -1.0;
59: *emax = -1.0;
60: }
61: return(0);
62: }
64: /*@
65: KSPComputeEigenvalues - Computes the extreme eigenvalues for the
66: preconditioned operator. Called after or during KSPSolve().
68: Not Collective
70: Input Parameter:
71: + ksp - iterative context obtained from KSPCreate()
72: - n - size of arrays r and c. The number of eigenvalues computed (neig) will, in
73: general, be less than this.
75: Output Parameters:
76: + r - real part of computed eigenvalues, provided by user with a dimension of at least n
77: . c - complex part of computed eigenvalues, provided by user with a dimension of at least n
78: - neig - actual number of eigenvalues computed (will be less than or equal to n)
80: Options Database Keys:
81: + -ksp_compute_eigenvalues - Prints eigenvalues to stdout
82: - -ksp_plot_eigenvalues - Plots eigenvalues in an x-window display
84: Notes:
85: The number of eigenvalues estimated depends on the size of the Krylov space
86: generated during the KSPSolve() ; for example, with
87: CG it corresponds to the number of CG iterations, for GMRES it is the number
88: of GMRES iterations SINCE the last restart. Any extra space in r[] and c[]
89: will be ignored.
91: KSPComputeEigenvalues() does not usually provide accurate estimates; it is
92: intended only for assistance in understanding the convergence of iterative
93: methods, not for eigenanalysis. For accurate computation of eigenvalues we recommend using
94: the excellent package SLEPc.
96: One must call KSPSetComputeEigenvalues() before calling KSPSetUp()
97: in order for this routine to work correctly.
99: Many users may just want to use the monitoring routine
100: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
101: to print the singular values at each iteration of the linear solve.
103: Level: advanced
105: .keywords: compute, extreme, singular, values
107: .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues(), KSP
108: @*/
109: PetscErrorCode KSPComputeEigenvalues(KSP ksp,PetscInt n,PetscReal r[],PetscReal c[],PetscInt *neig)
110: {
117: if (n<0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Requested < 0 Eigenvalues");
119: if (!ksp->calc_sings) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Eigenvalues not requested before KSPSetUp()");
121: if (n && ksp->ops->computeeigenvalues) {
122: (*ksp->ops->computeeigenvalues)(ksp,n,r,c,neig);
123: } else {
124: *neig = 0;
125: }
126: return(0);
127: }
129: /*@
130: KSPComputeRitz - Computes the Ritz or harmonic Ritz pairs associated to the
131: smallest or largest in modulus, for the preconditioned operator.
132: Called after KSPSolve().
134: Not Collective
136: Input Parameter:
137: + ksp - iterative context obtained from KSPCreate()
138: . ritz - PETSC_TRUE or PETSC_FALSE for ritz pairs or harmonic Ritz pairs, respectively
139: . small - PETSC_TRUE or PETSC_FALSE for smallest or largest (harmonic) Ritz values, respectively
140: . nrit - number of (harmonic) Ritz pairs to compute
142: Output Parameters:
143: + nrit - actual number of computed (harmonic) Ritz pairs
144: . S - multidimensional vector with Ritz vectors
145: . tetar - real part of the Ritz values
146: . tetai - imaginary part of the Ritz values
148: Notes:
149: -For GMRES, the (harmonic) Ritz pairs are computed from the Hessenberg matrix obtained during
150: the last complete cycle, or obtained at the end of the solution if the method is stopped before
151: a restart. Then, the number of actual (harmonic) Ritz pairs computed is less or equal to the restart
152: parameter for GMRES if a complete cycle has been performed or less or equal to the number of GMRES
153: iterations.
154: -Moreover, for real matrices, the (harmonic) Ritz pairs are possibly complex-valued. In such a case,
155: the routine selects the complex (harmonic) Ritz value and its conjugate, and two successive columns of S
156: are equal to the real and the imaginary parts of the associated vectors.
157: -the (harmonic) Ritz pairs are given in order of increasing (harmonic) Ritz values in modulus
158: -this is currently not implemented when PETSc is built with complex numbers
160: One must call KSPSetComputeRitz() before calling KSPSetUp()
161: in order for this routine to work correctly.
163: Level: advanced
165: .keywords: compute, ritz, values
167: .seealso: KSPSetComputeRitz(), KSP
168: @*/
169: PetscErrorCode KSPComputeRitz(KSP ksp,PetscBool ritz,PetscBool small,PetscInt *nrit,Vec S[],PetscReal tetar[],PetscReal tetai[])
170: {
175: if (!ksp->calc_ritz) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Ritz pairs not requested before KSPSetUp()");
176: if (ksp->ops->computeritz) {(*ksp->ops->computeritz)(ksp,ritz,small,nrit,S,tetar,tetai);}
177: return(0);
178: }
179: /*@
180: KSPSetUpOnBlocks - Sets up the preconditioner for each block in
181: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
182: methods.
184: Collective on KSP
186: Input Parameter:
187: . ksp - the KSP context
189: Notes:
190: KSPSetUpOnBlocks() is a routine that the user can optinally call for
191: more precise profiling (via -log_view) of the setup phase for these
192: block preconditioners. If the user does not call KSPSetUpOnBlocks(),
193: it will automatically be called from within KSPSolve().
195: Calling KSPSetUpOnBlocks() is the same as calling PCSetUpOnBlocks()
196: on the PC context within the KSP context.
198: Level: advanced
200: .keywords: setup, blocks
202: .seealso: PCSetUpOnBlocks(), KSPSetUp(), PCSetUp(), KSP
203: @*/
204: PetscErrorCode KSPSetUpOnBlocks(KSP ksp)
205: {
206: PC pc;
208: PCFailedReason pcreason;
212: KSPGetPC(ksp,&pc);
213: PCSetUpOnBlocks(pc);
214: PCGetSetUpFailedReason(pc,&pcreason);
215: if (pcreason) {
216: ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;
217: }
218: return(0);
219: }
221: /*@
222: KSPSetReusePreconditioner - reuse the current preconditioner, do not construct a new one even if the operator changes
224: Collective on KSP
226: Input Parameters:
227: + ksp - iterative context obtained from KSPCreate()
228: - flag - PETSC_TRUE to reuse the current preconditioner
230: Level: intermediate
232: .keywords: setup
234: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), PCSetReusePreconditioner(), KSP
235: @*/
236: PetscErrorCode KSPSetReusePreconditioner(KSP ksp,PetscBool flag)
237: {
238: PC pc;
243: KSPGetPC(ksp,&pc);
244: PCSetReusePreconditioner(pc,flag);
245: return(0);
246: }
248: /*@
249: KSPSetSkipPCSetFromOptions - prevents KSPSetFromOptions() from call PCSetFromOptions(). This is used if the same PC is shared by more than one KSP so its options are not resetable for each KSP
251: Collective on KSP
253: Input Parameters:
254: + ksp - iterative context obtained from KSPCreate()
255: - flag - PETSC_TRUE to skip calling the PCSetFromOptions()
257: Level: intermediate
259: .keywords: setup
261: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), PCSetReusePreconditioner(), KSP
262: @*/
263: PetscErrorCode KSPSetSkipPCSetFromOptions(KSP ksp,PetscBool flag)
264: {
267: ksp->skippcsetfromoptions = flag;
268: return(0);
269: }
271: /*@
272: KSPSetUp - Sets up the internal data structures for the
273: later use of an iterative solver.
275: Collective on KSP
277: Input Parameter:
278: . ksp - iterative context obtained from KSPCreate()
280: Level: developer
282: .keywords: setup
284: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), KSP
285: @*/
286: PetscErrorCode KSPSetUp(KSP ksp)
287: {
289: Mat A,B;
290: Mat mat,pmat;
291: MatNullSpace nullsp;
292: PCFailedReason pcreason;
293:
297: /* reset the convergence flag from the previous solves */
298: ksp->reason = KSP_CONVERGED_ITERATING;
300: if (!((PetscObject)ksp)->type_name) {
301: KSPSetType(ksp,KSPGMRES);
302: }
303: KSPSetUpNorms_Private(ksp,PETSC_TRUE,&ksp->normtype,&ksp->pc_side);
305: if (ksp->dmActive && !ksp->setupstage) {
306: /* first time in so build matrix and vector data structures using DM */
307: if (!ksp->vec_rhs) {DMCreateGlobalVector(ksp->dm,&ksp->vec_rhs);}
308: if (!ksp->vec_sol) {DMCreateGlobalVector(ksp->dm,&ksp->vec_sol);}
309: DMCreateMatrix(ksp->dm,&A);
310: KSPSetOperators(ksp,A,A);
311: PetscObjectDereference((PetscObject)A);
312: }
314: if (ksp->dmActive) {
315: DMKSP kdm;
316: DMGetDMKSP(ksp->dm,&kdm);
318: if (kdm->ops->computeinitialguess && ksp->setupstage != KSP_SETUP_NEWRHS) {
319: /* only computes initial guess the first time through */
320: (*kdm->ops->computeinitialguess)(ksp,ksp->vec_sol,kdm->initialguessctx);
321: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
322: }
323: if (kdm->ops->computerhs) {
324: (*kdm->ops->computerhs)(ksp,ksp->vec_rhs,kdm->rhsctx);
325: }
327: if (ksp->setupstage != KSP_SETUP_NEWRHS) {
328: if (kdm->ops->computeoperators) {
329: KSPGetOperators(ksp,&A,&B);
330: (*kdm->ops->computeoperators)(ksp,A,B,kdm->operatorsctx);
331: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You called KSPSetDM() but did not use DMKSPSetComputeOperators() or KSPSetDMActive(ksp,PETSC_FALSE);");
332: }
333: }
335: if (ksp->setupstage == KSP_SETUP_NEWRHS) return(0);
336: PetscLogEventBegin(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);
338: switch (ksp->setupstage) {
339: case KSP_SETUP_NEW:
340: (*ksp->ops->setup)(ksp);
341: break;
342: case KSP_SETUP_NEWMATRIX: { /* This should be replaced with a more general mechanism */
343: if (ksp->setupnewmatrix) {
344: (*ksp->ops->setup)(ksp);
345: }
346: } break;
347: default: break;
348: }
350: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
351: PCGetOperators(ksp->pc,&mat,&pmat);
352: /* scale the matrix if requested */
353: if (ksp->dscale) {
354: PetscScalar *xx;
355: PetscInt i,n;
356: PetscBool zeroflag = PETSC_FALSE;
357: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
358: if (!ksp->diagonal) { /* allocate vector to hold diagonal */
359: MatCreateVecs(pmat,&ksp->diagonal,0);
360: }
361: MatGetDiagonal(pmat,ksp->diagonal);
362: VecGetLocalSize(ksp->diagonal,&n);
363: VecGetArray(ksp->diagonal,&xx);
364: for (i=0; i<n; i++) {
365: if (xx[i] != 0.0) xx[i] = 1.0/PetscSqrtReal(PetscAbsScalar(xx[i]));
366: else {
367: xx[i] = 1.0;
368: zeroflag = PETSC_TRUE;
369: }
370: }
371: VecRestoreArray(ksp->diagonal,&xx);
372: if (zeroflag) {
373: PetscInfo(ksp,"Zero detected in diagonal of matrix, using 1 at those locations\n");
374: }
375: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
376: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
377: ksp->dscalefix2 = PETSC_FALSE;
378: }
379: PetscLogEventEnd(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);
380: PCSetErrorIfFailure(ksp->pc,ksp->errorifnotconverged);
381: PCSetUp(ksp->pc);
382: PCGetSetUpFailedReason(ksp->pc,&pcreason);
383: if (pcreason) {
384: ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;
385: }
387: MatGetNullSpace(mat,&nullsp);
388: if (nullsp) {
389: PetscBool test = PETSC_FALSE;
390: PetscOptionsGetBool(((PetscObject)ksp)->options,((PetscObject)ksp)->prefix,"-ksp_test_null_space",&test,NULL);
391: if (test) {
392: MatNullSpaceTest(nullsp,mat,NULL);
393: }
394: }
395: ksp->setupstage = KSP_SETUP_NEWRHS;
396: return(0);
397: }
399: /*@
400: KSPReasonView - Displays the reason a KSP solve converged or diverged to a viewer
402: Collective on KSP
404: Parameter:
405: + ksp - iterative context obtained from KSPCreate()
406: - viewer - the viewer to display the reason
409: Options Database Keys:
410: . -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations
412: Level: beginner
414: .keywords: solve, linear system
416: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
417: KSPSolveTranspose(), KSPGetIterationNumber(), KSP
418: @*/
419: PetscErrorCode KSPReasonView(KSP ksp,PetscViewer viewer)
420: {
422: PetscBool isAscii;
425: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
426: if (isAscii) {
427: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
428: if (ksp->reason > 0) {
429: if (((PetscObject) ksp)->prefix) {
430: PetscViewerASCIIPrintf(viewer,"Linear %s solve converged due to %s iterations %D\n",((PetscObject) ksp)->prefix,KSPConvergedReasons[ksp->reason],ksp->its);
431: } else {
432: PetscViewerASCIIPrintf(viewer,"Linear solve converged due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
433: }
434: } else {
435: if (((PetscObject) ksp)->prefix) {
436: PetscViewerASCIIPrintf(viewer,"Linear %s solve did not converge due to %s iterations %D\n",((PetscObject) ksp)->prefix,KSPConvergedReasons[ksp->reason],ksp->its);
437: } else {
438: PetscViewerASCIIPrintf(viewer,"Linear solve did not converge due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
439: }
440: if (ksp->reason == KSP_DIVERGED_PCSETUP_FAILED) {
441: PCFailedReason reason;
442: PCGetSetUpFailedReason(ksp->pc,&reason);
443: PetscViewerASCIIPrintf(viewer," PCSETUP_FAILED due to %s \n",PCFailedReasons[reason]);
444: }
445: }
446: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
447: }
448: return(0);
449: }
451: #if defined(PETSC_HAVE_THREADSAFETY)
452: #define KSPReasonViewFromOptions KSPReasonViewFromOptionsUnsafe
453: #else
454: #endif
455: /*@C
456: KSPReasonViewFromOptions - Processes command line options to determine if/how a KSPReason is to be viewed.
458: Collective on KSP
460: Input Parameters:
461: . ksp - the KSP object
463: Level: intermediate
465: @*/
466: PetscErrorCode KSPReasonViewFromOptions(KSP ksp)
467: {
468: PetscErrorCode ierr;
469: PetscViewer viewer;
470: PetscBool flg;
471: PetscViewerFormat format;
474: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_converged_reason",&viewer,&format,&flg);
475: if (flg) {
476: PetscViewerPushFormat(viewer,format);
477: KSPReasonView(ksp,viewer);
478: PetscViewerPopFormat(viewer);
479: PetscViewerDestroy(&viewer);
480: }
481: return(0);
482: }
484: #include <petscdraw.h>
485: /*@C
486: KSPSolve - Solves linear system.
488: Collective on KSP
490: Parameter:
491: + ksp - iterative context obtained from KSPCreate()
492: . b - the right hand side vector
493: - x - the solution (this may be the same vector as b, then b will be overwritten with answer)
495: Options Database Keys:
496: + -ksp_compute_eigenvalues - compute preconditioned operators eigenvalues
497: . -ksp_plot_eigenvalues - plot the computed eigenvalues in an X-window
498: . -ksp_plot_eigencontours - plot the computed eigenvalues in an X-window with contours
499: . -ksp_compute_eigenvalues_explicitly - compute the eigenvalues by forming the dense operator and using LAPACK
500: . -ksp_plot_eigenvalues_explicitly - plot the explicitly computing eigenvalues
501: . -ksp_view_mat binary - save matrix to the default binary viewer
502: . -ksp_view_pmat binary - save matrix used to build preconditioner to the default binary viewer
503: . -ksp_view_rhs binary - save right hand side vector to the default binary viewer
504: . -ksp_view_solution binary - save computed solution vector to the default binary viewer
505: (can be read later with src/ksp/examples/tutorials/ex10.c for testing solvers)
506: . -ksp_view_mat_explicit - for matrix-free operators, computes the matrix entries and views them
507: . -ksp_view_preconditioned_operator_explicit - computes the product of the preconditioner and matrix as an explicit matrix and views it
508: . -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations
509: . -ksp_final_residual - print 2-norm of true linear system residual at the end of the solution process
510: - -ksp_view - print the ksp data structure at the end of the system solution
512: Notes:
514: If one uses KSPSetDM() then x or b need not be passed. Use KSPGetSolution() to access the solution in this case.
516: The operator is specified with KSPSetOperators().
518: Call KSPGetConvergedReason() to determine if the solver converged or failed and
519: why. The number of iterations can be obtained from KSPGetIterationNumber().
521: If you provide a matrix that has a MatSetNullSpace() and MatSetTransposeNullSpace() this will use that information to solve singular systems
522: in the least squares sense with a norm minimizing solution.
523: $
524: $ A x = b where b = b_p + b_t where b_t is not in the range of A (and hence by the fundamental theorem of linear algebra is in the nullspace(A') see MatSetNullSpace()
525: $
526: $ KSP first removes b_t producing the linear system A x = b_p (which has multiple solutions) and solves this to find the ||x|| minimizing solution (and hence
527: $ it finds the solution x orthogonal to the nullspace(A). The algorithm is simply in each iteration of the Krylov method we remove the nullspace(A) from the search
528: $ direction thus the solution which is a linear combination of the search directions has no component in the nullspace(A).
529: $
530: $ We recommend always using GMRES for such singular systems.
531: $ If nullspace(A) = nullspace(A') (note symmetric matrices always satisfy this property) then both left and right preconditioning will work
532: $ If nullspace(A) != nullspace(A') then left preconditioning will work but right preconditioning may not work (or it may).
534: Developer Note: The reason we cannot always solve nullspace(A) != nullspace(A') systems with right preconditioning is because we need to remove at each iteration
535: the nullspace(AB) from the search direction. While we know the nullspace(A) the nullspace(AB) equals B^-1 times the nullspace(A) but except for trivial preconditioners
536: such as diagonal scaling we cannot apply the inverse of the preconditioner to a vector and thus cannot compute the nullspace(AB).
539: If using a direct method (e.g., via the KSP solver
540: KSPPREONLY and a preconditioner such as PCLU/PCILU),
541: then its=1. See KSPSetTolerances() and KSPConvergedDefault()
542: for more details.
544: Understanding Convergence:
545: The routines KSPMonitorSet(), KSPComputeEigenvalues(), and
546: KSPComputeEigenvaluesExplicitly() provide information on additional
547: options to monitor convergence and print eigenvalue information.
549: Level: beginner
551: .keywords: solve, linear system
553: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
554: KSPSolveTranspose(), KSPGetIterationNumber(), MatNullSpaceCreate(), MatSetNullSpace(), MatSetTransposeNullSpace(), KSP
555: @*/
556: PetscErrorCode KSPSolve(KSP ksp,Vec b,Vec x)
557: {
558: PetscErrorCode ierr;
559: PetscBool flag1,flag2,flag3,flg = PETSC_FALSE,inXisinB=PETSC_FALSE,guess_zero;
560: Mat mat,pmat;
561: MPI_Comm comm;
562: MatNullSpace nullsp;
563: Vec btmp,vec_rhs=0;
569: comm = PetscObjectComm((PetscObject)ksp);
570: if (x && x == b) {
571: if (!ksp->guess_zero) SETERRQ(comm,PETSC_ERR_ARG_INCOMP,"Cannot use x == b with nonzero initial guess");
572: VecDuplicate(b,&x);
573: inXisinB = PETSC_TRUE;
574: }
575: if (b) {
576: PetscObjectReference((PetscObject)b);
577: VecDestroy(&ksp->vec_rhs);
578: ksp->vec_rhs = b;
579: }
580: if (x) {
581: PetscObjectReference((PetscObject)x);
582: VecDestroy(&ksp->vec_sol);
583: ksp->vec_sol = x;
584: }
585: KSPViewFromOptions(ksp,NULL,"-ksp_view_pre");
587: if (ksp->presolve) {
588: (*ksp->presolve)(ksp,ksp->vec_rhs,ksp->vec_sol,ksp->prectx);
589: }
590: PetscLogEventBegin(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
592: /* reset the residual history list if requested */
593: if (ksp->res_hist_reset) ksp->res_hist_len = 0;
594: ksp->transpose_solve = PETSC_FALSE;
596: if (ksp->guess) {
597: PetscObjectState ostate,state;
599: KSPGuessSetUp(ksp->guess);
600: PetscObjectStateGet((PetscObject)ksp->vec_sol,&ostate);
601: KSPGuessFormGuess(ksp->guess,ksp->vec_rhs,ksp->vec_sol);
602: PetscObjectStateGet((PetscObject)ksp->vec_sol,&state);
603: if (state != ostate) {
604: ksp->guess_zero = PETSC_FALSE;
605: } else {
606: PetscInfo(ksp,"Using zero initial guess since the KSPGuess object did not change the vector\n");
607: ksp->guess_zero = PETSC_TRUE;
608: }
609: }
611: /* KSPSetUp() scales the matrix if needed */
612: KSPSetUp(ksp);
613: KSPSetUpOnBlocks(ksp);
615: VecLocked(ksp->vec_sol,3);
617: PCGetOperators(ksp->pc,&mat,&pmat);
618: /* diagonal scale RHS if called for */
619: if (ksp->dscale) {
620: VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
621: /* second time in, but matrix was scaled back to original */
622: if (ksp->dscalefix && ksp->dscalefix2) {
623: Mat mat,pmat;
625: PCGetOperators(ksp->pc,&mat,&pmat);
626: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
627: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
628: }
630: /* scale initial guess */
631: if (!ksp->guess_zero) {
632: if (!ksp->truediagonal) {
633: VecDuplicate(ksp->diagonal,&ksp->truediagonal);
634: VecCopy(ksp->diagonal,ksp->truediagonal);
635: VecReciprocal(ksp->truediagonal);
636: }
637: VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->truediagonal);
638: }
639: }
640: PCPreSolve(ksp->pc,ksp);
642: if (ksp->guess_zero) { VecSet(ksp->vec_sol,0.0);}
643: if (ksp->guess_knoll) { /* The Knoll trick is independent on the KSPGuess specified */
644: PCApply(ksp->pc,ksp->vec_rhs,ksp->vec_sol);
645: KSP_RemoveNullSpace(ksp,ksp->vec_sol);
646: ksp->guess_zero = PETSC_FALSE;
647: }
649: /* can we mark the initial guess as zero for this solve? */
650: guess_zero = ksp->guess_zero;
651: if (!ksp->guess_zero) {
652: PetscReal norm;
654: VecNormAvailable(ksp->vec_sol,NORM_2,&flg,&norm);
655: if (flg && !norm) ksp->guess_zero = PETSC_TRUE;
656: }
657: MatGetTransposeNullSpace(pmat,&nullsp);
658: if (nullsp) {
659: VecDuplicate(ksp->vec_rhs,&btmp);
660: VecCopy(ksp->vec_rhs,btmp);
661: MatNullSpaceRemove(nullsp,btmp);
662: vec_rhs = ksp->vec_rhs;
663: ksp->vec_rhs = btmp;
664: }
665: VecLockPush(ksp->vec_rhs);
666: if (ksp->reason == KSP_DIVERGED_PCSETUP_FAILED) {
667: VecSetInf(ksp->vec_sol);
668: }
669: (*ksp->ops->solve)(ksp);
670:
671: VecLockPop(ksp->vec_rhs);
672: if (nullsp) {
673: ksp->vec_rhs = vec_rhs;
674: VecDestroy(&btmp);
675: }
677: ksp->guess_zero = guess_zero;
680: if (!ksp->reason) SETERRQ(comm,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
681: ksp->totalits += ksp->its;
683: KSPReasonViewFromOptions(ksp);
684: PCPostSolve(ksp->pc,ksp);
686: /* diagonal scale solution if called for */
687: if (ksp->dscale) {
688: VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->diagonal);
689: /* unscale right hand side and matrix */
690: if (ksp->dscalefix) {
691: Mat mat,pmat;
693: VecReciprocal(ksp->diagonal);
694: VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
695: PCGetOperators(ksp->pc,&mat,&pmat);
696: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
697: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
698: VecReciprocal(ksp->diagonal);
699: ksp->dscalefix2 = PETSC_TRUE;
700: }
701: }
702: PetscLogEventEnd(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
703: if (ksp->postsolve) {
704: (*ksp->postsolve)(ksp,ksp->vec_rhs,ksp->vec_sol,ksp->postctx);
705: }
706: if (ksp->guess) {
707: KSPGuessUpdate(ksp->guess,ksp->vec_rhs,ksp->vec_sol);
708: }
710: PCGetOperators(ksp->pc,&mat,&pmat);
711: MatViewFromOptions(mat,(PetscObject)ksp,"-ksp_view_mat");
712: MatViewFromOptions(pmat,(PetscObject)ksp,"-ksp_view_pmat");
713: VecViewFromOptions(ksp->vec_rhs,(PetscObject)ksp,"-ksp_view_rhs");
715: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_compute_eigenvalues",NULL,NULL,&flag1);
716: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_plot_eigenvalues",NULL,NULL,&flag2);
717: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_plot_eigencontours",NULL,NULL,&flag3);
718: if (flag1 || flag2 || flag3) {
719: PetscInt nits,n,i,neig;
720: PetscReal *r,*c;
722: KSPGetIterationNumber(ksp,&nits);
723: n = nits+2;
725: if (!nits) {
726: PetscPrintf(comm,"Zero iterations in solver, cannot approximate any eigenvalues\n");
727: } else {
728: PetscMPIInt rank;
729: MPI_Comm_rank(comm,&rank);
730: PetscMalloc2(n,&r,n,&c);
731: KSPComputeEigenvalues(ksp,n,r,c,&neig);
732: if (flag1) {
733: PetscPrintf(comm,"Iteratively computed eigenvalues\n");
734: for (i=0; i<neig; i++) {
735: if (c[i] >= 0.0) {
736: PetscPrintf(comm,"%g + %gi\n",(double)r[i],(double)c[i]);
737: } else {
738: PetscPrintf(comm,"%g - %gi\n",(double)r[i],-(double)c[i]);
739: }
740: }
741: }
742: if (flag2 && !rank) {
743: PetscDraw draw;
744: PetscDrawSP drawsp;
746: if (!ksp->eigviewer) {
747: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Iteratively Computed Eigenvalues",PETSC_DECIDE,PETSC_DECIDE,400,400,&ksp->eigviewer);
748: }
749: PetscViewerDrawGetDraw(ksp->eigviewer,0,&draw);
750: PetscDrawSPCreate(draw,1,&drawsp);
751: for (i=0; i<neig; i++) {
752: PetscDrawSPAddPoint(drawsp,r+i,c+i);
753: }
754: PetscDrawSPDraw(drawsp,PETSC_TRUE);
755: PetscDrawSPSave(drawsp);
756: PetscDrawSPDestroy(&drawsp);
757: }
758: if (flag3 && !rank) {
759: KSPPlotEigenContours_Private(ksp,neig,r,c);
760: }
761: PetscFree2(r,c);
762: }
763: }
765: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_compute_singularvalues",NULL,NULL,&flag1);
766: if (flag1) {
767: PetscInt nits;
769: KSPGetIterationNumber(ksp,&nits);
770: if (!nits) {
771: PetscPrintf(comm,"Zero iterations in solver, cannot approximate any singular values\n");
772: } else {
773: PetscReal emax,emin;
775: KSPComputeExtremeSingularValues(ksp,&emax,&emin);
776: PetscPrintf(comm,"Iteratively computed extreme singular values: max %g min %g max/min %g\n",(double)emax,(double)emin,(double)(emax/emin));
777: }
778: }
780: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_compute_eigenvalues_explicitly",NULL,NULL,&flag1);
781: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_plot_eigenvalues_explicitly",NULL,NULL,&flag2);
782: if (flag1 || flag2) {
783: PetscInt n,i;
784: PetscReal *r,*c;
785: PetscMPIInt rank;
786: MPI_Comm_rank(comm,&rank);
787: VecGetSize(ksp->vec_sol,&n);
788: PetscMalloc2(n,&r,n,&c);
789: KSPComputeEigenvaluesExplicitly(ksp,n,r,c);
790: if (flag1) {
791: PetscPrintf(comm,"Explicitly computed eigenvalues\n");
792: for (i=0; i<n; i++) {
793: if (c[i] >= 0.0) {
794: PetscPrintf(comm,"%g + %gi\n",(double)r[i],(double)c[i]);
795: } else {
796: PetscPrintf(comm,"%g - %gi\n",(double)r[i],-(double)c[i]);
797: }
798: }
799: }
800: if (flag2 && !rank) {
801: PetscDraw draw;
802: PetscDrawSP drawsp;
804: if (!ksp->eigviewer) {
805: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Explicitly Computed Eigenvalues",PETSC_DECIDE,PETSC_DECIDE,400,400,&ksp->eigviewer);
806: }
807: PetscViewerDrawGetDraw(ksp->eigviewer,0,&draw);
808: PetscDrawSPCreate(draw,1,&drawsp);
809: PetscDrawSPReset(drawsp);
810: for (i=0; i<n; i++) {
811: PetscDrawSPAddPoint(drawsp,r+i,c+i);
812: }
813: PetscDrawSPDraw(drawsp,PETSC_TRUE);
814: PetscDrawSPSave(drawsp);
815: PetscDrawSPDestroy(&drawsp);
816: }
817: PetscFree2(r,c);
818: }
820: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_view_mat_explicit",NULL,NULL,&flag2);
821: if (flag2) {
822: Mat A,B;
823: PCGetOperators(ksp->pc,&A,NULL);
824: MatComputeExplicitOperator(A,&B);
825: MatViewFromOptions(B,(PetscObject)ksp,"-ksp_view_mat_explicit");
826: MatDestroy(&B);
827: }
828: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_view_preconditioned_operator_explicit",NULL,NULL,&flag2);
829: if (flag2) {
830: Mat B;
831: KSPComputeExplicitOperator(ksp,&B);
832: MatViewFromOptions(B,(PetscObject)ksp,"-ksp_view_preconditioned_operator_explicit");
833: MatDestroy(&B);
834: }
835: KSPViewFromOptions(ksp,NULL,"-ksp_view");
837: PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_final_residual",NULL,NULL,&flg);
838: if (flg) {
839: Mat A;
840: Vec t;
841: PetscReal norm;
842: if (ksp->dscale && !ksp->dscalefix) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot compute final scale with -ksp_diagonal_scale except also with -ksp_diagonal_scale_fix");
843: PCGetOperators(ksp->pc,&A,NULL);
844: VecDuplicate(ksp->vec_rhs,&t);
845: KSP_MatMult(ksp,A,ksp->vec_sol,t);
846: VecAYPX(t, -1.0, ksp->vec_rhs);
847: VecNorm(t,NORM_2,&norm);
848: VecDestroy(&t);
849: PetscPrintf(comm,"KSP final norm of residual %g\n",(double)norm);
850: }
851: VecViewFromOptions(ksp->vec_sol,(PetscObject)ksp,"-ksp_view_solution");
853: if (inXisinB) {
854: VecCopy(x,b);
855: VecDestroy(&x);
856: }
857: PetscObjectSAWsBlock((PetscObject)ksp);
858: if (ksp->errorifnotconverged && ksp->reason < 0) SETERRQ(comm,PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged");
859: return(0);
860: }
862: /*@
863: KSPSolveTranspose - Solves the transpose of a linear system.
865: Collective on KSP
867: Input Parameter:
868: + ksp - iterative context obtained from KSPCreate()
869: . b - right hand side vector
870: - x - solution vector
872: Notes: For complex numbers this solve the non-Hermitian transpose system.
874: This currently does NOT correctly use the null space of the operator and its transpose for solving singular systems.
876: Developer Notes: We need to implement a KSPSolveHermitianTranspose()
878: Level: developer
880: .keywords: solve, linear system
882: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
883: KSPSolve(), KSP
884: @*/
886: PetscErrorCode KSPSolveTranspose(KSP ksp,Vec b,Vec x)
887: {
889: PetscBool inXisinB=PETSC_FALSE;
890: Vec vec_rhs = 0,btmp;
891: Mat mat,pmat;
892: MatNullSpace nullsp;
898: if (x == b) {
899: VecDuplicate(b,&x);
900: inXisinB = PETSC_TRUE;
901: }
902: PetscObjectReference((PetscObject)b);
903: PetscObjectReference((PetscObject)x);
904: VecDestroy(&ksp->vec_rhs);
905: VecDestroy(&ksp->vec_sol);
907: ksp->vec_rhs = b;
908: ksp->vec_sol = x;
909: ksp->transpose_solve = PETSC_TRUE;
911: if (ksp->guess) {
912: PetscObjectState ostate,state;
914: KSPGuessSetUp(ksp->guess);
915: PetscObjectStateGet((PetscObject)ksp->vec_sol,&ostate);
916: KSPGuessFormGuess(ksp->guess,ksp->vec_rhs,ksp->vec_sol);
917: PetscObjectStateGet((PetscObject)ksp->vec_sol,&state);
918: if (state != ostate) {
919: ksp->guess_zero = PETSC_FALSE;
920: } else {
921: PetscInfo(ksp,"Using zero initial guess since the KSPGuess object did not change the vector\n");
922: ksp->guess_zero = PETSC_TRUE;
923: }
924: }
926: KSPSetUp(ksp);
927: KSPSetUpOnBlocks(ksp);
928: if (ksp->guess_zero) { VecSet(ksp->vec_sol,0.0);}
930: PCGetOperators(ksp->pc,&mat,&pmat);
931: MatGetNullSpace(pmat,&nullsp);
932: if (nullsp) {
933: VecDuplicate(ksp->vec_rhs,&btmp);
934: VecCopy(ksp->vec_rhs,btmp);
935: MatNullSpaceRemove(nullsp,btmp);
936: vec_rhs = ksp->vec_rhs;
937: ksp->vec_rhs = btmp;
938: }
940: (*ksp->ops->solve)(ksp);
941: if (nullsp) {
942: ksp->vec_rhs = vec_rhs;
943: VecDestroy(&btmp);
944: }
945: if (!ksp->reason) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
946: KSPReasonViewFromOptions(ksp);
948: if (ksp->guess) {
949: KSPGuessUpdate(ksp->guess,ksp->vec_rhs,ksp->vec_sol);
950: }
952: MatViewFromOptions(mat,(PetscObject)ksp,"-ksp_view_mat");
953: MatViewFromOptions(pmat,(PetscObject)ksp,"-ksp_view_pmat");
954: VecViewFromOptions(ksp->vec_rhs,(PetscObject)ksp,"-ksp_view_rhs");
955: VecViewFromOptions(ksp->vec_sol,(PetscObject)ksp,"-ksp_view_solution");
957: if (inXisinB) {
958: VecCopy(x,b);
959: VecDestroy(&x);
960: }
961: if (ksp->errorifnotconverged && ksp->reason < 0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged");
962: return(0);
963: }
965: /*@
966: KSPReset - Resets a KSP context to the kspsetupcalled = 0 state and removes any allocated Vecs and Mats
968: Collective on KSP
970: Input Parameter:
971: . ksp - iterative context obtained from KSPCreate()
973: Level: beginner
975: .keywords: destroy
977: .seealso: KSPCreate(), KSPSetUp(), KSPSolve(), KSP
978: @*/
979: PetscErrorCode KSPReset(KSP ksp)
980: {
985: if (!ksp) return(0);
986: if (ksp->ops->reset) {
987: (*ksp->ops->reset)(ksp);
988: }
989: if (ksp->pc) {PCReset(ksp->pc);}
990: if (ksp->guess) {
991: KSPGuess guess = ksp->guess;
992: if (guess->ops->reset) { (*guess->ops->reset)(guess); }
993: }
994: VecDestroyVecs(ksp->nwork,&ksp->work);
995: VecDestroy(&ksp->vec_rhs);
996: VecDestroy(&ksp->vec_sol);
997: VecDestroy(&ksp->diagonal);
998: VecDestroy(&ksp->truediagonal);
1000: ksp->setupstage = KSP_SETUP_NEW;
1001: return(0);
1002: }
1004: /*@
1005: KSPDestroy - Destroys KSP context.
1007: Collective on KSP
1009: Input Parameter:
1010: . ksp - iterative context obtained from KSPCreate()
1012: Level: beginner
1014: .keywords: destroy
1016: .seealso: KSPCreate(), KSPSetUp(), KSPSolve(), KSP
1017: @*/
1018: PetscErrorCode KSPDestroy(KSP *ksp)
1019: {
1021: PC pc;
1024: if (!*ksp) return(0);
1026: if (--((PetscObject)(*ksp))->refct > 0) {*ksp = 0; return(0);}
1028: PetscObjectSAWsViewOff((PetscObject)*ksp);
1030: /*
1031: Avoid a cascading call to PCReset(ksp->pc) from the following call:
1032: PCReset() shouldn't be called from KSPDestroy() as it is unprotected by pc's
1033: refcount (and may be shared, e.g., by other ksps).
1034: */
1035: pc = (*ksp)->pc;
1036: (*ksp)->pc = NULL;
1037: KSPReset((*ksp));
1038: (*ksp)->pc = pc;
1039: if ((*ksp)->ops->destroy) {(*(*ksp)->ops->destroy)(*ksp);}
1041: KSPGuessDestroy(&(*ksp)->guess);
1042: DMDestroy(&(*ksp)->dm);
1043: PCDestroy(&(*ksp)->pc);
1044: PetscFree((*ksp)->res_hist_alloc);
1045: if ((*ksp)->convergeddestroy) {
1046: (*(*ksp)->convergeddestroy)((*ksp)->cnvP);
1047: }
1048: KSPMonitorCancel((*ksp));
1049: PetscViewerDestroy(&(*ksp)->eigviewer);
1050: PetscHeaderDestroy(ksp);
1051: return(0);
1052: }
1054: /*@
1055: KSPSetPCSide - Sets the preconditioning side.
1057: Logically Collective on KSP
1059: Input Parameter:
1060: . ksp - iterative context obtained from KSPCreate()
1062: Output Parameter:
1063: . side - the preconditioning side, where side is one of
1064: .vb
1065: PC_LEFT - left preconditioning (default)
1066: PC_RIGHT - right preconditioning
1067: PC_SYMMETRIC - symmetric preconditioning
1068: .ve
1070: Options Database Keys:
1071: . -ksp_pc_side <right,left,symmetric>
1073: Notes:
1074: Left preconditioning is used by default for most Krylov methods except KSPFGMRES which only supports right preconditioning.
1076: For methods changing the side of the preconditioner changes the norm type that is used, see KSPSetNormType().
1078: Symmetric preconditioning is currently available only for the KSPQCG method. Note, however, that
1079: symmetric preconditioning can be emulated by using either right or left
1080: preconditioning and a pre or post processing step.
1082: Setting the PC side often affects the default norm type. See KSPSetNormType() for details.
1084: Level: intermediate
1086: .keywords: set, right, left, symmetric, side, preconditioner, flag
1088: .seealso: KSPGetPCSide(), KSPSetNormType(), KSPGetNormType(), KSP
1089: @*/
1090: PetscErrorCode KSPSetPCSide(KSP ksp,PCSide side)
1091: {
1095: ksp->pc_side = ksp->pc_side_set = side;
1096: return(0);
1097: }
1099: /*@
1100: KSPGetPCSide - Gets the preconditioning side.
1102: Not Collective
1104: Input Parameter:
1105: . ksp - iterative context obtained from KSPCreate()
1107: Output Parameter:
1108: . side - the preconditioning side, where side is one of
1109: .vb
1110: PC_LEFT - left preconditioning (default)
1111: PC_RIGHT - right preconditioning
1112: PC_SYMMETRIC - symmetric preconditioning
1113: .ve
1115: Level: intermediate
1117: .keywords: get, right, left, symmetric, side, preconditioner, flag
1119: .seealso: KSPSetPCSide(), KSP
1120: @*/
1121: PetscErrorCode KSPGetPCSide(KSP ksp,PCSide *side)
1122: {
1128: KSPSetUpNorms_Private(ksp,PETSC_TRUE,&ksp->normtype,&ksp->pc_side);
1129: *side = ksp->pc_side;
1130: return(0);
1131: }
1133: /*@
1134: KSPGetTolerances - Gets the relative, absolute, divergence, and maximum
1135: iteration tolerances used by the default KSP convergence tests.
1137: Not Collective
1139: Input Parameter:
1140: . ksp - the Krylov subspace context
1142: Output Parameters:
1143: + rtol - the relative convergence tolerance
1144: . abstol - the absolute convergence tolerance
1145: . dtol - the divergence tolerance
1146: - maxits - maximum number of iterations
1148: Notes:
1149: The user can specify NULL for any parameter that is not needed.
1151: Level: intermediate
1153: .keywords: get, tolerance, absolute, relative, divergence, convergence,
1154: maximum, iterations
1156: .seealso: KSPSetTolerances(), KSP
1157: @*/
1158: PetscErrorCode KSPGetTolerances(KSP ksp,PetscReal *rtol,PetscReal *abstol,PetscReal *dtol,PetscInt *maxits)
1159: {
1162: if (abstol) *abstol = ksp->abstol;
1163: if (rtol) *rtol = ksp->rtol;
1164: if (dtol) *dtol = ksp->divtol;
1165: if (maxits) *maxits = ksp->max_it;
1166: return(0);
1167: }
1169: /*@
1170: KSPSetTolerances - Sets the relative, absolute, divergence, and maximum
1171: iteration tolerances used by the default KSP convergence testers.
1173: Logically Collective on KSP
1175: Input Parameters:
1176: + ksp - the Krylov subspace context
1177: . rtol - the relative convergence tolerance, relative decrease in the (possibly preconditioned) residual norm
1178: . abstol - the absolute convergence tolerance absolute size of the (possibly preconditioned) residual norm
1179: . dtol - the divergence tolerance, amount (possibly preconditioned) residual norm can increase before KSPConvergedDefault() concludes that the method is diverging
1180: - maxits - maximum number of iterations to use
1182: Options Database Keys:
1183: + -ksp_atol <abstol> - Sets abstol
1184: . -ksp_rtol <rtol> - Sets rtol
1185: . -ksp_divtol <dtol> - Sets dtol
1186: - -ksp_max_it <maxits> - Sets maxits
1188: Notes:
1189: Use PETSC_DEFAULT to retain the default value of any of the tolerances.
1191: See KSPConvergedDefault() for details how these parameters are used in the default convergence test. See also KSPSetConvergenceTest()
1192: for setting user-defined stopping criteria.
1194: Level: intermediate
1196: .keywords: set, tolerance, absolute, relative, divergence,
1197: convergence, maximum, iterations
1199: .seealso: KSPGetTolerances(), KSPConvergedDefault(), KSPSetConvergenceTest(), KSP
1200: @*/
1201: PetscErrorCode KSPSetTolerances(KSP ksp,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt maxits)
1202: {
1210: if (rtol != PETSC_DEFAULT) {
1211: if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %g must be non-negative and less than 1.0",(double)rtol);
1212: ksp->rtol = rtol;
1213: }
1214: if (abstol != PETSC_DEFAULT) {
1215: if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol);
1216: ksp->abstol = abstol;
1217: }
1218: if (dtol != PETSC_DEFAULT) {
1219: if (dtol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Divergence tolerance %g must be larger than 1.0",(double)dtol);
1220: ksp->divtol = dtol;
1221: }
1222: if (maxits != PETSC_DEFAULT) {
1223: if (maxits < 0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxits);
1224: ksp->max_it = maxits;
1225: }
1226: return(0);
1227: }
1229: /*@
1230: KSPSetInitialGuessNonzero - Tells the iterative solver that the
1231: initial guess is nonzero; otherwise KSP assumes the initial guess
1232: is to be zero (and thus zeros it out before solving).
1234: Logically Collective on KSP
1236: Input Parameters:
1237: + ksp - iterative context obtained from KSPCreate()
1238: - flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero
1240: Options database keys:
1241: . -ksp_initial_guess_nonzero : use nonzero initial guess; this takes an optional truth value (0/1/no/yes/true/false)
1243: Level: beginner
1245: Notes:
1246: If this is not called the X vector is zeroed in the call to KSPSolve().
1248: .keywords: set, initial guess, nonzero
1250: .seealso: KSPGetInitialGuessNonzero(), KSPSetGuessType(), KSPGuessType, KSP
1251: @*/
1252: PetscErrorCode KSPSetInitialGuessNonzero(KSP ksp,PetscBool flg)
1253: {
1257: ksp->guess_zero = (PetscBool) !(int)flg;
1258: return(0);
1259: }
1261: /*@
1262: KSPGetInitialGuessNonzero - Determines whether the KSP solver is using
1263: a zero initial guess.
1265: Not Collective
1267: Input Parameter:
1268: . ksp - iterative context obtained from KSPCreate()
1270: Output Parameter:
1271: . flag - PETSC_TRUE if guess is nonzero, else PETSC_FALSE
1273: Level: intermediate
1275: .keywords: set, initial guess, nonzero
1277: .seealso: KSPSetInitialGuessNonzero(), KSP
1278: @*/
1279: PetscErrorCode KSPGetInitialGuessNonzero(KSP ksp,PetscBool *flag)
1280: {
1284: if (ksp->guess_zero) *flag = PETSC_FALSE;
1285: else *flag = PETSC_TRUE;
1286: return(0);
1287: }
1289: /*@
1290: KSPSetErrorIfNotConverged - Causes KSPSolve() to generate an error if the solver has not converged.
1292: Logically Collective on KSP
1294: Input Parameters:
1295: + ksp - iterative context obtained from KSPCreate()
1296: - flg - PETSC_TRUE indicates you want the error generated
1298: Options database keys:
1299: . -ksp_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false)
1301: Level: intermediate
1303: Notes:
1304: Normally PETSc continues if a linear solver fails to converge, you can call KSPGetConvergedReason() after a KSPSolve()
1305: to determine if it has converged.
1308: .seealso: KSPGetErrorIfNotConverged(), KSP
1309: @*/
1310: PetscErrorCode KSPSetErrorIfNotConverged(KSP ksp,PetscBool flg)
1311: {
1315: ksp->errorifnotconverged = flg;
1316: return(0);
1317: }
1319: /*@
1320: KSPGetErrorIfNotConverged - Will KSPSolve() generate an error if the solver does not converge?
1322: Not Collective
1324: Input Parameter:
1325: . ksp - iterative context obtained from KSPCreate()
1327: Output Parameter:
1328: . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE
1330: Level: intermediate
1332: .seealso: KSPSetErrorIfNotConverged(), KSP
1333: @*/
1334: PetscErrorCode KSPGetErrorIfNotConverged(KSP ksp,PetscBool *flag)
1335: {
1339: *flag = ksp->errorifnotconverged;
1340: return(0);
1341: }
1343: /*@
1344: KSPSetInitialGuessKnoll - Tells the iterative solver to use PCApply(pc,b,..) to compute the initial guess (The Knoll trick)
1346: Logically Collective on KSP
1348: Input Parameters:
1349: + ksp - iterative context obtained from KSPCreate()
1350: - flg - PETSC_TRUE or PETSC_FALSE
1352: Level: advanced
1354: Developer Note: the Knoll trick is not currently implemented using the KSPGuess class
1356: .keywords: set, initial guess, nonzero
1358: .seealso: KSPGetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero(), KSP
1359: @*/
1360: PetscErrorCode KSPSetInitialGuessKnoll(KSP ksp,PetscBool flg)
1361: {
1365: ksp->guess_knoll = flg;
1366: return(0);
1367: }
1369: /*@
1370: KSPGetInitialGuessKnoll - Determines whether the KSP solver is using the Knoll trick (using PCApply(pc,b,...) to compute
1371: the initial guess
1373: Not Collective
1375: Input Parameter:
1376: . ksp - iterative context obtained from KSPCreate()
1378: Output Parameter:
1379: . flag - PETSC_TRUE if using Knoll trick, else PETSC_FALSE
1381: Level: advanced
1383: .keywords: set, initial guess, nonzero
1385: .seealso: KSPSetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero(), KSP
1386: @*/
1387: PetscErrorCode KSPGetInitialGuessKnoll(KSP ksp,PetscBool *flag)
1388: {
1392: *flag = ksp->guess_knoll;
1393: return(0);
1394: }
1396: /*@
1397: KSPGetComputeSingularValues - Gets the flag indicating whether the extreme singular
1398: values will be calculated via a Lanczos or Arnoldi process as the linear
1399: system is solved.
1401: Not Collective
1403: Input Parameter:
1404: . ksp - iterative context obtained from KSPCreate()
1406: Output Parameter:
1407: . flg - PETSC_TRUE or PETSC_FALSE
1409: Options Database Key:
1410: . -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()
1412: Notes:
1413: Currently this option is not valid for all iterative methods.
1415: Many users may just want to use the monitoring routine
1416: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1417: to print the singular values at each iteration of the linear solve.
1419: Level: advanced
1421: .keywords: set, compute, singular values
1423: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue(), KSP
1424: @*/
1425: PetscErrorCode KSPGetComputeSingularValues(KSP ksp,PetscBool *flg)
1426: {
1430: *flg = ksp->calc_sings;
1431: return(0);
1432: }
1434: /*@
1435: KSPSetComputeSingularValues - Sets a flag so that the extreme singular
1436: values will be calculated via a Lanczos or Arnoldi process as the linear
1437: system is solved.
1439: Logically Collective on KSP
1441: Input Parameters:
1442: + ksp - iterative context obtained from KSPCreate()
1443: - flg - PETSC_TRUE or PETSC_FALSE
1445: Options Database Key:
1446: . -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()
1448: Notes:
1449: Currently this option is not valid for all iterative methods.
1451: Many users may just want to use the monitoring routine
1452: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1453: to print the singular values at each iteration of the linear solve.
1455: Level: advanced
1457: .keywords: set, compute, singular values
1459: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue(), KSP
1460: @*/
1461: PetscErrorCode KSPSetComputeSingularValues(KSP ksp,PetscBool flg)
1462: {
1466: ksp->calc_sings = flg;
1467: return(0);
1468: }
1470: /*@
1471: KSPGetComputeEigenvalues - Gets the flag indicating that the extreme eigenvalues
1472: values will be calculated via a Lanczos or Arnoldi process as the linear
1473: system is solved.
1475: Not Collective
1477: Input Parameter:
1478: . ksp - iterative context obtained from KSPCreate()
1480: Output Parameter:
1481: . flg - PETSC_TRUE or PETSC_FALSE
1483: Notes:
1484: Currently this option is not valid for all iterative methods.
1486: Level: advanced
1488: .keywords: set, compute, eigenvalues
1490: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly(), KSP
1491: @*/
1492: PetscErrorCode KSPGetComputeEigenvalues(KSP ksp,PetscBool *flg)
1493: {
1497: *flg = ksp->calc_sings;
1498: return(0);
1499: }
1501: /*@
1502: KSPSetComputeEigenvalues - Sets a flag so that the extreme eigenvalues
1503: values will be calculated via a Lanczos or Arnoldi process as the linear
1504: system is solved.
1506: Logically Collective on KSP
1508: Input Parameters:
1509: + ksp - iterative context obtained from KSPCreate()
1510: - flg - PETSC_TRUE or PETSC_FALSE
1512: Notes:
1513: Currently this option is not valid for all iterative methods.
1515: Level: advanced
1517: .keywords: set, compute, eigenvalues
1519: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly(), KSP
1520: @*/
1521: PetscErrorCode KSPSetComputeEigenvalues(KSP ksp,PetscBool flg)
1522: {
1526: ksp->calc_sings = flg;
1527: return(0);
1528: }
1530: /*@
1531: KSPSetComputeRitz - Sets a flag so that the Ritz or harmonic Ritz pairs
1532: will be calculated via a Lanczos or Arnoldi process as the linear
1533: system is solved.
1535: Logically Collective on KSP
1537: Input Parameters:
1538: + ksp - iterative context obtained from KSPCreate()
1539: - flg - PETSC_TRUE or PETSC_FALSE
1541: Notes:
1542: Currently this option is only valid for the GMRES method.
1544: Level: advanced
1546: .keywords: set, compute, ritz
1548: .seealso: KSPComputeRitz(), KSP
1549: @*/
1550: PetscErrorCode KSPSetComputeRitz(KSP ksp, PetscBool flg)
1551: {
1555: ksp->calc_ritz = flg;
1556: return(0);
1557: }
1559: /*@
1560: KSPGetRhs - Gets the right-hand-side vector for the linear system to
1561: be solved.
1563: Not Collective
1565: Input Parameter:
1566: . ksp - iterative context obtained from KSPCreate()
1568: Output Parameter:
1569: . r - right-hand-side vector
1571: Level: developer
1573: .keywords: get, right-hand-side, rhs
1575: .seealso: KSPGetSolution(), KSPSolve(), KSP
1576: @*/
1577: PetscErrorCode KSPGetRhs(KSP ksp,Vec *r)
1578: {
1582: *r = ksp->vec_rhs;
1583: return(0);
1584: }
1586: /*@
1587: KSPGetSolution - Gets the location of the solution for the
1588: linear system to be solved. Note that this may not be where the solution
1589: is stored during the iterative process; see KSPBuildSolution().
1591: Not Collective
1593: Input Parameters:
1594: . ksp - iterative context obtained from KSPCreate()
1596: Output Parameters:
1597: . v - solution vector
1599: Level: developer
1601: .keywords: get, solution
1603: .seealso: KSPGetRhs(), KSPBuildSolution(), KSPSolve(), KSP
1604: @*/
1605: PetscErrorCode KSPGetSolution(KSP ksp,Vec *v)
1606: {
1610: *v = ksp->vec_sol;
1611: return(0);
1612: }
1614: /*@
1615: KSPSetPC - Sets the preconditioner to be used to calculate the
1616: application of the preconditioner on a vector.
1618: Collective on KSP
1620: Input Parameters:
1621: + ksp - iterative context obtained from KSPCreate()
1622: - pc - the preconditioner object
1624: Notes:
1625: Use KSPGetPC() to retrieve the preconditioner context (for example,
1626: to free it at the end of the computations).
1628: Level: developer
1630: .keywords: set, precondition, Binv
1632: .seealso: KSPGetPC(), KSP
1633: @*/
1634: PetscErrorCode KSPSetPC(KSP ksp,PC pc)
1635: {
1642: PetscObjectReference((PetscObject)pc);
1643: PCDestroy(&ksp->pc);
1644: ksp->pc = pc;
1645: PetscLogObjectParent((PetscObject)ksp,(PetscObject)ksp->pc);
1646: return(0);
1647: }
1649: /*@
1650: KSPGetPC - Returns a pointer to the preconditioner context
1651: set with KSPSetPC().
1653: Not Collective
1655: Input Parameters:
1656: . ksp - iterative context obtained from KSPCreate()
1658: Output Parameter:
1659: . pc - preconditioner context
1661: Level: developer
1663: .keywords: get, preconditioner, Binv
1665: .seealso: KSPSetPC(), KSP
1666: @*/
1667: PetscErrorCode KSPGetPC(KSP ksp,PC *pc)
1668: {
1674: if (!ksp->pc) {
1675: PCCreate(PetscObjectComm((PetscObject)ksp),&ksp->pc);
1676: PetscObjectIncrementTabLevel((PetscObject)ksp->pc,(PetscObject)ksp,0);
1677: PetscLogObjectParent((PetscObject)ksp,(PetscObject)ksp->pc);
1678: }
1679: *pc = ksp->pc;
1680: return(0);
1681: }
1683: /*@
1684: KSPMonitor - runs the user provided monitor routines, if they exist
1686: Collective on KSP
1688: Input Parameters:
1689: + ksp - iterative context obtained from KSPCreate()
1690: . it - iteration number
1691: - rnorm - relative norm of the residual
1693: Notes:
1694: This routine is called by the KSP implementations.
1695: It does not typically need to be called by the user.
1697: Level: developer
1699: .seealso: KSPMonitorSet()
1700: @*/
1701: PetscErrorCode KSPMonitor(KSP ksp,PetscInt it,PetscReal rnorm)
1702: {
1703: PetscInt i, n = ksp->numbermonitors;
1707: for (i=0; i<n; i++) {
1708: (*ksp->monitor[i])(ksp,it,rnorm,ksp->monitorcontext[i]);
1709: }
1710: return(0);
1711: }
1713: /*
1715: Checks if two monitors are identical; if they are then it destroys the new one
1716: */
1717: PetscErrorCode PetscMonitorCompare(PetscErrorCode (*nmon)(void),void *nmctx,PetscErrorCode (*nmdestroy)(void**),PetscErrorCode (*mon)(void),void *mctx,PetscErrorCode (*mdestroy)(void**),PetscBool *identical)
1718: {
1719: *identical = PETSC_FALSE;
1720: if (nmon == mon && nmdestroy == mdestroy) {
1721: if (nmctx == mctx) *identical = PETSC_TRUE;
1722: else if (nmdestroy == (PetscErrorCode (*)(void**)) PetscViewerAndFormatDestroy) {
1723: PetscViewerAndFormat *old = (PetscViewerAndFormat*)mctx, *newo = (PetscViewerAndFormat*)nmctx;
1724: if (old->viewer == newo->viewer && old->format == newo->format) *identical = PETSC_TRUE;
1725: }
1726: if (*identical) {
1727: if (mdestroy) {
1729: (*mdestroy)(&nmctx);
1730: }
1731: }
1732: }
1733: return(0);
1734: }
1736: /*@C
1737: KSPMonitorSet - Sets an ADDITIONAL function to be called at every iteration to monitor
1738: the residual/error etc.
1740: Logically Collective on KSP
1742: Input Parameters:
1743: + ksp - iterative context obtained from KSPCreate()
1744: . monitor - pointer to function (if this is NULL, it turns off monitoring
1745: . mctx - [optional] context for private data for the
1746: monitor routine (use NULL if no context is desired)
1747: - monitordestroy - [optional] routine that frees monitor context
1748: (may be NULL)
1750: Calling Sequence of monitor:
1751: $ monitor (KSP ksp, int it, PetscReal rnorm, void *mctx)
1753: + ksp - iterative context obtained from KSPCreate()
1754: . it - iteration number
1755: . rnorm - (estimated) 2-norm of (preconditioned) residual
1756: - mctx - optional monitoring context, as set by KSPMonitorSet()
1758: Options Database Keys:
1759: + -ksp_monitor - sets KSPMonitorDefault()
1760: . -ksp_monitor_true_residual - sets KSPMonitorTrueResidualNorm()
1761: . -ksp_monitor_max - sets KSPMonitorTrueResidualMaxNorm()
1762: . -ksp_monitor_lg_residualnorm - sets line graph monitor,
1763: uses KSPMonitorLGResidualNormCreate()
1764: . -ksp_monitor_lg_true_residualnorm - sets line graph monitor,
1765: uses KSPMonitorLGResidualNormCreate()
1766: . -ksp_monitor_singular_value - sets KSPMonitorSingularValue()
1767: - -ksp_monitor_cancel - cancels all monitors that have
1768: been hardwired into a code by
1769: calls to KSPMonitorSet(), but
1770: does not cancel those set via
1771: the options database.
1773: Notes:
1774: The default is to do nothing. To print the residual, or preconditioned
1775: residual if KSPSetNormType(ksp,KSP_NORM_PRECONDITIONED) was called, use
1776: KSPMonitorDefault() as the monitoring routine, with a ASCII viewer as the
1777: context.
1779: Several different monitoring routines may be set by calling
1780: KSPMonitorSet() multiple times; all will be called in the
1781: order in which they were set.
1783: Fortran notes: Only a single monitor function can be set for each KSP object
1785: Level: beginner
1787: .keywords: set, monitor
1789: .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSPMonitorCancel(), KSP
1790: @*/
1791: PetscErrorCode KSPMonitorSet(KSP ksp,PetscErrorCode (*monitor)(KSP,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
1792: {
1793: PetscInt i;
1795: PetscBool identical;
1799: for (i=0; i<ksp->numbermonitors;i++) {
1800: PetscMonitorCompare((PetscErrorCode (*)(void))monitor,mctx,monitordestroy,(PetscErrorCode (*)(void))ksp->monitor[i],ksp->monitorcontext[i],ksp->monitordestroy[i],&identical);
1801: if (identical) return(0);
1802: }
1803: if (ksp->numbermonitors >= MAXKSPMONITORS) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Too many KSP monitors set");
1804: ksp->monitor[ksp->numbermonitors] = monitor;
1805: ksp->monitordestroy[ksp->numbermonitors] = monitordestroy;
1806: ksp->monitorcontext[ksp->numbermonitors++] = (void*)mctx;
1807: return(0);
1808: }
1810: /*@
1811: KSPMonitorCancel - Clears all monitors for a KSP object.
1813: Logically Collective on KSP
1815: Input Parameters:
1816: . ksp - iterative context obtained from KSPCreate()
1818: Options Database Key:
1819: . -ksp_monitor_cancel - Cancels all monitors that have
1820: been hardwired into a code by calls to KSPMonitorSet(),
1821: but does not cancel those set via the options database.
1823: Level: intermediate
1825: .keywords: set, monitor
1827: .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSPMonitorSet(), KSP
1828: @*/
1829: PetscErrorCode KSPMonitorCancel(KSP ksp)
1830: {
1832: PetscInt i;
1836: for (i=0; i<ksp->numbermonitors; i++) {
1837: if (ksp->monitordestroy[i]) {
1838: (*ksp->monitordestroy[i])(&ksp->monitorcontext[i]);
1839: }
1840: }
1841: ksp->numbermonitors = 0;
1842: return(0);
1843: }
1845: /*@C
1846: KSPGetMonitorContext - Gets the monitoring context, as set by
1847: KSPMonitorSet() for the FIRST monitor only.
1849: Not Collective
1851: Input Parameter:
1852: . ksp - iterative context obtained from KSPCreate()
1854: Output Parameter:
1855: . ctx - monitoring context
1857: Level: intermediate
1859: .keywords: get, monitor, context
1861: .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSP
1862: @*/
1863: PetscErrorCode KSPGetMonitorContext(KSP ksp,void **ctx)
1864: {
1867: *ctx = (ksp->monitorcontext[0]);
1868: return(0);
1869: }
1871: /*@
1872: KSPSetResidualHistory - Sets the array used to hold the residual history.
1873: If set, this array will contain the residual norms computed at each
1874: iteration of the solver.
1876: Not Collective
1878: Input Parameters:
1879: + ksp - iterative context obtained from KSPCreate()
1880: . a - array to hold history
1881: . na - size of a
1882: - reset - PETSC_TRUE indicates the history counter is reset to zero
1883: for each new linear solve
1885: Level: advanced
1887: Notes: The array is NOT freed by PETSc so the user needs to keep track of
1888: it and destroy once the KSP object is destroyed.
1890: If 'a' is NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
1891: default array of length 10000 is allocated.
1893: .keywords: set, residual, history, norm
1895: .seealso: KSPGetResidualHistory(), KSP
1897: @*/
1898: PetscErrorCode KSPSetResidualHistory(KSP ksp,PetscReal a[],PetscInt na,PetscBool reset)
1899: {
1905: PetscFree(ksp->res_hist_alloc);
1906: if (na != PETSC_DECIDE && na != PETSC_DEFAULT && a) {
1907: ksp->res_hist = a;
1908: ksp->res_hist_max = na;
1909: } else {
1910: if (na != PETSC_DECIDE && na != PETSC_DEFAULT) ksp->res_hist_max = na;
1911: else ksp->res_hist_max = 10000; /* like default ksp->max_it */
1912: PetscCalloc1(ksp->res_hist_max,&ksp->res_hist_alloc);
1914: ksp->res_hist = ksp->res_hist_alloc;
1915: }
1916: ksp->res_hist_len = 0;
1917: ksp->res_hist_reset = reset;
1918: return(0);
1919: }
1921: /*@C
1922: KSPGetResidualHistory - Gets the array used to hold the residual history
1923: and the number of residuals it contains.
1925: Not Collective
1927: Input Parameter:
1928: . ksp - iterative context obtained from KSPCreate()
1930: Output Parameters:
1931: + a - pointer to array to hold history (or NULL)
1932: - na - number of used entries in a (or NULL)
1934: Level: advanced
1936: Notes:
1937: Can only be called after a KSPSetResidualHistory() otherwise a and na are set to zero
1939: The Fortran version of this routine has a calling sequence
1940: $ call KSPGetResidualHistory(KSP ksp, integer na, integer ierr)
1941: note that you have passed a Fortran array into KSPSetResidualHistory() and you need
1942: to access the residual values from this Fortran array you provided. Only the na (number of
1943: residual norms currently held) is set.
1945: .keywords: get, residual, history, norm
1947: .seealso: KSPGetResidualHistory(), KSP
1949: @*/
1950: PetscErrorCode KSPGetResidualHistory(KSP ksp,PetscReal *a[],PetscInt *na)
1951: {
1954: if (a) *a = ksp->res_hist;
1955: if (na) *na = ksp->res_hist_len;
1956: return(0);
1957: }
1959: /*@C
1960: KSPSetConvergenceTest - Sets the function to be used to determine
1961: convergence.
1963: Logically Collective on KSP
1965: Input Parameters:
1966: + ksp - iterative context obtained from KSPCreate()
1967: . converge - pointer to int function
1968: . cctx - context for private data for the convergence routine (may be null)
1969: - destroy - a routine for destroying the context (may be null)
1971: Calling sequence of converge:
1972: $ converge (KSP ksp, int it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)
1974: + ksp - iterative context obtained from KSPCreate()
1975: . it - iteration number
1976: . rnorm - (estimated) 2-norm of (preconditioned) residual
1977: . reason - the reason why it has converged or diverged
1978: - cctx - optional convergence context, as set by KSPSetConvergenceTest()
1981: Notes:
1982: Must be called after the KSP type has been set so put this after
1983: a call to KSPSetType(), or KSPSetFromOptions().
1985: The default convergence test, KSPConvergedDefault(), aborts if the
1986: residual grows to more than 10000 times the initial residual.
1988: The default is a combination of relative and absolute tolerances.
1989: The residual value that is tested may be an approximation; routines
1990: that need exact values should compute them.
1992: In the default PETSc convergence test, the precise values of reason
1993: are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
1995: Level: advanced
1997: .keywords: set, convergence, test, context
1999: .seealso: KSPConvergedDefault(), KSPGetConvergenceContext(), KSPSetTolerances(), KSP
2000: @*/
2001: PetscErrorCode KSPSetConvergenceTest(KSP ksp,PetscErrorCode (*converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
2002: {
2007: if (ksp->convergeddestroy) {
2008: (*ksp->convergeddestroy)(ksp->cnvP);
2009: }
2010: ksp->converged = converge;
2011: ksp->convergeddestroy = destroy;
2012: ksp->cnvP = (void*)cctx;
2013: return(0);
2014: }
2016: /*@C
2017: KSPGetConvergenceContext - Gets the convergence context set with
2018: KSPSetConvergenceTest().
2020: Not Collective
2022: Input Parameter:
2023: . ksp - iterative context obtained from KSPCreate()
2025: Output Parameter:
2026: . ctx - monitoring context
2028: Level: advanced
2030: .keywords: get, convergence, test, context
2032: .seealso: KSPConvergedDefault(), KSPSetConvergenceTest(), KSP
2033: @*/
2034: PetscErrorCode KSPGetConvergenceContext(KSP ksp,void **ctx)
2035: {
2038: *ctx = ksp->cnvP;
2039: return(0);
2040: }
2042: /*@C
2043: KSPBuildSolution - Builds the approximate solution in a vector provided.
2044: This routine is NOT commonly needed (see KSPSolve()).
2046: Collective on KSP
2048: Input Parameter:
2049: . ctx - iterative context obtained from KSPCreate()
2051: Output Parameter:
2052: Provide exactly one of
2053: + v - location to stash solution.
2054: - V - the solution is returned in this location. This vector is created
2055: internally. This vector should NOT be destroyed by the user with
2056: VecDestroy().
2058: Notes:
2059: This routine can be used in one of two ways
2060: .vb
2061: KSPBuildSolution(ksp,NULL,&V);
2062: or
2063: KSPBuildSolution(ksp,v,NULL); or KSPBuildSolution(ksp,v,&v);
2064: .ve
2065: In the first case an internal vector is allocated to store the solution
2066: (the user cannot destroy this vector). In the second case the solution
2067: is generated in the vector that the user provides. Note that for certain
2068: methods, such as KSPCG, the second case requires a copy of the solution,
2069: while in the first case the call is essentially free since it simply
2070: returns the vector where the solution already is stored. For some methods
2071: like GMRES this is a reasonably expensive operation and should only be
2072: used in truly needed.
2074: Level: advanced
2076: .keywords: build, solution
2078: .seealso: KSPGetSolution(), KSPBuildResidual(), KSP
2079: @*/
2080: PetscErrorCode KSPBuildSolution(KSP ksp,Vec v,Vec *V)
2081: {
2086: if (!V && !v) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONG,"Must provide either v or V");
2087: if (!V) V = &v;
2088: (*ksp->ops->buildsolution)(ksp,v,V);
2089: return(0);
2090: }
2092: /*@C
2093: KSPBuildResidual - Builds the residual in a vector provided.
2095: Collective on KSP
2097: Input Parameter:
2098: . ksp - iterative context obtained from KSPCreate()
2100: Output Parameters:
2101: + v - optional location to stash residual. If v is not provided,
2102: then a location is generated.
2103: . t - work vector. If not provided then one is generated.
2104: - V - the residual
2106: Notes:
2107: Regardless of whether or not v is provided, the residual is
2108: returned in V.
2110: Level: advanced
2112: .keywords: KSP, build, residual
2114: .seealso: KSPBuildSolution()
2115: @*/
2116: PetscErrorCode KSPBuildResidual(KSP ksp,Vec t,Vec v,Vec *V)
2117: {
2119: PetscBool flag = PETSC_FALSE;
2120: Vec w = v,tt = t;
2124: if (!w) {
2125: VecDuplicate(ksp->vec_rhs,&w);
2126: PetscLogObjectParent((PetscObject)ksp,(PetscObject)w);
2127: }
2128: if (!tt) {
2129: VecDuplicate(ksp->vec_sol,&tt); flag = PETSC_TRUE;
2130: PetscLogObjectParent((PetscObject)ksp,(PetscObject)tt);
2131: }
2132: (*ksp->ops->buildresidual)(ksp,tt,w,V);
2133: if (flag) {VecDestroy(&tt);}
2134: return(0);
2135: }
2137: /*@
2138: KSPSetDiagonalScale - Tells KSP to symmetrically diagonally scale the system
2139: before solving. This actually CHANGES the matrix (and right hand side).
2141: Logically Collective on KSP
2143: Input Parameter:
2144: + ksp - the KSP context
2145: - scale - PETSC_TRUE or PETSC_FALSE
2147: Options Database Key:
2148: + -ksp_diagonal_scale -
2149: - -ksp_diagonal_scale_fix - scale the matrix back AFTER the solve
2152: Notes: Scales the matrix by D^(-1/2) A D^(-1/2) [D^(1/2) x ] = D^(-1/2) b
2153: where D_{ii} is 1/abs(A_{ii}) unless A_{ii} is zero and then it is 1.
2155: BE CAREFUL with this routine: it actually scales the matrix and right
2156: hand side that define the system. After the system is solved the matrix
2157: and right hand side remain scaled unless you use KSPSetDiagonalScaleFix()
2159: This should NOT be used within the SNES solves if you are using a line
2160: search.
2162: If you use this with the PCType Eisenstat preconditioner than you can
2163: use the PCEisenstatSetNoDiagonalScaling() option, or -pc_eisenstat_no_diagonal_scaling
2164: to save some unneeded, redundant flops.
2166: Level: intermediate
2168: .keywords: set, options, prefix, database
2170: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScaleFix(), KSP
2171: @*/
2172: PetscErrorCode KSPSetDiagonalScale(KSP ksp,PetscBool scale)
2173: {
2177: ksp->dscale = scale;
2178: return(0);
2179: }
2181: /*@
2182: KSPGetDiagonalScale - Checks if KSP solver scales the matrix and
2183: right hand side
2185: Not Collective
2187: Input Parameter:
2188: . ksp - the KSP context
2190: Output Parameter:
2191: . scale - PETSC_TRUE or PETSC_FALSE
2193: Notes:
2194: BE CAREFUL with this routine: it actually scales the matrix and right
2195: hand side that define the system. After the system is solved the matrix
2196: and right hand side remain scaled unless you use KSPSetDiagonalScaleFix()
2198: Level: intermediate
2200: .keywords: set, options, prefix, database
2202: .seealso: KSPSetDiagonalScale(), KSPSetDiagonalScaleFix(), KSP
2203: @*/
2204: PetscErrorCode KSPGetDiagonalScale(KSP ksp,PetscBool *scale)
2205: {
2209: *scale = ksp->dscale;
2210: return(0);
2211: }
2213: /*@
2214: KSPSetDiagonalScaleFix - Tells KSP to diagonally scale the system
2215: back after solving.
2217: Logically Collective on KSP
2219: Input Parameter:
2220: + ksp - the KSP context
2221: - fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
2222: rescale (default)
2224: Notes:
2225: Must be called after KSPSetDiagonalScale()
2227: Using this will slow things down, because it rescales the matrix before and
2228: after each linear solve. This is intended mainly for testing to allow one
2229: to easily get back the original system to make sure the solution computed is
2230: accurate enough.
2232: Level: intermediate
2234: .keywords: set, options, prefix, database
2236: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPGetDiagonalScaleFix(), KSP
2237: @*/
2238: PetscErrorCode KSPSetDiagonalScaleFix(KSP ksp,PetscBool fix)
2239: {
2243: ksp->dscalefix = fix;
2244: return(0);
2245: }
2247: /*@
2248: KSPGetDiagonalScaleFix - Determines if KSP diagonally scales the system
2249: back after solving.
2251: Not Collective
2253: Input Parameter:
2254: . ksp - the KSP context
2256: Output Parameter:
2257: . fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
2258: rescale (default)
2260: Notes:
2261: Must be called after KSPSetDiagonalScale()
2263: If PETSC_TRUE will slow things down, because it rescales the matrix before and
2264: after each linear solve. This is intended mainly for testing to allow one
2265: to easily get back the original system to make sure the solution computed is
2266: accurate enough.
2268: Level: intermediate
2270: .keywords: set, options, prefix, database
2272: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPSetDiagonalScaleFix(), KSP
2273: @*/
2274: PetscErrorCode KSPGetDiagonalScaleFix(KSP ksp,PetscBool *fix)
2275: {
2279: *fix = ksp->dscalefix;
2280: return(0);
2281: }
2283: /*@C
2284: KSPSetComputeOperators - set routine to compute the linear operators
2286: Logically Collective
2288: Input Arguments:
2289: + ksp - the KSP context
2290: . func - function to compute the operators
2291: - ctx - optional context
2293: Calling sequence of func:
2294: $ func(KSP ksp,Mat A,Mat B,void *ctx)
2296: + ksp - the KSP context
2297: . A - the linear operator
2298: . B - preconditioning matrix
2299: - ctx - optional user-provided context
2301: Notes: The user provided func() will be called automatically at the very next call to KSPSolve(). It will not be called at future KSPSolve() calls
2302: unless either KSPSetComputeOperators() or KSPSetOperators() is called before that KSPSolve() is called.
2304: To reuse the same preconditioner for the next KSPSolve() and not compute a new one based on the most recently computed matrix call KSPSetReusePreconditioner()
2306: Level: beginner
2308: .seealso: KSPSetOperators(), KSPSetComputeRHS(), DMKSPSetComputeOperators(), KSPSetComputeInitialGuess()
2309: @*/
2310: PetscErrorCode KSPSetComputeOperators(KSP ksp,PetscErrorCode (*func)(KSP,Mat,Mat,void*),void *ctx)
2311: {
2313: DM dm;
2317: KSPGetDM(ksp,&dm);
2318: DMKSPSetComputeOperators(dm,func,ctx);
2319: if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX;
2320: return(0);
2321: }
2323: /*@C
2324: KSPSetComputeRHS - set routine to compute the right hand side of the linear system
2326: Logically Collective
2328: Input Arguments:
2329: + ksp - the KSP context
2330: . func - function to compute the right hand side
2331: - ctx - optional context
2333: Calling sequence of func:
2334: $ func(KSP ksp,Vec b,void *ctx)
2336: + ksp - the KSP context
2337: . b - right hand side of linear system
2338: - ctx - optional user-provided context
2340: Notes: The routine you provide will be called EACH you call KSPSolve() to prepare the new right hand side for that solve
2342: Level: beginner
2344: .seealso: KSPSolve(), DMKSPSetComputeRHS(), KSPSetComputeOperators()
2345: @*/
2346: PetscErrorCode KSPSetComputeRHS(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2347: {
2349: DM dm;
2353: KSPGetDM(ksp,&dm);
2354: DMKSPSetComputeRHS(dm,func,ctx);
2355: return(0);
2356: }
2358: /*@C
2359: KSPSetComputeInitialGuess - set routine to compute the initial guess of the linear system
2361: Logically Collective
2363: Input Arguments:
2364: + ksp - the KSP context
2365: . func - function to compute the initial guess
2366: - ctx - optional context
2368: Calling sequence of func:
2369: $ func(KSP ksp,Vec x,void *ctx)
2371: + ksp - the KSP context
2372: . x - solution vector
2373: - ctx - optional user-provided context
2375: Level: beginner
2377: .seealso: KSPSolve(), KSPSetComputeRHS(), KSPSetComputeOperators(), DMKSPSetComputeInitialGuess()
2378: @*/
2379: PetscErrorCode KSPSetComputeInitialGuess(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2380: {
2382: DM dm;
2386: KSPGetDM(ksp,&dm);
2387: DMKSPSetComputeInitialGuess(dm,func,ctx);
2388: return(0);
2389: }