Actual source code: petsc-kspimpl.h

petsc-3.3-p3 2012-08-29
  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;

122: PETSC_STATIC_INLINE PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm)
123: {
125:   if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) {
126:     ksp->res_hist[ksp->res_hist_len++] = norm;
127:   }
128:   return(0);
129: }

131: PETSC_EXTERN PetscErrorCode KSPDefaultDestroy(KSP);
132: PETSC_EXTERN PetscErrorCode KSPGetVecs(KSP,PetscInt,Vec**,PetscInt,Vec**);
133: PETSC_EXTERN PetscErrorCode KSPDefaultGetWork(KSP,PetscInt);
134: PETSC_EXTERN PetscErrorCode KSPSetUpNorms_Private(KSP,KSPNormType*,PCSide*);

136: /* Context for resolution-dependent KSP callback information */
137: typedef struct _n_KSPDM *KSPDM;
138: struct _n_KSPDM {
139:   PetscErrorCode (*computeoperators)(KSP,Mat,Mat,MatStructure*,void*);
140:   PetscErrorCode (*computerhs)(KSP,Vec,void*);
141:   void *operatorsctx;
142:   void *rhsctx;

144:   /* This context/destroy pair allows implementation-specific routines such as DMDA local functions. */
145:   PetscErrorCode (*destroy)(KSPDM);
146:   void *data;

148:   /* This is NOT reference counted. The KSP that originally created this context is cached here to implement copy-on-write.
149:    * Fields in the KSPDM should only be written if the KSP matches originalksp.
150:    */
151:   DM originaldm;

153:   void (*fortran_func_pointers[2])(void); /* Store our own function pointers so they are associated with the KSPDM instead of the DM */
154: };
155: PETSC_EXTERN PetscErrorCode DMKSPGetContext(DM,KSPDM*);
156: PETSC_EXTERN PetscErrorCode DMKSPGetContextWrite(DM,KSPDM*);
157: PETSC_EXTERN PetscErrorCode DMKSPCopyContext(DM,DM);
158: PETSC_EXTERN PetscErrorCode DMKSPDuplicateContext(DM,DM);

160: /*
161:        These allow the various Krylov methods to apply to either the linear system or its transpose.
162: */
163: #define KSP_RemoveNullSpace(ksp,y) ((ksp->nullsp && ksp->pc_side == PC_LEFT) ? MatNullSpaceRemove(ksp->nullsp,y,PETSC_NULL) : 0)

165: #define KSP_MatMult(ksp,A,x,y)          (!ksp->transpose_solve) ? MatMult(A,x,y)                                                            : MatMultTranspose(A,x,y) 
166: #define KSP_MatMultTranspose(ksp,A,x,y) (!ksp->transpose_solve) ? MatMultTranspose(A,x,y)                                                   : MatMult(A,x,y) 
167: #define KSP_PCApply(ksp,x,y)            (!ksp->transpose_solve) ? (PCApply(ksp->pc,x,y) || KSP_RemoveNullSpace(ksp,y))                      : PCApplyTranspose(ksp->pc,x,y) 
168: #define KSP_PCApplyTranspose(ksp,x,y)   (!ksp->transpose_solve) ? PCApplyTranspose(ksp->pc,x,y)                                             : (PCApply(ksp->pc,x,y) || KSP_RemoveNullSpace(ksp,y)) 
169: #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)
170: #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)

172: PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization, KSP_SetUp, KSP_Solve;


175: #endif