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