Actual source code: kspimpl.h
petsc-3.7.6 2017-04-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 KSPMatRegisterAll(void);
12: typedef struct _KSPOps *KSPOps;
14: struct _KSPOps {
15: PetscErrorCode (*buildsolution)(KSP,Vec,Vec*); /* Returns a pointer to the solution, or
16: calculates the solution in a
17: user-provided area. */
18: PetscErrorCode (*buildresidual)(KSP,Vec,Vec,Vec*); /* Returns a pointer to the residual, or
19: calculates the residual in a
20: user-provided area. */
21: PetscErrorCode (*solve)(KSP); /* actual solver */
22: PetscErrorCode (*setup)(KSP);
23: PetscErrorCode (*setfromoptions)(PetscOptionItems*,KSP);
24: PetscErrorCode (*publishoptions)(KSP);
25: PetscErrorCode (*computeextremesingularvalues)(KSP,PetscReal*,PetscReal*);
26: PetscErrorCode (*computeeigenvalues)(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt *);
27: PetscErrorCode (*computeritz)(KSP,PetscBool,PetscBool,PetscInt*,Vec[],PetscReal*,PetscReal*);
28: PetscErrorCode (*destroy)(KSP);
29: PetscErrorCode (*view)(KSP,PetscViewer);
30: PetscErrorCode (*reset)(KSP);
31: PetscErrorCode (*load)(KSP,PetscViewer);
32: };
34: typedef struct {PetscInt model,curl,maxl;Mat mat; KSP ksp;}* KSPGuessFischer;
36: /*
37: Maximum number of monitors you can run with a single KSP
38: */
39: #define MAXKSPMONITORS 5
40: typedef enum {KSP_SETUP_NEW, KSP_SETUP_NEWMATRIX, KSP_SETUP_NEWRHS} KSPSetUpStage;
42: /*
43: Defines the KSP data structure.
44: */
45: struct _p_KSP {
46: PETSCHEADER(struct _KSPOps);
47: DM dm;
48: PetscBool dmAuto; /* DM was created automatically by KSP */
49: PetscBool dmActive; /* KSP should use DM for computing operators */
50: /*------------------------- User parameters--------------------------*/
51: PetscInt max_it; /* maximum number of iterations */
52: KSPFischerGuess guess;
53: PetscBool guess_zero, /* flag for whether initial guess is 0 */
54: calc_sings, /* calculate extreme Singular Values */
55: calc_ritz, /* calculate (harmonic) Ritz pairs */
56: guess_knoll; /* use initial guess of PCApply(ksp->B,b */
57: PCSide pc_side; /* flag for left, right, or symmetric preconditioning */
58: PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */
59: PetscReal rtol, /* relative tolerance */
60: abstol, /* absolute tolerance */
61: ttol, /* (not set by user) */
62: divtol; /* divergence tolerance */
63: PetscReal rnorm0; /* initial residual norm (used for divergence testing) */
64: PetscReal rnorm; /* current residual norm */
65: KSPConvergedReason reason;
66: PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */
68: Vec vec_sol,vec_rhs; /* pointer to where user has stashed
69: the solution and rhs, these are
70: never touched by the code, only
71: passed back to the user */
72: PetscReal *res_hist; /* If !0 stores residual at iterations*/
73: PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
74: PetscInt res_hist_len; /* current size of residual history array */
75: PetscInt res_hist_max; /* actual amount of data in residual_history */
76: PetscBool res_hist_reset; /* reset history to size zero for each new solve */
78: PetscInt chknorm; /* only compute/check norm if iterations is great than this */
79: PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the
80: MPI_Allreduce() for computing the inner products for the next iteration. */
81: /* --------User (or default) routines (most return -1 on error) --------*/
82: PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */
83: PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void**); /* */
84: void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */
85: PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */
87: PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
88: PetscErrorCode (*convergeddestroy)(void*);
89: void *cnvP;
91: void *user; /* optional user-defined context */
93: PC pc;
95: void *data; /* holder for misc stuff associated
96: with a particular iterative solver */
98: /* ----------------Default work-area management -------------------- */
99: PetscInt nwork;
100: Vec *work;
102: KSPSetUpStage setupstage;
104: PetscInt its; /* number of iterations so far computed in THIS linear solve*/
105: PetscInt totalits; /* number of iterations used by this KSP object since it was created */
107: PetscBool transpose_solve; /* solve transpose system instead */
109: KSPNormType normtype; /* type of norm used for convergence tests */
111: PCSide pc_side_set; /* PC type set explicitly by user */
112: KSPNormType normtype_set; /* Norm type set explicitly by user */
114: /* Allow diagonally scaling the matrix before computing the preconditioner or using
115: the Krylov method. Note this is NOT just Jacobi preconditioning */
117: PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */
118: PetscBool dscalefix; /* unscale system after solve */
119: PetscBool dscalefix2; /* system has been unscaled */
120: Vec diagonal; /* 1/sqrt(diag of matrix) */
121: Vec truediagonal;
123: PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */
125: PetscViewer eigviewer; /* Viewer where computed eigenvalues are displayed */
127: PetscErrorCode (*presolve)(KSP,Vec,Vec,void*);
128: PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*);
129: void *prectx,*postctx;
130: };
132: typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
133: PetscReal coef;
134: PetscReal bnrm;
135: } KSPDynTolCtx;
137: typedef struct {
138: PetscBool initialrtol; /* default relative residual decrease is computing from initial residual, not rhs */
139: PetscBool mininitialrtol; /* default relative residual decrease is computing from min of initial residual and rhs */
140: Vec work;
141: } KSPConvergedDefaultCtx;
145: PETSC_STATIC_INLINE PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm)
146: {
150: PetscObjectSAWsTakeAccess((PetscObject)ksp);
151: if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) {
152: ksp->res_hist[ksp->res_hist_len++] = norm;
153: }
154: PetscObjectSAWsGrantAccess((PetscObject)ksp);
155: return(0);
156: }
158: PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP,KSPNormType*,PCSide*);
160: PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP,PetscInt,const PetscReal*,const PetscReal*);
162: typedef struct _p_DMKSP *DMKSP;
163: typedef struct _DMKSPOps *DMKSPOps;
164: struct _DMKSPOps {
165: PetscErrorCode (*computeoperators)(KSP,Mat,Mat,void*);
166: PetscErrorCode (*computerhs)(KSP,Vec,void*);
167: PetscErrorCode (*computeinitialguess)(KSP,Vec,void*);
168: PetscErrorCode (*destroy)(DMKSP*);
169: PetscErrorCode (*duplicate)(DMKSP,DMKSP);
170: };
172: struct _p_DMKSP {
173: PETSCHEADER(struct _DMKSPOps);
174: void *operatorsctx;
175: void *rhsctx;
176: void *initialguessctx;
177: void *data;
179: /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
180: * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user
181: * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by
182: * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates
183: * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes
184: * to the original level will no longer propagate to that level.
185: */
186: DM originaldm;
188: void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */
189: };
190: PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM,DMKSP*);
191: PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM,DMKSP*);
192: PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM,DM);
194: /*
195: These allow the various Krylov methods to apply to either the linear system or its transpose.
196: */
199: PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpace(KSP ksp,Vec y)
200: {
203: if (ksp->pc_side == PC_LEFT) {
204: Mat A;
205: MatNullSpace nullsp;
206: PCGetOperators(ksp->pc,&A,NULL);
207: MatGetNullSpace(A,&nullsp);
208: if (nullsp) {
209: MatNullSpaceRemove(nullsp,y);
210: }
211: }
212: return(0);
213: }
217: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y)
218: {
221: if (!ksp->transpose_solve) {MatMult(A,x,y);}
222: else {MatMultTranspose(A,x,y);}
223: return(0);
224: }
228: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y)
229: {
232: if (!ksp->transpose_solve) {MatMultTranspose(A,x,y);}
233: else {MatMult(A,x,y);}
234: return(0);
235: }
239: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y)
240: {
243: if (!ksp->transpose_solve) {
244: PCApply(ksp->pc,x,y);
245: KSP_RemoveNullSpace(ksp,y);
246: } else {
247: PCApplyTranspose(ksp->pc,x,y);
248: }
249: return(0);
250: }
254: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y)
255: {
258: if (!ksp->transpose_solve) {
259: PCApplyTranspose(ksp->pc,x,y);
260: } else {
261: PCApply(ksp->pc,x,y);
262: KSP_RemoveNullSpace(ksp,y);
263: }
264: return(0);
265: }
269: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w)
270: {
273: if (!ksp->transpose_solve) {
274: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
275: KSP_RemoveNullSpace(ksp,y);
276: } else {
277: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
278: }
279: return(0);
280: }
284: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w)
285: {
288: if (!ksp->transpose_solve) {
289: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
290: KSP_RemoveNullSpace(ksp,y);
291: } else {
292: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
293: }
294: return(0);
295: }
297: PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization, KSP_SetUp, KSP_Solve;
298: 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;
300: PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*);
301: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*);
303: /*
304: Either generate an error or mark as diverged when a scalar from an inner product is Nan or Inf
305: */
306: #define KSPCheckDot(ksp,beta) \
307: if (PetscIsInfOrNanScalar(beta)) { \
308: if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");\
309: else {\
311: PCFailedReason pcreason;\
312: PetscInt sendbuf,pcreason_max; \
313: PCGetSetUpFailedReason(ksp->pc,&pcreason);\
314: sendbuf = (PetscInt)pcreason; \
315: MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPIU_MAX,PetscObjectComm((PetscObject)ksp)); \
316: if (pcreason_max) {\
317: ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;\
318: VecSetInf(ksp->vec_sol);\
319: } else {\
320: ksp->reason = KSP_DIVERGED_NANORINF;\
321: }\
322: return(0);\
323: }\
324: }
326: /*
327: Either generate an error or mark as diverged when a real from a norm is Nan or Inf
328: */
329: #define KSPCheckNorm(ksp,beta) \
330: if (PetscIsInfOrNanReal(beta)) { \
331: if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf norm");\
332: else {\
334: PCFailedReason pcreason;\
335: PetscInt sendbuf,pcreason_max; \
336: PCGetSetUpFailedReason(ksp->pc,&pcreason);\
337: sendbuf = (PetscInt)pcreason; \
338: MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPIU_MAX,PetscObjectComm((PetscObject)ksp)); \
339: if (pcreason_max) {\
340: ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;\
341: VecSetInf(ksp->vec_sol);\
342: } else {\
343: ksp->reason = KSP_DIVERGED_NANORINF;\
344: }\
345: return(0);\
346: }\
347: }
349: #endif