Actual source code: kspimpl.h
2: #ifndef _KSPIMPL_H
3: #define _KSPIMPL_H
5: #include <petscksp.h>
7: typedef struct _KSPOps *KSPOps;
9: struct _KSPOps {
10: PetscErrorCode (*buildsolution)(KSP,Vec,Vec*); /* Returns a pointer to the solution, or
11: calculates the solution in a
12: user-provided area. */
13: PetscErrorCode (*buildresidual)(KSP,Vec,Vec,Vec*); /* Returns a pointer to the residual, or
14: calculates the residual in a
15: user-provided area. */
16: PetscErrorCode (*solve)(KSP); /* actual solver */
17: PetscErrorCode (*setup)(KSP);
18: PetscErrorCode (*setfromoptions)(KSP);
19: PetscErrorCode (*publishoptions)(KSP);
20: PetscErrorCode (*computeextremesingularvalues)(KSP,PetscReal*,PetscReal*);
21: PetscErrorCode (*computeeigenvalues)(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt *);
22: PetscErrorCode (*destroy)(KSP);
23: PetscErrorCode (*view)(KSP,PetscViewer);
24: PetscErrorCode (*reset)(KSP);
25: };
27: typedef struct {PetscInt model,curl,maxl;Mat mat; KSP ksp;}* KSPGuessFischer;
29: /*
30: Maximum number of monitors you can run with a single KSP
31: */
32: #define MAXKSPMONITORS 5
33: typedef enum {KSP_SETUP_NEW, KSP_SETUP_NEWMATRIX, KSP_SETUP_NEWRHS} KSPSetUpStage;
35: /*
36: Defines the KSP data structure.
37: */
38: struct _p_KSP {
39: PETSCHEADER(struct _KSPOps);
40: DM dm;
41: PetscBool dmActive;
42: /*------------------------- User parameters--------------------------*/
43: PetscInt max_it; /* maximum number of iterations */
44: KSPFischerGuess guess;
45: PetscBool guess_zero, /* flag for whether initial guess is 0 */
46: calc_sings, /* calculate extreme Singular Values */
47: guess_knoll; /* use initial guess of PCApply(ksp->B,b */
48: PCSide pc_side; /* flag for left, right, or symmetric preconditioning */
49: PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */
50: PetscReal rtol, /* relative tolerance */
51: abstol, /* absolute tolerance */
52: ttol, /* (not set by user) */
53: divtol; /* divergence tolerance */
54: PetscReal rnorm0; /* initial residual norm (used for divergence testing) */
55: PetscReal rnorm; /* current residual norm */
56: KSPConvergedReason reason;
57: PetscBool printreason; /* prints converged reason after solve */
58: PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */
60: Vec vec_sol,vec_rhs; /* pointer to where user has stashed
61: the solution and rhs, these are
62: never touched by the code, only
63: passed back to the user */
64: PetscReal *res_hist; /* If !0 stores residual at iterations*/
65: PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
66: PetscInt res_hist_len; /* current size of residual history array */
67: PetscInt res_hist_max; /* actual amount of data in residual_history */
68: PetscBool res_hist_reset; /* reset history to size zero for each new solve */
70: PetscInt chknorm; /* only compute/check norm if iterations is great than this */
71: PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the
72: MPI_Allreduce() for computing the inner products for the next iteration. */
73: /* --------User (or default) routines (most return -1 on error) --------*/
74: PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */
75: PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void**); /* */
76: void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */
77: PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */
79: PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
80: PetscErrorCode (*convergeddestroy)(void*);
81: void *cnvP;
83: void *user; /* optional user-defined context */
85: PC pc;
87: void *data; /* holder for misc stuff associated
88: with a particular iterative solver */
90: /* ----------------Default work-area management -------------------- */
91: PetscInt nwork;
92: Vec *work;
94: KSPSetUpStage setupstage;
96: PetscInt its; /* number of iterations so far computed */
98: PetscBool transpose_solve; /* solve transpose system instead */
100: KSPNormType normtype; /* type of norm used for convergence tests */
102: /* Allow diagonally scaling the matrix before computing the preconditioner or using
103: the Krylov method. Note this is NOT just Jacobi preconditioning */
105: PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */
106: PetscBool dscalefix; /* unscale system after solve */
107: PetscBool dscalefix2; /* system has been unscaled */
108: Vec diagonal; /* 1/sqrt(diag of matrix) */
109: Vec truediagonal;
111: MatNullSpace nullsp; /* Null space of the operator, removed from Krylov space */
112: };
114: typedef struct {
115: PetscBool initialrtol; /* default relative residual decrease is computing from initial residual, not rhs */
116: PetscBool mininitialrtol; /* default relative residual decrease is computing from min of initial residual and rhs */
117: Vec work;
118: } KSPDefaultConvergedCtx;
120: #define KSPLogResidualHistory(ksp,norm) \
121: {if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) \
122: ksp->res_hist[ksp->res_hist_len++] = norm;}
129: /*
130: These allow the various Krylov methods to apply to either the linear system or its transpose.
131: */
132: #define KSP_RemoveNullSpace(ksp,y) ((ksp->nullsp && ksp->pc_side == PC_LEFT) ? MatNullSpaceRemove(ksp->nullsp,y,PETSC_NULL) : 0)
134: #define KSP_MatMult(ksp,A,x,y) (!ksp->transpose_solve) ? MatMult(A,x,y) : MatMultTranspose(A,x,y)
135: #define KSP_MatMultTranspose(ksp,A,x,y) (!ksp->transpose_solve) ? MatMultTranspose(A,x,y) : MatMult(A,x,y)
136: #define KSP_PCApply(ksp,x,y) (!ksp->transpose_solve) ? (PCApply(ksp->pc,x,y) || KSP_RemoveNullSpace(ksp,y)) : PCApplyTranspose(ksp->pc,x,y)
137: #define KSP_PCApplyTranspose(ksp,x,y) (!ksp->transpose_solve) ? PCApplyTranspose(ksp->pc,x,y) : (PCApply(ksp->pc,x,y) || KSP_RemoveNullSpace(ksp,y))
138: #define KSP_PCApplyBAorAB(ksp,x,y,w) (!ksp->transpose_solve) ? (PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w) || KSP_RemoveNullSpace(ksp,y)) : PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w)
139: #define KSP_PCApplyBAorABTranspose(ksp,x,y,w) (!ksp->transpose_solve) ? (PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w) || KSP_RemoveNullSpace(ksp,y)) : PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w)
144: #endif