Actual source code: kspimpl.h
petsc-master 2019-10-18
2: #ifndef _KSPIMPL_H
3: #define _KSPIMPL_H
5: #include <petscksp.h>
6: #include <petsc/private/petscimpl.h>
8: PETSC_EXTERN PetscBool KSPRegisterAllCalled;
9: PETSC_EXTERN PetscErrorCode KSPRegisterAll(void);
10: PETSC_EXTERN PetscErrorCode KSPGuessRegisterAll(void);
11: PETSC_EXTERN PetscErrorCode KSPMatRegisterAll(void);
13: typedef struct _KSPOps *KSPOps;
15: struct _KSPOps {
16: PetscErrorCode (*buildsolution)(KSP,Vec,Vec*); /* Returns a pointer to the solution, or
17: calculates the solution in a
18: user-provided area. */
19: PetscErrorCode (*buildresidual)(KSP,Vec,Vec,Vec*); /* Returns a pointer to the residual, or
20: calculates the residual in a
21: user-provided area. */
22: PetscErrorCode (*solve)(KSP); /* actual solver */
23: PetscErrorCode (*setup)(KSP);
24: PetscErrorCode (*setfromoptions)(PetscOptionItems*,KSP);
25: PetscErrorCode (*publishoptions)(KSP);
26: PetscErrorCode (*computeextremesingularvalues)(KSP,PetscReal*,PetscReal*);
27: PetscErrorCode (*computeeigenvalues)(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt *);
28: PetscErrorCode (*computeritz)(KSP,PetscBool,PetscBool,PetscInt*,Vec[],PetscReal*,PetscReal*);
29: PetscErrorCode (*destroy)(KSP);
30: PetscErrorCode (*view)(KSP,PetscViewer);
31: PetscErrorCode (*reset)(KSP);
32: PetscErrorCode (*load)(KSP,PetscViewer);
33: };
35: typedef struct _KSPGuessOps *KSPGuessOps;
37: struct _KSPGuessOps {
38: PetscErrorCode (*formguess)(KSPGuess,Vec,Vec); /* Form initial guess */
39: PetscErrorCode (*update)(KSPGuess,Vec,Vec); /* Update database */
40: PetscErrorCode (*setfromoptions)(KSPGuess);
41: PetscErrorCode (*setup)(KSPGuess);
42: PetscErrorCode (*destroy)(KSPGuess);
43: PetscErrorCode (*view)(KSPGuess,PetscViewer);
44: PetscErrorCode (*reset)(KSPGuess);
45: };
47: /*
48: Defines the KSPGuess data structure.
49: */
50: struct _p_KSPGuess {
51: PETSCHEADER(struct _KSPGuessOps);
52: KSP ksp; /* the parent KSP */
53: Mat A; /* the current linear operator */
54: PetscObjectState omatstate; /* previous linear operator state */
55: void *data; /* pointer to the specific implementation */
56: };
58: PETSC_EXTERN PetscErrorCode KSPGuessCreate_Fischer(KSPGuess);
59: PETSC_EXTERN PetscErrorCode KSPGuessCreate_POD(KSPGuess);
61: /*
62: Maximum number of monitors you can run with a single KSP
63: */
64: #define MAXKSPMONITORS 5
65: typedef enum {KSP_SETUP_NEW, KSP_SETUP_NEWMATRIX, KSP_SETUP_NEWRHS} KSPSetUpStage;
67: /*
68: Defines the KSP data structure.
69: */
70: struct _p_KSP {
71: PETSCHEADER(struct _KSPOps);
72: DM dm;
73: PetscBool dmAuto; /* DM was created automatically by KSP */
74: PetscBool dmActive; /* KSP should use DM for computing operators */
75: /*------------------------- User parameters--------------------------*/
76: PetscInt max_it; /* maximum number of iterations */
77: KSPGuess guess;
78: PetscBool guess_zero, /* flag for whether initial guess is 0 */
79: calc_sings, /* calculate extreme Singular Values */
80: calc_ritz, /* calculate (harmonic) Ritz pairs */
81: guess_knoll; /* use initial guess of PCApply(ksp->B,b */
82: PCSide pc_side; /* flag for left, right, or symmetric preconditioning */
83: PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */
84: PetscReal rtol, /* relative tolerance */
85: abstol, /* absolute tolerance */
86: ttol, /* (not set by user) */
87: divtol; /* divergence tolerance */
88: PetscReal rnorm0; /* initial residual norm (used for divergence testing) */
89: PetscReal rnorm; /* current residual norm */
90: KSPConvergedReason reason;
91: PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */
93: Vec vec_sol,vec_rhs; /* pointer to where user has stashed
94: the solution and rhs, these are
95: never touched by the code, only
96: passed back to the user */
97: PetscReal *res_hist; /* If !0 stores residual at iterations*/
98: PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
99: PetscInt res_hist_len; /* current size of residual history array */
100: PetscInt res_hist_max; /* actual amount of data in residual_history */
101: PetscBool res_hist_reset; /* reset history to size zero for each new solve */
103: PetscInt chknorm; /* only compute/check norm if iterations is great than this */
104: PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the
105: MPI_Allreduce() for computing the inner products for the next iteration. */
106: /* --------User (or default) routines (most return -1 on error) --------*/
107: PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */
108: PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void**); /* */
109: void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */
110: PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */
112: PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
113: PetscErrorCode (*convergeddestroy)(void*);
114: void *cnvP;
116: void *user; /* optional user-defined context */
118: PC pc;
120: void *data; /* holder for misc stuff associated
121: with a particular iterative solver */
123: PetscBool view, viewPre, viewReason, viewMat, viewPMat, viewRhs, viewSol, viewMatExp, viewEV, viewSV, viewEVExp, viewFinalRes, viewPOpExp, viewDScale;
124: PetscViewer viewer, viewerPre, viewerReason, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp, viewerDScale;
125: PetscViewerFormat format, formatPre, formatReason, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp, formatDScale;
127: /* ----------------Default work-area management -------------------- */
128: PetscInt nwork;
129: Vec *work;
131: KSPSetUpStage setupstage;
132: PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */
134: PetscInt its; /* number of iterations so far computed in THIS linear solve*/
135: PetscInt totalits; /* number of iterations used by this KSP object since it was created */
137: PetscBool transpose_solve; /* solve transpose system instead */
139: KSPNormType normtype; /* type of norm used for convergence tests */
141: PCSide pc_side_set; /* PC type set explicitly by user */
142: KSPNormType normtype_set; /* Norm type set explicitly by user */
144: /* Allow diagonally scaling the matrix before computing the preconditioner or using
145: the Krylov method. Note this is NOT just Jacobi preconditioning */
147: PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */
148: PetscBool dscalefix; /* unscale system after solve */
149: PetscBool dscalefix2; /* system has been unscaled */
150: Vec diagonal; /* 1/sqrt(diag of matrix) */
151: Vec truediagonal;
153: PetscInt setfromoptionscalled;
154: PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */
156: PetscViewer eigviewer; /* Viewer where computed eigenvalues are displayed */
158: PetscErrorCode (*presolve)(KSP,Vec,Vec,void*);
159: PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*);
160: void *prectx,*postctx;
161: };
163: typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
164: PetscReal coef;
165: PetscReal bnrm;
166: } KSPDynTolCtx;
168: typedef struct {
169: PetscBool initialrtol; /* default relative residual decrease is computing from initial residual, not rhs */
170: PetscBool mininitialrtol; /* default relative residual decrease is computing from min of initial residual and rhs */
171: Vec work;
172: } KSPConvergedDefaultCtx;
174: PETSC_STATIC_INLINE PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm)
175: {
179: PetscObjectSAWsTakeAccess((PetscObject)ksp);
180: if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) {
181: ksp->res_hist[ksp->res_hist_len++] = norm;
182: }
183: PetscObjectSAWsGrantAccess((PetscObject)ksp);
184: return(0);
185: }
187: PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP,PetscBool,KSPNormType*,PCSide*);
189: PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP,PetscInt,const PetscReal*,const PetscReal*);
191: typedef struct _p_DMKSP *DMKSP;
192: typedef struct _DMKSPOps *DMKSPOps;
193: struct _DMKSPOps {
194: PetscErrorCode (*computeoperators)(KSP,Mat,Mat,void*);
195: PetscErrorCode (*computerhs)(KSP,Vec,void*);
196: PetscErrorCode (*computeinitialguess)(KSP,Vec,void*);
197: PetscErrorCode (*destroy)(DMKSP*);
198: PetscErrorCode (*duplicate)(DMKSP,DMKSP);
199: };
201: struct _p_DMKSP {
202: PETSCHEADER(struct _DMKSPOps);
203: void *operatorsctx;
204: void *rhsctx;
205: void *initialguessctx;
206: void *data;
208: /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
209: * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user
210: * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by
211: * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates
212: * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes
213: * to the original level will no longer propagate to that level.
214: */
215: DM originaldm;
217: void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */
218: };
219: PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM,DMKSP*);
220: PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM,DMKSP*);
221: PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM,DM);
223: /*
224: These allow the various Krylov methods to apply to either the linear system or its transpose.
225: */
226: PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpace(KSP ksp,Vec y)
227: {
230: if (ksp->pc_side == PC_LEFT) {
231: Mat A;
232: MatNullSpace nullsp;
233: PCGetOperators(ksp->pc,&A,NULL);
234: MatGetNullSpace(A,&nullsp);
235: if (nullsp) {
236: MatNullSpaceRemove(nullsp,y);
237: }
238: }
239: return(0);
240: }
242: PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp,Vec y)
243: {
246: if (ksp->pc_side == PC_LEFT) {
247: Mat A;
248: MatNullSpace nullsp;
249: PCGetOperators(ksp->pc,&A,NULL);
250: MatGetTransposeNullSpace(A,&nullsp);
251: if (nullsp) {
252: MatNullSpaceRemove(nullsp,y);
253: }
254: }
255: return(0);
256: }
258: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y)
259: {
262: if (!ksp->transpose_solve) {MatMult(A,x,y);}
263: else {MatMultTranspose(A,x,y);}
264: return(0);
265: }
267: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y)
268: {
271: if (!ksp->transpose_solve) {MatMultTranspose(A,x,y);}
272: else {MatMult(A,x,y);}
273: return(0);
274: }
276: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y)
277: {
280: if (!ksp->transpose_solve) {
281: PCApply(ksp->pc,x,y);
282: KSP_RemoveNullSpace(ksp,y);
283: } else {
284: PCApplyTranspose(ksp->pc,x,y);
285: KSP_RemoveNullSpaceTranspose(ksp,y);
286: }
287: return(0);
288: }
290: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y)
291: {
294: if (!ksp->transpose_solve) {
295: PCApplyTranspose(ksp->pc,x,y);
296: KSP_RemoveNullSpaceTranspose(ksp,y);
297: } else {
298: PCApply(ksp->pc,x,y);
299: KSP_RemoveNullSpace(ksp,y);
300: }
301: return(0);
302: }
304: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w)
305: {
308: if (!ksp->transpose_solve) {
309: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
310: KSP_RemoveNullSpace(ksp,y);
311: } else {
312: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
313: KSP_RemoveNullSpaceTranspose(ksp,y);
314: }
315: return(0);
316: }
318: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w)
319: {
322: if (!ksp->transpose_solve) {
323: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
324: } else {
325: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
326: }
327: return(0);
328: }
330: PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization;
331: PETSC_EXTERN PetscLogEvent KSP_SetUp;
332: PETSC_EXTERN PetscLogEvent KSP_Solve;
333: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0;
334: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1;
335: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2;
336: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3;
337: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4;
338: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S;
339: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L;
340: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U;
341: PETSC_EXTERN PetscLogEvent KSP_SolveTranspose;
343: PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*);
344: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*);
346: /*MC
347: KSPCheckDot - Checks if the result of a dot product used by the corresponding KSP contains Inf or NaN. These indicate that the previous
348: application of the preconditioner generated an error
350: Collective on ksp
352: Input Parameter:
353: . ksp - the linear solver (KSP) context.
355: Output Parameter:
356: . beta - the result of the inner product
358: Level: developer
360: Developer Note:
361: this is used to manage returning from KSP solvers whose preconditioners have failed in some way
363: .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckNorm(), KSPCheckSolve()
364: M*/
365: #define KSPCheckDot(ksp,beta) do { \
366: if (PetscIsInfOrNanScalar(beta)) { \
367: if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");\
368: else {\
370: PCFailedReason pcreason;\
371: PetscInt sendbuf,pcreason_max; \
372: PCGetFailedReason(ksp->pc,&pcreason);\
373: sendbuf = (PetscInt)pcreason; \
374: MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp)); \
375: if (pcreason_max) {\
376: ksp->reason = KSP_DIVERGED_PC_FAILED;\
377: VecSetInf(ksp->vec_sol);\
378: } else {\
379: ksp->reason = KSP_DIVERGED_NANORINF;\
380: }\
381: return(0);\
382: }\
383: } } while (0)
385: /*MC
386: KSPCheckNorm - Checks if the result of a norm used by the corresponding KSP contains Inf or NaN. These indicate that the previous
387: application of the preconditioner generated an error
389: Collective on ksp
391: Input Parameter:
392: . ksp - the linear solver (KSP) context.
394: Output Parameter:
395: . beta - the result of the norm
397: Level: developer
399: Developer Note:
400: this is used to manage returning from KSP solvers whose preconditioners have failed in some way
402: .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckDot(), KSPCheckSolve()
403: M*/
404: #define KSPCheckNorm(ksp,beta) do { \
405: if (PetscIsInfOrNanReal(beta)) { \
406: if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf norm");\
407: else {\
409: PCFailedReason pcreason;\
410: PetscInt sendbuf,pcreason_max; \
411: PCGetFailedReason(ksp->pc,&pcreason);\
412: sendbuf = (PetscInt)pcreason; \
413: MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp)); \
414: if (pcreason_max) {\
415: ksp->reason = KSP_DIVERGED_PC_FAILED;\
416: VecSetInf(ksp->vec_sol);\
417: } else {\
418: ksp->reason = KSP_DIVERGED_NANORINF;\
419: }\
420: return(0);\
421: }\
422: } } while (0)
424: #endif