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