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