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 PetscObjectParameterDeclare(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 PetscObjectParameterDeclare(PetscReal, rtol); /* relative tolerance */ 97 PetscObjectParameterDeclare(PetscReal, abstol); /* absolute tolerance */ 98 PetscObjectParameterDeclare(PetscReal, ttol); /* (not set by user) */ 99 PetscObjectParameterDeclare(PetscReal, 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 PetscCount res_hist_len; /* current entry count of residual history array */ 112 PetscCount res_hist_max; /* total entry count 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 PetscCount err_hist_len; /* current entry count of error history array */ 117 PetscCount err_hist_max; /* total entry count 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 KSPMonitorFn *monitor[MAXKSPMONITORS]; 128 PetscCtxDestroyFn *monitordestroy[MAXKSPMONITORS]; 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 PetscViewer convergedreasonviewer; 134 PetscViewerFormat convergedreasonformat; 135 KSPConvergedReasonViewFn *reasonview[MAXKSPREASONVIEWS]; /* KSP converged reason view */ 136 PetscCtxDestroyFn *reasonviewdestroy[MAXKSPREASONVIEWS]; /* optional destroy routine */ 137 void *reasonviewcontext[MAXKSPREASONVIEWS]; /* viewer context */ 138 PetscInt numberreasonviews; /* current number of reason viewers */ 139 140 KSPConvergenceTestFn *converged; 141 PetscCtxDestroyFn *convergeddestroy; 142 void *cnvP; 143 144 void *ctx; /* optional user-defined context */ 145 146 PC pc; 147 148 void *data; /* holder for misc stuff associated with a particular iterative solver */ 149 150 PetscBool view, viewPre, viewRate, viewMat, viewPMat, viewRhs, viewSol, viewMatExp, viewEV, viewSV, viewEVExp, viewFinalRes, viewPOpExp, viewDScale; 151 PetscViewer viewer, viewerPre, viewerRate, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp, viewerDScale; 152 PetscViewerFormat format, formatPre, formatRate, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp, formatDScale; 153 154 /* ----------------Default work-area management -------------------- */ 155 PetscInt nwork; 156 Vec *work; 157 158 KSPSetUpStage setupstage; 159 PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */ 160 161 PetscInt its; /* number of iterations so far computed in THIS linear solve*/ 162 PetscInt totalits; /* number of iterations used by this KSP object since it was created */ 163 164 PetscBool transpose_solve; /* solve transpose system instead */ 165 struct { 166 Mat AT, BT; 167 PetscBool use_explicittranspose; /* transpose the system explicitly in KSPSolveTranspose */ 168 PetscBool reuse_transpose; /* reuse the previous transposed system */ 169 } transpose; 170 171 KSPNormType normtype; /* type of norm used for convergence tests */ 172 173 PCSide pc_side_set; /* PC type set explicitly by user */ 174 KSPNormType normtype_set; /* Norm type set explicitly by user */ 175 176 /* Allow diagonally scaling the matrix before computing the preconditioner or using 177 the Krylov method. Note this is NOT just Jacobi preconditioning */ 178 179 PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */ 180 PetscBool dscalefix; /* unscale system after solve */ 181 PetscBool dscalefix2; /* system has been unscaled */ 182 Vec diagonal; /* 1/sqrt(diag of matrix) */ 183 Vec truediagonal; 184 185 /* Allow declaring convergence when negative curvature is detected */ 186 PetscBool converged_neg_curve; 187 188 PetscInt setfromoptionscalled; 189 PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */ 190 191 PetscErrorCode (*presolve)(KSP, Vec, Vec, void *); 192 PetscErrorCode (*postsolve)(KSP, Vec, Vec, void *); 193 void *prectx, *postctx; 194 195 PetscInt nestlevel; /* how many levels of nesting does the KSP have */ 196 }; 197 198 typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */ 199 PetscReal coef; 200 PetscReal bnrm; 201 } KSPDynTolCtx; 202 203 typedef struct { 204 PetscBool initialrtol; /* default relative residual decrease is computed from initial residual, not rhs */ 205 PetscBool mininitialrtol; /* default relative residual decrease is computed from min of initial residual and rhs */ 206 PetscBool convmaxits; /* if true, the convergence test returns KSP_CONVERGED_ITS if the maximum number of iterations is reached */ 207 Vec work; 208 } KSPConvergedDefaultCtx; 209 210 static inline PetscErrorCode KSPLogResidualHistory(KSP ksp, PetscReal norm) 211 { 212 PetscFunctionBegin; 213 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 214 if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) ksp->res_hist[ksp->res_hist_len++] = norm; 215 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 216 PetscFunctionReturn(PETSC_SUCCESS); 217 } 218 219 static inline PetscErrorCode KSPLogErrorHistory(KSP ksp) 220 { 221 DM dm; 222 223 PetscFunctionBegin; 224 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp)); 225 PetscCall(KSPGetDM(ksp, &dm)); 226 if (dm && ksp->err_hist && ksp->err_hist_max > ksp->err_hist_len) { 227 PetscSimplePointFn *exactSol; 228 void *exactCtx; 229 PetscDS ds; 230 Vec u; 231 PetscReal error; 232 PetscInt Nf; 233 234 PetscCall(KSPBuildSolution(ksp, NULL, &u)); 235 /* TODO Was needed to correct for Newton solution, but I just need to set a solution */ 236 //PetscCall(VecScale(u, -1.0)); 237 /* TODO Case when I have a solution */ 238 if (0) { 239 PetscCall(DMGetDS(dm, &ds)); 240 PetscCall(PetscDSGetNumFields(ds, &Nf)); 241 PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle number of fields %" PetscInt_FMT " > 1 right now", Nf); 242 PetscCall(PetscDSGetExactSolution(ds, 0, &exactSol, &exactCtx)); 243 PetscCall(DMComputeL2FieldDiff(dm, 0.0, &exactSol, &exactCtx, u, &error)); 244 } else { 245 /* The null solution A 0 = 0 */ 246 PetscCall(VecNorm(u, NORM_2, &error)); 247 } 248 ksp->err_hist[ksp->err_hist_len++] = error; 249 } 250 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp)); 251 PetscFunctionReturn(PETSC_SUCCESS); 252 } 253 254 static inline PetscScalar KSPNoisyHash_Private(PetscInt xx) 255 { 256 unsigned int x = (unsigned int)xx; 257 x = ((x >> 16) ^ x) * 0x45d9f3b; 258 x = ((x >> 16) ^ x) * 0x45d9f3b; 259 x = ((x >> 16) ^ x); 260 return (PetscScalar)(((PetscInt64)x - 2147483648) * 5.e-10); /* center around zero, scaled about -1. to 1.*/ 261 } 262 263 static inline PetscErrorCode KSPSetNoisy_Private(Mat A, Vec v) 264 { 265 PetscScalar *a; 266 PetscInt n, istart; 267 MatNullSpace nullsp = NULL; 268 269 PetscFunctionBegin; 270 if (A) PetscCall(MatGetNullSpace(A, &nullsp)); 271 PetscCall(VecGetOwnershipRange(v, &istart, NULL)); 272 PetscCall(VecGetLocalSize(v, &n)); 273 PetscCall(VecGetArrayWrite(v, &a)); 274 for (PetscInt i = 0; i < n; ++i) a[i] = KSPNoisyHash_Private(i + istart); 275 PetscCall(VecRestoreArrayWrite(v, &a)); 276 if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, v)); 277 PetscFunctionReturn(PETSC_SUCCESS); 278 } 279 280 PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP, PetscBool, KSPNormType *, PCSide *); 281 282 PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP, PetscInt, const PetscReal *, const PetscReal *); 283 284 typedef struct _p_DMKSP *DMKSP; 285 typedef struct _DMKSPOps *DMKSPOps; 286 struct _DMKSPOps { 287 KSPComputeOperatorsFn *computeoperators; 288 KSPComputeRHSFn *computerhs; 289 KSPComputeInitialGuessFn *computeinitialguess; 290 PetscErrorCode (*destroy)(DMKSP *); 291 PetscErrorCode (*duplicate)(DMKSP, DMKSP); 292 }; 293 294 /*S 295 DMKSP - Object held by a `DM` that contains all the callback functions and their contexts needed by a `KSP` 296 297 Level: developer 298 299 Notes: 300 Users provides callback functions and their contexts to `KSP` using, for example, `KSPSetComputeRHS()`. These values are stored 301 in a `DMKSP` that is contained in the `DM` associated with the `KSP`. If no `DM` was provided by 302 the user with `KSPSetDM()` it is automatically created by `KSPGetDM()` with `DMShellCreate()`. 303 304 Users very rarely need to worked directly with the `DMKSP` object, rather they work with the `KSP` and the `DM` they created 305 306 Multiple `DM` can share a single `DMKSP`, often each `DM` is associated with 307 a grid refinement level. `DMGetDMKSP()` returns the `DMKSP` associated with a `DM`. `DMGetDMKSPWrite()` returns a unique 308 `DMKSP` that is only associated with the current `DM`, making a copy of the shared `DMKSP` if needed (copy-on-write). 309 310 Developer Notes: 311 It is rather subtle why `DMKSP`, `DMSNES`, and `DMTS` are needed instead of simply storing the user callback functions and contexts in `DM` or `KSP`, `SNES`, or `TS`. 312 It is to support composable solvers such as geometric multigrid. We want, by default, the same callback functions and contexts for all the levels in the computation, 313 but we need to also support different callbacks and contexts on each level. The copy-on-write approach of `DMGetDMKSPWrite()` makes this possible. 314 315 The `originaldm` inside the `DMKSP` is NOT reference counted (to prevent a reference count loop between a `DM` and a `DMKSP`). 316 The `DM` on which this context was first created is cached here to implement one-way 317 copy-on-write. When `DMGetDMKSPWrite()` sees a request using a different `DM`, it makes a copy of the `TSDM`. Thus, if a user 318 only interacts directly with one level, e.g., using `TSSetIFunction()`, then coarse levels of a multilevel item 319 integrator are built, then the user changes the routine with another call to `TSSetIFunction()`, it automatically 320 propagates to all the levels. If instead, they get out a specific level and set the function on that level, 321 subsequent changes to the original level will no longer propagate to that level. 322 323 .seealso: [](ch_ts), `KSP`, `KSPCreate()`, `DM`, `DMGetDMKSPWrite()`, `DMGetDMKSP()`, `DMSNES`, `DMTS`, `DMKSPSetComputeOperators()`, `DMKSPGetComputeOperators()`, 324 `DMKSPSetComputeRHS()`, `DMKSPSetComputeInitialGuess()` 325 S*/ 326 struct _p_DMKSP { 327 PETSCHEADER(struct _DMKSPOps); 328 void *operatorsctx; 329 void *rhsctx; 330 void *initialguessctx; 331 void *data; 332 333 /* See developer note for `DMKSP` above */ 334 DM originaldm; 335 336 void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */ 337 }; 338 PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM, DMKSP *); 339 PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM, DMKSP *); 340 PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM, DM); 341 342 /* 343 These allow the various Krylov methods to apply to either the linear system or its transpose. 344 */ 345 static inline PetscErrorCode KSP_RemoveNullSpace(KSP ksp, Vec y) 346 { 347 PetscFunctionBegin; 348 if (ksp->pc_side == PC_LEFT) { 349 Mat A; 350 MatNullSpace nullsp; 351 352 PetscCall(PCGetOperators(ksp->pc, &A, NULL)); 353 PetscCall(MatGetNullSpace(A, &nullsp)); 354 if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y)); 355 } 356 PetscFunctionReturn(PETSC_SUCCESS); 357 } 358 359 static inline PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp, Vec y) 360 { 361 PetscFunctionBegin; 362 if (ksp->pc_side == PC_LEFT) { 363 Mat A; 364 MatNullSpace nullsp; 365 366 PetscCall(PCGetOperators(ksp->pc, &A, NULL)); 367 PetscCall(MatGetTransposeNullSpace(A, &nullsp)); 368 if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, y)); 369 } 370 PetscFunctionReturn(PETSC_SUCCESS); 371 } 372 373 static inline PetscErrorCode KSP_MatMult(KSP ksp, Mat A, Vec x, Vec y) 374 { 375 PetscFunctionBegin; 376 if (ksp->transpose_solve) PetscCall(MatMultTranspose(A, x, y)); 377 else PetscCall(MatMult(A, x, y)); 378 PetscFunctionReturn(PETSC_SUCCESS); 379 } 380 381 static inline PetscErrorCode KSP_MatMultTranspose(KSP ksp, Mat A, Vec x, Vec y) 382 { 383 PetscFunctionBegin; 384 if (ksp->transpose_solve) PetscCall(MatMult(A, x, y)); 385 else PetscCall(MatMultTranspose(A, x, y)); 386 PetscFunctionReturn(PETSC_SUCCESS); 387 } 388 389 static inline PetscErrorCode KSP_MatMultHermitianTranspose(KSP ksp, Mat A, Vec x, Vec y) 390 { 391 PetscFunctionBegin; 392 if (!ksp->transpose_solve) PetscCall(MatMultHermitianTranspose(A, x, y)); 393 else { 394 Vec w; 395 396 PetscCall(VecDuplicate(x, &w)); 397 PetscCall(VecCopy(x, w)); 398 PetscCall(VecConjugate(w)); 399 PetscCall(MatMult(A, w, y)); 400 PetscCall(VecDestroy(&w)); 401 PetscCall(VecConjugate(y)); 402 } 403 PetscFunctionReturn(PETSC_SUCCESS); 404 } 405 406 static inline PetscErrorCode KSP_PCApply(KSP ksp, Vec x, Vec y) 407 { 408 PetscFunctionBegin; 409 if (ksp->transpose_solve) { 410 PetscCall(PCApplyTranspose(ksp->pc, x, y)); 411 PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y)); 412 } else { 413 PetscCall(PCApply(ksp->pc, x, y)); 414 PetscCall(KSP_RemoveNullSpace(ksp, y)); 415 } 416 PetscFunctionReturn(PETSC_SUCCESS); 417 } 418 419 static inline PetscErrorCode KSP_PCApplyTranspose(KSP ksp, Vec x, Vec y) 420 { 421 PetscFunctionBegin; 422 if (ksp->transpose_solve) { 423 PetscCall(PCApply(ksp->pc, x, y)); 424 PetscCall(KSP_RemoveNullSpace(ksp, y)); 425 } else { 426 PetscCall(PCApplyTranspose(ksp->pc, x, y)); 427 PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y)); 428 } 429 PetscFunctionReturn(PETSC_SUCCESS); 430 } 431 432 static inline PetscErrorCode KSP_PCApplyHermitianTranspose(KSP ksp, Vec x, Vec y) 433 { 434 PetscFunctionBegin; 435 PetscCall(VecConjugate(x)); 436 PetscCall(KSP_PCApplyTranspose(ksp, x, y)); 437 PetscCall(VecConjugate(x)); 438 PetscCall(VecConjugate(y)); 439 PetscFunctionReturn(PETSC_SUCCESS); 440 } 441 442 static inline PetscErrorCode KSP_PCMatApply(KSP ksp, Mat X, Mat Y) 443 { 444 PetscFunctionBegin; 445 if (ksp->transpose_solve) PetscCall(PCMatApplyTranspose(ksp->pc, X, Y)); 446 else PetscCall(PCMatApply(ksp->pc, X, Y)); 447 PetscFunctionReturn(PETSC_SUCCESS); 448 } 449 450 static inline PetscErrorCode KSP_PCMatApplyTranspose(KSP ksp, Mat X, Mat Y) 451 { 452 PetscFunctionBegin; 453 if (!ksp->transpose_solve) PetscCall(PCMatApplyTranspose(ksp->pc, X, Y)); 454 else PetscCall(PCMatApply(ksp->pc, X, Y)); 455 PetscFunctionReturn(PETSC_SUCCESS); 456 } 457 458 static inline PetscErrorCode KSP_PCApplyBAorAB(KSP ksp, Vec x, Vec y, Vec w) 459 { 460 PetscFunctionBegin; 461 if (ksp->transpose_solve) { 462 PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w)); 463 PetscCall(KSP_RemoveNullSpaceTranspose(ksp, y)); 464 } else { 465 PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w)); 466 PetscCall(KSP_RemoveNullSpace(ksp, y)); 467 } 468 PetscFunctionReturn(PETSC_SUCCESS); 469 } 470 471 static inline PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp, Vec x, Vec y, Vec w) 472 { 473 PetscFunctionBegin; 474 if (ksp->transpose_solve) PetscCall(PCApplyBAorAB(ksp->pc, ksp->pc_side, x, y, w)); 475 else PetscCall(PCApplyBAorABTranspose(ksp->pc, ksp->pc_side, x, y, w)); 476 PetscFunctionReturn(PETSC_SUCCESS); 477 } 478 479 PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization; 480 PETSC_EXTERN PetscLogEvent KSP_SetUp; 481 PETSC_EXTERN PetscLogEvent KSP_Solve; 482 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0; 483 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1; 484 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2; 485 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3; 486 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4; 487 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S; 488 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L; 489 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U; 490 PETSC_EXTERN PetscLogEvent KSP_SolveTranspose; 491 PETSC_EXTERN PetscLogEvent KSP_MatSolve; 492 PETSC_EXTERN PetscLogEvent KSP_MatSolveTranspose; 493 494 PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat, IS, IS, IS, IS, MatReuse, Mat *, MatSchurComplementAinvType, MatReuse, Mat *); 495 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC, PetscBool *); 496 497 /*MC 498 KSPCheckDot - Checks if the result of a dot product used by the corresponding `KSP` contains infinity or NaN. These indicate that the previous 499 application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`. 500 501 Collective 502 503 Input Parameter: 504 . ksp - the linear solver `KSP` context. 505 506 Output Parameter: 507 . beta - the result of the inner product 508 509 Level: developer 510 511 Developer Notes: 512 Used to manage returning from `KSP` solvers collectively whose preconditioners have failed, possibly only a subset of MPI processes, in some way 513 514 It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products. 515 516 .seealso: `PCFailedReason`, `KSPConvergedReason`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckNorm()`, `KSPCheckSolve()`, 517 `KSPSetErrorIfNotConverged()` 518 M*/ 519 #define KSPCheckDot(ksp, beta) \ 520 do { \ 521 if (PetscIsInfOrNanScalar(beta)) { \ 522 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to infinity or NaN inner product"); \ 523 { \ 524 PCFailedReason pcreason; \ 525 PetscCall(PCReduceFailedReason(ksp->pc)); \ 526 PetscCall(PCGetFailedReason(ksp->pc, &pcreason)); \ 527 PetscCall(VecFlag(ksp->vec_sol, pcreason)); \ 528 if (pcreason) { \ 529 ksp->reason = KSP_DIVERGED_PC_FAILED; \ 530 } else { \ 531 ksp->reason = KSP_DIVERGED_NANORINF; \ 532 } \ 533 PetscFunctionReturn(PETSC_SUCCESS); \ 534 } \ 535 } \ 536 } while (0) 537 538 /*MC 539 KSPCheckNorm - Checks if the result of a norm used by the corresponding `KSP` contains `inf` or `NaN`. These indicate that the previous 540 application of the preconditioner generated an error. Sets a `KSPConvergedReason` and returns if the `PC` set a `PCFailedReason`. 541 542 Collective 543 544 Input Parameter: 545 . ksp - the linear solver `KSP` context. 546 547 Output Parameter: 548 . beta - the result of the norm 549 550 Level: developer 551 552 Developer Notes: 553 Used to manage returning from `KSP` solvers collectively whose preconditioners have failed, possibly only a subset of MPI processes, in some way. 554 555 It uses the fact that `KSP` piggy-backs the collectivity of certain error conditions on the results of norms and inner products. 556 557 .seealso: `PCFailedReason`, `KSPConvergedReason`, `KSP`, `KSPCreate()`, `KSPSetType()`, `KSP`, `KSPCheckDot()`, `KSPCheckSolve()`, 558 `KSPSetErrorIfNotConverged()` 559 M*/ 560 #define KSPCheckNorm(ksp, beta) \ 561 do { \ 562 if (PetscIsInfOrNanReal(beta)) { \ 563 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to infinity or NaN norm"); \ 564 { \ 565 PCFailedReason pcreason; \ 566 PetscCall(PCReduceFailedReason(ksp->pc)); \ 567 PetscCall(PCGetFailedReason(ksp->pc, &pcreason)); \ 568 PetscCall(VecFlag(ksp->vec_sol, pcreason)); \ 569 if (pcreason) { \ 570 ksp->reason = KSP_DIVERGED_PC_FAILED; \ 571 } else { \ 572 ksp->reason = KSP_DIVERGED_NANORINF; \ 573 } \ 574 ksp->rnorm = beta; \ 575 PetscFunctionReturn(PETSC_SUCCESS); \ 576 } \ 577 } \ 578 } while (0) 579 580 PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]); 581 PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP, PetscInt, PetscReal *); 582