Actual source code: kspimpl.h
petsc-3.8.4 2018-03-24
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: void *data; /* pointer to the specific implementation */
55: };
57: PETSC_EXTERN PetscErrorCode KSPGuessCreate_Fischer(KSPGuess);
58: PETSC_EXTERN PetscErrorCode KSPGuessCreate_POD(KSPGuess);
60: /*
61: Maximum number of monitors you can run with a single KSP
62: */
63: #define MAXKSPMONITORS 5
64: typedef enum {KSP_SETUP_NEW, KSP_SETUP_NEWMATRIX, KSP_SETUP_NEWRHS} KSPSetUpStage;
66: /*
67: Defines the KSP data structure.
68: */
69: struct _p_KSP {
70: PETSCHEADER(struct _KSPOps);
71: DM dm;
72: PetscBool dmAuto; /* DM was created automatically by KSP */
73: PetscBool dmActive; /* KSP should use DM for computing operators */
74: /*------------------------- User parameters--------------------------*/
75: PetscInt max_it; /* maximum number of iterations */
76: KSPGuess guess;
77: PetscBool guess_zero, /* flag for whether initial guess is 0 */
78: calc_sings, /* calculate extreme Singular Values */
79: calc_ritz, /* calculate (harmonic) Ritz pairs */
80: guess_knoll; /* use initial guess of PCApply(ksp->B,b */
81: PCSide pc_side; /* flag for left, right, or symmetric preconditioning */
82: PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */
83: PetscReal rtol, /* relative tolerance */
84: abstol, /* absolute tolerance */
85: ttol, /* (not set by user) */
86: divtol; /* divergence tolerance */
87: PetscReal rnorm0; /* initial residual norm (used for divergence testing) */
88: PetscReal rnorm; /* current residual norm */
89: KSPConvergedReason reason;
90: PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */
92: Vec vec_sol,vec_rhs; /* pointer to where user has stashed
93: the solution and rhs, these are
94: never touched by the code, only
95: passed back to the user */
96: PetscReal *res_hist; /* If !0 stores residual at iterations*/
97: PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
98: PetscInt res_hist_len; /* current size of residual history array */
99: PetscInt res_hist_max; /* actual amount of data in residual_history */
100: PetscBool res_hist_reset; /* reset history to size zero for each new solve */
102: PetscInt chknorm; /* only compute/check norm if iterations is great than this */
103: PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the
104: MPI_Allreduce() for computing the inner products for the next iteration. */
105: /* --------User (or default) routines (most return -1 on error) --------*/
106: PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */
107: PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void**); /* */
108: void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */
109: PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */
111: PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
112: PetscErrorCode (*convergeddestroy)(void*);
113: void *cnvP;
115: void *user; /* optional user-defined context */
117: PC pc;
119: void *data; /* holder for misc stuff associated
120: with a particular iterative solver */
122: /* ----------------Default work-area management -------------------- */
123: PetscInt nwork;
124: Vec *work;
126: KSPSetUpStage setupstage;
127: PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */
129: PetscInt its; /* number of iterations so far computed in THIS linear solve*/
130: PetscInt totalits; /* number of iterations used by this KSP object since it was created */
132: PetscBool transpose_solve; /* solve transpose system instead */
134: KSPNormType normtype; /* type of norm used for convergence tests */
136: PCSide pc_side_set; /* PC type set explicitly by user */
137: KSPNormType normtype_set; /* Norm type set explicitly by user */
139: /* Allow diagonally scaling the matrix before computing the preconditioner or using
140: the Krylov method. Note this is NOT just Jacobi preconditioning */
142: PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */
143: PetscBool dscalefix; /* unscale system after solve */
144: PetscBool dscalefix2; /* system has been unscaled */
145: Vec diagonal; /* 1/sqrt(diag of matrix) */
146: Vec truediagonal;
148: PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */
150: PetscViewer eigviewer; /* Viewer where computed eigenvalues are displayed */
152: PetscErrorCode (*presolve)(KSP,Vec,Vec,void*);
153: PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*);
154: void *prectx,*postctx;
155: };
157: typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
158: PetscReal coef;
159: PetscReal bnrm;
160: } KSPDynTolCtx;
162: typedef struct {
163: PetscBool initialrtol; /* default relative residual decrease is computing from initial residual, not rhs */
164: PetscBool mininitialrtol; /* default relative residual decrease is computing from min of initial residual and rhs */
165: Vec work;
166: } KSPConvergedDefaultCtx;
168: PETSC_STATIC_INLINE PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm)
169: {
173: PetscObjectSAWsTakeAccess((PetscObject)ksp);
174: if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) {
175: ksp->res_hist[ksp->res_hist_len++] = norm;
176: }
177: PetscObjectSAWsGrantAccess((PetscObject)ksp);
178: return(0);
179: }
181: PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP,PetscBool,KSPNormType*,PCSide*);
183: PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP,PetscInt,const PetscReal*,const PetscReal*);
185: typedef struct _p_DMKSP *DMKSP;
186: typedef struct _DMKSPOps *DMKSPOps;
187: struct _DMKSPOps {
188: PetscErrorCode (*computeoperators)(KSP,Mat,Mat,void*);
189: PetscErrorCode (*computerhs)(KSP,Vec,void*);
190: PetscErrorCode (*computeinitialguess)(KSP,Vec,void*);
191: PetscErrorCode (*destroy)(DMKSP*);
192: PetscErrorCode (*duplicate)(DMKSP,DMKSP);
193: };
195: struct _p_DMKSP {
196: PETSCHEADER(struct _DMKSPOps);
197: void *operatorsctx;
198: void *rhsctx;
199: void *initialguessctx;
200: void *data;
202: /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
203: * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user
204: * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by
205: * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates
206: * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes
207: * to the original level will no longer propagate to that level.
208: */
209: DM originaldm;
211: void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */
212: };
213: PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM,DMKSP*);
214: PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM,DMKSP*);
215: PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM,DM);
217: /*
218: These allow the various Krylov methods to apply to either the linear system or its transpose.
219: */
220: PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpace(KSP ksp,Vec y)
221: {
224: if (ksp->pc_side == PC_LEFT) {
225: Mat A;
226: MatNullSpace nullsp;
227: PCGetOperators(ksp->pc,&A,NULL);
228: MatGetNullSpace(A,&nullsp);
229: if (nullsp) {
230: MatNullSpaceRemove(nullsp,y);
231: }
232: }
233: return(0);
234: }
236: PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp,Vec y)
237: {
240: if (ksp->pc_side == PC_LEFT) {
241: Mat A;
242: MatNullSpace nullsp;
243: PCGetOperators(ksp->pc,&A,NULL);
244: MatGetTransposeNullSpace(A,&nullsp);
245: if (nullsp) {
246: MatNullSpaceRemove(nullsp,y);
247: }
248: }
249: return(0);
250: }
252: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y)
253: {
256: if (!ksp->transpose_solve) {MatMult(A,x,y);}
257: else {MatMultTranspose(A,x,y);}
258: return(0);
259: }
261: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y)
262: {
265: if (!ksp->transpose_solve) {MatMultTranspose(A,x,y);}
266: else {MatMult(A,x,y);}
267: return(0);
268: }
270: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y)
271: {
274: if (!ksp->transpose_solve) {
275: PCApply(ksp->pc,x,y);
276: KSP_RemoveNullSpace(ksp,y);
277: } else {
278: PCApplyTranspose(ksp->pc,x,y);
279: KSP_RemoveNullSpaceTranspose(ksp,y);
280: }
281: return(0);
282: }
284: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y)
285: {
288: if (!ksp->transpose_solve) {
289: PCApplyTranspose(ksp->pc,x,y);
290: KSP_RemoveNullSpaceTranspose(ksp,y);
291: } else {
292: PCApply(ksp->pc,x,y);
293: KSP_RemoveNullSpace(ksp,y);
294: }
295: return(0);
296: }
298: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w)
299: {
302: if (!ksp->transpose_solve) {
303: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
304: KSP_RemoveNullSpace(ksp,y);
305: } else {
306: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
307: KSP_RemoveNullSpaceTranspose(ksp,y);
308: }
309: return(0);
310: }
312: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w)
313: {
316: if (!ksp->transpose_solve) {
317: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
318: } else {
319: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
320: }
321: return(0);
322: }
324: PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization, KSP_SetUp, KSP_Solve;
325: PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0,KSP_Solve_FS_1,KSP_Solve_FS_2,KSP_Solve_FS_3,KSP_Solve_FS_4,KSP_Solve_FS_S,KSP_Solve_FS_L,KSP_Solve_FS_U;
327: PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*);
328: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*);
330: /*
331: Either generate an error or mark as diverged when a scalar from an inner product is Nan or Inf
332: */
333: #define KSPCheckDot(ksp,beta) \
334: if (PetscIsInfOrNanScalar(beta)) { \
335: if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");\
336: else {\
338: PCFailedReason pcreason;\
339: PetscInt sendbuf,pcreason_max; \
340: PCGetSetUpFailedReason(ksp->pc,&pcreason);\
341: sendbuf = (PetscInt)pcreason; \
342: MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp)); \
343: if (pcreason_max) {\
344: ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;\
345: VecSetInf(ksp->vec_sol);\
346: } else {\
347: ksp->reason = KSP_DIVERGED_NANORINF;\
348: }\
349: return(0);\
350: }\
351: }
353: /*
354: Either generate an error or mark as diverged when a real from a norm is Nan or Inf
355: */
356: #define KSPCheckNorm(ksp,beta) \
357: if (PetscIsInfOrNanReal(beta)) { \
358: if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf norm");\
359: else {\
361: PCFailedReason pcreason;\
362: PetscInt sendbuf,pcreason_max; \
363: PCGetSetUpFailedReason(ksp->pc,&pcreason);\
364: sendbuf = (PetscInt)pcreason; \
365: MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp)); \
366: if (pcreason_max) {\
367: ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;\
368: VecSetInf(ksp->vec_sol);\
369: } else {\
370: ksp->reason = KSP_DIVERGED_NANORINF;\
371: }\
372: return(0);\
373: }\
374: }
376: #endif