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