Actual source code: kspimpl.h
petsc-3.6.2 2015-10-02
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)(PetscOptions*,KSP);
24: PetscErrorCode (*publishoptions)(KSP);
25: PetscErrorCode (*computeextremesingularvalues)(KSP,PetscReal*,PetscReal*);
26: PetscErrorCode (*computeeigenvalues)(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt *);
27: PetscErrorCode (*destroy)(KSP);
28: PetscErrorCode (*view)(KSP,PetscViewer);
29: PetscErrorCode (*reset)(KSP);
30: PetscErrorCode (*load)(KSP,PetscViewer);
31: };
33: typedef struct {PetscInt model,curl,maxl;Mat mat; KSP ksp;}* KSPGuessFischer;
35: /*
36: Maximum number of monitors you can run with a single KSP
37: */
38: #define MAXKSPMONITORS 5
39: typedef enum {KSP_SETUP_NEW, KSP_SETUP_NEWMATRIX, KSP_SETUP_NEWRHS} KSPSetUpStage;
41: /*
42: Defines the KSP data structure.
43: */
44: struct _p_KSP {
45: PETSCHEADER(struct _KSPOps);
46: DM dm;
47: PetscBool dmAuto; /* DM was created automatically by KSP */
48: PetscBool dmActive; /* KSP should use DM for computing operators */
49: /*------------------------- User parameters--------------------------*/
50: PetscInt max_it; /* maximum number of iterations */
51: KSPFischerGuess guess;
52: PetscBool guess_zero, /* flag for whether initial guess is 0 */
53: calc_sings, /* calculate extreme Singular Values */
54: guess_knoll; /* use initial guess of PCApply(ksp->B,b */
55: PCSide pc_side; /* flag for left, right, or symmetric preconditioning */
56: PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */
57: PetscReal rtol, /* relative tolerance */
58: abstol, /* absolute tolerance */
59: ttol, /* (not set by user) */
60: divtol; /* divergence tolerance */
61: PetscReal rnorm0; /* initial residual norm (used for divergence testing) */
62: PetscReal rnorm; /* current residual norm */
63: KSPConvergedReason reason;
64: PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */
66: Vec vec_sol,vec_rhs; /* pointer to where user has stashed
67: the solution and rhs, these are
68: never touched by the code, only
69: passed back to the user */
70: PetscReal *res_hist; /* If !0 stores residual at iterations*/
71: PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
72: PetscInt res_hist_len; /* current size of residual history array */
73: PetscInt res_hist_max; /* actual amount of data in residual_history */
74: PetscBool res_hist_reset; /* reset history to size zero for each new solve */
76: PetscInt chknorm; /* only compute/check norm if iterations is great than this */
77: PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the
78: MPI_Allreduce() for computing the inner products for the next iteration. */
79: /* --------User (or default) routines (most return -1 on error) --------*/
80: PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */
81: PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void**); /* */
82: void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */
83: PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */
85: PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
86: PetscErrorCode (*convergeddestroy)(void*);
87: void *cnvP;
89: void *user; /* optional user-defined context */
91: PC pc;
93: void *data; /* holder for misc stuff associated
94: with a particular iterative solver */
96: /* ----------------Default work-area management -------------------- */
97: PetscInt nwork;
98: Vec *work;
100: KSPSetUpStage setupstage;
102: PetscInt its; /* number of iterations so far computed in THIS linear solve*/
103: PetscInt totalits; /* number of iterations used by this KSP object since it was created */
105: PetscBool transpose_solve; /* solve transpose system instead */
107: KSPNormType normtype; /* type of norm used for convergence tests */
109: PCSide pc_side_set; /* PC type set explicitly by user */
110: KSPNormType normtype_set; /* Norm type set explicitly by user */
112: /* Allow diagonally scaling the matrix before computing the preconditioner or using
113: the Krylov method. Note this is NOT just Jacobi preconditioning */
115: PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */
116: PetscBool dscalefix; /* unscale system after solve */
117: PetscBool dscalefix2; /* system has been unscaled */
118: Vec diagonal; /* 1/sqrt(diag of matrix) */
119: Vec truediagonal;
121: PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */
123: PetscViewer eigviewer; /* Viewer where computed eigenvalues are displayed */
125: PetscErrorCode (*presolve)(KSP,Vec,Vec,void*);
126: PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*);
127: void *prectx,*postctx;
128: };
130: typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
131: PetscReal coef;
132: PetscReal bnrm;
133: } KSPDynTolCtx;
135: typedef struct {
136: PetscBool initialrtol; /* default relative residual decrease is computing from initial residual, not rhs */
137: PetscBool mininitialrtol; /* default relative residual decrease is computing from min of initial residual and rhs */
138: Vec work;
139: } KSPConvergedDefaultCtx;
143: PETSC_STATIC_INLINE PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm)
144: {
148: PetscObjectSAWsTakeAccess((PetscObject)ksp);
149: if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) {
150: ksp->res_hist[ksp->res_hist_len++] = norm;
151: }
152: PetscObjectSAWsGrantAccess((PetscObject)ksp);
153: return(0);
154: }
156: PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP,KSPNormType*,PCSide*);
158: PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP,PetscInt,const PetscReal*,const PetscReal*);
160: typedef struct _p_DMKSP *DMKSP;
161: typedef struct _DMKSPOps *DMKSPOps;
162: struct _DMKSPOps {
163: PetscErrorCode (*computeoperators)(KSP,Mat,Mat,void*);
164: PetscErrorCode (*computerhs)(KSP,Vec,void*);
165: PetscErrorCode (*computeinitialguess)(KSP,Vec,void*);
166: PetscErrorCode (*destroy)(DMKSP*);
167: PetscErrorCode (*duplicate)(DMKSP,DMKSP);
168: };
170: struct _p_DMKSP {
171: PETSCHEADER(struct _DMKSPOps);
172: void *operatorsctx;
173: void *rhsctx;
174: void *initialguessctx;
175: void *data;
177: /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
178: * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user
179: * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by
180: * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates
181: * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes
182: * to the original level will no longer propagate to that level.
183: */
184: DM originaldm;
186: void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */
187: };
188: PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM,DMKSP*);
189: PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM,DMKSP*);
190: PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM,DM);
192: /*
193: These allow the various Krylov methods to apply to either the linear system or its transpose.
194: */
197: PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpace(KSP ksp,Vec y)
198: {
201: if (ksp->pc_side == PC_LEFT) {
202: Mat A;
203: MatNullSpace nullsp;
204: PCGetOperators(ksp->pc,&A,NULL);
205: MatGetNullSpace(A,&nullsp);
206: if (nullsp) {
207: MatNullSpaceRemove(nullsp,y);
208: }
209: }
210: return(0);
211: }
215: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y)
216: {
219: if (!ksp->transpose_solve) {MatMult(A,x,y);}
220: else {MatMultTranspose(A,x,y);}
221: return(0);
222: }
226: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y)
227: {
230: if (!ksp->transpose_solve) {MatMultTranspose(A,x,y);}
231: else {MatMult(A,x,y);}
232: return(0);
233: }
237: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y)
238: {
241: if (!ksp->transpose_solve) {
242: PCApply(ksp->pc,x,y);
243: KSP_RemoveNullSpace(ksp,y);
244: } else {
245: PCApplyTranspose(ksp->pc,x,y);
246: }
247: return(0);
248: }
252: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y)
253: {
256: if (!ksp->transpose_solve) {
257: PCApplyTranspose(ksp->pc,x,y);
258: } else {
259: PCApply(ksp->pc,x,y);
260: KSP_RemoveNullSpace(ksp,y);
261: }
262: return(0);
263: }
267: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w)
268: {
271: if (!ksp->transpose_solve) {
272: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
273: KSP_RemoveNullSpace(ksp,y);
274: } else {
275: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
276: }
277: return(0);
278: }
282: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w)
283: {
286: if (!ksp->transpose_solve) {
287: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
288: KSP_RemoveNullSpace(ksp,y);
289: } else {
290: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
291: }
292: return(0);
293: }
295: PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization, KSP_SetUp, KSP_Solve;
297: PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*);
299: /*
300: Either generate an error or mark as diverged when a scalar from an inner product is Nan or Inf
301: */
302: #define KSPCheckDot(ksp,beta) \
303: if (PetscIsInfOrNanScalar(beta)) { \
304: if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");\
305: else {\
306: ksp->reason = KSP_DIVERGED_NANORINF;\
307: return(0);\
308: }\
309: }
311: /*
312: Either generate an error or mark as diverged when a real from a norm is Nan or Inf
313: */
314: #define KSPCheckNorm(ksp,beta) \
315: if (PetscIsInfOrNanReal(beta)) { \
316: if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf norm");\
317: else {\
318: ksp->reason = KSP_DIVERGED_NANORINF;\
319: return(0);\
320: }\
321: }
323: #endif