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: };
26: typedef struct {PetscInt model,curl,maxl;Mat mat; KSP ksp;}* KSPGuessFischer;
28: /*
29: Maximum number of monitors you can run with a single KSP
30: */
31: #define MAXKSPMONITORS 5
33: /*
34: Defines the KSP data structure.
35: */
36: struct _p_KSP {
37: PETSCHEADER(struct _KSPOps);
38: /*------------------------- User parameters--------------------------*/
39: PetscInt max_it; /* maximum number of iterations */
40: KSPFischerGuess guess;
41: PetscTruth guess_zero, /* flag for whether initial guess is 0 */
42: calc_sings, /* calculate extreme Singular Values */
43: guess_knoll; /* use initial guess of PCApply(ksp->B,b */
44: PCSide pc_side; /* flag for left, right, or symmetric
45: preconditioning */
46: PetscReal rtol, /* relative tolerance */
47: abstol, /* absolute tolerance */
48: ttol, /* (not set by user) */
49: divtol; /* divergence tolerance */
50: PetscReal rnorm0; /* initial residual norm (used for divergence testing) */
51: PetscReal rnorm; /* current residual norm */
52: KSPConvergedReason reason;
53: PetscTruth printreason; /* prints converged reason after solve */
55: Vec vec_sol,vec_rhs; /* pointer to where user has stashed
56: the solution and rhs, these are
57: never touched by the code, only
58: passed back to the user */
59: PetscReal *res_hist; /* If !0 stores residual at iterations*/
60: PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */
61: PetscInt res_hist_len; /* current size of residual history array */
62: PetscInt res_hist_max; /* actual amount of data in residual_history */
63: PetscTruth res_hist_reset; /* reset history to size zero for each new solve */
65: PetscInt chknorm; /* only compute/check norm if iterations is great than this */
66: PetscTruth lagnorm; /* Lag the residual norm calculation so that it is computed as part of the
67: MPI_Allreduce() for computing the inner products for the next iteration. */
68: /* --------User (or default) routines (most return -1 on error) --------*/
69: PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */
70: PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void*); /* */
71: void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */
72: PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */
74: PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
75: PetscErrorCode (*convergeddestroy)(void*);
76: void *cnvP;
78: PC pc;
80: void *data; /* holder for misc stuff associated
81: with a particular iterative solver */
83: /* ----------------Default work-area management -------------------- */
84: PetscInt nwork;
85: Vec *work;
87: PetscInt setupcalled;
89: PetscInt its; /* number of iterations so far computed */
91: PetscTruth transpose_solve; /* solve transpose system instead */
93: KSPNormType normtype; /* type of norm used for convergence tests */
95: /* Allow diagonally scaling the matrix before computing the preconditioner or using
96: the Krylov method. Note this is NOT just Jacobi preconditioning */
98: PetscTruth dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */
99: PetscTruth dscalefix; /* unscale system after solve */
100: PetscTruth dscalefix2; /* system has been unscaled */
101: Vec diagonal; /* 1/sqrt(diag of matrix) */
102: Vec truediagonal;
104: MatNullSpace nullsp; /* Null space of the operator, removed from Krylov space */
105: };
107: #define KSPLogResidualHistory(ksp,norm) \
108: {if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) \
109: ksp->res_hist[ksp->res_hist_len++] = norm;}
111: #define KSPMonitor(ksp,it,rnorm) \
112: { PetscErrorCode _ierr; PetscInt _i,_im = ksp->numbermonitors; \
113: for (_i=0; _i<_im; _i++) {\
114: _(*ksp->monitor[_i])(ksp,it,rnorm,ksp->monitorcontext[_i]);CHKERRQ(_ierr); \
115: } \
116: }
118: EXTERN PetscErrorCode KSPDefaultDestroy(KSP);
119: EXTERN PetscErrorCode KSPGetVecs(KSP,PetscInt,Vec**,PetscInt,Vec**);
120: EXTERN PetscErrorCode KSPDefaultGetWork(KSP,PetscInt);
121: EXTERN PetscErrorCode KSPDefaultFreeWork(KSP);
123: /*
124: These allow the various Krylov methods to apply to either the linear system or its transpose.
125: */
126: #define KSP_RemoveNullSpace(ksp,y) ((ksp->nullsp && ksp->pc_side == PC_LEFT) ? MatNullSpaceRemove(ksp->nullsp,y,PETSC_NULL) : 0)
128: #define KSP_MatMult(ksp,A,x,y) (!ksp->transpose_solve) ? MatMult(A,x,y) : MatMultTranspose(A,x,y)
129: #define KSP_MatMultTranspose(ksp,A,x,y) (!ksp->transpose_solve) ? MatMultTranspose(A,x,y) : MatMult(A,x,y)
130: #define KSP_PCApply(ksp,x,y) (!ksp->transpose_solve) ? (PCApply(ksp->pc,x,y) || KSP_RemoveNullSpace(ksp,y)) : PCApplyTranspose(ksp->pc,x,y)
131: #define KSP_PCApplyTranspose(ksp,x,y) (!ksp->transpose_solve) ? PCApplyTranspose(ksp->pc,x,y) : (PCApply(ksp->pc,x,y) || KSP_RemoveNullSpace(ksp,y))
132: #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)
133: #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)
138: #endif