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