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