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