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_MatMultHermitianTranspose(KSP ksp,Mat A,Vec x,Vec y) 361 { 362 PetscErrorCode ierr; 363 PetscFunctionBegin; 364 if (!ksp->transpose_solve) {ierr = MatMultHermitianTranspose(A,x,y);CHKERRQ(ierr);} 365 else { 366 Vec w; 367 ierr = VecDuplicate(x,&w);CHKERRQ(ierr); 368 ierr = VecCopy(x,w);CHKERRQ(ierr); 369 ierr = VecConjugate(w);CHKERRQ(ierr); 370 ierr = MatMult(A,w,y);CHKERRQ(ierr); 371 ierr = VecDestroy(&w);CHKERRQ(ierr); 372 ierr = VecConjugate(y);CHKERRQ(ierr); 373 } 374 PetscFunctionReturn(0); 375 } 376 377 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y) 378 { 379 PetscErrorCode ierr; 380 PetscFunctionBegin; 381 if (!ksp->transpose_solve) { 382 ierr = PCApply(ksp->pc,x,y);CHKERRQ(ierr); 383 ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr); 384 } else { 385 ierr = PCApplyTranspose(ksp->pc,x,y);CHKERRQ(ierr); 386 ierr = KSP_RemoveNullSpaceTranspose(ksp,y);CHKERRQ(ierr); 387 } 388 PetscFunctionReturn(0); 389 } 390 391 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y) 392 { 393 PetscErrorCode ierr; 394 PetscFunctionBegin; 395 if (!ksp->transpose_solve) { 396 ierr = PCApplyTranspose(ksp->pc,x,y);CHKERRQ(ierr); 397 ierr = KSP_RemoveNullSpaceTranspose(ksp,y);CHKERRQ(ierr); 398 } else { 399 ierr = PCApply(ksp->pc,x,y);CHKERRQ(ierr); 400 ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr); 401 } 402 PetscFunctionReturn(0); 403 } 404 405 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyHermitianTranspose(KSP ksp,Vec x,Vec y) 406 { 407 PetscErrorCode ierr; 408 PetscFunctionBegin; 409 ierr = VecConjugate(x);CHKERRQ(ierr); 410 ierr = KSP_PCApplyTranspose(ksp,x,y);CHKERRQ(ierr); 411 ierr = VecConjugate(x);CHKERRQ(ierr); 412 ierr = VecConjugate(y);CHKERRQ(ierr); 413 PetscFunctionReturn(0); 414 } 415 416 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w) 417 { 418 PetscErrorCode ierr; 419 PetscFunctionBegin; 420 if (!ksp->transpose_solve) { 421 ierr = PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 422 ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr); 423 } else { 424 ierr = PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 425 ierr = KSP_RemoveNullSpaceTranspose(ksp,y);CHKERRQ(ierr); 426 } 427 PetscFunctionReturn(0); 428 } 429 430 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w) 431 { 432 PetscErrorCode ierr; 433 PetscFunctionBegin; 434 if (!ksp->transpose_solve) { 435 ierr = PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 436 } else { 437 ierr = PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 438 } 439 PetscFunctionReturn(0); 440 } 441 442 PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization; 443 PETSC_EXTERN PetscLogEvent KSP_SetUp; 444 PETSC_EXTERN PetscLogEvent KSP_Solve; 445 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0; 446 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1; 447 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2; 448 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3; 449 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4; 450 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S; 451 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L; 452 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U; 453 PETSC_EXTERN PetscLogEvent KSP_SolveTranspose; 454 PETSC_EXTERN PetscLogEvent KSP_MatSolve; 455 456 PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*); 457 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*); 458 459 /*MC 460 KSPCheckDot - Checks if the result of a dot product used by the corresponding KSP contains Inf or NaN. These indicate that the previous 461 application of the preconditioner generated an error 462 463 Collective on ksp 464 465 Input Parameter: 466 . ksp - the linear solver (KSP) context. 467 468 Output Parameter: 469 . beta - the result of the inner product 470 471 Level: developer 472 473 Developer Note: 474 this is used to manage returning from KSP solvers whose preconditioners have failed in some way 475 476 .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckNorm(), KSPCheckSolve() 477 M*/ 478 #define KSPCheckDot(ksp,beta) do { \ 479 if (PetscIsInfOrNanScalar(beta)) { \ 480 if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");\ 481 else {\ 482 PetscErrorCode ierr;\ 483 PCFailedReason pcreason;\ 484 PetscInt sendbuf,recvbuf; \ 485 ierr = PCGetFailedReasonRank(ksp->pc,&pcreason);CHKERRQ(ierr);\ 486 sendbuf = (PetscInt)pcreason; \ 487 ierr = MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));CHKERRMPI(ierr);\ 488 if (recvbuf) { \ 489 ierr = PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf);CHKERRQ(ierr); \ 490 ksp->reason = KSP_DIVERGED_PC_FAILED;\ 491 ierr = VecSetInf(ksp->vec_sol);CHKERRQ(ierr);\ 492 } else {\ 493 ksp->reason = KSP_DIVERGED_NANORINF;\ 494 }\ 495 PetscFunctionReturn(0);\ 496 }\ 497 } } while (0) 498 499 /*MC 500 KSPCheckNorm - Checks if the result of a norm used by the corresponding KSP contains Inf or NaN. These indicate that the previous 501 application of the preconditioner generated an error 502 503 Collective on ksp 504 505 Input Parameter: 506 . ksp - the linear solver (KSP) context. 507 508 Output Parameter: 509 . beta - the result of the norm 510 511 Level: developer 512 513 Developer Note: 514 this is used to manage returning from KSP solvers whose preconditioners have failed in some way 515 516 .seealso: KSPCreate(), KSPSetType(), KSP, KSPCheckDot(), KSPCheckSolve() 517 M*/ 518 #define KSPCheckNorm(ksp,beta) do { \ 519 if (PetscIsInfOrNanReal(beta)) { \ 520 if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf norm");\ 521 else {\ 522 PetscErrorCode ierr;\ 523 PCFailedReason pcreason;\ 524 PetscInt sendbuf,recvbuf; \ 525 ierr = PCGetFailedReasonRank(ksp->pc,&pcreason);CHKERRQ(ierr);\ 526 sendbuf = (PetscInt)pcreason; \ 527 ierr = MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));CHKERRMPI(ierr);\ 528 if (recvbuf) { \ 529 ierr = PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf);CHKERRQ(ierr); \ 530 ksp->reason = KSP_DIVERGED_PC_FAILED; \ 531 ierr = VecSetInf(ksp->vec_sol);CHKERRQ(ierr);\ 532 } else {\ 533 ierr = PCSetFailedReason(ksp->pc,PC_NOERROR);CHKERRQ(ierr); \ 534 ksp->reason = KSP_DIVERGED_NANORINF;\ 535 }\ 536 PetscFunctionReturn(0);\ 537 }\ 538 } } while (0) 539 540 #endif 541 542 PETSC_INTERN PetscErrorCode KSPMonitorMakeKey_Internal(const char[], PetscViewerType, PetscViewerFormat, char[]); 543 PETSC_INTERN PetscErrorCode KSPMonitorRange_Private(KSP,PetscInt,PetscReal*); 544