1 2 #ifndef _KSPIMPL_H 3 #define _KSPIMPL_H 4 5 #include <petscksp.h> 6 #include <petscds.h> 7 #include <petsc/private/petscimpl.h> 8 9 /* SUBMANSEC = KSP */ 10 11 PETSC_EXTERN PetscBool KSPRegisterAllCalled; 12 PETSC_EXTERN PetscBool KSPMonitorRegisterAllCalled; 13 PETSC_EXTERN PetscErrorCode KSPRegisterAll(void); 14 PETSC_EXTERN PetscErrorCode KSPMonitorRegisterAll(void); 15 PETSC_EXTERN PetscErrorCode KSPGuessRegisterAll(void); 16 PETSC_EXTERN PetscErrorCode KSPMatRegisterAll(void); 17 18 typedef struct _KSPOps *KSPOps; 19 20 struct _KSPOps { 21 PetscErrorCode (*buildsolution)(KSP, Vec, Vec *); /* Returns a pointer to the solution, or 22 calculates the solution in a 23 user-provided area. */ 24 PetscErrorCode (*buildresidual)(KSP, Vec, Vec, Vec *); /* Returns a pointer to the residual, or 25 calculates the residual in a 26 user-provided area. */ 27 PetscErrorCode (*solve)(KSP); /* actual solver */ 28 PetscErrorCode (*matsolve)(KSP, Mat, Mat); /* multiple dense RHS solver */ 29 PetscErrorCode (*setup)(KSP); 30 PetscErrorCode (*setfromoptions)(KSP, PetscOptionItems *); 31 PetscErrorCode (*publishoptions)(KSP); 32 PetscErrorCode (*computeextremesingularvalues)(KSP, PetscReal *, PetscReal *); 33 PetscErrorCode (*computeeigenvalues)(KSP, PetscInt, PetscReal *, PetscReal *, PetscInt *); 34 PetscErrorCode (*computeritz)(KSP, PetscBool, PetscBool, PetscInt *, Vec[], PetscReal *, PetscReal *); 35 PetscErrorCode (*destroy)(KSP); 36 PetscErrorCode (*view)(KSP, PetscViewer); 37 PetscErrorCode (*reset)(KSP); 38 PetscErrorCode (*load)(KSP, PetscViewer); 39 }; 40 41 typedef struct _KSPGuessOps *KSPGuessOps; 42 43 struct _KSPGuessOps { 44 PetscErrorCode (*formguess)(KSPGuess, Vec, Vec); /* Form initial guess */ 45 PetscErrorCode (*update)(KSPGuess, Vec, Vec); /* Update database */ 46 PetscErrorCode (*setfromoptions)(KSPGuess); 47 PetscErrorCode (*settolerance)(KSPGuess, PetscReal); 48 PetscErrorCode (*setup)(KSPGuess); 49 PetscErrorCode (*destroy)(KSPGuess); 50 PetscErrorCode (*view)(KSPGuess, PetscViewer); 51 PetscErrorCode (*reset)(KSPGuess); 52 }; 53 54 /* 55 Defines the KSPGuess data structure. 56 */ 57 struct _p_KSPGuess { 58 PETSCHEADER(struct _KSPGuessOps); 59 KSP ksp; /* the parent KSP */ 60 Mat A; /* the current linear operator */ 61 PetscObjectState omatstate; /* previous linear operator state */ 62 void *data; /* pointer to the specific implementation */ 63 }; 64 65 PETSC_EXTERN PetscErrorCode KSPGuessCreate_Fischer(KSPGuess); 66 PETSC_EXTERN PetscErrorCode KSPGuessCreate_POD(KSPGuess); 67 68 /* 69 Maximum number of monitors you can run with a single KSP 70 */ 71 #define MAXKSPMONITORS 5 72 #define MAXKSPREASONVIEWS 5 73 typedef enum { 74 KSP_SETUP_NEW, 75 KSP_SETUP_NEWMATRIX, 76 KSP_SETUP_NEWRHS 77 } KSPSetUpStage; 78 79 /* 80 Defines the KSP data structure. 81 */ 82 struct _p_KSP { 83 PETSCHEADER(struct _KSPOps); 84 DM dm; 85 PetscBool dmAuto; /* DM was created automatically by KSP */ 86 PetscBool dmActive; /* KSP should use DM for computing operators */ 87 /*------------------------- User parameters--------------------------*/ 88 PetscInt max_it; /* maximum number of iterations */ 89 KSPGuess guess; 90 PetscBool guess_zero, /* flag for whether initial guess is 0 */ 91 calc_sings, /* calculate extreme Singular Values */ 92 calc_ritz, /* calculate (harmonic) Ritz pairs */ 93 guess_knoll; /* use initial guess of PCApply(ksp->B,b */ 94 PCSide pc_side; /* flag for left, right, or symmetric preconditioning */ 95 PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */ 96 PetscReal rtol, /* relative tolerance */ 97 abstol, /* absolute tolerance */ 98 ttol, /* (not set by user) */ 99 divtol; /* divergence tolerance */ 100 PetscReal rnorm0; /* initial residual norm (used for divergence testing) */ 101 PetscReal rnorm; /* current residual norm */ 102 KSPConvergedReason reason; 103 PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */ 104 105 Vec vec_sol, vec_rhs; /* pointer to where user has stashed 106 the solution and rhs, these are 107 never touched by the code, only 108 passed back to the user */ 109 PetscReal *res_hist; /* If !0 stores residual each at iteration */ 110 PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */ 111 size_t res_hist_len; /* current size of residual history array */ 112 size_t res_hist_max; /* actual amount of storage in residual history */ 113 PetscBool res_hist_reset; /* reset history to length zero for each new solve */ 114 PetscReal *err_hist; /* If !0 stores error at each iteration */ 115 PetscReal *err_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */ 116 size_t err_hist_len; /* current size of error history array */ 117 size_t err_hist_max; /* actual amount of storage in error history */ 118 PetscBool err_hist_reset; /* reset history to length zero for each new solve */ 119 120 PetscInt chknorm; /* only compute/check norm if iterations is great than this */ 121 PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the 122 MPI_Allreduce() for computing the inner products for the next iteration. */ 123 124 PetscInt nmax; /* maximum number of right-hand sides to be handled simultaneously */ 125 126 /* --------User (or default) routines (most return -1 on error) --------*/ 127 PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP, PetscInt, PetscReal, void *); /* returns control to user after */ 128 PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void **); /* */ 129 void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */ 130 PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */ 131 PetscBool pauseFinal; /* Pause all drawing monitor at the final iterate */ 132 133 PetscErrorCode (*reasonview[MAXKSPREASONVIEWS])(KSP, void *); /* KSP converged reason view */ 134 PetscErrorCode (*reasonviewdestroy[MAXKSPREASONVIEWS])(void **); /* Optional destroy routine */ 135 void *reasonviewcontext[MAXKSPREASONVIEWS]; /* User context */ 136 PetscInt numberreasonviews; /* Number if reason viewers */ 137 138 PetscErrorCode (*converged)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *); 139 PetscErrorCode (*convergeddestroy)(void *); 140 void *cnvP; 141 142 void *user; /* optional user-defined context */ 143 144 PC pc; 145 146 void *data; /* holder for misc stuff associated 147 with a particular iterative solver */ 148 149 PetscBool view, viewPre, viewRate, viewMat, viewPMat, viewRhs, viewSol, viewMatExp, viewEV, viewSV, viewEVExp, viewFinalRes, viewPOpExp, viewDScale; 150 PetscViewer viewer, viewerPre, viewerRate, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp, viewerDScale; 151 PetscViewerFormat format, formatPre, formatRate, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp, formatDScale; 152 153 /* ----------------Default work-area management -------------------- */ 154 PetscInt nwork; 155 Vec *work; 156 157 KSPSetUpStage setupstage; 158 PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */ 159 160 PetscInt its; /* number of iterations so far computed in THIS linear solve*/ 161 PetscInt totalits; /* number of iterations used by this KSP object since it was created */ 162 163 PetscBool transpose_solve; /* solve transpose system instead */ 164 struct { 165 Mat AT, BT; 166 PetscBool use_explicittranspose; /* transpose the system explicitly in KSPSolveTranspose */ 167 PetscBool reuse_transpose; /* reuse the previous transposed system */ 168 } transpose; 169 170 KSPNormType normtype; /* type of norm used for convergence tests */ 171 172 PCSide pc_side_set; /* PC type set explicitly by user */ 173 KSPNormType normtype_set; /* Norm type set explicitly by user */ 174 175 /* Allow diagonally scaling the matrix before computing the preconditioner or using 176 the Krylov method. Note this is NOT just Jacobi preconditioning */ 177 178 PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */ 179 PetscBool dscalefix; /* unscale system after solve */ 180 PetscBool dscalefix2; /* system has been unscaled */ 181 Vec diagonal; /* 1/sqrt(diag of matrix) */ 182 Vec truediagonal; 183 184 /* Allow declaring convergence when negative curvature is detected */ 185 PetscBool converged_neg_curve; 186 187 PetscInt setfromoptionscalled; 188 PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */ 189 190 PetscErrorCode (*presolve)(KSP, Vec, Vec, void *); 191 PetscErrorCode (*postsolve)(KSP, Vec, Vec, void *); 192 void *prectx, *postctx; 193 }; 194 195 typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */ 196 PetscReal coef; 197 PetscReal bnrm; 198 } KSPDynTolCtx; 199 200 typedef struct { 201 PetscBool initialrtol; /* default relative residual decrease is computed from initial residual, not rhs */ 202 PetscBool mininitialrtol; /* default relative residual decrease is computed from min of initial residual and rhs */ 203 PetscBool convmaxits; /* if true, the convergence test returns KSP_CONVERGED_ITS if the maximum number of iterations is reached */ 204 Vec work; 205 } KSPConvergedDefaultCtx; 206 207 static inline PetscErrorCode KSPLogResidualHistory(KSP ksp, PetscReal norm) 208 { 209 PetscFunctionBegin; 210 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 211 if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) ksp->res_hist[ksp->res_hist_len++] = norm; 212 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 213 PetscFunctionReturn(PETSC_SUCCESS); 214 } 215 216 static inline PetscErrorCode KSPLogErrorHistory(KSP ksp) 217 { 218 DM dm; 219 220 PetscFunctionBegin; 221 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 222 PetscCall(KSPGetDM(ksp, &dm)); 223 if (dm && ksp->err_hist && ksp->err_hist_max > ksp->err_hist_len) { 224 PetscSimplePointFunc exactSol; 225 void *exactCtx; 226 PetscDS ds; 227 Vec u; 228 PetscReal error; 229 PetscInt Nf; 230 231 PetscCall(KSPBuildSolution(ksp, NULL, &u)); 232 /* TODO Was needed to correct for Newton solution, but I just need to set a solution */ 233 //PetscCall(VecScale(u, -1.0)); 234 /* TODO Case when I have a solution */ 235 if (0) { 236 PetscCall(DMGetDS(dm, &ds)); 237 PetscCall(PetscDSGetNumFields(ds, &Nf)); 238 PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle number of fields %" PetscInt_FMT " > 1 right now", Nf); 239 PetscCall(PetscDSGetExactSolution(ds, 0, &exactSol, &exactCtx)); 240 PetscCall(DMComputeL2FieldDiff(dm, 0.0, &exactSol, &exactCtx, u, &error)); 241 } else { 242 /* The null solution A 0 = 0 */ 243 PetscCall(VecNorm(u, NORM_2, &error)); 244 } 245 ksp->err_hist[ksp->err_hist_len++] = error; 246 } 247 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 248 PetscFunctionReturn(PETSC_SUCCESS); 249 } 250 251 static inline PetscScalar KSPNoisyHash_Private(PetscInt xx) 252 { 253 unsigned int x = (unsigned int)xx; 254 x = ((x >> 16) ^ x) * 0x45d9f3b; 255 x = ((x >> 16) ^ x) * 0x45d9f3b; 256 x = ((x >> 16) ^ x); 257 return (PetscScalar)((PetscInt64)x - 2147483648) * 5.e-10; /* center around zero, scaled about -1. to 1.*/ 258 } 259 260 static inline PetscErrorCode KSPSetNoisy_Private(Vec v) 261 { 262 PetscScalar *a; 263 PetscInt n, istart; 264 265 PetscFunctionBegin; 266 PetscCall(VecGetOwnershipRange(v, &istart, NULL)); 267 PetscCall(VecGetLocalSize(v, &n)); 268 PetscCall(VecGetArrayWrite(v, &a)); 269 for (PetscInt i = 0; i < n; ++i) a[i] = KSPNoisyHash_Private(i + istart); 270 PetscCall(VecRestoreArrayWrite(v, &a)); 271 PetscFunctionReturn(PETSC_SUCCESS); 272 } 273 274 PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP, PetscBool, KSPNormType *, PCSide *); 275 276 PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP, PetscInt, const PetscReal *, const PetscReal *); 277 278 typedef struct _p_DMKSP *DMKSP; 279 typedef struct _DMKSPOps *DMKSPOps; 280 struct _DMKSPOps { 281 PetscErrorCode (*computeoperators)(KSP, Mat, Mat, void *); 282 PetscErrorCode (*computerhs)(KSP, Vec, void *); 283 PetscErrorCode (*computeinitialguess)(KSP, Vec, void *); 284 PetscErrorCode (*destroy)(DMKSP *); 285 PetscErrorCode (*duplicate)(DMKSP, DMKSP); 286 }; 287 288 struct _p_DMKSP { 289 PETSCHEADER(struct _DMKSPOps); 290 void *operatorsctx; 291 void *rhsctx; 292 void *initialguessctx; 293 void *data; 294 295 /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way 296 * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user 297 * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by 298 * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates 299 * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes 300 * to the original level will no longer propagate to that level. 301 */ 302 DM originaldm; 303 304 void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */ 305 }; 306 PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM, DMKSP *); 307 PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM, DMKSP *); 308 PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM, DM); 309 310 /* 311 These allow the various Krylov methods to apply to either the linear system or its transpose. 312 */ 313 static inline PetscErrorCode KSP_RemoveNullSpace(KSP ksp, Vec y) 314 { 315 PetscFunctionBegin; 316 if (ksp->pc_side == PC_LEFT) { 317 Mat A; 318 MatNullSpace nullsp; 319 320 PetscCall(PCGetOperators(ksp->pc, &A, NULL)); 321 PetscCall(MatGetNullSpace(A, &nullsp)); 322 if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y)); 323 } 324 PetscFunctionReturn(PETSC_SUCCESS); 325 } 326 327 static inline PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp, Vec y) 328 { 329 PetscFunctionBegin; 330 if (ksp->pc_side == PC_LEFT) { 331 Mat A; 332 MatNullSpace nullsp; 333 334 PetscCall(PCGetOperators(ksp->pc, &A, NULL)); 335 PetscCall(MatGetTransposeNullSpace(A, &nullsp)); 336 if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y)); 337 } 338 PetscFunctionReturn(PETSC_SUCCESS); 339 } 340 341 static inline PetscErrorCode KSP_MatMult(KSP ksp, Mat A, Vec x, Vec y) 342 { 343 PetscFunctionBegin; 344 if (ksp->transpose_solve) PetscCall(MatMultTranspose(A, x, y)); 345 else PetscCall(MatMult(A, x, y)); 346 PetscFunctionReturn(PETSC_SUCCESS); 347 } 348 349 static inline PetscErrorCode KSP_MatMultTranspose(KSP ksp, Mat A, Vec x, Vec y) 350 { 351 PetscFunctionBegin; 352 if (ksp->transpose_solve) PetscCall(MatMult(A, x, y)); 353 else PetscCall(MatMultTranspose(A, x, y)); 354 PetscFunctionReturn(PETSC_SUCCESS); 355 } 356 357 static inline PetscErrorCode KSP_MatMultHermitianTranspose(KSP ksp, Mat A, Vec x, Vec y) 358 { 359 PetscFunctionBegin; 360 if (!ksp->transpose_solve) PetscCall(MatMultHermitianTranspose(A, x, y)); 361 else { 362 Vec w; 363 364 PetscCall(VecDuplicate(x, &w)); 365 PetscCall(VecCopy(x, w)); 366 PetscCall(VecConjugate(w)); 367 PetscCall(MatMult(A, w, y)); 368 PetscCall(VecDestroy(&w)); 369 PetscCall(VecConjugate(y)); 370 } 371 PetscFunctionReturn(PETSC_SUCCESS); 372 } 373 374 static inline PetscErrorCode KSP_PCApply(KSP ksp, Vec x, Vec y) 375 { 376 PetscFunctionBegin; 377 if (ksp->transpose_solve) { 378 PetscCall(PCApplyTranspose(ksp->pc, x, y)); 379 PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y)); 380 } else { 381 PetscCall(PCApply(ksp->pc, x, y)); 382 PetscCall(KSP_RemoveNullSpace(ksp, y)); 383 } 384 PetscFunctionReturn(PETSC_SUCCESS); 385 } 386 387 static inline PetscErrorCode KSP_PCApplyTranspose(KSP ksp, Vec x, Vec y) 388 { 389 PetscFunctionBegin; 390 if (ksp->transpose_solve) { 391 PetscCall(PCApply(ksp->pc, x, y)); 392 PetscCall(KSP_RemoveNullSpace(ksp, y)); 393 } else { 394 PetscCall(PCApplyTranspose(ksp->pc, x, y)); 395 PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y)); 396 } 397 PetscFunctionReturn(PETSC_SUCCESS); 398 } 399 400 static inline PetscErrorCode KSP_PCApplyHermitianTranspose(KSP ksp, Vec x, Vec y) 401 { 402 PetscFunctionBegin; 403 PetscCall(VecConjugate(x)); 404 PetscCall(KSP_PCApplyTranspose(ksp, x, y)); 405 PetscCall(VecConjugate(x)); 406 PetscCall(VecConjugate(y)); 407 PetscFunctionReturn(PETSC_SUCCESS); 408 } 409 410 static inline PetscErrorCode KSP_PCApplyBAorAB(KSP ksp, Vec x, Vec y, Vec w) 411 { 412 PetscFunctionBegin; 413 if (ksp->transpose_solve) { 414 PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w)); 415 PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y)); 416 } else { 417 PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w)); 418 PetscCall(KSP_RemoveNullSpace(ksp, y)); 419 } 420 PetscFunctionReturn(PETSC_SUCCESS); 421 } 422 423 static inline PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp, Vec x, Vec y, Vec w) 424 { 425 PetscFunctionBegin; 426 if (ksp->transpose_solve) PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w)); 427 else PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w)); 428 PetscFunctionReturn(PETSC_SUCCESS); 429 } 430 431 PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization; 432 PETSC_EXTERN PetscLogEvent KSP_SetUp; 433 PETSC_EXTERN PetscLogEvent KSP_Solve; 434 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0; 435 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1; 436 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2; 437 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3; 438 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4; 439 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S; 440 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L; 441 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U; 442 PETSC_EXTERN PetscLogEvent KSP_SolveTranspose; 443 PETSC_EXTERN PetscLogEvent KSP_MatSolve; 444 445 PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat, IS, IS, IS, IS, MatReuse, Mat *, MatSchurComplementAinvType, MatReuse, Mat *); 446 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC, PetscBool *); 447 448 /*MC 449 KSPCheckDot - Checks if the result of a dot product used by the corresponding `KSP` contains Inf or NaN. These indicate that the previous 450 application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`. 451 452 Collective 453 454 Input Parameter: 455 . ksp - the linear solver `KSP` context. 456 457 Output Parameter: 458 . beta - the result of the inner product 459 460 Level: developer 461 462 Developer Notes: 463 Used to manage returning from `KSP` solvers whose preconditioners have failed, possibly only a subset of MPI ranks, in some way 464 465 It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products. 466 467 .seealso: `PCFailedReason`, `KSPConvergedReason`, `PCGetFailedReasonRank()`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckNorm()`, `KSPCheckSolve()` 468 M*/ 469 #define KSPCheckDot(ksp, beta) \ 470 do { \ 471 if (PetscIsInfOrNanScalar(beta)) { \ 472 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to Nan or Inf inner product"); \ 473 { \ 474 PCFailedReason pcreason; \ 475 PetscInt sendbuf, recvbuf; \ 476 PetscCall(PCGetFailedReasonRank(ksp->pc, &pcreason)); \ 477 sendbuf = (PetscInt)pcreason; \ 478 PetscCallMPI(MPI_Allreduce(&sendbuf, &recvbuf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ksp))); \ 479 if (recvbuf) { \ 480 PetscCall(PCSetFailedReason(ksp->pc, (PCFailedReason)recvbuf)); \ 481 ksp->reason = KSP_DIVERGED_PC_FAILED; \ 482 PetscCall(VecSetInf(ksp->vec_sol)); \ 483 } else { \ 484 ksp->reason = KSP_DIVERGED_NANORINF; \ 485 } \ 486 PetscFunctionReturn(PETSC_SUCCESS); \ 487 } \ 488 } \ 489 } while (0) 490 491 /*MC 492 KSPCheckNorm - Checks if the result of a norm used by the corresponding `KSP` contains `inf` or `NaN`. These indicate that the previous 493 application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`. 494 495 Collective 496 497 Input Parameter: 498 . ksp - the linear solver `KSP` context. 499 500 Output Parameter: 501 . beta - the result of the norm 502 503 Level: developer 504 505 Developer Notes: 506 Used to manage returning from `KSP` solvers whose preconditioners have failed, possibly only a subset of MPI ranks, in some way. 507 508 It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products. 509 510 .seealso: `PCFailedReason`, `KSPConvergedReason`, `PCGetFailedReasonRank()`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckDot()`, `KSPCheckSolve()` 511 M*/ 512 #define KSPCheckNorm(ksp, beta) \ 513 do { \ 514 if (PetscIsInfOrNanReal(beta)) { \ 515 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to Nan or Inf norm"); \ 516 { \ 517 PCFailedReason pcreason; \ 518 PetscInt sendbuf, recvbuf; \ 519 PetscCall(PCGetFailedReasonRank(ksp->pc, &pcreason)); \ 520 sendbuf = (PetscInt)pcreason; \ 521 PetscCallMPI(MPI_Allreduce(&sendbuf, &recvbuf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ksp))); \ 522 if (recvbuf) { \ 523 PetscCall(PCSetFailedReason(ksp->pc, (PCFailedReason)recvbuf)); \ 524 ksp->reason = KSP_DIVERGED_PC_FAILED; \ 525 PetscCall(VecSetInf(ksp->vec_sol)); \ 526 ksp->rnorm = beta; \ 527 } else { \ 528 PetscCall(PCSetFailedReason(ksp->pc, PC_NOERROR)); \ 529 ksp->reason = KSP_DIVERGED_NANORINF; \ 530 ksp->rnorm = beta; \ 531 } \ 532 PetscFunctionReturn(PETSC_SUCCESS); \ 533 } \ 534 } \ 535 } while (0) 536 537 #endif 538 539 PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]); 540 PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP, PetscInt, PetscReal *); 541