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