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 PetscInt setfromoptionscalled; 185 PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */ 186 187 PetscErrorCode (*presolve)(KSP, Vec, Vec, void *); 188 PetscErrorCode (*postsolve)(KSP, Vec, Vec, void *); 189 void *prectx, *postctx; 190 }; 191 192 typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */ 193 PetscReal coef; 194 PetscReal bnrm; 195 } KSPDynTolCtx; 196 197 typedef struct { 198 PetscBool initialrtol; /* default relative residual decrease is computed from initial residual, not rhs */ 199 PetscBool mininitialrtol; /* default relative residual decrease is computed from min of initial residual and rhs */ 200 PetscBool convmaxits; /* if true, the convergence test returns KSP_CONVERGED_ITS if the maximum number of iterations is reached */ 201 Vec work; 202 } KSPConvergedDefaultCtx; 203 204 static inline PetscErrorCode KSPLogResidualHistory(KSP ksp, PetscReal norm) 205 { 206 PetscFunctionBegin; 207 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 208 if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) ksp->res_hist[ksp->res_hist_len++] = norm; 209 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 210 PetscFunctionReturn(PETSC_SUCCESS); 211 } 212 213 static inline PetscErrorCode KSPLogErrorHistory(KSP ksp) 214 { 215 DM dm; 216 217 PetscFunctionBegin; 218 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 219 PetscCall(KSPGetDM(ksp, &dm)); 220 if (dm && ksp->err_hist && ksp->err_hist_max > ksp->err_hist_len) { 221 PetscSimplePointFunc exactSol; 222 void *exactCtx; 223 PetscDS ds; 224 Vec u; 225 PetscReal error; 226 PetscInt Nf; 227 228 PetscCall(KSPBuildSolution(ksp, NULL, &u)); 229 /* TODO Was needed to correct for Newton solution, but I just need to set a solution */ 230 //PetscCall(VecScale(u, -1.0)); 231 /* TODO Case when I have a solution */ 232 if (0) { 233 PetscCall(DMGetDS(dm, &ds)); 234 PetscCall(PetscDSGetNumFields(ds, &Nf)); 235 PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle number of fields %" PetscInt_FMT " > 1 right now", Nf); 236 PetscCall(PetscDSGetExactSolution(ds, 0, &exactSol, &exactCtx)); 237 PetscCall(DMComputeL2FieldDiff(dm, 0.0, &exactSol, &exactCtx, u, &error)); 238 } else { 239 /* The null solution A 0 = 0 */ 240 PetscCall(VecNorm(u, NORM_2, &error)); 241 } 242 ksp->err_hist[ksp->err_hist_len++] = error; 243 } 244 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 245 PetscFunctionReturn(PETSC_SUCCESS); 246 } 247 248 static inline PetscScalar KSPNoisyHash_Private(PetscInt xx) 249 { 250 unsigned int x = (unsigned int)xx; 251 x = ((x >> 16) ^ x) * 0x45d9f3b; 252 x = ((x >> 16) ^ x) * 0x45d9f3b; 253 x = ((x >> 16) ^ x); 254 return (PetscScalar)((PetscInt64)x - 2147483648) * 5.e-10; /* center around zero, scaled about -1. to 1.*/ 255 } 256 257 static inline PetscErrorCode KSPSetNoisy_Private(Vec v) 258 { 259 PetscScalar *a; 260 PetscInt n, istart; 261 262 PetscFunctionBegin; 263 PetscCall(VecGetOwnershipRange(v, &istart, NULL)); 264 PetscCall(VecGetLocalSize(v, &n)); 265 PetscCall(VecGetArrayWrite(v, &a)); 266 for (PetscInt i = 0; i < n; ++i) a[i] = KSPNoisyHash_Private(i + istart); 267 PetscCall(VecRestoreArrayWrite(v, &a)); 268 PetscFunctionReturn(PETSC_SUCCESS); 269 } 270 271 PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP, PetscBool, KSPNormType *, PCSide *); 272 273 PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP, PetscInt, const PetscReal *, const PetscReal *); 274 275 typedef struct _p_DMKSP *DMKSP; 276 typedef struct _DMKSPOps *DMKSPOps; 277 struct _DMKSPOps { 278 PetscErrorCode (*computeoperators)(KSP, Mat, Mat, void *); 279 PetscErrorCode (*computerhs)(KSP, Vec, void *); 280 PetscErrorCode (*computeinitialguess)(KSP, Vec, void *); 281 PetscErrorCode (*destroy)(DMKSP *); 282 PetscErrorCode (*duplicate)(DMKSP, DMKSP); 283 }; 284 285 struct _p_DMKSP { 286 PETSCHEADER(struct _DMKSPOps); 287 void *operatorsctx; 288 void *rhsctx; 289 void *initialguessctx; 290 void *data; 291 292 /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way 293 * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user 294 * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by 295 * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates 296 * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes 297 * to the original level will no longer propagate to that level. 298 */ 299 DM originaldm; 300 301 void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */ 302 }; 303 PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM, DMKSP *); 304 PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM, DMKSP *); 305 PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM, DM); 306 307 /* 308 These allow the various Krylov methods to apply to either the linear system or its transpose. 309 */ 310 static inline PetscErrorCode KSP_RemoveNullSpace(KSP ksp, Vec y) 311 { 312 PetscFunctionBegin; 313 if (ksp->pc_side == PC_LEFT) { 314 Mat A; 315 MatNullSpace nullsp; 316 317 PetscCall(PCGetOperators(ksp->pc, &A, NULL)); 318 PetscCall(MatGetNullSpace(A, &nullsp)); 319 if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y)); 320 } 321 PetscFunctionReturn(PETSC_SUCCESS); 322 } 323 324 static inline PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp, Vec y) 325 { 326 PetscFunctionBegin; 327 if (ksp->pc_side == PC_LEFT) { 328 Mat A; 329 MatNullSpace nullsp; 330 331 PetscCall(PCGetOperators(ksp->pc, &A, NULL)); 332 PetscCall(MatGetTransposeNullSpace(A, &nullsp)); 333 if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y)); 334 } 335 PetscFunctionReturn(PETSC_SUCCESS); 336 } 337 338 static inline PetscErrorCode KSP_MatMult(KSP ksp, Mat A, Vec x, Vec y) 339 { 340 PetscFunctionBegin; 341 if (ksp->transpose_solve) PetscCall(MatMultTranspose(A, x, y)); 342 else PetscCall(MatMult(A, x, y)); 343 PetscFunctionReturn(PETSC_SUCCESS); 344 } 345 346 static inline PetscErrorCode KSP_MatMultTranspose(KSP ksp, Mat A, Vec x, Vec y) 347 { 348 PetscFunctionBegin; 349 if (ksp->transpose_solve) PetscCall(MatMult(A, x, y)); 350 else PetscCall(MatMultTranspose(A, x, y)); 351 PetscFunctionReturn(PETSC_SUCCESS); 352 } 353 354 static inline PetscErrorCode KSP_MatMultHermitianTranspose(KSP ksp, Mat A, Vec x, Vec y) 355 { 356 PetscFunctionBegin; 357 if (!ksp->transpose_solve) PetscCall(MatMultHermitianTranspose(A, x, y)); 358 else { 359 Vec w; 360 361 PetscCall(VecDuplicate(x, &w)); 362 PetscCall(VecCopy(x, w)); 363 PetscCall(VecConjugate(w)); 364 PetscCall(MatMult(A, w, y)); 365 PetscCall(VecDestroy(&w)); 366 PetscCall(VecConjugate(y)); 367 } 368 PetscFunctionReturn(PETSC_SUCCESS); 369 } 370 371 static inline PetscErrorCode KSP_PCApply(KSP ksp, Vec x, Vec y) 372 { 373 PetscFunctionBegin; 374 if (ksp->transpose_solve) { 375 PetscCall(PCApplyTranspose(ksp->pc, x, y)); 376 PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y)); 377 } else { 378 PetscCall(PCApply(ksp->pc, x, y)); 379 PetscCall(KSP_RemoveNullSpace(ksp, y)); 380 } 381 PetscFunctionReturn(PETSC_SUCCESS); 382 } 383 384 static inline PetscErrorCode KSP_PCApplyTranspose(KSP ksp, Vec x, Vec y) 385 { 386 PetscFunctionBegin; 387 if (ksp->transpose_solve) { 388 PetscCall(PCApply(ksp->pc, x, y)); 389 PetscCall(KSP_RemoveNullSpace(ksp, y)); 390 } else { 391 PetscCall(PCApplyTranspose(ksp->pc, x, y)); 392 PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y)); 393 } 394 PetscFunctionReturn(PETSC_SUCCESS); 395 } 396 397 static inline PetscErrorCode KSP_PCApplyHermitianTranspose(KSP ksp, Vec x, Vec y) 398 { 399 PetscFunctionBegin; 400 PetscCall(VecConjugate(x)); 401 PetscCall(KSP_PCApplyTranspose(ksp, x, y)); 402 PetscCall(VecConjugate(x)); 403 PetscCall(VecConjugate(y)); 404 PetscFunctionReturn(PETSC_SUCCESS); 405 } 406 407 static inline PetscErrorCode KSP_PCApplyBAorAB(KSP ksp, Vec x, Vec y, Vec w) 408 { 409 PetscFunctionBegin; 410 if (ksp->transpose_solve) { 411 PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w)); 412 PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y)); 413 } else { 414 PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w)); 415 PetscCall(KSP_RemoveNullSpace(ksp, y)); 416 } 417 PetscFunctionReturn(PETSC_SUCCESS); 418 } 419 420 static inline PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp, Vec x, Vec y, Vec w) 421 { 422 PetscFunctionBegin; 423 if (ksp->transpose_solve) PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w)); 424 else PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w)); 425 PetscFunctionReturn(PETSC_SUCCESS); 426 } 427 428 PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization; 429 PETSC_EXTERN PetscLogEvent KSP_SetUp; 430 PETSC_EXTERN PetscLogEvent KSP_Solve; 431 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0; 432 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1; 433 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2; 434 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3; 435 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4; 436 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S; 437 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L; 438 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U; 439 PETSC_EXTERN PetscLogEvent KSP_SolveTranspose; 440 PETSC_EXTERN PetscLogEvent KSP_MatSolve; 441 442 PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat, IS, IS, IS, IS, MatReuse, Mat *, MatSchurComplementAinvType, MatReuse, Mat *); 443 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC, PetscBool *); 444 445 /*MC 446 KSPCheckDot - Checks if the result of a dot product used by the corresponding `KSP` contains Inf or NaN. These indicate that the previous 447 application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`. 448 449 Collective 450 451 Input Parameter: 452 . ksp - the linear solver `KSP` context. 453 454 Output Parameter: 455 . beta - the result of the inner product 456 457 Level: developer 458 459 Developer Notes: 460 Used to manage returning from `KSP` solvers whose preconditioners have failed, possibly only a subset of MPI ranks, in some way 461 462 It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products. 463 464 .seealso: `PCFailedReason`, `KSPConvergedReason`, `PCGetFailedReasonRank()`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckNorm()`, `KSPCheckSolve()` 465 M*/ 466 #define KSPCheckDot(ksp, beta) \ 467 do { \ 468 if (PetscIsInfOrNanScalar(beta)) { \ 469 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to Nan or Inf inner product"); \ 470 { \ 471 PCFailedReason pcreason; \ 472 PetscInt sendbuf, recvbuf; \ 473 PetscCall(PCGetFailedReasonRank(ksp->pc, &pcreason)); \ 474 sendbuf = (PetscInt)pcreason; \ 475 PetscCallMPI(MPI_Allreduce(&sendbuf, &recvbuf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ksp))); \ 476 if (recvbuf) { \ 477 PetscCall(PCSetFailedReason(ksp->pc, (PCFailedReason)recvbuf)); \ 478 ksp->reason = KSP_DIVERGED_PC_FAILED; \ 479 PetscCall(VecSetInf(ksp->vec_sol)); \ 480 } else { \ 481 ksp->reason = KSP_DIVERGED_NANORINF; \ 482 } \ 483 PetscFunctionReturn(PETSC_SUCCESS); \ 484 } \ 485 } \ 486 } while (0) 487 488 /*MC 489 KSPCheckNorm - Checks if the result of a norm used by the corresponding `KSP` contains `inf` or `NaN`. These indicate that the previous 490 application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`. 491 492 Collective 493 494 Input Parameter: 495 . ksp - the linear solver `KSP` context. 496 497 Output Parameter: 498 . beta - the result of the norm 499 500 Level: developer 501 502 Developer Notes: 503 Used to manage returning from `KSP` solvers whose preconditioners have failed, possibly only a subset of MPI ranks, in some way. 504 505 It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products. 506 507 .seealso: `PCFailedReason`, `KSPConvergedReason`, `PCGetFailedReasonRank()`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckDot()`, `KSPCheckSolve()` 508 M*/ 509 #define KSPCheckNorm(ksp, beta) \ 510 do { \ 511 if (PetscIsInfOrNanReal(beta)) { \ 512 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to Nan or Inf norm"); \ 513 { \ 514 PCFailedReason pcreason; \ 515 PetscInt sendbuf, recvbuf; \ 516 PetscCall(PCGetFailedReasonRank(ksp->pc, &pcreason)); \ 517 sendbuf = (PetscInt)pcreason; \ 518 PetscCallMPI(MPI_Allreduce(&sendbuf, &recvbuf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ksp))); \ 519 if (recvbuf) { \ 520 PetscCall(PCSetFailedReason(ksp->pc, (PCFailedReason)recvbuf)); \ 521 ksp->reason = KSP_DIVERGED_PC_FAILED; \ 522 PetscCall(VecSetInf(ksp->vec_sol)); \ 523 ksp->rnorm = beta; \ 524 } else { \ 525 PetscCall(PCSetFailedReason(ksp->pc, PC_NOERROR)); \ 526 ksp->reason = KSP_DIVERGED_NANORINF; \ 527 ksp->rnorm = beta; \ 528 } \ 529 PetscFunctionReturn(PETSC_SUCCESS); \ 530 } \ 531 } \ 532 } while (0) 533 534 #endif 535 536 PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]); 537 PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP, PetscInt, PetscReal *); 538