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