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 {KSP_SETUP_NEW, KSP_SETUP_NEWMATRIX, KSP_SETUP_NEWRHS} KSPSetUpStage; 74 75 /* 76 Defines the KSP data structure. 77 */ 78 struct _p_KSP { 79 PETSCHEADER(struct _KSPOps); 80 DM dm; 81 PetscBool dmAuto; /* DM was created automatically by KSP */ 82 PetscBool dmActive; /* KSP should use DM for computing operators */ 83 /*------------------------- User parameters--------------------------*/ 84 PetscInt max_it; /* maximum number of iterations */ 85 KSPGuess guess; 86 PetscBool guess_zero, /* flag for whether initial guess is 0 */ 87 calc_sings, /* calculate extreme Singular Values */ 88 calc_ritz, /* calculate (harmonic) Ritz pairs */ 89 guess_knoll; /* use initial guess of PCApply(ksp->B,b */ 90 PCSide pc_side; /* flag for left, right, or symmetric preconditioning */ 91 PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */ 92 PetscReal rtol, /* relative tolerance */ 93 abstol, /* absolute tolerance */ 94 ttol, /* (not set by user) */ 95 divtol; /* divergence tolerance */ 96 PetscReal rnorm0; /* initial residual norm (used for divergence testing) */ 97 PetscReal rnorm; /* current residual norm */ 98 KSPConvergedReason reason; 99 PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */ 100 101 Vec vec_sol,vec_rhs; /* pointer to where user has stashed 102 the solution and rhs, these are 103 never touched by the code, only 104 passed back to the user */ 105 PetscReal *res_hist; /* If !0 stores residual each at iteration */ 106 PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */ 107 size_t res_hist_len; /* current size of residual history array */ 108 size_t res_hist_max; /* actual amount of storage in residual history */ 109 PetscBool res_hist_reset; /* reset history to length zero for each new solve */ 110 PetscReal *err_hist; /* If !0 stores error at each iteration */ 111 PetscReal *err_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */ 112 size_t err_hist_len; /* current size of error history array */ 113 size_t err_hist_max; /* actual amount of storage in error history */ 114 PetscBool err_hist_reset; /* reset history to length zero for each new solve */ 115 116 PetscInt chknorm; /* only compute/check norm if iterations is great than this */ 117 PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the 118 MPI_Allreduce() for computing the inner products for the next iteration. */ 119 120 PetscInt nmax; /* maximum number of right-hand sides to be handled simultaneously */ 121 122 /* --------User (or default) routines (most return -1 on error) --------*/ 123 PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */ 124 PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void**); /* */ 125 void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */ 126 PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */ 127 PetscBool pauseFinal; /* Pause all drawing monitor at the final iterate */ 128 129 PetscErrorCode (*reasonview[MAXKSPREASONVIEWS])(KSP,void*); /* KSP converged reason view */ 130 PetscErrorCode (*reasonviewdestroy[MAXKSPREASONVIEWS])(void**); /* Optional destroy routine */ 131 void *reasonviewcontext[MAXKSPREASONVIEWS]; /* User context */ 132 PetscInt numberreasonviews; /* Number if reason viewers */ 133 134 PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*); 135 PetscErrorCode (*convergeddestroy)(void*); 136 void *cnvP; 137 138 void *user; /* optional user-defined context */ 139 140 PC pc; 141 142 void *data; /* holder for misc stuff associated 143 with a particular iterative solver */ 144 145 PetscBool view, viewPre, viewRate, viewMat, viewPMat, viewRhs, viewSol, viewMatExp, viewEV, viewSV, viewEVExp, viewFinalRes, viewPOpExp, viewDScale; 146 PetscViewer viewer, viewerPre, viewerRate, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp, viewerDScale; 147 PetscViewerFormat format, formatPre, formatRate, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp, formatDScale; 148 149 /* ----------------Default work-area management -------------------- */ 150 PetscInt nwork; 151 Vec *work; 152 153 KSPSetUpStage setupstage; 154 PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */ 155 156 PetscInt its; /* number of iterations so far computed in THIS linear solve*/ 157 PetscInt totalits; /* number of iterations used by this KSP object since it was created */ 158 159 PetscBool transpose_solve; /* solve transpose system instead */ 160 struct { 161 Mat AT,BT; 162 PetscBool use_explicittranspose; /* transpose the system explicitly in KSPSolveTranspose */ 163 PetscBool reuse_transpose; /* reuse the previous transposed system */ 164 } transpose; 165 166 KSPNormType normtype; /* type of norm used for convergence tests */ 167 168 PCSide pc_side_set; /* PC type set explicitly by user */ 169 KSPNormType normtype_set; /* Norm type set explicitly by user */ 170 171 /* Allow diagonally scaling the matrix before computing the preconditioner or using 172 the Krylov method. Note this is NOT just Jacobi preconditioning */ 173 174 PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */ 175 PetscBool dscalefix; /* unscale system after solve */ 176 PetscBool dscalefix2; /* system has been unscaled */ 177 Vec diagonal; /* 1/sqrt(diag of matrix) */ 178 Vec truediagonal; 179 180 PetscInt setfromoptionscalled; 181 PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */ 182 183 PetscViewer eigviewer; /* Viewer where computed eigenvalues are displayed */ 184 185 PetscErrorCode (*presolve)(KSP,Vec,Vec,void*); 186 PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*); 187 void *prectx,*postctx; 188 }; 189 190 typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */ 191 PetscReal coef; 192 PetscReal bnrm; 193 } KSPDynTolCtx; 194 195 typedef struct { 196 PetscBool initialrtol; /* default relative residual decrease is computed from initial residual, not rhs */ 197 PetscBool mininitialrtol; /* default relative residual decrease is computed from min of initial residual and rhs */ 198 PetscBool convmaxits; /* if true, the convergence test returns KSP_CONVERGED_ITS if the maximum number of iterations is reached */ 199 Vec work; 200 } KSPConvergedDefaultCtx; 201 202 static inline PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm) 203 { 204 PetscFunctionBegin; 205 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 206 if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) { 207 ksp->res_hist[ksp->res_hist_len++] = norm; 208 } 209 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 210 PetscFunctionReturn(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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 on ksp 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) do { \ 467 if (PetscIsInfOrNanScalar(beta)) { \ 468 PetscCheck(!ksp->errorifnotconverged,PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");\ 469 {\ 470 PCFailedReason pcreason;\ 471 PetscInt sendbuf,recvbuf; \ 472 PetscCall(PCGetFailedReasonRank(ksp->pc,&pcreason));\ 473 sendbuf = (PetscInt)pcreason; \ 474 PetscCallMPI(MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp)));\ 475 if (recvbuf) { \ 476 PetscCall(PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf)); \ 477 ksp->reason = KSP_DIVERGED_PC_FAILED;\ 478 PetscCall(VecSetInf(ksp->vec_sol));\ 479 } else {\ 480 ksp->reason = KSP_DIVERGED_NANORINF;\ 481 }\ 482 PetscFunctionReturn(0);\ 483 }\ 484 } } while (0) 485 486 /*MC 487 KSPCheckNorm - Checks if the result of a norm used by the corresponding KSP contains Inf or NaN. These indicate that the previous 488 application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`. 489 490 Collective on ksp 491 492 Input Parameter: 493 . ksp - the linear solver (KSP) context. 494 495 Output Parameter: 496 . beta - the result of the norm 497 498 Level: developer 499 500 Developer Notes: 501 Used to manage returning from `KSP` solvers whose preconditioners have failed, possibly only a subset of MPI ranks, in some way. 502 503 It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products. 504 505 .seealso: `PCFailedReason`, `KSPConvergedReason`, `PCGetFailedReasonRank()`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckDot()`, `KSPCheckSolve()` 506 M*/ 507 #define KSPCheckNorm(ksp,beta) do { \ 508 if (PetscIsInfOrNanReal(beta)) { \ 509 PetscCheck(!ksp->errorifnotconverged,PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf norm");\ 510 {\ 511 PCFailedReason pcreason;\ 512 PetscInt sendbuf,recvbuf; \ 513 PetscCall(PCGetFailedReasonRank(ksp->pc,&pcreason));\ 514 sendbuf = (PetscInt)pcreason; \ 515 PetscCallMPI(MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp)));\ 516 if (recvbuf) { \ 517 PetscCall(PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf)); \ 518 ksp->reason = KSP_DIVERGED_PC_FAILED; \ 519 PetscCall(VecSetInf(ksp->vec_sol));\ 520 ksp->rnorm = beta; \ 521 } else {\ 522 PetscCall(PCSetFailedReason(ksp->pc,PC_NOERROR)); \ 523 ksp->reason = KSP_DIVERGED_NANORINF;\ 524 ksp->rnorm = beta; \ 525 } \ 526 PetscFunctionReturn(0);\ 527 }\ 528 } } while (0) 529 530 #endif 531 532 PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]); 533 PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP,PetscInt,PetscReal*); 534