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