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 (*settolerance)(KSPGuess,PetscReal);
46: PetscErrorCode (*setup)(KSPGuess);
47: PetscErrorCode (*destroy)(KSPGuess);
48: PetscErrorCode (*view)(KSPGuess,PetscViewer);
49: PetscErrorCode (*reset)(KSPGuess);
50: };
52: /*
53: Defines the KSPGuess data structure.
54: */
55: struct _p_KSPGuess {
56: PETSCHEADER(struct _KSPGuessOps);
57: KSP ksp; /* the parent KSP */
58: Mat A; /* the current linear operator */
59: PetscObjectState omatstate; /* previous linear operator state */
60: void *data; /* pointer to the specific implementation */
61: };
63: PETSC_EXTERN PetscErrorCode KSPGuessCreate_Fischer(KSPGuess);
64: PETSC_EXTERN PetscErrorCode KSPGuessCreate_POD(KSPGuess);
66: /*
67: Maximum number of monitors you can run with a single KSP
68: */
69: #define MAXKSPMONITORS 5
70: #define MAXKSPREASONVIEWS 5
71: typedef enum {KSP_SETUP_NEW, KSP_SETUP_NEWMATRIX, KSP_SETUP_NEWRHS} KSPSetUpStage;
73: /*
74: Defines the KSP data structure.
75: */
76: struct _p_KSP {
77: PETSCHEADER(struct _KSPOps);
78: DM dm;
79: PetscBool dmAuto; /* DM was created automatically by KSP */
80: PetscBool dmActive; /* KSP should use DM for computing operators */
81: /*------------------------- User parameters--------------------------*/
82: PetscInt max_it; /* maximum number of iterations */
83: KSPGuess guess;
84: PetscBool guess_zero, /* flag for whether initial guess is 0 */
85: calc_sings, /* calculate extreme Singular Values */
86: calc_ritz, /* calculate (harmonic) Ritz pairs */
87: guess_knoll; /* use initial guess of PCApply(ksp->B,b */
88: PCSide pc_side; /* flag for left, right, or symmetric preconditioning */
89: PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */
90: PetscReal rtol, /* relative tolerance */
91: abstol, /* absolute tolerance */
92: ttol, /* (not set by user) */
93: divtol; /* divergence tolerance */
94: PetscReal rnorm0; /* initial residual norm (used for divergence testing) */
95: PetscReal rnorm; /* current residual norm */
96: KSPConvergedReason reason;
97: PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */
99: Vec vec_sol,vec_rhs; /* pointer to where user has stashed
100: the solution and rhs, these are
101: never touched by the code, only
102: passed back to the user */
103: PetscReal *res_hist; /* If !0 stores residual each at iteration */
104: PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
105: size_t res_hist_len; /* current size of residual history array */
106: size_t res_hist_max; /* actual amount of storage in residual history */
107: PetscBool res_hist_reset; /* reset history to length zero for each new solve */
108: PetscReal *err_hist; /* If !0 stores error at each iteration */
109: PetscReal *err_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
110: size_t err_hist_len; /* current size of error history array */
111: size_t err_hist_max; /* actual amount of storage in error history */
112: PetscBool err_hist_reset; /* reset history to length zero for each new solve */
114: PetscInt chknorm; /* only compute/check norm if iterations is great than this */
115: PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the
116: MPI_Allreduce() for computing the inner products for the next iteration. */
118: PetscInt nmax; /* maximum number of right-hand sides to be handled simultaneously */
120: /* --------User (or default) routines (most return -1 on error) --------*/
121: PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */
122: PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void**); /* */
123: void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */
124: PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */
125: PetscBool pauseFinal; /* Pause all drawing monitor at the final iterate */
127: PetscErrorCode (*reasonview[MAXKSPREASONVIEWS])(KSP,void*); /* KSP converged reason view */
128: PetscErrorCode (*reasonviewdestroy[MAXKSPREASONVIEWS])(void**); /* Optional destroy routine */
129: void *reasonviewcontext[MAXKSPREASONVIEWS]; /* User context */
130: PetscInt numberreasonviews; /* Number if reason viewers */
132: PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
133: PetscErrorCode (*convergeddestroy)(void*);
134: void *cnvP;
136: void *user; /* optional user-defined context */
138: PC pc;
140: void *data; /* holder for misc stuff associated
141: with a particular iterative solver */
143: PetscBool view, viewPre, viewRate, viewMat, viewPMat, viewRhs, viewSol, viewMatExp, viewEV, viewSV, viewEVExp, viewFinalRes, viewPOpExp, viewDScale;
144: PetscViewer viewer, viewerPre, viewerRate, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp, viewerDScale;
145: PetscViewerFormat format, formatPre, formatRate, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp, formatDScale;
147: /* ----------------Default work-area management -------------------- */
148: PetscInt nwork;
149: Vec *work;
151: KSPSetUpStage setupstage;
152: PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */
154: PetscInt its; /* number of iterations so far computed in THIS linear solve*/
155: PetscInt totalits; /* number of iterations used by this KSP object since it was created */
157: PetscBool transpose_solve; /* solve transpose system instead */
158: struct {
159: Mat AT,BT;
160: PetscBool use_explicittranspose; /* transpose the system explicitly in KSPSolveTranspose */
161: PetscBool reuse_transpose; /* reuse the previous transposed system */
162: } transpose;
164: KSPNormType normtype; /* type of norm used for convergence tests */
166: PCSide pc_side_set; /* PC type set explicitly by user */
167: KSPNormType normtype_set; /* Norm type set explicitly by user */
169: /* Allow diagonally scaling the matrix before computing the preconditioner or using
170: the Krylov method. Note this is NOT just Jacobi preconditioning */
172: PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */
173: PetscBool dscalefix; /* unscale system after solve */
174: PetscBool dscalefix2; /* system has been unscaled */
175: Vec diagonal; /* 1/sqrt(diag of matrix) */
176: Vec truediagonal;
178: PetscInt setfromoptionscalled;
179: PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */
181: PetscViewer eigviewer; /* Viewer where computed eigenvalues are displayed */
183: PetscErrorCode (*presolve)(KSP,Vec,Vec,void*);
184: PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*);
185: void *prectx,*postctx;
186: };
188: typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
189: PetscReal coef;
190: PetscReal bnrm;
191: } KSPDynTolCtx;
193: typedef struct {
194: PetscBool initialrtol; /* default relative residual decrease is computed from initial residual, not rhs */
195: PetscBool mininitialrtol; /* default relative residual decrease is computed from min of initial residual and rhs */
196: PetscBool convmaxits; /* if true, the convergence test returns KSP_CONVERGED_ITS if the maximum number of iterations is reached */
197: Vec work;
198: } KSPConvergedDefaultCtx;
200: static inline PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm)
201: {
202: PetscObjectSAWsTakeAccess((PetscObject)ksp);
203: if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) {
204: ksp->res_hist[ksp->res_hist_len++] = norm;
205: }
206: PetscObjectSAWsGrantAccess((PetscObject)ksp);
207: return 0;
208: }
210: static inline PetscErrorCode KSPLogErrorHistory(KSP ksp)
211: {
212: DM dm;
214: PetscObjectSAWsTakeAccess((PetscObject) ksp);
215: KSPGetDM(ksp, &dm);
216: if (dm && ksp->err_hist && ksp->err_hist_max > ksp->err_hist_len) {
217: PetscSimplePointFunc exactSol;
218: void *exactCtx;
219: PetscDS ds;
220: Vec u;
221: PetscReal error;
222: PetscInt Nf;
224: KSPBuildSolution(ksp, NULL, &u);
225: /* TODO Was needed to correct for Newton solution, but I just need to set a solution */
226: //VecScale(u, -1.0);
227: /* TODO Case when I have a solution */
228: if (0) {
229: DMGetDS(dm, &ds);
230: PetscDSGetNumFields(ds, &Nf);
232: PetscDSGetExactSolution(ds, 0, &exactSol, &exactCtx);
233: DMComputeL2FieldDiff(dm, 0.0, &exactSol, &exactCtx, u, &error);
234: } else {
235: /* The null solution A 0 = 0 */
236: VecNorm(u, NORM_2, &error);
237: }
238: ksp->err_hist[ksp->err_hist_len++] = error;
239: }
240: PetscObjectSAWsGrantAccess((PetscObject) ksp);
241: return 0;
242: }
244: static inline PetscScalar KSPNoisyHash_Private(PetscInt xx)
245: {
246: unsigned int x = (unsigned int) xx;
247: x = ((x >> 16) ^ x) * 0x45d9f3b;
248: x = ((x >> 16) ^ x) * 0x45d9f3b;
249: x = ((x >> 16) ^ x);
250: return (PetscScalar)((PetscInt64)x-2147483648)*5.e-10; /* center around zero, scaled about -1. to 1.*/
251: }
253: static inline PetscErrorCode KSPSetNoisy_Private(Vec v)
254: {
255: PetscScalar *a;
256: PetscInt n, istart;
258: VecGetOwnershipRange(v, &istart, NULL);
259: VecGetLocalSize(v, &n);
260: VecGetArrayWrite(v, &a);
261: for (PetscInt i = 0; i < n; ++i) a[i] = KSPNoisyHash_Private(i+istart);
262: VecRestoreArrayWrite(v, &a);
263: return 0;
264: }
266: PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP,PetscBool,KSPNormType*,PCSide*);
268: PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP,PetscInt,const PetscReal*,const PetscReal*);
270: typedef struct _p_DMKSP *DMKSP;
271: typedef struct _DMKSPOps *DMKSPOps;
272: struct _DMKSPOps {
273: PetscErrorCode (*computeoperators)(KSP,Mat,Mat,void*);
274: PetscErrorCode (*computerhs)(KSP,Vec,void*);
275: PetscErrorCode (*computeinitialguess)(KSP,Vec,void*);
276: PetscErrorCode (*destroy)(DMKSP*);
277: PetscErrorCode (*duplicate)(DMKSP,DMKSP);
278: };
280: struct _p_DMKSP {
281: PETSCHEADER(struct _DMKSPOps);
282: void *operatorsctx;
283: void *rhsctx;
284: void *initialguessctx;
285: void *data;
287: /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
288: * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user
289: * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by
290: * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates
291: * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes
292: * to the original level will no longer propagate to that level.
293: */
294: DM originaldm;
296: void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */
297: };
298: PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM,DMKSP*);
299: PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM,DMKSP*);
300: PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM,DM);
302: /*
303: These allow the various Krylov methods to apply to either the linear system or its transpose.
304: */
305: static inline PetscErrorCode KSP_RemoveNullSpace(KSP ksp,Vec y)
306: {
307: if (ksp->pc_side == PC_LEFT) {
308: Mat A;
309: MatNullSpace nullsp;
311: PCGetOperators(ksp->pc,&A,NULL);
312: MatGetNullSpace(A,&nullsp);
313: if (nullsp) MatNullSpaceRemove(nullsp,y);
314: }
315: return 0;
316: }
318: static inline PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp,Vec y)
319: {
320: if (ksp->pc_side == PC_LEFT) {
321: Mat A;
322: MatNullSpace nullsp;
324: PCGetOperators(ksp->pc,&A,NULL);
325: MatGetTransposeNullSpace(A,&nullsp);
326: if (nullsp) MatNullSpaceRemove(nullsp,y);
327: }
328: return 0;
329: }
331: static inline PetscErrorCode KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y)
332: {
333: if (ksp->transpose_solve) MatMultTranspose(A,x,y);
334: else MatMult(A,x,y);
335: return 0;
336: }
338: static inline PetscErrorCode KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y)
339: {
340: if (ksp->transpose_solve) MatMult(A,x,y);
341: else MatMultTranspose(A,x,y);
342: return 0;
343: }
345: static inline PetscErrorCode KSP_MatMultHermitianTranspose(KSP ksp,Mat A,Vec x,Vec y)
346: {
347: if (!ksp->transpose_solve) MatMultHermitianTranspose(A,x,y);
348: else {
349: Vec w;
351: VecDuplicate(x,&w);
352: VecCopy(x,w);
353: VecConjugate(w);
354: MatMult(A,w,y);
355: VecDestroy(&w);
356: VecConjugate(y);
357: }
358: return 0;
359: }
361: static inline PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y)
362: {
363: if (ksp->transpose_solve) {
364: PCApplyTranspose(ksp->pc,x,y);
365: KSP_RemoveNullSpaceTranspose(ksp,y);
366: } else {
367: PCApply(ksp->pc,x,y);
368: KSP_RemoveNullSpace(ksp,y);
369: }
370: return 0;
371: }
373: static inline PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y)
374: {
375: if (ksp->transpose_solve) {
376: PCApply(ksp->pc,x,y);
377: KSP_RemoveNullSpace(ksp,y);
378: } else {
379: PCApplyTranspose(ksp->pc,x,y);
380: KSP_RemoveNullSpaceTranspose(ksp,y);
381: }
382: return 0;
383: }
385: static inline PetscErrorCode KSP_PCApplyHermitianTranspose(KSP ksp,Vec x,Vec y)
386: {
387: VecConjugate(x);
388: KSP_PCApplyTranspose(ksp,x,y);
389: VecConjugate(x);
390: VecConjugate(y);
391: return 0;
392: }
394: static inline PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w)
395: {
396: if (ksp->transpose_solve) {
397: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
398: KSP_RemoveNullSpaceTranspose(ksp,y);
399: } else {
400: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
401: KSP_RemoveNullSpace(ksp,y);
402: }
403: return 0;
404: }
406: static inline PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w)
407: {
408: if (ksp->transpose_solve) PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
409: else PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
410: return 0;
411: }
413: PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization;
414: PETSC_EXTERN PetscLogEvent KSP_SetUp;
415: PETSC_EXTERN PetscLogEvent KSP_Solve;
416: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0;
417: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1;
418: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2;
419: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3;
420: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4;
421: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S;
422: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L;
423: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U;
424: PETSC_EXTERN PetscLogEvent KSP_SolveTranspose;
425: PETSC_EXTERN PetscLogEvent KSP_MatSolve;
427: PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*);
428: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*);
430: /*MC
431: KSPCheckDot - Checks if the result of a dot product used by the corresponding KSP contains Inf or NaN. These indicate that the previous
432: application of the preconditioner generated an error
434: Collective on ksp
436: Input Parameter:
437: . ksp - the linear solver (KSP) context.
439: Output Parameter:
440: . beta - the result of the inner product
442: Level: developer
444: Developer Note:
445: this is used to manage returning from KSP solvers whose preconditioners have failed in some way
447: .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckNorm(), KSPCheckSolve()
448: M*/
449: #define KSPCheckDot(ksp,beta) do { \
450: if (PetscIsInfOrNanScalar(beta)) { \
452: else {\
453: PCFailedReason pcreason;\
454: PetscInt sendbuf,recvbuf; \
455: PCGetFailedReasonRank(ksp->pc,&pcreason);\
456: sendbuf = (PetscInt)pcreason; \
457: MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));\
458: if (recvbuf) { \
459: PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf); \
460: ksp->reason = KSP_DIVERGED_PC_FAILED;\
461: VecSetInf(ksp->vec_sol);\
462: } else {\
463: ksp->reason = KSP_DIVERGED_NANORINF;\
464: }\
465: return 0;\
466: }\
467: } } while (0)
469: /*MC
470: KSPCheckNorm - Checks if the result of a norm used by the corresponding KSP contains Inf or NaN. These indicate that the previous
471: application of the preconditioner generated an error
473: Collective on ksp
475: Input Parameter:
476: . ksp - the linear solver (KSP) context.
478: Output Parameter:
479: . beta - the result of the norm
481: Level: developer
483: Developer Note:
484: this is used to manage returning from KSP solvers whose preconditioners have failed in some way
486: .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckDot(), KSPCheckSolve()
487: M*/
488: #define KSPCheckNorm(ksp,beta) do { \
489: if (PetscIsInfOrNanReal(beta)) { \
491: else {\
492: PCFailedReason pcreason;\
493: PetscInt sendbuf,recvbuf; \
494: PCGetFailedReasonRank(ksp->pc,&pcreason);\
495: sendbuf = (PetscInt)pcreason; \
496: MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));\
497: if (recvbuf) { \
498: PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf); \
499: ksp->reason = KSP_DIVERGED_PC_FAILED; \
500: VecSetInf(ksp->vec_sol);\
501: ksp->rnorm = beta; \
502: } else {\
503: PCSetFailedReason(ksp->pc,PC_NOERROR); \
504: ksp->reason = KSP_DIVERGED_NANORINF;\
505: ksp->rnorm = beta; \
506: } \
507: return 0;\
508: }\
509: } } while (0)
511: #endif
513: PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]);
514: PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP,PetscInt,PetscReal*);