1 2 #ifndef _KSPIMPL_H 3 #define _KSPIMPL_H 4 5 #include <petscksp.h> 6 #include <petsc/private/petscimpl.h> 7 8 PETSC_EXTERN PetscBool KSPRegisterAllCalled; 9 PETSC_EXTERN PetscErrorCode KSPRegisterAll(void); 10 PETSC_EXTERN PetscErrorCode KSPMatRegisterAll(void); 11 12 typedef struct _KSPOps *KSPOps; 13 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)(PetscOptionItems*,KSP); 24 PetscErrorCode (*publishoptions)(KSP); 25 PetscErrorCode (*computeextremesingularvalues)(KSP,PetscReal*,PetscReal*); 26 PetscErrorCode (*computeeigenvalues)(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt *); 27 PetscErrorCode (*computeritz)(KSP,PetscBool,PetscBool,PetscInt*,Vec[],PetscReal*,PetscReal*); 28 PetscErrorCode (*destroy)(KSP); 29 PetscErrorCode (*view)(KSP,PetscViewer); 30 PetscErrorCode (*reset)(KSP); 31 PetscErrorCode (*load)(KSP,PetscViewer); 32 }; 33 34 typedef struct {PetscInt model,curl,maxl;Mat mat; KSP ksp;}* KSPGuessFischer; 35 36 /* 37 Maximum number of monitors you can run with a single KSP 38 */ 39 #define MAXKSPMONITORS 5 40 typedef enum {KSP_SETUP_NEW, KSP_SETUP_NEWMATRIX, KSP_SETUP_NEWRHS} KSPSetUpStage; 41 42 /* 43 Defines the KSP data structure. 44 */ 45 struct _p_KSP { 46 PETSCHEADER(struct _KSPOps); 47 DM dm; 48 PetscBool dmAuto; /* DM was created automatically by KSP */ 49 PetscBool dmActive; /* KSP should use DM for computing operators */ 50 /*------------------------- User parameters--------------------------*/ 51 PetscInt max_it; /* maximum number of iterations */ 52 KSPFischerGuess guess; 53 PetscBool guess_zero, /* flag for whether initial guess is 0 */ 54 calc_sings, /* calculate extreme Singular Values */ 55 calc_ritz, /* calculate (harmonic) Ritz pairs */ 56 guess_knoll; /* use initial guess of PCApply(ksp->B,b */ 57 PCSide pc_side; /* flag for left, right, or symmetric preconditioning */ 58 PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */ 59 PetscReal rtol, /* relative tolerance */ 60 abstol, /* absolute tolerance */ 61 ttol, /* (not set by user) */ 62 divtol; /* divergence tolerance */ 63 PetscReal rnorm0; /* initial residual norm (used for divergence testing) */ 64 PetscReal rnorm; /* current residual norm */ 65 KSPConvergedReason reason; 66 PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */ 67 68 Vec vec_sol,vec_rhs; /* pointer to where user has stashed 69 the solution and rhs, these are 70 never touched by the code, only 71 passed back to the user */ 72 PetscReal *res_hist; /* If !0 stores residual at iterations*/ 73 PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */ 74 PetscInt res_hist_len; /* current size of residual history array */ 75 PetscInt res_hist_max; /* actual amount of data in residual_history */ 76 PetscBool res_hist_reset; /* reset history to size zero for each new solve */ 77 78 PetscInt chknorm; /* only compute/check norm if iterations is great than this */ 79 PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the 80 MPI_Allreduce() for computing the inner products for the next iteration. */ 81 /* --------User (or default) routines (most return -1 on error) --------*/ 82 PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */ 83 PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void**); /* */ 84 void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */ 85 PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */ 86 87 PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*); 88 PetscErrorCode (*convergeddestroy)(void*); 89 void *cnvP; 90 91 void *user; /* optional user-defined context */ 92 93 PC pc; 94 95 void *data; /* holder for misc stuff associated 96 with a particular iterative solver */ 97 98 /* ----------------Default work-area management -------------------- */ 99 PetscInt nwork; 100 Vec *work; 101 102 KSPSetUpStage setupstage; 103 PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */ 104 105 PetscInt its; /* number of iterations so far computed in THIS linear solve*/ 106 PetscInt totalits; /* number of iterations used by this KSP object since it was created */ 107 108 PetscBool transpose_solve; /* solve transpose system instead */ 109 110 KSPNormType normtype; /* type of norm used for convergence tests */ 111 112 PCSide pc_side_set; /* PC type set explicitly by user */ 113 KSPNormType normtype_set; /* Norm type set explicitly by user */ 114 115 /* Allow diagonally scaling the matrix before computing the preconditioner or using 116 the Krylov method. Note this is NOT just Jacobi preconditioning */ 117 118 PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */ 119 PetscBool dscalefix; /* unscale system after solve */ 120 PetscBool dscalefix2; /* system has been unscaled */ 121 Vec diagonal; /* 1/sqrt(diag of matrix) */ 122 Vec truediagonal; 123 124 PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */ 125 126 PetscViewer eigviewer; /* Viewer where computed eigenvalues are displayed */ 127 128 PetscErrorCode (*presolve)(KSP,Vec,Vec,void*); 129 PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*); 130 void *prectx,*postctx; 131 }; 132 133 typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */ 134 PetscReal coef; 135 PetscReal bnrm; 136 } KSPDynTolCtx; 137 138 typedef struct { 139 PetscBool initialrtol; /* default relative residual decrease is computing from initial residual, not rhs */ 140 PetscBool mininitialrtol; /* default relative residual decrease is computing from min of initial residual and rhs */ 141 Vec work; 142 } KSPConvergedDefaultCtx; 143 144 PETSC_STATIC_INLINE PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm) 145 { 146 PetscErrorCode ierr; 147 148 PetscFunctionBegin; 149 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr); 150 if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) { 151 ksp->res_hist[ksp->res_hist_len++] = norm; 152 } 153 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr); 154 PetscFunctionReturn(0); 155 } 156 157 PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP,PetscBool,KSPNormType*,PCSide*); 158 159 PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP,PetscInt,const PetscReal*,const PetscReal*); 160 161 typedef struct _p_DMKSP *DMKSP; 162 typedef struct _DMKSPOps *DMKSPOps; 163 struct _DMKSPOps { 164 PetscErrorCode (*computeoperators)(KSP,Mat,Mat,void*); 165 PetscErrorCode (*computerhs)(KSP,Vec,void*); 166 PetscErrorCode (*computeinitialguess)(KSP,Vec,void*); 167 PetscErrorCode (*destroy)(DMKSP*); 168 PetscErrorCode (*duplicate)(DMKSP,DMKSP); 169 }; 170 171 struct _p_DMKSP { 172 PETSCHEADER(struct _DMKSPOps); 173 void *operatorsctx; 174 void *rhsctx; 175 void *initialguessctx; 176 void *data; 177 178 /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way 179 * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user 180 * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by 181 * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates 182 * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes 183 * to the original level will no longer propagate to that level. 184 */ 185 DM originaldm; 186 187 void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */ 188 }; 189 PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM,DMKSP*); 190 PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM,DMKSP*); 191 PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM,DM); 192 193 /* 194 These allow the various Krylov methods to apply to either the linear system or its transpose. 195 */ 196 PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpace(KSP ksp,Vec y) 197 { 198 PetscErrorCode ierr; 199 PetscFunctionBegin; 200 if (ksp->pc_side == PC_LEFT) { 201 Mat A; 202 MatNullSpace nullsp; 203 ierr = PCGetOperators(ksp->pc,&A,NULL);CHKERRQ(ierr); 204 ierr = MatGetNullSpace(A,&nullsp);CHKERRQ(ierr); 205 if (nullsp) { 206 ierr = MatNullSpaceRemove(nullsp,y);CHKERRQ(ierr); 207 } 208 } 209 PetscFunctionReturn(0); 210 } 211 212 PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp,Vec y) 213 { 214 PetscErrorCode ierr; 215 PetscFunctionBegin; 216 if (ksp->pc_side == PC_LEFT) { 217 Mat A; 218 MatNullSpace nullsp; 219 ierr = PCGetOperators(ksp->pc,&A,NULL);CHKERRQ(ierr); 220 ierr = MatGetTransposeNullSpace(A,&nullsp);CHKERRQ(ierr); 221 if (nullsp) { 222 ierr = MatNullSpaceRemove(nullsp,y);CHKERRQ(ierr); 223 } 224 } 225 PetscFunctionReturn(0); 226 } 227 228 PETSC_STATIC_INLINE PetscErrorCode KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y) 229 { 230 PetscErrorCode ierr; 231 PetscFunctionBegin; 232 if (!ksp->transpose_solve) {ierr = MatMult(A,x,y);CHKERRQ(ierr);} 233 else {ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);} 234 PetscFunctionReturn(0); 235 } 236 237 PETSC_STATIC_INLINE PetscErrorCode KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y) 238 { 239 PetscErrorCode ierr; 240 PetscFunctionBegin; 241 if (!ksp->transpose_solve) {ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);} 242 else {ierr = MatMult(A,x,y);CHKERRQ(ierr);} 243 PetscFunctionReturn(0); 244 } 245 246 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y) 247 { 248 PetscErrorCode ierr; 249 PetscFunctionBegin; 250 if (!ksp->transpose_solve) { 251 ierr = PCApply(ksp->pc,x,y);CHKERRQ(ierr); 252 ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr); 253 } else { 254 ierr = PCApplyTranspose(ksp->pc,x,y);CHKERRQ(ierr); 255 ierr = KSP_RemoveNullSpaceTranspose(ksp,y);CHKERRQ(ierr); 256 } 257 PetscFunctionReturn(0); 258 } 259 260 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y) 261 { 262 PetscErrorCode ierr; 263 PetscFunctionBegin; 264 if (!ksp->transpose_solve) { 265 ierr = PCApplyTranspose(ksp->pc,x,y);CHKERRQ(ierr); 266 ierr = KSP_RemoveNullSpaceTranspose(ksp,y);CHKERRQ(ierr); 267 } else { 268 ierr = PCApply(ksp->pc,x,y);CHKERRQ(ierr); 269 ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr); 270 } 271 PetscFunctionReturn(0); 272 } 273 274 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w) 275 { 276 PetscErrorCode ierr; 277 PetscFunctionBegin; 278 if (!ksp->transpose_solve) { 279 ierr = PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 280 ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr); 281 } else { 282 ierr = PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 283 ierr = KSP_RemoveNullSpaceTranspose(ksp,y);CHKERRQ(ierr); 284 } 285 PetscFunctionReturn(0); 286 } 287 288 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w) 289 { 290 PetscErrorCode ierr; 291 PetscFunctionBegin; 292 if (!ksp->transpose_solve) { 293 ierr = PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 294 } else { 295 ierr = PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 296 } 297 PetscFunctionReturn(0); 298 } 299 300 PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization, KSP_SetUp, KSP_Solve; 301 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0,KSP_Solve_FS_1,KSP_Solve_FS_2,KSP_Solve_FS_3,KSP_Solve_FS_4,KSP_Solve_FS_S,KSP_Solve_FS_L,KSP_Solve_FS_U; 302 303 PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*); 304 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*); 305 306 /* 307 Either generate an error or mark as diverged when a scalar from an inner product is Nan or Inf 308 */ 309 #define KSPCheckDot(ksp,beta) \ 310 if (PetscIsInfOrNanScalar(beta)) { \ 311 if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");\ 312 else {\ 313 PetscErrorCode ierr;\ 314 PCFailedReason pcreason;\ 315 PetscInt sendbuf,pcreason_max; \ 316 ierr = PCGetSetUpFailedReason(ksp->pc,&pcreason);CHKERRQ(ierr);\ 317 sendbuf = (PetscInt)pcreason; \ 318 ierr = MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr); \ 319 if (pcreason_max) {\ 320 ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;\ 321 ierr = VecSetInf(ksp->vec_sol);CHKERRQ(ierr);\ 322 } else {\ 323 ksp->reason = KSP_DIVERGED_NANORINF;\ 324 }\ 325 PetscFunctionReturn(0);\ 326 }\ 327 } 328 329 /* 330 Either generate an error or mark as diverged when a real from a norm is Nan or Inf 331 */ 332 #define KSPCheckNorm(ksp,beta) \ 333 if (PetscIsInfOrNanReal(beta)) { \ 334 if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf norm");\ 335 else {\ 336 PetscErrorCode ierr;\ 337 PCFailedReason pcreason;\ 338 PetscInt sendbuf,pcreason_max; \ 339 ierr = PCGetSetUpFailedReason(ksp->pc,&pcreason);CHKERRQ(ierr);\ 340 sendbuf = (PetscInt)pcreason; \ 341 ierr = MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr); \ 342 if (pcreason_max) {\ 343 ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;\ 344 ierr = VecSetInf(ksp->vec_sol);CHKERRQ(ierr);\ 345 } else {\ 346 ksp->reason = KSP_DIVERGED_NANORINF;\ 347 }\ 348 PetscFunctionReturn(0);\ 349 }\ 350 } 351 352 #endif 353