Actual source code: kspimpl.h

petsc-3.6.2 2015-10-02
Report Typos and Errors
  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