Actual source code: itfunc.c
petsc-3.3-p7 2013-05-11
2: /*
3: Interface KSP routines that the user calls.
4: */
6: #include <petsc-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 ig = PETSC_FALSE;
179: Mat A,B;
180: MatStructure stflg = SAME_NONZERO_PATTERN;
185: /* reset the convergence flag from the previous solves */
186: ksp->reason = KSP_CONVERGED_ITERATING;
188: if (!((PetscObject)ksp)->type_name){
189: KSPSetType(ksp,KSPGMRES);
190: }
191: KSPSetUpNorms_Private(ksp,&ksp->normtype,&ksp->pc_side);
193: if (ksp->dmActive && !ksp->setupstage) {
194: /* first time in so build matrix and vector data structures using DM */
195: if (!ksp->vec_rhs) {DMCreateGlobalVector(ksp->dm,&ksp->vec_rhs);}
196: if (!ksp->vec_sol) {DMCreateGlobalVector(ksp->dm,&ksp->vec_sol);}
197: DMCreateMatrix(ksp->dm,MATAIJ,&A);
198: KSPSetOperators(ksp,A,A,stflg);
199: PetscObjectDereference((PetscObject)A);
200: }
202: if (ksp->dmActive) {
203: KSPDM kdm;
204: DMKSPGetContext(ksp->dm,&kdm);
206: DMHasInitialGuess(ksp->dm,&ig);
207: if (ig && ksp->setupstage != KSP_SETUP_NEWRHS) {
208: /* only computes initial guess the first time through */
209: DMComputeInitialGuess(ksp->dm,ksp->vec_sol);
210: KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
211: }
212: if (kdm->computerhs) {
213: (*kdm->computerhs)(ksp,ksp->vec_rhs,kdm->rhsctx);
214: } else {
215: PetscBool irhs;
216: DMHasFunction(ksp->dm,&irhs);
217: if (irhs) {
218: DMComputeFunction(ksp->dm,PETSC_NULL,ksp->vec_rhs);
219: }
220: }
222: if (ksp->setupstage != KSP_SETUP_NEWRHS) {
223: KSPGetOperators(ksp,&A,&B,PETSC_NULL);
224: if (kdm->computeoperators) {
225: (*kdm->computeoperators)(ksp,A,B,&stflg,kdm->operatorsctx);
226: } else { /* Eventually remove this code path */
227: if (0) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_USER,"Must call KSPSetComputeOperators()");
228: DMComputeJacobian(ksp->dm,PETSC_NULL,A,B,&stflg);
229: }
230: KSPSetOperators(ksp,A,B,stflg);
231: }
232: }
234: if (ksp->setupstage == KSP_SETUP_NEWRHS) return(0);
235: PetscLogEventBegin(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);
237: switch (ksp->setupstage) {
238: case KSP_SETUP_NEW:
239: (*ksp->ops->setup)(ksp);
240: break;
241: case KSP_SETUP_NEWMATRIX: { /* This should be replaced with a more general mechanism */
242: KSPChebyshevSetNewMatrix(ksp);
243: } break;
244: default: break;
245: }
247: /* scale the matrix if requested */
248: if (ksp->dscale) {
249: Mat mat,pmat;
250: PetscScalar *xx;
251: PetscInt i,n;
252: PetscBool zeroflag = PETSC_FALSE;
253: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
254: PCGetOperators(ksp->pc,&mat,&pmat,PETSC_NULL);
255: if (!ksp->diagonal) { /* allocate vector to hold diagonal */
256: MatGetVecs(pmat,&ksp->diagonal,0);
257: }
258: MatGetDiagonal(pmat,ksp->diagonal);
259: VecGetLocalSize(ksp->diagonal,&n);
260: VecGetArray(ksp->diagonal,&xx);
261: for (i=0; i<n; i++) {
262: if (xx[i] != 0.0) xx[i] = 1.0/PetscSqrtReal(PetscAbsScalar(xx[i]));
263: else {
264: xx[i] = 1.0;
265: zeroflag = PETSC_TRUE;
266: }
267: }
268: VecRestoreArray(ksp->diagonal,&xx);
269: if (zeroflag) {
270: PetscInfo(ksp,"Zero detected in diagonal of matrix, using 1 at those locations\n");
271: }
272: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
273: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
274: ksp->dscalefix2 = PETSC_FALSE;
275: }
276: PetscLogEventEnd(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);
277: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
278: PCSetUp(ksp->pc);
279: if (ksp->nullsp) {
280: PetscBool test = PETSC_FALSE;
281: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_test_null_space",&test,PETSC_NULL);
282: if (test) {
283: Mat mat;
284: PCGetOperators(ksp->pc,&mat,PETSC_NULL,PETSC_NULL);
285: MatNullSpaceTest(ksp->nullsp,mat,PETSC_NULL);
286: }
287: }
288: ksp->setupstage = KSP_SETUP_NEWRHS;
289: return(0);
290: }
294: /*@
295: KSPSolve - Solves linear system.
297: Collective on KSP
299: Parameter:
300: + ksp - iterative context obtained from KSPCreate()
301: . b - the right hand side vector
302: - x - the solution (this may be the same vector as b, then b will be overwritten with answer)
304: Options Database Keys:
305: + -ksp_compute_eigenvalues - compute preconditioned operators eigenvalues
306: . -ksp_plot_eigenvalues - plot the computed eigenvalues in an X-window
307: . -ksp_compute_eigenvalues_explicitly - compute the eigenvalues by forming the dense operator and useing LAPACK
308: . -ksp_plot_eigenvalues_explicitly - plot the explicitly computing eigenvalues
309: . -ksp_view_binary - save matrix and right hand side that define linear system to the default binary viewer (can be
310: read later with src/ksp/examples/tutorials/ex10.c for testing solvers)
311: . -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations
312: . -ksp_final_residual - print 2-norm of true linear system residual at the end of the solution process
313: - -ksp_view - print the ksp data structure at the end of the system solution
315: Notes:
317: If one uses KSPSetDM() then x or b need not be passed. Use KSPGetSolution() to access the solution in this case.
319: The operator is specified with KSPSetOperators().
321: Call KSPGetConvergedReason() to determine if the solver converged or failed and
322: why. The number of iterations can be obtained from KSPGetIterationNumber().
323:
324: If using a direct method (e.g., via the KSP solver
325: KSPPREONLY and a preconditioner such as PCLU/PCILU),
326: then its=1. See KSPSetTolerances() and KSPDefaultConverged()
327: for more details.
329: Understanding Convergence:
330: The routines KSPMonitorSet(), KSPComputeEigenvalues(), and
331: KSPComputeEigenvaluesExplicitly() provide information on additional
332: options to monitor convergence and print eigenvalue information.
334: Level: beginner
336: .keywords: KSP, solve, linear system
338: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPDefaultConverged(),
339: KSPSolveTranspose(), KSPGetIterationNumber()
340: @*/
341: PetscErrorCode KSPSolve(KSP ksp,Vec b,Vec x)
342: {
344: PetscMPIInt rank;
345: PetscBool flag1,flag2,flg = PETSC_FALSE,inXisinB=PETSC_FALSE,guess_zero;
346: char view[10];
347: char filename[PETSC_MAX_PATH_LEN];
348: PetscViewer viewer;
349:
356: if (x && x == b) {
357: if (!ksp->guess_zero) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_INCOMP,"Cannot use x == b with nonzero initial guess");
358: VecDuplicate(b,&x);
359: inXisinB = PETSC_TRUE;
360: }
361: if (b) {
362: PetscObjectReference((PetscObject)b);
363: VecDestroy(&ksp->vec_rhs);
364: ksp->vec_rhs = b;
365: }
366: if (x) {
367: PetscObjectReference((PetscObject)x);
368: VecDestroy(&ksp->vec_sol);
369: ksp->vec_sol = x;
370: }
372: flag1 = flag2 = PETSC_FALSE;
373: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_view_binary",&flag1,PETSC_NULL);
374: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_view_binary_pre",&flag2,PETSC_NULL);
375: if (flag1 || flag2) {
376: Mat mat,premat;
377: PetscViewer viewer = PETSC_VIEWER_BINARY_(((PetscObject)ksp)->comm);
378: PCGetOperators(ksp->pc,&mat,&premat,PETSC_NULL);
379: if (flag1) {MatView(mat,viewer);}
380: if (flag2) {MatView(premat,viewer);}
381: VecView(ksp->vec_rhs,viewer);
382: }
383: PetscLogEventBegin(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
385: /* reset the residual history list if requested */
386: if (ksp->res_hist_reset) ksp->res_hist_len = 0;
388: PetscOptionsGetString(((PetscObject)ksp)->prefix,"-ksp_view_before",view,10,&flg);
389: if (flg) {
390: PetscViewer viewer;
391: PetscViewerASCIIGetStdout(((PetscObject)ksp)->comm,&viewer);
392: KSPView(ksp,viewer);
393: }
395: ksp->transpose_solve = PETSC_FALSE;
397: if (ksp->guess) {
398: KSPFischerGuessFormGuess(ksp->guess,ksp->vec_rhs,ksp->vec_sol);
399: ksp->guess_zero = PETSC_FALSE;
400: }
401: /* KSPSetUp() scales the matrix if needed */
402: KSPSetUp(ksp);
403: KSPSetUpOnBlocks(ksp);
405: /* diagonal scale RHS if called for */
406: if (ksp->dscale) {
407: VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
408: /* second time in, but matrix was scaled back to original */
409: if (ksp->dscalefix && ksp->dscalefix2) {
410: Mat mat,pmat;
412: PCGetOperators(ksp->pc,&mat,&pmat,PETSC_NULL);
413: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
414: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
415: }
417: /* scale initial guess */
418: if (!ksp->guess_zero) {
419: if (!ksp->truediagonal) {
420: VecDuplicate(ksp->diagonal,&ksp->truediagonal);
421: VecCopy(ksp->diagonal,ksp->truediagonal);
422: VecReciprocal(ksp->truediagonal);
423: }
424: VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->truediagonal);
425: }
426: }
427: PCPreSolve(ksp->pc,ksp);
429: if (ksp->guess_zero) { VecSet(ksp->vec_sol,0.0);}
430: if (ksp->guess_knoll) {
431: PCApply(ksp->pc,ksp->vec_rhs,ksp->vec_sol);
432: KSP_RemoveNullSpace(ksp,ksp->vec_sol);
433: ksp->guess_zero = PETSC_FALSE;
434: }
436: /* can we mark the initial guess as zero for this solve? */
437: guess_zero = ksp->guess_zero;
438: if (!ksp->guess_zero) {
439: PetscReal norm;
441: VecNormAvailable(ksp->vec_sol,NORM_2,&flg,&norm);
442: if (flg && !norm) {
443: ksp->guess_zero = PETSC_TRUE;
444: }
445: }
446: (*ksp->ops->solve)(ksp);
447: ksp->guess_zero = guess_zero;
449: if (!ksp->reason) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
450: if (ksp->printreason) {
451: PetscViewerASCIIAddTab(PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm),((PetscObject)ksp)->tablevel);
452: if (ksp->reason > 0) {
453: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm),"Linear solve converged due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
454: } else {
455: PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm),"Linear solve did not converge due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
456: }
457: PetscViewerASCIISubtractTab(PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm),((PetscObject)ksp)->tablevel);
458: }
459: PCPostSolve(ksp->pc,ksp);
461: /* diagonal scale solution if called for */
462: if (ksp->dscale) {
463: VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->diagonal);
464: /* unscale right hand side and matrix */
465: if (ksp->dscalefix) {
466: Mat mat,pmat;
468: VecReciprocal(ksp->diagonal);
469: VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
470: PCGetOperators(ksp->pc,&mat,&pmat,PETSC_NULL);
471: MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
472: if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
473: VecReciprocal(ksp->diagonal);
474: ksp->dscalefix2 = PETSC_TRUE;
475: }
476: }
477: PetscLogEventEnd(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
479: if (ksp->guess) {
480: KSPFischerGuessUpdate(ksp->guess,ksp->vec_sol);
481: }
483: MPI_Comm_rank(((PetscObject)ksp)->comm,&rank);
485: flag1 = PETSC_FALSE;
486: flag2 = PETSC_FALSE;
487: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_compute_eigenvalues",&flag1,PETSC_NULL);
488: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_plot_eigenvalues",&flag2,PETSC_NULL);
489: if (flag1 || flag2) {
490: PetscInt nits,n,i,neig;
491: PetscReal *r,*c;
492:
493: KSPGetIterationNumber(ksp,&nits);
494: n = nits+2;
496: if (!nits) {
497: PetscPrintf(((PetscObject)ksp)->comm,"Zero iterations in solver, cannot approximate any eigenvalues\n");
498: } else {
499: PetscMalloc(2*n*sizeof(PetscReal),&r);
500: c = r + n;
501: KSPComputeEigenvalues(ksp,n,r,c,&neig);
502: if (flag1) {
503: PetscPrintf(((PetscObject)ksp)->comm,"Iteratively computed eigenvalues\n");
504: for (i=0; i<neig; i++) {
505: if (c[i] >= 0.0) {PetscPrintf(((PetscObject)ksp)->comm,"%G + %Gi\n",r[i],c[i]);}
506: else {PetscPrintf(((PetscObject)ksp)->comm,"%G - %Gi\n",r[i],-c[i]);}
507: }
508: }
509: if (flag2 && !rank) {
510: PetscDraw draw;
511: PetscDrawSP drawsp;
513: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Iteratively Computed Eigenvalues",PETSC_DECIDE,PETSC_DECIDE,300,300,&viewer);
514: PetscViewerDrawGetDraw(viewer,0,&draw);
515: PetscDrawSPCreate(draw,1,&drawsp);
516: for (i=0; i<neig; i++) {
517: PetscDrawSPAddPoint(drawsp,r+i,c+i);
518: }
519: PetscDrawSPDraw(drawsp);
520: PetscDrawSPDestroy(&drawsp);
521: PetscViewerDestroy(&viewer);
522: }
523: PetscFree(r);
524: }
525: }
527: flag1 = PETSC_FALSE;
528: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_compute_singularvalues",&flag1,PETSC_NULL);
529: if (flag1) {
530: PetscInt nits;
532: KSPGetIterationNumber(ksp,&nits);
533: if (!nits) {
534: PetscPrintf(((PetscObject)ksp)->comm,"Zero iterations in solver, cannot approximate any singular values\n");
535: } else {
536: PetscReal emax,emin;
538: KSPComputeExtremeSingularValues(ksp,&emax,&emin);
539: PetscPrintf(((PetscObject)ksp)->comm,"Iteratively computed extreme singular values: max %G min %G max/min %G\n",emax,emin,emax/emin);
540: }
541: }
544: flag1 = PETSC_FALSE;
545: flag2 = PETSC_FALSE;
546: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_compute_eigenvalues_explicitly",&flag1,PETSC_NULL);
547: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_plot_eigenvalues_explicitly",&flag2,PETSC_NULL);
548: if (flag1 || flag2) {
549: PetscInt n,i;
550: PetscReal *r,*c;
551: VecGetSize(ksp->vec_sol,&n);
552: PetscMalloc2(n,PetscReal,&r,n,PetscReal,&c);
553: KSPComputeEigenvaluesExplicitly(ksp,n,r,c);
554: if (flag1) {
555: PetscPrintf(((PetscObject)ksp)->comm,"Explicitly computed eigenvalues\n");
556: for (i=0; i<n; i++) {
557: if (c[i] >= 0.0) {PetscPrintf(((PetscObject)ksp)->comm,"%G + %Gi\n",r[i],c[i]);}
558: else {PetscPrintf(((PetscObject)ksp)->comm,"%G - %Gi\n",r[i],-c[i]);}
559: }
560: }
561: if (flag2 && !rank) {
562: PetscDraw draw;
563: PetscDrawSP drawsp;
565: PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Explicitly Computed Eigenvalues",0,320,300,300,&viewer);
566: PetscViewerDrawGetDraw(viewer,0,&draw);
567: PetscDrawSPCreate(draw,1,&drawsp);
568: for (i=0; i<n; i++) {
569: PetscDrawSPAddPoint(drawsp,r+i,c+i);
570: }
571: PetscDrawSPDraw(drawsp);
572: PetscDrawSPDestroy(&drawsp);
573: PetscViewerDestroy(&viewer);
574: }
575: PetscFree2(r,c);
576: }
578: flag2 = PETSC_FALSE;
579: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_view_operator",&flag2,PETSC_NULL);
580: if (flag2) {
581: Mat A,B;
582: PetscViewer viewer;
584: PCGetOperators(ksp->pc,&A,PETSC_NULL,PETSC_NULL);
585: MatComputeExplicitOperator(A,&B);
586: PetscViewerASCIIGetStdout(((PetscObject)ksp)->comm,&viewer);
587: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
588: MatView(B,viewer);
589: PetscViewerPopFormat(viewer);
590: MatDestroy(&B);
591: }
592: flag2 = PETSC_FALSE;
593: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_view_operator_binary",&flag2,PETSC_NULL);
594: if (flag2) {
595: Mat A,B;
596: PCGetOperators(ksp->pc,&A,PETSC_NULL,PETSC_NULL);
597: MatComputeExplicitOperator(A,&B);
598: MatView(B,PETSC_VIEWER_BINARY_(((PetscObject)ksp)->comm));
599: MatDestroy(&B);
600: }
601: flag2 = PETSC_FALSE;
602: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_view_preconditioned_operator_binary",&flag2,PETSC_NULL);
603: if (flag2) {
604: Mat B;
605: KSPComputeExplicitOperator(ksp,&B);
606: MatView(B,PETSC_VIEWER_BINARY_(((PetscObject)ksp)->comm));
607: MatDestroy(&B);
608: }
609: flag2 = PETSC_FALSE;
610: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_view_preconditioner_binary",&flag2,PETSC_NULL);
611: if (flag2) {
612: Mat B;
613: PCComputeExplicitOperator(ksp->pc,&B);
614: MatView(B,PETSC_VIEWER_BINARY_(((PetscObject)ksp)->comm));
615: MatDestroy(&B);
616: }
617: PetscOptionsGetString(((PetscObject)ksp)->prefix,"-ksp_view",filename,PETSC_MAX_PATH_LEN,&flg);
618: if (flg && !PetscPreLoadingOn) {
619: PetscViewerASCIIOpen(((PetscObject)ksp)->comm,filename,&viewer);
620: KSPView(ksp,viewer);
621: PetscViewerDestroy(&viewer);
622: }
623: flg = PETSC_FALSE;
624: PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_final_residual",&flg,PETSC_NULL);
625: if (flg) {
626: Mat A;
627: Vec t;
628: PetscReal norm;
629: 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");
630: PCGetOperators(ksp->pc,&A,0,0);
631: VecDuplicate(ksp->vec_rhs,&t);
632: KSP_MatMult(ksp,A,ksp->vec_sol,t);
633: VecAYPX(t, -1.0, ksp->vec_rhs);
634: VecNorm(t,NORM_2,&norm);
635: VecDestroy(&t);
636: PetscPrintf(((PetscObject)ksp)->comm,"KSP final norm of residual %G\n",norm);
637: }
638: if (inXisinB) {
639: VecCopy(x,b);
640: VecDestroy(&x);
641: }
642: if (ksp->errorifnotconverged && ksp->reason < 0) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged");
643: return(0);
644: }
648: /*@
649: KSPSolveTranspose - Solves the transpose of a linear system.
651: Collective on KSP
653: Input Parameter:
654: + ksp - iterative context obtained from KSPCreate()
655: . b - right hand side vector
656: - x - solution vector
658: Notes: For complex numbers this solve the non-Hermitian transpose system.
660: Developer Notes: We need to implement a KSPSolveHermitianTranspose()
662: Level: developer
664: .keywords: KSP, solve, linear system
666: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPDefaultConverged(),
667: KSPSolve()
668: @*/
670: PetscErrorCode KSPSolveTranspose(KSP ksp,Vec b,Vec x)
671: {
673: PetscBool inXisinB=PETSC_FALSE;
679: if (x == b) {
680: VecDuplicate(b,&x);
681: inXisinB = PETSC_TRUE;
682: }
683: PetscObjectReference((PetscObject)b);
684: PetscObjectReference((PetscObject)x);
685: VecDestroy(&ksp->vec_rhs);
686: VecDestroy(&ksp->vec_sol);
687: ksp->vec_rhs = b;
688: ksp->vec_sol = x;
689: ksp->transpose_solve = PETSC_TRUE;
690: KSPSetUp(ksp);
691: if (ksp->guess_zero) { VecSet(ksp->vec_sol,0.0);}
692: (*ksp->ops->solve)(ksp);
693: if (!ksp->reason) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
694: if (ksp->printreason) {
695: if (ksp->reason > 0) {
696: PetscPrintf(((PetscObject)ksp)->comm,"Linear solve converged due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
697: } else {
698: PetscPrintf(((PetscObject)ksp)->comm,"Linear solve did not converge due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
699: }
700: }
701: if (inXisinB) {
702: VecCopy(x,b);
703: VecDestroy(&x);
704: }
705: if (ksp->errorifnotconverged && ksp->reason < 0) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged");
706: return(0);
707: }
711: /*@
712: KSPReset - Resets a KSP context to the kspsetupcalled = 0 state and removes any allocated Vecs and Mats
714: Collective on KSP
716: Input Parameter:
717: . ksp - iterative context obtained from KSPCreate()
719: Level: beginner
721: .keywords: KSP, destroy
723: .seealso: KSPCreate(), KSPSetUp(), KSPSolve()
724: @*/
725: PetscErrorCode KSPReset(KSP ksp)
726: {
731: if (!ksp) return(0);
732: if (ksp->ops->reset) {
733: (*ksp->ops->reset)(ksp);
734: }
735: if (ksp->pc) {PCReset(ksp->pc);}
736: KSPFischerGuessDestroy(&ksp->guess);
737: VecDestroyVecs(ksp->nwork,&ksp->work);
738: VecDestroy(&ksp->vec_rhs);
739: VecDestroy(&ksp->vec_sol);
740: VecDestroy(&ksp->diagonal);
741: VecDestroy(&ksp->truediagonal);
742: MatNullSpaceDestroy(&ksp->nullsp);
743: ksp->setupstage = KSP_SETUP_NEW;
744: return(0);
745: }
749: /*@
750: KSPDestroy - Destroys KSP context.
752: Collective on KSP
754: Input Parameter:
755: . ksp - iterative context obtained from KSPCreate()
757: Level: beginner
759: .keywords: KSP, destroy
761: .seealso: KSPCreate(), KSPSetUp(), KSPSolve()
762: @*/
763: PetscErrorCode KSPDestroy(KSP *ksp)
764: {
766: PC pc;
769: if (!*ksp) return(0);
771: if (--((PetscObject)(*ksp))->refct > 0) {*ksp = 0; return(0);}
772:
773: /*
774: Avoid a cascading call to PCReset(ksp->pc) from the following call:
775: PCReset() shouldn't be called from KSPDestroy() as it is unprotected by pc's
776: refcount (and may be shared, e.g., by other ksps).
777: */
778: pc = (*ksp)->pc;
779: (*ksp)->pc = PETSC_NULL;
780: KSPReset((*ksp));
781: (*ksp)->pc = pc;
782: PetscObjectDepublish((*ksp));
783: if ((*ksp)->ops->destroy) {(*(*ksp)->ops->destroy)(*ksp);}
785: DMDestroy(&(*ksp)->dm);
786: PCDestroy(&(*ksp)->pc);
787: PetscFree((*ksp)->res_hist_alloc);
788: if ((*ksp)->convergeddestroy) {
789: (*(*ksp)->convergeddestroy)((*ksp)->cnvP);
790: }
791: KSPMonitorCancel((*ksp));
792: PetscHeaderDestroy(ksp);
793: return(0);
794: }
798: /*@
799: KSPSetPCSide - Sets the preconditioning side.
801: Logically Collective on KSP
803: Input Parameter:
804: . ksp - iterative context obtained from KSPCreate()
806: Output Parameter:
807: . side - the preconditioning side, where side is one of
808: .vb
809: PC_LEFT - left preconditioning (default)
810: PC_RIGHT - right preconditioning
811: PC_SYMMETRIC - symmetric preconditioning
812: .ve
814: Options Database Keys:
815: . -ksp_pc_side <right,left,symmetric>
817: Notes:
818: Left preconditioning is used by default for most Krylov methods except KSPFGMRES which only supports right preconditioning.
819: Symmetric preconditioning is currently available only for the KSPQCG method. Note, however, that
820: symmetric preconditioning can be emulated by using either right or left
821: preconditioning and a pre or post processing step.
823: Level: intermediate
825: .keywords: KSP, set, right, left, symmetric, side, preconditioner, flag
827: .seealso: KSPGetPCSide()
828: @*/
829: PetscErrorCode KSPSetPCSide(KSP ksp,PCSide side)
830: {
834: ksp->pc_side = side;
835: return(0);
836: }
840: /*@
841: KSPGetPCSide - Gets the preconditioning side.
843: Not Collective
845: Input Parameter:
846: . ksp - iterative context obtained from KSPCreate()
848: Output Parameter:
849: . side - the preconditioning side, where side is one of
850: .vb
851: PC_LEFT - left preconditioning (default)
852: PC_RIGHT - right preconditioning
853: PC_SYMMETRIC - symmetric preconditioning
854: .ve
856: Level: intermediate
858: .keywords: KSP, get, right, left, symmetric, side, preconditioner, flag
860: .seealso: KSPSetPCSide()
861: @*/
862: PetscErrorCode KSPGetPCSide(KSP ksp,PCSide *side)
863: {
869: KSPSetUpNorms_Private(ksp,&ksp->normtype,&ksp->pc_side);
870: *side = ksp->pc_side;
871: return(0);
872: }
876: /*@
877: KSPGetTolerances - Gets the relative, absolute, divergence, and maximum
878: iteration tolerances used by the default KSP convergence tests.
880: Not Collective
882: Input Parameter:
883: . ksp - the Krylov subspace context
884:
885: Output Parameters:
886: + rtol - the relative convergence tolerance
887: . abstol - the absolute convergence tolerance
888: . dtol - the divergence tolerance
889: - maxits - maximum number of iterations
891: Notes:
892: The user can specify PETSC_NULL for any parameter that is not needed.
894: Level: intermediate
896: .keywords: KSP, get, tolerance, absolute, relative, divergence, convergence,
897: maximum, iterations
899: .seealso: KSPSetTolerances()
900: @*/
901: PetscErrorCode KSPGetTolerances(KSP ksp,PetscReal *rtol,PetscReal *abstol,PetscReal *dtol,PetscInt *maxits)
902: {
905: if (abstol) *abstol = ksp->abstol;
906: if (rtol) *rtol = ksp->rtol;
907: if (dtol) *dtol = ksp->divtol;
908: if (maxits) *maxits = ksp->max_it;
909: return(0);
910: }
914: /*@
915: KSPSetTolerances - Sets the relative, absolute, divergence, and maximum
916: iteration tolerances used by the default KSP convergence testers.
918: Logically Collective on KSP
920: Input Parameters:
921: + ksp - the Krylov subspace context
922: . rtol - the relative convergence tolerance
923: (relative decrease in the residual norm)
924: . abstol - the absolute convergence tolerance
925: (absolute size of the residual norm)
926: . dtol - the divergence tolerance
927: (amount residual can increase before KSPDefaultConverged()
928: concludes that the method is diverging)
929: - maxits - maximum number of iterations to use
931: Options Database Keys:
932: + -ksp_atol <abstol> - Sets abstol
933: . -ksp_rtol <rtol> - Sets rtol
934: . -ksp_divtol <dtol> - Sets dtol
935: - -ksp_max_it <maxits> - Sets maxits
937: Notes:
938: Use PETSC_DEFAULT to retain the default value of any of the tolerances.
940: See KSPDefaultConverged() for details on the use of these parameters
941: in the default convergence test. See also KSPSetConvergenceTest()
942: for setting user-defined stopping criteria.
944: Level: intermediate
946: .keywords: KSP, set, tolerance, absolute, relative, divergence,
947: convergence, maximum, iterations
949: .seealso: KSPGetTolerances(), KSPDefaultConverged(), KSPSetConvergenceTest()
950: @*/
951: PetscErrorCode KSPSetTolerances(KSP ksp,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt maxits)
952: {
960: if (rtol != PETSC_DEFAULT) {
961: 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;
962: }
963: if (abstol != PETSC_DEFAULT) {
964: if (abstol < 0.0) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %G must be non-negative",abstol);
965: ksp->abstol = abstol;
966: }
967: if (dtol != PETSC_DEFAULT) {
968: if (dtol < 0.0) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Divergence tolerance %G must be larger than 1.0",dtol);
969: ksp->divtol = dtol;
970: }
971: if (maxits != PETSC_DEFAULT) {
972: if (maxits < 0) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxits);
973: ksp->max_it = maxits;
974: }
975: return(0);
976: }
980: /*@
981: KSPSetInitialGuessNonzero - Tells the iterative solver that the
982: initial guess is nonzero; otherwise KSP assumes the initial guess
983: is to be zero (and thus zeros it out before solving).
985: Logically Collective on KSP
987: Input Parameters:
988: + ksp - iterative context obtained from KSPCreate()
989: - flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero
991: Options database keys:
992: . -ksp_initial_guess_nonzero : use nonzero initial guess; this takes an optional truth value (0/1/no/yes/true/false)
994: Level: beginner
996: Notes:
997: If this is not called the X vector is zeroed in the call to KSPSolve().
999: .keywords: KSP, set, initial guess, nonzero
1001: .seealso: KSPGetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll()
1002: @*/
1003: PetscErrorCode KSPSetInitialGuessNonzero(KSP ksp,PetscBool flg)
1004: {
1008: ksp->guess_zero = (PetscBool)!(int)flg;
1009: return(0);
1010: }
1014: /*@
1015: KSPGetInitialGuessNonzero - Determines whether the KSP solver is using
1016: a zero initial guess.
1018: Not Collective
1020: Input Parameter:
1021: . ksp - iterative context obtained from KSPCreate()
1023: Output Parameter:
1024: . flag - PETSC_TRUE if guess is nonzero, else PETSC_FALSE
1026: Level: intermediate
1028: .keywords: KSP, set, initial guess, nonzero
1030: .seealso: KSPSetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll()
1031: @*/
1032: PetscErrorCode KSPGetInitialGuessNonzero(KSP ksp,PetscBool *flag)
1033: {
1037: if (ksp->guess_zero) *flag = PETSC_FALSE;
1038: else *flag = PETSC_TRUE;
1039: return(0);
1040: }
1044: /*@
1045: KSPSetErrorIfNotConverged - Causes KSPSolve() to generate an error if the solver has not converged.
1047: Logically Collective on KSP
1049: Input Parameters:
1050: + ksp - iterative context obtained from KSPCreate()
1051: - flg - PETSC_TRUE indicates you want the error generated
1053: Options database keys:
1054: . -ksp_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false)
1056: Level: intermediate
1058: Notes:
1059: Normally PETSc continues if a linear solver fails to converge, you can call KSPGetConvergedReason() after a KSPSolve()
1060: to determine if it has converged.
1062: .keywords: KSP, set, initial guess, nonzero
1064: .seealso: KSPGetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll(), KSPGetErrorIfNotConverged()
1065: @*/
1066: PetscErrorCode KSPSetErrorIfNotConverged(KSP ksp,PetscBool flg)
1067: {
1071: ksp->errorifnotconverged = flg;
1072: return(0);
1073: }
1077: /*@
1078: KSPGetErrorIfNotConverged - Will KSPSolve() generate an error if the solver does not converge?
1080: Not Collective
1082: Input Parameter:
1083: . ksp - iterative context obtained from KSPCreate()
1085: Output Parameter:
1086: . flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE
1088: Level: intermediate
1090: .keywords: KSP, set, initial guess, nonzero
1092: .seealso: KSPSetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll(), KSPSetErrorIfNotConverged()
1093: @*/
1094: PetscErrorCode KSPGetErrorIfNotConverged(KSP ksp,PetscBool *flag)
1095: {
1099: *flag = ksp->errorifnotconverged;
1100: return(0);
1101: }
1105: /*@
1106: KSPSetInitialGuessKnoll - Tells the iterative solver to use PCApply(pc,b,..) to compute the initial guess (The Knoll trick)
1108: Logically Collective on KSP
1110: Input Parameters:
1111: + ksp - iterative context obtained from KSPCreate()
1112: - flg - PETSC_TRUE or PETSC_FALSE
1114: Level: advanced
1117: .keywords: KSP, set, initial guess, nonzero
1119: .seealso: KSPGetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero()
1120: @*/
1121: PetscErrorCode KSPSetInitialGuessKnoll(KSP ksp,PetscBool flg)
1122: {
1126: ksp->guess_knoll = flg;
1127: return(0);
1128: }
1132: /*@
1133: KSPGetInitialGuessKnoll - Determines whether the KSP solver is using the Knoll trick (using PCApply(pc,b,...) to compute
1134: the initial guess
1136: Not Collective
1138: Input Parameter:
1139: . ksp - iterative context obtained from KSPCreate()
1141: Output Parameter:
1142: . flag - PETSC_TRUE if using Knoll trick, else PETSC_FALSE
1144: Level: advanced
1146: .keywords: KSP, set, initial guess, nonzero
1148: .seealso: KSPSetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero()
1149: @*/
1150: PetscErrorCode KSPGetInitialGuessKnoll(KSP ksp,PetscBool *flag)
1151: {
1155: *flag = ksp->guess_knoll;
1156: return(0);
1157: }
1161: /*@
1162: KSPGetComputeSingularValues - Gets the flag indicating whether the extreme singular
1163: values will be calculated via a Lanczos or Arnoldi process as the linear
1164: system is solved.
1166: Not Collective
1168: Input Parameter:
1169: . ksp - iterative context obtained from KSPCreate()
1171: Output Parameter:
1172: . flg - PETSC_TRUE or PETSC_FALSE
1174: Options Database Key:
1175: . -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()
1177: Notes:
1178: Currently this option is not valid for all iterative methods.
1180: Many users may just want to use the monitoring routine
1181: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1182: to print the singular values at each iteration of the linear solve.
1184: Level: advanced
1186: .keywords: KSP, set, compute, singular values
1188: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue()
1189: @*/
1190: PetscErrorCode KSPGetComputeSingularValues(KSP ksp,PetscBool *flg)
1191: {
1195: *flg = ksp->calc_sings;
1196: return(0);
1197: }
1201: /*@
1202: KSPSetComputeSingularValues - Sets a flag so that the extreme singular
1203: values will be calculated via a Lanczos or Arnoldi process as the linear
1204: system is solved.
1206: Logically Collective on KSP
1208: Input Parameters:
1209: + ksp - iterative context obtained from KSPCreate()
1210: - flg - PETSC_TRUE or PETSC_FALSE
1212: Options Database Key:
1213: . -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()
1215: Notes:
1216: Currently this option is not valid for all iterative methods.
1218: Many users may just want to use the monitoring routine
1219: KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1220: to print the singular values at each iteration of the linear solve.
1222: Level: advanced
1224: .keywords: KSP, set, compute, singular values
1226: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue()
1227: @*/
1228: PetscErrorCode KSPSetComputeSingularValues(KSP ksp,PetscBool flg)
1229: {
1233: ksp->calc_sings = flg;
1234: return(0);
1235: }
1239: /*@
1240: KSPGetComputeEigenvalues - Gets the flag indicating that the extreme eigenvalues
1241: values will be calculated via a Lanczos or Arnoldi process as the linear
1242: system is solved.
1244: Not Collective
1246: Input Parameter:
1247: . ksp - iterative context obtained from KSPCreate()
1249: Output Parameter:
1250: . flg - PETSC_TRUE or PETSC_FALSE
1252: Notes:
1253: Currently this option is not valid for all iterative methods.
1255: Level: advanced
1257: .keywords: KSP, set, compute, eigenvalues
1259: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly()
1260: @*/
1261: PetscErrorCode KSPGetComputeEigenvalues(KSP ksp,PetscBool *flg)
1262: {
1266: *flg = ksp->calc_sings;
1267: return(0);
1268: }
1272: /*@
1273: KSPSetComputeEigenvalues - Sets a flag so that the extreme eigenvalues
1274: values will be calculated via a Lanczos or Arnoldi process as the linear
1275: system is solved.
1277: Logically Collective on KSP
1279: Input Parameters:
1280: + ksp - iterative context obtained from KSPCreate()
1281: - flg - PETSC_TRUE or PETSC_FALSE
1283: Notes:
1284: Currently this option is not valid for all iterative methods.
1286: Level: advanced
1288: .keywords: KSP, set, compute, eigenvalues
1290: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly()
1291: @*/
1292: PetscErrorCode KSPSetComputeEigenvalues(KSP ksp,PetscBool flg)
1293: {
1297: ksp->calc_sings = flg;
1298: return(0);
1299: }
1303: /*@
1304: KSPGetRhs - Gets the right-hand-side vector for the linear system to
1305: be solved.
1307: Not Collective
1309: Input Parameter:
1310: . ksp - iterative context obtained from KSPCreate()
1312: Output Parameter:
1313: . r - right-hand-side vector
1315: Level: developer
1317: .keywords: KSP, get, right-hand-side, rhs
1319: .seealso: KSPGetSolution(), KSPSolve()
1320: @*/
1321: PetscErrorCode KSPGetRhs(KSP ksp,Vec *r)
1322: {
1326: *r = ksp->vec_rhs;
1327: return(0);
1328: }
1332: /*@
1333: KSPGetSolution - Gets the location of the solution for the
1334: linear system to be solved. Note that this may not be where the solution
1335: is stored during the iterative process; see KSPBuildSolution().
1337: Not Collective
1339: Input Parameters:
1340: . ksp - iterative context obtained from KSPCreate()
1342: Output Parameters:
1343: . v - solution vector
1345: Level: developer
1347: .keywords: KSP, get, solution
1349: .seealso: KSPGetRhs(), KSPBuildSolution(), KSPSolve()
1350: @*/
1351: PetscErrorCode KSPGetSolution(KSP ksp,Vec *v)
1352: {
1356: *v = ksp->vec_sol;
1357: return(0);
1358: }
1362: /*@
1363: KSPSetPC - Sets the preconditioner to be used to calculate the
1364: application of the preconditioner on a vector.
1366: Collective on KSP
1368: Input Parameters:
1369: + ksp - iterative context obtained from KSPCreate()
1370: - pc - the preconditioner object
1372: Notes:
1373: Use KSPGetPC() to retrieve the preconditioner context (for example,
1374: to free it at the end of the computations).
1376: Level: developer
1378: .keywords: KSP, set, precondition, Binv
1380: .seealso: KSPGetPC()
1381: @*/
1382: PetscErrorCode KSPSetPC(KSP ksp,PC pc)
1383: {
1390: PetscObjectReference((PetscObject)pc);
1391: PCDestroy(&ksp->pc);
1392: ksp->pc = pc;
1393: PetscLogObjectParent(ksp,ksp->pc);
1394: return(0);
1395: }
1399: /*@
1400: KSPGetPC - Returns a pointer to the preconditioner context
1401: set with KSPSetPC().
1403: Not Collective
1405: Input Parameters:
1406: . ksp - iterative context obtained from KSPCreate()
1408: Output Parameter:
1409: . pc - preconditioner context
1411: Level: developer
1413: .keywords: KSP, get, preconditioner, Binv
1415: .seealso: KSPSetPC()
1416: @*/
1417: PetscErrorCode KSPGetPC(KSP ksp,PC *pc)
1418: {
1424: if (!ksp->pc) {
1425: PCCreate(((PetscObject)ksp)->comm,&ksp->pc);
1426: PetscObjectIncrementTabLevel((PetscObject)ksp->pc,(PetscObject)ksp,0);
1427: PetscLogObjectParent(ksp,ksp->pc);
1428: }
1429: *pc = ksp->pc;
1430: return(0);
1431: }
1435: /*@
1436: KSPMonitor - runs the user provided monitor routines, if they exist
1438: Collective on KSP
1440: Input Parameters:
1441: + ksp - iterative context obtained from KSPCreate()
1442: . it - iteration number
1443: - rnorm - relative norm of the residual
1445: Notes:
1446: This routine is called by the KSP implementations.
1447: It does not typically need to be called by the user.
1449: Level: developer
1451: .seealso: KSPMonitorSet()
1452: @*/
1453: PetscErrorCode KSPMonitor(KSP ksp,PetscInt it,PetscReal rnorm)
1454: {
1455: PetscInt i, n = ksp->numbermonitors;
1459: for (i=0; i<n; i++) {
1460: (*ksp->monitor[i])(ksp,it,rnorm,ksp->monitorcontext[i]);
1461: }
1462: return(0);
1463: }
1467: /*@C
1468: KSPMonitorSet - Sets an ADDITIONAL function to be called at every iteration to monitor
1469: the residual/error etc.
1470:
1471: Logically Collective on KSP
1473: Input Parameters:
1474: + ksp - iterative context obtained from KSPCreate()
1475: . monitor - pointer to function (if this is PETSC_NULL, it turns off monitoring
1476: . mctx - [optional] context for private data for the
1477: monitor routine (use PETSC_NULL if no context is desired)
1478: - monitordestroy - [optional] routine that frees monitor context
1479: (may be PETSC_NULL)
1481: Calling Sequence of monitor:
1482: $ monitor (KSP ksp, int it, PetscReal rnorm, void *mctx)
1484: + ksp - iterative context obtained from KSPCreate()
1485: . it - iteration number
1486: . rnorm - (estimated) 2-norm of (preconditioned) residual
1487: - mctx - optional monitoring context, as set by KSPMonitorSet()
1489: Options Database Keys:
1490: + -ksp_monitor - sets KSPMonitorDefault()
1491: . -ksp_monitor_true_residual - sets KSPMonitorTrueResidualNorm()
1492: . -ksp_monitor_draw - sets line graph monitor,
1493: uses KSPMonitorLGCreate()
1494: . -ksp_monitor_draw_true_residual - sets line graph monitor,
1495: uses KSPMonitorLGCreate()
1496: . -ksp_monitor_singular_value - sets KSPMonitorSingularValue()
1497: - -ksp_monitor_cancel - cancels all monitors that have
1498: been hardwired into a code by
1499: calls to KSPMonitorSet(), but
1500: does not cancel those set via
1501: the options database.
1503: Notes:
1504: The default is to do nothing. To print the residual, or preconditioned
1505: residual if KSPSetNormType(ksp,KSP_NORM_PRECONDITIONED) was called, use
1506: KSPMonitorDefault() as the monitoring routine, with a null monitoring
1507: context.
1509: Several different monitoring routines may be set by calling
1510: KSPMonitorSet() multiple times; all will be called in the
1511: order in which they were set.
1513: Fortran notes: Only a single monitor function can be set for each KSP object
1515: Level: beginner
1517: .keywords: KSP, set, monitor
1519: .seealso: KSPMonitorDefault(), KSPMonitorLGCreate(), KSPMonitorCancel()
1520: @*/
1521: PetscErrorCode KSPMonitorSet(KSP ksp,PetscErrorCode (*monitor)(KSP,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
1522: {
1523: PetscInt i;
1528: if (ksp->numbermonitors >= MAXKSPMONITORS) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Too many KSP monitors set");
1529: for (i=0; i<ksp->numbermonitors;i++) {
1530: if (monitor == ksp->monitor[i] && monitordestroy == ksp->monitordestroy[i] && mctx == ksp->monitorcontext[i]) {
1531: if (monitordestroy) {
1532: (*monitordestroy)(&mctx);
1533: }
1534: return(0);
1535: }
1536: }
1537: ksp->monitor[ksp->numbermonitors] = monitor;
1538: ksp->monitordestroy[ksp->numbermonitors] = monitordestroy;
1539: ksp->monitorcontext[ksp->numbermonitors++] = (void*)mctx;
1540: return(0);
1541: }
1545: /*@
1546: KSPMonitorCancel - Clears all monitors for a KSP object.
1548: Logically Collective on KSP
1550: Input Parameters:
1551: . ksp - iterative context obtained from KSPCreate()
1553: Options Database Key:
1554: . -ksp_monitor_cancel - Cancels all monitors that have
1555: been hardwired into a code by calls to KSPMonitorSet(),
1556: but does not cancel those set via the options database.
1558: Level: intermediate
1560: .keywords: KSP, set, monitor
1562: .seealso: KSPMonitorDefault(), KSPMonitorLGCreate(), KSPMonitorSet()
1563: @*/
1564: PetscErrorCode KSPMonitorCancel(KSP ksp)
1565: {
1567: PetscInt i;
1571: for (i=0; i<ksp->numbermonitors; i++) {
1572: if (ksp->monitordestroy[i]) {
1573: (*ksp->monitordestroy[i])(&ksp->monitorcontext[i]);
1574: }
1575: }
1576: ksp->numbermonitors = 0;
1577: return(0);
1578: }
1582: /*@C
1583: KSPGetMonitorContext - Gets the monitoring context, as set by
1584: KSPMonitorSet() for the FIRST monitor only.
1586: Not Collective
1588: Input Parameter:
1589: . ksp - iterative context obtained from KSPCreate()
1591: Output Parameter:
1592: . ctx - monitoring context
1594: Level: intermediate
1596: .keywords: KSP, get, monitor, context
1598: .seealso: KSPMonitorDefault(), KSPMonitorLGCreate()
1599: @*/
1600: PetscErrorCode KSPGetMonitorContext(KSP ksp,void **ctx)
1601: {
1604: *ctx = (ksp->monitorcontext[0]);
1605: return(0);
1606: }
1610: /*@
1611: KSPSetResidualHistory - Sets the array used to hold the residual history.
1612: If set, this array will contain the residual norms computed at each
1613: iteration of the solver.
1615: Not Collective
1617: Input Parameters:
1618: + ksp - iterative context obtained from KSPCreate()
1619: . a - array to hold history
1620: . na - size of a
1621: - reset - PETSC_TRUE indicates the history counter is reset to zero
1622: for each new linear solve
1624: Level: advanced
1626: Notes: The array is NOT freed by PETSc so the user needs to keep track of
1627: it and destroy once the KSP object is destroyed.
1629: If 'a' is PETSC_NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
1630: default array of length 10000 is allocated.
1632: .keywords: KSP, set, residual, history, norm
1634: .seealso: KSPGetResidualHistory()
1636: @*/
1637: PetscErrorCode KSPSetResidualHistory(KSP ksp,PetscReal a[],PetscInt na,PetscBool reset)
1638: {
1644: PetscFree(ksp->res_hist_alloc);
1645: if (na != PETSC_DECIDE && na != PETSC_DEFAULT && a) {
1646: ksp->res_hist = a;
1647: ksp->res_hist_max = na;
1648: } else {
1649: if (na != PETSC_DECIDE && na != PETSC_DEFAULT)
1650: ksp->res_hist_max = na;
1651: else
1652: ksp->res_hist_max = 10000; /* like default ksp->max_it */
1653: PetscMalloc(ksp->res_hist_max*sizeof(PetscReal),&ksp->res_hist_alloc);
1654: ksp->res_hist = ksp->res_hist_alloc;
1655: }
1656: ksp->res_hist_len = 0;
1657: ksp->res_hist_reset = reset;
1659: return(0);
1660: }
1664: /*@C
1665: KSPGetResidualHistory - Gets the array used to hold the residual history
1666: and the number of residuals it contains.
1668: Not Collective
1670: Input Parameter:
1671: . ksp - iterative context obtained from KSPCreate()
1673: Output Parameters:
1674: + a - pointer to array to hold history (or PETSC_NULL)
1675: - na - number of used entries in a (or PETSC_NULL)
1677: Level: advanced
1679: Notes:
1680: Can only be called after a KSPSetResidualHistory() otherwise a and na are set to zero
1682: The Fortran version of this routine has a calling sequence
1683: $ call KSPGetResidualHistory(KSP ksp, integer na, integer ierr)
1684: note that you have passed a Fortran array into KSPSetResidualHistory() and you need
1685: to access the residual values from this Fortran array you provided. Only the na (number of
1686: residual norms currently held) is set.
1688: .keywords: KSP, get, residual, history, norm
1690: .seealso: KSPGetResidualHistory()
1692: @*/
1693: PetscErrorCode KSPGetResidualHistory(KSP ksp,PetscReal *a[],PetscInt *na)
1694: {
1697: if (a) *a = ksp->res_hist;
1698: if (na) *na = ksp->res_hist_len;
1699: return(0);
1700: }
1704: /*@C
1705: KSPSetConvergenceTest - Sets the function to be used to determine
1706: convergence.
1708: Logically Collective on KSP
1710: Input Parameters:
1711: + ksp - iterative context obtained from KSPCreate()
1712: . converge - pointer to int function
1713: . cctx - context for private data for the convergence routine (may be null)
1714: - destroy - a routine for destroying the context (may be null)
1716: Calling sequence of converge:
1717: $ converge (KSP ksp, int it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)
1719: + ksp - iterative context obtained from KSPCreate()
1720: . it - iteration number
1721: . rnorm - (estimated) 2-norm of (preconditioned) residual
1722: . reason - the reason why it has converged or diverged
1723: - cctx - optional convergence context, as set by KSPSetConvergenceTest()
1726: Notes:
1727: Must be called after the KSP type has been set so put this after
1728: a call to KSPSetType(), or KSPSetFromOptions().
1730: The default convergence test, KSPDefaultConverged(), aborts if the
1731: residual grows to more than 10000 times the initial residual.
1733: The default is a combination of relative and absolute tolerances.
1734: The residual value that is tested may be an approximation; routines
1735: that need exact values should compute them.
1737: In the default PETSc convergence test, the precise values of reason
1738: are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
1740: Level: advanced
1742: .keywords: KSP, set, convergence, test, context
1744: .seealso: KSPDefaultConverged(), KSPGetConvergenceContext(), KSPSetTolerances()
1745: @*/
1746: PetscErrorCode KSPSetConvergenceTest(KSP ksp,PetscErrorCode (*converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
1747: {
1752: if (ksp->convergeddestroy) {
1753: (*ksp->convergeddestroy)(ksp->cnvP);
1754: }
1755: ksp->converged = converge;
1756: ksp->convergeddestroy = destroy;
1757: ksp->cnvP = (void*)cctx;
1758: return(0);
1759: }
1763: /*@C
1764: KSPGetConvergenceContext - Gets the convergence context set with
1765: KSPSetConvergenceTest().
1767: Not Collective
1769: Input Parameter:
1770: . ksp - iterative context obtained from KSPCreate()
1772: Output Parameter:
1773: . ctx - monitoring context
1775: Level: advanced
1777: .keywords: KSP, get, convergence, test, context
1779: .seealso: KSPDefaultConverged(), KSPSetConvergenceTest()
1780: @*/
1781: PetscErrorCode KSPGetConvergenceContext(KSP ksp,void **ctx)
1782: {
1785: *ctx = ksp->cnvP;
1786: return(0);
1787: }
1791: /*@C
1792: KSPBuildSolution - Builds the approximate solution in a vector provided.
1793: This routine is NOT commonly needed (see KSPSolve()).
1795: Collective on KSP
1797: Input Parameter:
1798: . ctx - iterative context obtained from KSPCreate()
1800: Output Parameter:
1801: Provide exactly one of
1802: + v - location to stash solution.
1803: - V - the solution is returned in this location. This vector is created
1804: internally. This vector should NOT be destroyed by the user with
1805: VecDestroy().
1807: Notes:
1808: This routine can be used in one of two ways
1809: .vb
1810: KSPBuildSolution(ksp,PETSC_NULL,&V);
1811: or
1812: KSPBuildSolution(ksp,v,PETSC_NULL); or KSPBuildSolution(ksp,v,&v);
1813: .ve
1814: In the first case an internal vector is allocated to store the solution
1815: (the user cannot destroy this vector). In the second case the solution
1816: is generated in the vector that the user provides. Note that for certain
1817: methods, such as KSPCG, the second case requires a copy of the solution,
1818: while in the first case the call is essentially free since it simply
1819: returns the vector where the solution already is stored. For some methods
1820: like GMRES this is a reasonably expensive operation and should only be
1821: used in truly needed.
1823: Level: advanced
1825: .keywords: KSP, build, solution
1827: .seealso: KSPGetSolution(), KSPBuildResidual()
1828: @*/
1829: PetscErrorCode KSPBuildSolution(KSP ksp,Vec v,Vec *V)
1830: {
1835: if (!V && !v) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_WRONG,"Must provide either v or V");
1836: if (!V) V = &v;
1837: (*ksp->ops->buildsolution)(ksp,v,V);
1838: return(0);
1839: }
1843: /*@C
1844: KSPBuildResidual - Builds the residual in a vector provided.
1846: Collective on KSP
1848: Input Parameter:
1849: . ksp - iterative context obtained from KSPCreate()
1851: Output Parameters:
1852: + v - optional location to stash residual. If v is not provided,
1853: then a location is generated.
1854: . t - work vector. If not provided then one is generated.
1855: - V - the residual
1857: Notes:
1858: Regardless of whether or not v is provided, the residual is
1859: returned in V.
1861: Level: advanced
1863: .keywords: KSP, build, residual
1865: .seealso: KSPBuildSolution()
1866: @*/
1867: PetscErrorCode KSPBuildResidual(KSP ksp,Vec t,Vec v,Vec *V)
1868: {
1870: PetscBool flag = PETSC_FALSE;
1871: Vec w = v,tt = t;
1875: if (!w) {
1876: VecDuplicate(ksp->vec_rhs,&w);
1877: PetscLogObjectParent((PetscObject)ksp,w);
1878: }
1879: if (!tt) {
1880: VecDuplicate(ksp->vec_sol,&tt); flag = PETSC_TRUE;
1881: PetscLogObjectParent((PetscObject)ksp,tt);
1882: }
1883: (*ksp->ops->buildresidual)(ksp,tt,w,V);
1884: if (flag) {VecDestroy(&tt);}
1885: return(0);
1886: }
1890: /*@
1891: KSPSetDiagonalScale - Tells KSP to symmetrically diagonally scale the system
1892: before solving. This actually CHANGES the matrix (and right hand side).
1894: Logically Collective on KSP
1896: Input Parameter:
1897: + ksp - the KSP context
1898: - scale - PETSC_TRUE or PETSC_FALSE
1900: Options Database Key:
1901: + -ksp_diagonal_scale -
1902: - -ksp_diagonal_scale_fix - scale the matrix back AFTER the solve
1905: BE CAREFUL with this routine: it actually scales the matrix and right
1906: hand side that define the system. After the system is solved the matrix
1907: and right hand side remain scaled unless you use KSPSetDiagonalScaleFix()
1909: This should NOT be used within the SNES solves if you are using a line
1910: search.
1911:
1912: If you use this with the PCType Eisenstat preconditioner than you can
1913: use the PCEisenstatNoDiagonalScaling() option, or -pc_eisenstat_no_diagonal_scaling
1914: to save some unneeded, redundant flops.
1916: Level: intermediate
1918: .keywords: KSP, set, options, prefix, database
1920: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScaleFix()
1921: @*/
1922: PetscErrorCode KSPSetDiagonalScale(KSP ksp,PetscBool scale)
1923: {
1927: ksp->dscale = scale;
1928: return(0);
1929: }
1933: /*@
1934: KSPGetDiagonalScale - Checks if KSP solver scales the matrix and
1935: right hand side
1937: Not Collective
1939: Input Parameter:
1940: . ksp - the KSP context
1942: Output Parameter:
1943: . scale - PETSC_TRUE or PETSC_FALSE
1945: Notes:
1946: BE CAREFUL with this routine: it actually scales the matrix and right
1947: hand side that define the system. After the system is solved the matrix
1948: and right hand side remain scaled unless you use KSPSetDiagonalScaleFix()
1950: Level: intermediate
1952: .keywords: KSP, set, options, prefix, database
1954: .seealso: KSPSetDiagonalScale(), KSPSetDiagonalScaleFix()
1955: @*/
1956: PetscErrorCode KSPGetDiagonalScale(KSP ksp,PetscBool *scale)
1957: {
1961: *scale = ksp->dscale;
1962: return(0);
1963: }
1967: /*@
1968: KSPSetDiagonalScaleFix - Tells KSP to diagonally scale the system
1969: back after solving.
1971: Logically Collective on KSP
1973: Input Parameter:
1974: + ksp - the KSP context
1975: - fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
1976: rescale (default)
1978: Notes:
1979: Must be called after KSPSetDiagonalScale()
1981: Using this will slow things down, because it rescales the matrix before and
1982: after each linear solve. This is intended mainly for testing to allow one
1983: to easily get back the original system to make sure the solution computed is
1984: accurate enough.
1986: Level: intermediate
1988: .keywords: KSP, set, options, prefix, database
1990: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPGetDiagonalScaleFix()
1991: @*/
1992: PetscErrorCode KSPSetDiagonalScaleFix(KSP ksp,PetscBool fix)
1993: {
1997: ksp->dscalefix = fix;
1998: return(0);
1999: }
2003: /*@
2004: KSPGetDiagonalScaleFix - Determines if KSP diagonally scales the system
2005: back after solving.
2007: Not Collective
2009: Input Parameter:
2010: . ksp - the KSP context
2012: Output Parameter:
2013: . fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
2014: rescale (default)
2016: Notes:
2017: Must be called after KSPSetDiagonalScale()
2019: If PETSC_TRUE will slow things down, because it rescales the matrix before and
2020: after each linear solve. This is intended mainly for testing to allow one
2021: to easily get back the original system to make sure the solution computed is
2022: accurate enough.
2024: Level: intermediate
2026: .keywords: KSP, set, options, prefix, database
2028: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPSetDiagonalScaleFix()
2029: @*/
2030: PetscErrorCode KSPGetDiagonalScaleFix(KSP ksp,PetscBool *fix)
2031: {
2035: *fix = ksp->dscalefix;
2036: return(0);
2037: }
2041: /*@C
2042: KSPSetComputeOperators - set routine to compute the linear operators
2044: Logically Collective
2046: Input Arguments:
2047: + ksp - the KSP context
2048: . func - function to compute the operators
2049: - ctx - optional context
2051: Calling sequence of func:
2052: $ func(KSP ksp,Mat *A,Mat *B,MatStructure *mstruct,void *ctx)
2054: + ksp - the KSP context
2055: . A - the linear operator
2056: . B - preconditioning matrix
2057: . mstruct - flag indicating structure, same as in KSPSetOperators(), one of SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER
2058: - ctx - optional user-provided context
2060: Level: beginner
2062: .seealso: KSPSetOperators()
2063: @*/
2064: PetscErrorCode KSPSetComputeOperators(KSP ksp,PetscErrorCode (*func)(KSP,Mat,Mat,MatStructure*,void*),void *ctx)
2065: {
2067: DM dm;
2071: KSPGetDM(ksp,&dm);
2072: DMKSPSetComputeOperators(dm,func,ctx);
2073: if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX;
2074: return(0);
2075: }
2079: /*@C
2080: KSPSetComputeRHS - set routine to compute the right hand side of the linear system
2082: Logically Collective
2084: Input Arguments:
2085: + ksp - the KSP context
2086: . func - function to compute the right hand side
2087: - ctx - optional context
2089: Calling sequence of func:
2090: $ func(KSP ksp,Vec b,void *ctx)
2092: + ksp - the KSP context
2093: . b - right hand side of linear system
2094: - ctx - optional user-provided context
2096: Level: beginner
2098: .seealso: KSPSolve()
2099: @*/
2100: PetscErrorCode KSPSetComputeRHS(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2101: {
2103: DM dm;
2107: KSPGetDM(ksp,&dm);
2108: DMKSPSetComputeRHS(dm,func,ctx);
2109: return(0);
2110: }