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