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 PetscBool view, viewPre, viewReason, viewMat, viewPMat, viewRhs, viewSol, viewMatExp, viewEV, viewSV, viewEVExp, viewFinalRes, viewPOpExp, viewDScale; 124 PetscViewer viewer, viewerPre, viewerReason, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp, viewerDScale; 125 PetscViewerFormat format, formatPre, formatReason, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp, formatDScale; 126 127 /* ----------------Default work-area management -------------------- */ 128 PetscInt nwork; 129 Vec *work; 130 131 KSPSetUpStage setupstage; 132 PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */ 133 134 PetscInt its; /* number of iterations so far computed in THIS linear solve*/ 135 PetscInt totalits; /* number of iterations used by this KSP object since it was created */ 136 137 PetscBool transpose_solve; /* solve transpose system instead */ 138 139 KSPNormType normtype; /* type of norm used for convergence tests */ 140 141 PCSide pc_side_set; /* PC type set explicitly by user */ 142 KSPNormType normtype_set; /* Norm type set explicitly by user */ 143 144 /* Allow diagonally scaling the matrix before computing the preconditioner or using 145 the Krylov method. Note this is NOT just Jacobi preconditioning */ 146 147 PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */ 148 PetscBool dscalefix; /* unscale system after solve */ 149 PetscBool dscalefix2; /* system has been unscaled */ 150 Vec diagonal; /* 1/sqrt(diag of matrix) */ 151 Vec truediagonal; 152 153 PetscInt setfromoptionscalled; 154 PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */ 155 156 PetscViewer eigviewer; /* Viewer where computed eigenvalues are displayed */ 157 158 PetscErrorCode (*presolve)(KSP,Vec,Vec,void*); 159 PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*); 160 void *prectx,*postctx; 161 }; 162 163 typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */ 164 PetscReal coef; 165 PetscReal bnrm; 166 } KSPDynTolCtx; 167 168 typedef struct { 169 PetscBool initialrtol; /* default relative residual decrease is computing from initial residual, not rhs */ 170 PetscBool mininitialrtol; /* default relative residual decrease is computing from min of initial residual and rhs */ 171 Vec work; 172 } KSPConvergedDefaultCtx; 173 174 PETSC_STATIC_INLINE PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm) 175 { 176 PetscErrorCode ierr; 177 178 PetscFunctionBegin; 179 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr); 180 if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) { 181 ksp->res_hist[ksp->res_hist_len++] = norm; 182 } 183 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr); 184 PetscFunctionReturn(0); 185 } 186 187 PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP,PetscBool,KSPNormType*,PCSide*); 188 189 PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP,PetscInt,const PetscReal*,const PetscReal*); 190 191 typedef struct _p_DMKSP *DMKSP; 192 typedef struct _DMKSPOps *DMKSPOps; 193 struct _DMKSPOps { 194 PetscErrorCode (*computeoperators)(KSP,Mat,Mat,void*); 195 PetscErrorCode (*computerhs)(KSP,Vec,void*); 196 PetscErrorCode (*computeinitialguess)(KSP,Vec,void*); 197 PetscErrorCode (*destroy)(DMKSP*); 198 PetscErrorCode (*duplicate)(DMKSP,DMKSP); 199 }; 200 201 struct _p_DMKSP { 202 PETSCHEADER(struct _DMKSPOps); 203 void *operatorsctx; 204 void *rhsctx; 205 void *initialguessctx; 206 void *data; 207 208 /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way 209 * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user 210 * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by 211 * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates 212 * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes 213 * to the original level will no longer propagate to that level. 214 */ 215 DM originaldm; 216 217 void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */ 218 }; 219 PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM,DMKSP*); 220 PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM,DMKSP*); 221 PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM,DM); 222 223 /* 224 These allow the various Krylov methods to apply to either the linear system or its transpose. 225 */ 226 PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpace(KSP ksp,Vec y) 227 { 228 PetscErrorCode ierr; 229 PetscFunctionBegin; 230 if (ksp->pc_side == PC_LEFT) { 231 Mat A; 232 MatNullSpace nullsp; 233 ierr = PCGetOperators(ksp->pc,&A,NULL);CHKERRQ(ierr); 234 ierr = MatGetNullSpace(A,&nullsp);CHKERRQ(ierr); 235 if (nullsp) { 236 ierr = MatNullSpaceRemove(nullsp,y);CHKERRQ(ierr); 237 } 238 } 239 PetscFunctionReturn(0); 240 } 241 242 PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp,Vec y) 243 { 244 PetscErrorCode ierr; 245 PetscFunctionBegin; 246 if (ksp->pc_side == PC_LEFT) { 247 Mat A; 248 MatNullSpace nullsp; 249 ierr = PCGetOperators(ksp->pc,&A,NULL);CHKERRQ(ierr); 250 ierr = MatGetTransposeNullSpace(A,&nullsp);CHKERRQ(ierr); 251 if (nullsp) { 252 ierr = MatNullSpaceRemove(nullsp,y);CHKERRQ(ierr); 253 } 254 } 255 PetscFunctionReturn(0); 256 } 257 258 PETSC_STATIC_INLINE PetscErrorCode KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y) 259 { 260 PetscErrorCode ierr; 261 PetscFunctionBegin; 262 if (!ksp->transpose_solve) {ierr = MatMult(A,x,y);CHKERRQ(ierr);} 263 else {ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);} 264 PetscFunctionReturn(0); 265 } 266 267 PETSC_STATIC_INLINE PetscErrorCode KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y) 268 { 269 PetscErrorCode ierr; 270 PetscFunctionBegin; 271 if (!ksp->transpose_solve) {ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);} 272 else {ierr = MatMult(A,x,y);CHKERRQ(ierr);} 273 PetscFunctionReturn(0); 274 } 275 276 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y) 277 { 278 PetscErrorCode ierr; 279 PetscFunctionBegin; 280 if (!ksp->transpose_solve) { 281 ierr = PCApply(ksp->pc,x,y);CHKERRQ(ierr); 282 ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr); 283 } else { 284 ierr = PCApplyTranspose(ksp->pc,x,y);CHKERRQ(ierr); 285 ierr = KSP_RemoveNullSpaceTranspose(ksp,y);CHKERRQ(ierr); 286 } 287 PetscFunctionReturn(0); 288 } 289 290 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y) 291 { 292 PetscErrorCode ierr; 293 PetscFunctionBegin; 294 if (!ksp->transpose_solve) { 295 ierr = PCApplyTranspose(ksp->pc,x,y);CHKERRQ(ierr); 296 ierr = KSP_RemoveNullSpaceTranspose(ksp,y);CHKERRQ(ierr); 297 } else { 298 ierr = PCApply(ksp->pc,x,y);CHKERRQ(ierr); 299 ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr); 300 } 301 PetscFunctionReturn(0); 302 } 303 304 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w) 305 { 306 PetscErrorCode ierr; 307 PetscFunctionBegin; 308 if (!ksp->transpose_solve) { 309 ierr = PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 310 ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr); 311 } else { 312 ierr = PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 313 ierr = KSP_RemoveNullSpaceTranspose(ksp,y);CHKERRQ(ierr); 314 } 315 PetscFunctionReturn(0); 316 } 317 318 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w) 319 { 320 PetscErrorCode ierr; 321 PetscFunctionBegin; 322 if (!ksp->transpose_solve) { 323 ierr = PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 324 } else { 325 ierr = PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 326 } 327 PetscFunctionReturn(0); 328 } 329 330 PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization; 331 PETSC_EXTERN PetscLogEvent KSP_SetUp; 332 PETSC_EXTERN PetscLogEvent KSP_Solve; 333 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0; 334 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1; 335 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2; 336 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3; 337 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4; 338 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S; 339 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L; 340 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U; 341 342 PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*); 343 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*); 344 345 /*MC 346 KSPCheckDot - Checks if the result of a dot product used by the corresponding KSP contains Inf or NaN. These indicate that the previous 347 application of the preconditioner generated an error 348 349 Collective on KSP 350 351 Input Parameter: 352 . ksp - the linear solver (KSP) context. 353 354 Output Parameter: 355 . beta - the result of the inner product 356 357 Level: developer 358 359 Developer Note: 360 this is used to manage returning from KSP solvers whose preconditioners have failed in some way 361 362 .keywords: KSP, PC, divergence, convergence 363 364 .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckNorm(), KSPCheckSolve() 365 M*/ 366 #define KSPCheckDot(ksp,beta) do { \ 367 if (PetscIsInfOrNanScalar(beta)) { \ 368 if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");\ 369 else {\ 370 PetscErrorCode ierr;\ 371 PCFailedReason pcreason;\ 372 PetscInt sendbuf,pcreason_max; \ 373 ierr = PCGetFailedReason(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_PC_FAILED;\ 378 ierr = VecSetInf(ksp->vec_sol);CHKERRQ(ierr);\ 379 } else {\ 380 ksp->reason = KSP_DIVERGED_NANORINF;\ 381 }\ 382 PetscFunctionReturn(0);\ 383 }\ 384 } } while (0) 385 386 /*MC 387 KSPCheckNorm - Checks if the result of a norm used by the corresponding KSP contains Inf or NaN. These indicate that the previous 388 application of the preconditioner generated an error 389 390 Collective on KSP 391 392 Input Parameter: 393 . ksp - the linear solver (KSP) context. 394 395 Output Parameter: 396 . beta - the result of the norm 397 398 Level: developer 399 400 Developer Note: 401 this is used to manage returning from KSP solvers whose preconditioners have failed in some way 402 403 .keywords: KSP, PC, divergence, convergence 404 405 .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckDot(), KSPCheckSolve() 406 M*/ 407 #define KSPCheckNorm(ksp,beta) do { \ 408 if (PetscIsInfOrNanReal(beta)) { \ 409 if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf norm");\ 410 else {\ 411 PetscErrorCode ierr;\ 412 PCFailedReason pcreason;\ 413 PetscInt sendbuf,pcreason_max; \ 414 ierr = PCGetFailedReason(ksp->pc,&pcreason);CHKERRQ(ierr);\ 415 sendbuf = (PetscInt)pcreason; \ 416 ierr = MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr); \ 417 if (pcreason_max) {\ 418 ksp->reason = KSP_DIVERGED_PC_FAILED;\ 419 ierr = VecSetInf(ksp->vec_sol);CHKERRQ(ierr);\ 420 } else {\ 421 ksp->reason = KSP_DIVERGED_NANORINF;\ 422 }\ 423 PetscFunctionReturn(0);\ 424 }\ 425 } } while (0) 426 427 #endif 428