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)(PetscOptionItems*,KSP); 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 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 Note: 460 this is used to manage returning from KSP solvers whose preconditioners have failed in some way 461 462 .seealso: `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckNorm()`, `KSPCheckSolve()` 463 M*/ 464 #define KSPCheckDot(ksp,beta) do { \ 465 if (PetscIsInfOrNanScalar(beta)) { \ 466 PetscCheck(!ksp->errorifnotconverged,PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");\ 467 {\ 468 PCFailedReason pcreason;\ 469 PetscInt sendbuf,recvbuf; \ 470 PetscCall(PCGetFailedReasonRank(ksp->pc,&pcreason));\ 471 sendbuf = (PetscInt)pcreason; \ 472 PetscCallMPI(MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp)));\ 473 if (recvbuf) { \ 474 PetscCall(PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf)); \ 475 ksp->reason = KSP_DIVERGED_PC_FAILED;\ 476 PetscCall(VecSetInf(ksp->vec_sol));\ 477 } else {\ 478 ksp->reason = KSP_DIVERGED_NANORINF;\ 479 }\ 480 PetscFunctionReturn(0);\ 481 }\ 482 } } while (0) 483 484 /*MC 485 KSPCheckNorm - Checks if the result of a norm used by the corresponding KSP contains Inf or NaN. These indicate that the previous 486 application of the preconditioner generated an error 487 488 Collective on ksp 489 490 Input Parameter: 491 . ksp - the linear solver (KSP) context. 492 493 Output Parameter: 494 . beta - the result of the norm 495 496 Level: developer 497 498 Developer Note: 499 this is used to manage returning from KSP solvers whose preconditioners have failed in some way 500 501 .seealso: `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckDot()`, `KSPCheckSolve()` 502 M*/ 503 #define KSPCheckNorm(ksp,beta) do { \ 504 if (PetscIsInfOrNanReal(beta)) { \ 505 PetscCheck(!ksp->errorifnotconverged,PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf norm");\ 506 {\ 507 PCFailedReason pcreason;\ 508 PetscInt sendbuf,recvbuf; \ 509 PetscCall(PCGetFailedReasonRank(ksp->pc,&pcreason));\ 510 sendbuf = (PetscInt)pcreason; \ 511 PetscCallMPI(MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp)));\ 512 if (recvbuf) { \ 513 PetscCall(PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf)); \ 514 ksp->reason = KSP_DIVERGED_PC_FAILED; \ 515 PetscCall(VecSetInf(ksp->vec_sol));\ 516 ksp->rnorm = beta; \ 517 } else {\ 518 PetscCall(PCSetFailedReason(ksp->pc,PC_NOERROR)); \ 519 ksp->reason = KSP_DIVERGED_NANORINF;\ 520 ksp->rnorm = beta; \ 521 } \ 522 PetscFunctionReturn(0);\ 523 }\ 524 } } while (0) 525 526 #endif 527 528 PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]); 529 PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP,PetscInt,PetscReal*); 530