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