Actual source code: petsc-kspimpl.h
petsc-3.5.4 2015-05-23
2: #ifndef _KSPIMPL_H
3: #define _KSPIMPL_H
5: #include <petscksp.h>
6: #include <petsc-private/petscimpl.h>
8: typedef struct _KSPOps *KSPOps;
10: struct _KSPOps {
11: PetscErrorCode (*buildsolution)(KSP,Vec,Vec*); /* Returns a pointer to the solution, or
12: calculates the solution in a
13: user-provided area. */
14: PetscErrorCode (*buildresidual)(KSP,Vec,Vec,Vec*); /* Returns a pointer to the residual, or
15: calculates the residual in a
16: user-provided area. */
17: PetscErrorCode (*solve)(KSP); /* actual solver */
18: PetscErrorCode (*setup)(KSP);
19: PetscErrorCode (*setfromoptions)(KSP);
20: PetscErrorCode (*publishoptions)(KSP);
21: PetscErrorCode (*computeextremesingularvalues)(KSP,PetscReal*,PetscReal*);
22: PetscErrorCode (*computeeigenvalues)(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt *);
23: PetscErrorCode (*destroy)(KSP);
24: PetscErrorCode (*view)(KSP,PetscViewer);
25: PetscErrorCode (*reset)(KSP);
26: PetscErrorCode (*load)(KSP,PetscViewer);
27: };
29: typedef struct {PetscInt model,curl,maxl;Mat mat; KSP ksp;}* KSPGuessFischer;
31: /*
32: Maximum number of monitors you can run with a single KSP
33: */
34: #define MAXKSPMONITORS 5
35: typedef enum {KSP_SETUP_NEW, KSP_SETUP_NEWMATRIX, KSP_SETUP_NEWRHS} KSPSetUpStage;
37: /*
38: Defines the KSP data structure.
39: */
40: struct _p_KSP {
41: PETSCHEADER(struct _KSPOps);
42: DM dm;
43: PetscBool dmAuto; /* DM was created automatically by KSP */
44: PetscBool dmActive; /* KSP should use DM for computing operators */
45: /*------------------------- User parameters--------------------------*/
46: PetscInt max_it; /* maximum number of iterations */
47: KSPFischerGuess guess;
48: PetscBool guess_zero, /* flag for whether initial guess is 0 */
49: calc_sings, /* calculate extreme Singular Values */
50: guess_knoll; /* use initial guess of PCApply(ksp->B,b */
51: PCSide pc_side; /* flag for left, right, or symmetric preconditioning */
52: PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */
53: PetscReal rtol, /* relative tolerance */
54: abstol, /* absolute tolerance */
55: ttol, /* (not set by user) */
56: divtol; /* divergence tolerance */
57: PetscReal rnorm0; /* initial residual norm (used for divergence testing) */
58: PetscReal rnorm; /* current residual norm */
59: KSPConvergedReason reason;
60: PetscBool printreason; /* prints converged reason after solve */
61: PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */
63: Vec vec_sol,vec_rhs; /* pointer to where user has stashed
64: the solution and rhs, these are
65: never touched by the code, only
66: passed back to the user */
67: PetscReal *res_hist; /* If !0 stores residual at iterations*/
68: PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
69: PetscInt res_hist_len; /* current size of residual history array */
70: PetscInt res_hist_max; /* actual amount of data in residual_history */
71: PetscBool res_hist_reset; /* reset history to size zero for each new solve */
73: PetscInt chknorm; /* only compute/check norm if iterations is great than this */
74: PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the
75: MPI_Allreduce() for computing the inner products for the next iteration. */
76: /* --------User (or default) routines (most return -1 on error) --------*/
77: PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */
78: PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void**); /* */
79: void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */
80: PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */
82: PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
83: PetscErrorCode (*convergeddestroy)(void*);
84: void *cnvP;
86: void *user; /* optional user-defined context */
88: PC pc;
90: void *data; /* holder for misc stuff associated
91: with a particular iterative solver */
93: /* ----------------Default work-area management -------------------- */
94: PetscInt nwork;
95: Vec *work;
97: KSPSetUpStage setupstage;
99: PetscInt its; /* number of iterations so far computed */
101: PetscBool transpose_solve; /* solve transpose system instead */
103: KSPNormType normtype; /* type of norm used for convergence tests */
105: PCSide pc_side_set; /* PC type set explicitly by user */
106: KSPNormType normtype_set; /* Norm type set explicitly by user */
108: /* Allow diagonally scaling the matrix before computing the preconditioner or using
109: the Krylov method. Note this is NOT just Jacobi preconditioning */
111: PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */
112: PetscBool dscalefix; /* unscale system after solve */
113: PetscBool dscalefix2; /* system has been unscaled */
114: Vec diagonal; /* 1/sqrt(diag of matrix) */
115: Vec truediagonal;
117: MatNullSpace nullsp; /* Null space of the operator, removed from Krylov space */
119: PetscViewer eigviewer; /* Viewer where computed eigenvalues are displayed */
121: PetscErrorCode (*presolve)(KSP,Vec,Vec,void*);
122: PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*);
123: void *prectx,*postctx;
124: };
126: typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
127: PetscReal coef;
128: PetscReal bnrm;
129: } KSPDynTolCtx;
131: typedef struct {
132: PetscBool initialrtol; /* default relative residual decrease is computing from initial residual, not rhs */
133: PetscBool mininitialrtol; /* default relative residual decrease is computing from min of initial residual and rhs */
134: Vec work;
135: } KSPConvergedDefaultCtx;
139: PETSC_STATIC_INLINE PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm)
140: {
144: PetscObjectSAWsTakeAccess((PetscObject)ksp);
145: if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) {
146: ksp->res_hist[ksp->res_hist_len++] = norm;
147: }
148: PetscObjectSAWsGrantAccess((PetscObject)ksp);
149: return(0);
150: }
152: PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP,KSPNormType*,PCSide*);
154: PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP,PetscInt,const PetscReal*,const PetscReal*);
156: typedef struct _p_DMKSP *DMKSP;
157: typedef struct _DMKSPOps *DMKSPOps;
158: struct _DMKSPOps {
159: PetscErrorCode (*computeoperators)(KSP,Mat,Mat,void*);
160: PetscErrorCode (*computerhs)(KSP,Vec,void*);
161: PetscErrorCode (*computeinitialguess)(KSP,Vec,void*);
162: PetscErrorCode (*destroy)(DMKSP*);
163: PetscErrorCode (*duplicate)(DMKSP,DMKSP);
164: };
166: struct _p_DMKSP {
167: PETSCHEADER(struct _DMKSPOps);
168: void *operatorsctx;
169: void *rhsctx;
170: void *initialguessctx;
171: void *data;
173: /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way
174: * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user
175: * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by
176: * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates
177: * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes
178: * to the original level will no longer propagate to that level.
179: */
180: DM originaldm;
182: void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */
183: };
184: PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM,DMKSP*);
185: PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM,DMKSP*);
186: PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM,DM);
188: /*
189: These allow the various Krylov methods to apply to either the linear system or its transpose.
190: */
193: PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpace(KSP ksp,Vec y)
194: {
197: if (ksp->nullsp && ksp->pc_side == PC_LEFT) {MatNullSpaceRemove(ksp->nullsp,y);}
198: return(0);
199: }
203: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y)
204: {
207: if (!ksp->transpose_solve) {MatMult(A,x,y);}
208: else {MatMultTranspose(A,x,y);}
209: return(0);
210: }
214: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y)
215: {
218: if (!ksp->transpose_solve) {MatMultTranspose(A,x,y);}
219: else {MatMult(A,x,y);}
220: return(0);
221: }
225: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y)
226: {
229: if (!ksp->transpose_solve) {
230: PCApply(ksp->pc,x,y);
231: KSP_RemoveNullSpace(ksp,y);
232: } else {
233: PCApplyTranspose(ksp->pc,x,y);
234: }
235: return(0);
236: }
240: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y)
241: {
244: if (!ksp->transpose_solve) {
245: PCApplyTranspose(ksp->pc,x,y);
246: } else {
247: PCApply(ksp->pc,x,y);
248: KSP_RemoveNullSpace(ksp,y);
249: }
250: return(0);
251: }
255: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w)
256: {
259: if (!ksp->transpose_solve) {
260: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
261: KSP_RemoveNullSpace(ksp,y);
262: } else {
263: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
264: }
265: return(0);
266: }
270: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w)
271: {
274: if (!ksp->transpose_solve) {
275: PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
276: KSP_RemoveNullSpace(ksp,y);
277: } else {
278: PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
279: }
280: return(0);
281: }
283: PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization, KSP_SetUp, KSP_Solve;
285: PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*);
287: #endif