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 PetscInt min_it; /* minimum number of iterations */ 90 KSPGuess guess; 91 PetscBool guess_zero, /* flag for whether initial guess is 0 */ 92 guess_not_read, /* guess is not read, does not need to be zeroed */ 93 calc_sings, /* calculate extreme Singular Values */ 94 calc_ritz, /* calculate (harmonic) Ritz pairs */ 95 guess_knoll; /* use initial guess of PCApply(ksp->B,b */ 96 PCSide pc_side; /* flag for left, right, or symmetric preconditioning */ 97 PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */ 98 PetscReal rtol, /* relative tolerance */ 99 abstol, /* absolute tolerance */ 100 ttol, /* (not set by user) */ 101 divtol; /* divergence tolerance */ 102 PetscReal rnorm0; /* initial residual norm (used for divergence testing) */ 103 PetscReal rnorm; /* current residual norm */ 104 KSPConvergedReason reason; 105 PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */ 106 107 Vec vec_sol, vec_rhs; /* pointer to where user has stashed 108 the solution and rhs, these are 109 never touched by the code, only 110 passed back to the user */ 111 PetscReal *res_hist; /* If !0 stores residual each at iteration */ 112 PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */ 113 size_t res_hist_len; /* current size of residual history array */ 114 size_t res_hist_max; /* actual amount of storage in residual history */ 115 PetscBool res_hist_reset; /* reset history to length zero for each new solve */ 116 PetscReal *err_hist; /* If !0 stores error at each iteration */ 117 PetscReal *err_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */ 118 size_t err_hist_len; /* current size of error history array */ 119 size_t err_hist_max; /* actual amount of storage in error history */ 120 PetscBool err_hist_reset; /* reset history to length zero for each new solve */ 121 122 PetscInt chknorm; /* only compute/check norm if iterations is great than this */ 123 PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the 124 MPI_Allreduce() for computing the inner products for the next iteration. */ 125 126 PetscInt nmax; /* maximum number of right-hand sides to be handled simultaneously */ 127 128 /* --------User (or default) routines (most return -1 on error) --------*/ 129 PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP, PetscInt, PetscReal, void *); /* returns control to user after */ 130 PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void **); /* */ 131 void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */ 132 PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */ 133 PetscBool pauseFinal; /* Pause all drawing monitor at the final iterate */ 134 135 PetscErrorCode (*reasonview[MAXKSPREASONVIEWS])(KSP, void *); /* KSP converged reason view */ 136 PetscErrorCode (*reasonviewdestroy[MAXKSPREASONVIEWS])(void **); /* Optional destroy routine */ 137 void *reasonviewcontext[MAXKSPREASONVIEWS]; /* User context */ 138 PetscInt numberreasonviews; /* Number if reason viewers */ 139 140 PetscErrorCode (*converged)(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *); 141 PetscErrorCode (*convergeddestroy)(void *); 142 void *cnvP; 143 144 void *user; /* optional user-defined context */ 145 146 PC pc; 147 148 void *data; /* holder for misc stuff associated 149 with a particular iterative solver */ 150 151 PetscBool view, viewPre, viewRate, viewMat, viewPMat, viewRhs, viewSol, viewMatExp, viewEV, viewSV, viewEVExp, viewFinalRes, viewPOpExp, viewDScale; 152 PetscViewer viewer, viewerPre, viewerRate, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp, viewerDScale; 153 PetscViewerFormat format, formatPre, formatRate, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp, formatDScale; 154 155 /* ----------------Default work-area management -------------------- */ 156 PetscInt nwork; 157 Vec *work; 158 159 KSPSetUpStage setupstage; 160 PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */ 161 162 PetscInt its; /* number of iterations so far computed in THIS linear solve*/ 163 PetscInt totalits; /* number of iterations used by this KSP object since it was created */ 164 165 PetscBool transpose_solve; /* solve transpose system instead */ 166 struct { 167 Mat AT, BT; 168 PetscBool use_explicittranspose; /* transpose the system explicitly in KSPSolveTranspose */ 169 PetscBool reuse_transpose; /* reuse the previous transposed system */ 170 } transpose; 171 172 KSPNormType normtype; /* type of norm used for convergence tests */ 173 174 PCSide pc_side_set; /* PC type set explicitly by user */ 175 KSPNormType normtype_set; /* Norm type set explicitly by user */ 176 177 /* Allow diagonally scaling the matrix before computing the preconditioner or using 178 the Krylov method. Note this is NOT just Jacobi preconditioning */ 179 180 PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */ 181 PetscBool dscalefix; /* unscale system after solve */ 182 PetscBool dscalefix2; /* system has been unscaled */ 183 Vec diagonal; /* 1/sqrt(diag of matrix) */ 184 Vec truediagonal; 185 186 /* Allow declaring convergence when negative curvature is detected */ 187 PetscBool converged_neg_curve; 188 189 PetscInt setfromoptionscalled; 190 PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */ 191 192 PetscErrorCode (*presolve)(KSP, Vec, Vec, void *); 193 PetscErrorCode (*postsolve)(KSP, Vec, Vec, void *); 194 void *prectx, *postctx; 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 PetscInt sendbuf, recvbuf; \ 503 PetscCall(PCGetFailedReasonRank(ksp->pc, &pcreason)); \ 504 sendbuf = (PetscInt)pcreason; \ 505 PetscCall(MPIU_Allreduce(&sendbuf, &recvbuf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ksp))); \ 506 if (recvbuf) { \ 507 PetscCall(PCSetFailedReason(ksp->pc, (PCFailedReason)recvbuf)); \ 508 ksp->reason = KSP_DIVERGED_PC_FAILED; \ 509 PetscCall(VecSetInf(ksp->vec_sol)); \ 510 } else { \ 511 ksp->reason = KSP_DIVERGED_NANORINF; \ 512 } \ 513 PetscFunctionReturn(PETSC_SUCCESS); \ 514 } \ 515 } \ 516 } while (0) 517 518 /*MC 519 KSPCheckNorm - Checks if the result of a norm used by the corresponding `KSP` contains `inf` or `NaN`. These indicate that the previous 520 application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`. 521 522 Collective 523 524 Input Parameter: 525 . ksp - the linear solver `KSP` context. 526 527 Output Parameter: 528 . beta - the result of the norm 529 530 Level: developer 531 532 Developer Notes: 533 Used to manage returning from `KSP` solvers whose preconditioners have failed, possibly only a subset of MPI ranks, in some way. 534 535 It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products. 536 537 .seealso: `PCFailedReason`, `KSPConvergedReason`, `PCGetFailedReasonRank()`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckDot()`, `KSPCheckSolve()` 538 M*/ 539 #define KSPCheckNorm(ksp, beta) \ 540 do { \ 541 if (PetscIsInfOrNanReal(beta)) { \ 542 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to Nan or Inf norm"); \ 543 { \ 544 PCFailedReason pcreason; \ 545 PetscInt sendbuf, recvbuf; \ 546 PetscCall(PCGetFailedReasonRank(ksp->pc, &pcreason)); \ 547 sendbuf = (PetscInt)pcreason; \ 548 PetscCall(MPIU_Allreduce(&sendbuf, &recvbuf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ksp))); \ 549 if (recvbuf) { \ 550 PetscCall(PCSetFailedReason(ksp->pc, (PCFailedReason)recvbuf)); \ 551 ksp->reason = KSP_DIVERGED_PC_FAILED; \ 552 PetscCall(VecSetInf(ksp->vec_sol)); \ 553 ksp->rnorm = beta; \ 554 } else { \ 555 PetscCall(PCSetFailedReason(ksp->pc, PC_NOERROR)); \ 556 ksp->reason = KSP_DIVERGED_NANORINF; \ 557 ksp->rnorm = beta; \ 558 } \ 559 PetscFunctionReturn(PETSC_SUCCESS); \ 560 } \ 561 } \ 562 } while (0) 563 564 #endif 565 566 PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]); 567 PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP, PetscInt, PetscReal *); 568