Actual source code: kspimpl.h
2: #ifndef _KSPIMPL_H
3: #define _KSPIMPL_H
5: #include <petscksp.h>
6: #include <petscds.h>
7: #include <petsc/private/petscimpl.h>
9: PETSC_EXTERN PetscBool KSPRegisterAllCalled;
10: PETSC_EXTERN PetscBool KSPMonitorRegisterAllCalled;
11: PETSC_EXTERN PetscErrorCode KSPRegisterAll(void);
12: PETSC_EXTERN PetscErrorCode KSPMonitorRegisterAll(void);
13: PETSC_EXTERN PetscErrorCode KSPGuessRegisterAll(void);
14: PETSC_EXTERN PetscErrorCode KSPMatRegisterAll(void);
16: typedef struct _KSPOps *KSPOps;
18: struct _KSPOps {
19: PetscErrorCode (*buildsolution)(KSP,Vec,Vec*); /* Returns a pointer to the solution, or
20: calculates the solution in a
21: user-provided area. */
22: PetscErrorCode (*buildresidual)(KSP,Vec,Vec,Vec*); /* Returns a pointer to the residual, or
23: calculates the residual in a
24: user-provided area. */
25: PetscErrorCode (*solve)(KSP); /* actual solver */
26: PetscErrorCode (*matsolve)(KSP,Mat,Mat); /* multiple dense RHS solver */
27: PetscErrorCode (*setup)(KSP);
28: PetscErrorCode (*setfromoptions)(PetscOptionItems*,KSP);
29: PetscErrorCode (*publishoptions)(KSP);
30: PetscErrorCode (*computeextremesingularvalues)(KSP,PetscReal*,PetscReal*);
31: PetscErrorCode (*computeeigenvalues)(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt *);
32: PetscErrorCode (*computeritz)(KSP,PetscBool,PetscBool,PetscInt*,Vec[],PetscReal*,PetscReal*);
33: PetscErrorCode (*destroy)(KSP);
34: PetscErrorCode (*view)(KSP,PetscViewer);
35: PetscErrorCode (*reset)(KSP);
36: PetscErrorCode (*load)(KSP,PetscViewer);
37: };
39: typedef struct _KSPGuessOps *KSPGuessOps;
41: struct _KSPGuessOps {
42: PetscErrorCode (*formguess)(KSPGuess,Vec,Vec); /* Form initial guess */
43: PetscErrorCode (*update)(KSPGuess,Vec,Vec); /* Update database */
44: PetscErrorCode (*setfromoptions)(KSPGuess);
45: PetscErrorCode (*setup)(KSPGuess);
46: PetscErrorCode (*destroy)(KSPGuess);
47: PetscErrorCode (*view)(KSPGuess,PetscViewer);
48: PetscErrorCode (*reset)(KSPGuess);
49: };
51: /*
52: Defines the KSPGuess data structure.
53: */
54: struct _p_KSPGuess {
55: PETSCHEADER(struct _KSPGuessOps);
56: KSP ksp; /* the parent KSP */
57: Mat A; /* the current linear operator */
58: PetscObjectState omatstate; /* previous linear operator state */
59: void *data; /* pointer to the specific implementation */
60: };
62: PETSC_EXTERN PetscErrorCode KSPGuessCreate_Fischer(KSPGuess);
63: PETSC_EXTERN PetscErrorCode KSPGuessCreate_POD(KSPGuess);
65: /*
66: Maximum number of monitors you can run with a single KSP
67: */
68: #define MAXKSPMONITORS 5
69: #define MAXKSPREASONVIEWS 5
70: typedef enum {KSP_SETUP_NEW, KSP_SETUP_NEWMATRIX, KSP_SETUP_NEWRHS} KSPSetUpStage;
72: /*
73: Defines the KSP data structure.
74: */
75: struct _p_KSP {
76: PETSCHEADER(struct _KSPOps);
77: DM dm;
78: PetscBool dmAuto; /* DM was created automatically by KSP */
79: PetscBool dmActive; /* KSP should use DM for computing operators */
80: /*------------------------- User parameters--------------------------*/
81: PetscInt max_it; /* maximum number of iterations */
82: KSPGuess guess;
83: PetscBool guess_zero, /* flag for whether initial guess is 0 */
84: calc_sings, /* calculate extreme Singular Values */
85: calc_ritz, /* calculate (harmonic) Ritz pairs */
86: guess_knoll; /* use initial guess of PCApply(ksp->B,b */
87: PCSide pc_side; /* flag for left, right, or symmetric preconditioning */
88: PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */
89: PetscReal rtol, /* relative tolerance */
90: abstol, /* absolute tolerance */
91: ttol, /* (not set by user) */
92: divtol; /* divergence tolerance */
93: PetscReal rnorm0; /* initial residual norm (used for divergence testing) */
94: PetscReal rnorm; /* current residual norm */
95: KSPConvergedReason reason;
96: PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */
98: Vec vec_sol,vec_rhs; /* pointer to where user has stashed
99: the solution and rhs, these are
100: never touched by the code, only
101: passed back to the user */
102: PetscReal *res_hist; /* If !0 stores residual each at iteration */
103: PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
104: PetscInt res_hist_len; /* current size of residual history array */
105: PetscInt res_hist_max; /* actual amount of storage in residual history */
106: PetscBool res_hist_reset; /* reset history to length zero for each new solve */
107: PetscReal *err_hist; /* If !0 stores error at each iteration */
108: PetscReal *err_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
109: PetscInt err_hist_len; /* current size of error history array */
110: PetscInt err_hist_max; /* actual amount of storage in error history */
111: PetscBool err_hist_reset; /* reset history to length zero for each new solve */
113: PetscInt chknorm; /* only compute/check norm if iterations is great than this */
114: PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the
115: MPI_Allreduce() for computing the inner products for the next iteration. */
117: PetscInt nmax; /* maximum number of right-hand sides to be handled simultaneously */
119: /* --------User (or default) routines (most return -1 on error) --------*/
120: PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */
121: PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void**); /* */
122: void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */
123: PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */
124: PetscBool pauseFinal; /* Pause all drawing monitor at the final iterate */
126: PetscErrorCode (*reasonview[MAXKSPREASONVIEWS])(KSP,void*); /* KSP converged reason view */
127: PetscErrorCode (*reasonviewdestroy[MAXKSPREASONVIEWS])(void**); /* Optional destroy routine */
128: void *reasonviewcontext[MAXKSPREASONVIEWS]; /* User context */
129: PetscInt numberreasonviews; /* Number if reason viewers */
131: PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
132: PetscErrorCode (*convergeddestroy)(void*);
133: void *cnvP;
135: void *user; /* optional user-defined context */
137: PC pc;
139: void *data; /* holder for misc stuff associated
140: with a particular iterative solver */
142: PetscBool view, viewPre, viewRate, viewMat, viewPMat, viewRhs, viewSol, viewMatExp, viewEV, viewSV, viewEVExp, viewFinalRes, viewPOpExp, viewDScale;
143: PetscViewer viewer, viewerPre, viewerRate, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp, viewerDScale;
144: PetscViewerFormat format, formatPre, formatRate, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp, formatDScale;
146: /* ----------------Default work-area management -------------------- */
147: PetscInt nwork;
148: Vec *work;
150: KSPSetUpStage setupstage;
151: PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */
153: PetscInt its; /* number of iterations so far computed in THIS linear solve*/
154: PetscInt totalits; /* number of iterations used by this KSP object since it was created */
156: PetscBool transpose_solve; /* solve transpose system instead */
157: struct {
158: Mat AT,BT;
159: PetscBool use_explicittranspose; /* transpose the system explicitly in KSPSolveTranspose */
160: PetscBool reuse_transpose; /* reuse the previous transposed system */
161: } transpose;
163: KSPNormType normtype; /* type of norm used for convergence tests */
165: PCSide pc_side_set; /* PC type set explicitly by user */
166: KSPNormType normtype_set; /* Norm type set explicitly by user */
168: /* Allow diagonally scaling the matrix before computing the preconditioner or using
169: the Krylov method. Note this is NOT just Jacobi preconditioning */
171: PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */
172: PetscBool dscalefix; /* unscale system after solve */
173: PetscBool dscalefix2; /* system has been unscaled */
174: Vec diagonal; /* 1/sqrt(diag of matrix) */
175: Vec truediagonal;
177: PetscInt setfromoptionscalled;
178: PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */
180: PetscViewer eigviewer; /* Viewer where computed eigenvalues are displayed */
182: PetscErrorCode (*presolve)(KSP,Vec,Vec,void*);
183: PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*);
184: void *prectx,*postctx;
185: };
187: typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
188: PetscReal coef;
189: PetscReal bnrm;
190: } KSPDynTolCtx;
192: typedef struct {
193: PetscBool initialrtol; /* default relative residual decrease is computed from initial residual, not rhs */
194: PetscBool mininitialrtol; /* default relative residual decrease is computed from min of initial residual and rhs */
195: PetscBool convmaxits; /* if true, the convergence test returns KSP_CONVERGED_ITS if the maximum number of iterations is reached */
196: Vec work;
197: } KSPConvergedDefaultCtx;
199: PETSC_STATIC_INLINE PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm)
200: {
204: PetscObjectSAWsTakeAccess((PetscObject)ksp);
205: if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) {
206: ksp->res_hist[ksp->res_hist_len++] = norm;
207: }
208: PetscObjectSAWsGrantAccess((PetscObject)ksp);
209: return(0);
210: }
212: PETSC_STATIC_INLINE PetscErrorCode KSPLogErrorHistory(KSP ksp)
213: {
214: DM dm;
218: PetscObjectSAWsTakeAccess((PetscObject) ksp);
219: KSPGetDM(ksp, &dm);
220: if (dm && ksp->err_hist && ksp->err_hist_max > ksp->err_hist_len) {
221: PetscSimplePointFunc exactSol;
222: void *exactCtx;
223: PetscDS ds;
224: Vec u;
225: PetscReal error;
226: PetscInt Nf;
228: KSPBuildSolution(ksp, NULL, &u);
229: /* TODO Was needed to correct for Newton solution, but I just need to set a solution */
230: //VecScale(u, -1.0);
231: /* TODO Case when I have a solution */
232: if (0) {
233: DMGetDS(dm, &ds);
234: PetscDSGetNumFields(ds, &Nf);
235: if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle number of fields %D > 1 right now", Nf);
236: PetscDSGetExactSolution(ds, 0, &exactSol, &exactCtx);
237: DMComputeL2FieldDiff(dm, 0.0, &exactSol, &exactCtx, u, &error);
238: } else {
239: /* The null solution A 0 = 0 */
240: VecNorm(u, NORM_2, &error);
241: }
242: ksp->err_hist[ksp->err_hist_len++] = error;
243: }
244: PetscObjectSAWsGrantAccess((PetscObject) ksp);
245: return(0);
246: }
248: PETSC_STATIC_INLINE PetscScalar KSPNoisyHash_Private(PetscInt xx)
249: {
250: unsigned int x = xx;
251: x = ((x >> 16) ^ x) * 0x45d9f3b;
252: x = ((x >> 16) ^ x) * 0x45d9f3b;
253: x = ((x >> 16) ^ x);
254: return (PetscScalar)((PetscInt64)x-2147483648)*5.e-10; /* center around zero, scaled about -1. to 1.*/
255: }
257: PETSC_STATIC_INLINE PetscErrorCode KSPSetNoisy_Private(Vec v)
258: {
259: PetscScalar *a;
260: PetscInt n, istart, i;
263: VecGetOwnershipRange(v, &istart, NULL);
264: VecGetLocalSize(v, &n);
265: VecGetArrayWrite(v, &a);
266: for (i = 0; i < n; ++i) a[i] = KSPNoisyHash_Private(i+istart);
267: VecRestoreArrayWrite(v, &a);
268: return(0);
269: }
271: PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP,PetscBool,KSPNormType*,PCSide*);
273: PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP,PetscInt,const PetscReal*,const PetscReal*);
275: typedef struct _p_DMKSP *DMKSP;
276: typedef struct _DMKSPOps *DMKSPOps;
277: struct _DMKSPOps {
278: PetscErrorCode (*computeoperators)(KSP,Mat,Mat,void*);
279: PetscErrorCode (*computerhs)(KSP,Vec,void*);
280: PetscErrorCode (*computeinitialguess)(KSP,Vec,void*);
281: PetscErrorCode (*destroy)(DMKSP*);
282: PetscErrorCode (*duplicate)(DMKSP,DMKSP);
283: };
285: struct _p_DMKSP {
286: PETSCHEADER(struct _DMKSPOps);
287: void *operatorsctx;
288: void *rhsctx;
289: void *initialguessctx;
290: void *data;
292: /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
293: * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user
294: * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by
295: * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates
296: * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes
297: * to the original level will no longer propagate to that level.
298: */
299: DM originaldm;
301: void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */
302: };
303: PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM,DMKSP*);
304: PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM,DMKSP*);
305: PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM,DM);
307: /*
308: These allow the various Krylov methods to apply to either the linear system or its transpose.
309: */
310: PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpace(KSP ksp,Vec y)
311: {
314: if (ksp->pc_side == PC_LEFT) {
315: Mat A;
316: MatNullSpace nullsp;
317: PCGetOperators(ksp->pc,&A,NULL);
318: MatGetNullSpace(A,&nullsp);
319: if (nullsp) {
320: MatNullSpaceRemove(nullsp,y);
321: }
322: }
323: return(0);
324: }
326: PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp,Vec y)
327: {
330: if (ksp->pc_side == PC_LEFT) {
331: Mat A;
332: MatNullSpace nullsp;
333: PCGetOperators(ksp->pc,&A,NULL);
334: MatGetTransposeNullSpace(A,&nullsp);
335: if (nullsp) {
336: MatNullSpaceRemove(nullsp,y);
337: }
338: }
339: return(0);
340: }
342: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y)
343: {
346: if (!ksp->transpose_solve) {MatMult(A,x,y);}
347: else {MatMultTranspose(A,x,y);}
348: return(0);
349: }
351: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y)
352: {
355: if (!ksp->transpose_solve) {MatMultTranspose(A,x,y);}
356: else {MatMult(A,x,y);}
357: return(0);
358: }
360: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y)
361: {
364: if (!ksp->transpose_solve) {
365: PCApply(ksp->pc,x,y);
366: KSP_RemoveNullSpace(ksp,y);
367: } else {
368: PCApplyTranspose(ksp->pc,x,y);
369: KSP_RemoveNullSpaceTranspose(ksp,y);
370: }
371: return(0);
372: }
374: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y)
375: {
378: if (!ksp->transpose_solve) {
379: PCApplyTranspose(ksp->pc,x,y);
380: KSP_RemoveNullSpaceTranspose(ksp,y);
381: } else {
382: PCApply(ksp->pc,x,y);
383: KSP_RemoveNullSpace(ksp,y);
384: }
385: return(0);
386: }
388: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w)
389: {
392: if (!ksp->transpose_solve) {
393: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
394: KSP_RemoveNullSpace(ksp,y);
395: } else {
396: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
397: KSP_RemoveNullSpaceTranspose(ksp,y);
398: }
399: return(0);
400: }
402: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w)
403: {
406: if (!ksp->transpose_solve) {
407: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
408: } else {
409: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
410: }
411: return(0);
412: }
414: PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization;
415: PETSC_EXTERN PetscLogEvent KSP_SetUp;
416: PETSC_EXTERN PetscLogEvent KSP_Solve;
417: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0;
418: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1;
419: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2;
420: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3;
421: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4;
422: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S;
423: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L;
424: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U;
425: PETSC_EXTERN PetscLogEvent KSP_SolveTranspose;
426: PETSC_EXTERN PetscLogEvent KSP_MatSolve;
428: PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*);
429: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*);
431: /*MC
432: KSPCheckDot - Checks if the result of a dot product used by the corresponding KSP contains Inf or NaN. These indicate that the previous
433: application of the preconditioner generated an error
435: Collective on ksp
437: Input Parameter:
438: . ksp - the linear solver (KSP) context.
440: Output Parameter:
441: . beta - the result of the inner product
443: Level: developer
445: Developer Note:
446: this is used to manage returning from KSP solvers whose preconditioners have failed in some way
448: .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckNorm(), KSPCheckSolve()
449: M*/
450: #define KSPCheckDot(ksp,beta) do { \
451: if (PetscIsInfOrNanScalar(beta)) { \
452: if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");\
453: else {\
455: PCFailedReason pcreason;\
456: PetscInt sendbuf,recvbuf; \
457: PCGetFailedReasonRank(ksp->pc,&pcreason);\
458: sendbuf = (PetscInt)pcreason; \
459: MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));\
460: if (recvbuf) { \
461: PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf); \
462: ksp->reason = KSP_DIVERGED_PC_FAILED;\
463: VecSetInf(ksp->vec_sol);\
464: } else {\
465: ksp->reason = KSP_DIVERGED_NANORINF;\
466: }\
467: return(0);\
468: }\
469: } } while (0)
471: /*MC
472: KSPCheckNorm - Checks if the result of a norm used by the corresponding KSP contains Inf or NaN. These indicate that the previous
473: application of the preconditioner generated an error
475: Collective on ksp
477: Input Parameter:
478: . ksp - the linear solver (KSP) context.
480: Output Parameter:
481: . beta - the result of the norm
483: Level: developer
485: Developer Note:
486: this is used to manage returning from KSP solvers whose preconditioners have failed in some way
488: .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckDot(), KSPCheckSolve()
489: M*/
490: #define KSPCheckNorm(ksp,beta) do { \
491: if (PetscIsInfOrNanReal(beta)) { \
492: if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf norm");\
493: else {\
495: PCFailedReason pcreason;\
496: PetscInt sendbuf,recvbuf; \
497: PCGetFailedReasonRank(ksp->pc,&pcreason);\
498: sendbuf = (PetscInt)pcreason; \
499: MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));\
500: if (recvbuf) { \
501: PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf); \
502: ksp->reason = KSP_DIVERGED_PC_FAILED; \
503: VecSetInf(ksp->vec_sol);\
504: } else {\
505: PCSetFailedReason(ksp->pc,PC_NOERROR); \
506: ksp->reason = KSP_DIVERGED_NANORINF;\
507: }\
508: return(0);\
509: }\
510: } } while (0)
512: #endif
514: PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]);
515: PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP,PetscInt,PetscReal*);