1 2 #ifndef _KSPIMPL_H 3 #define _KSPIMPL_H 4 5 #include <petscksp.h> 6 #include <petsc/private/petscimpl.h> 7 8 PETSC_EXTERN PetscBool KSPRegisterAllCalled; 9 PETSC_EXTERN PetscErrorCode KSPRegisterAll(void); 10 PETSC_EXTERN PetscErrorCode KSPGuessRegisterAll(void); 11 PETSC_EXTERN PetscErrorCode KSPMatRegisterAll(void); 12 13 typedef struct _KSPOps *KSPOps; 14 15 struct _KSPOps { 16 PetscErrorCode (*buildsolution)(KSP,Vec,Vec*); /* Returns a pointer to the solution, or 17 calculates the solution in a 18 user-provided area. */ 19 PetscErrorCode (*buildresidual)(KSP,Vec,Vec,Vec*); /* Returns a pointer to the residual, or 20 calculates the residual in a 21 user-provided area. */ 22 PetscErrorCode (*solve)(KSP); /* actual solver */ 23 PetscErrorCode (*setup)(KSP); 24 PetscErrorCode (*setfromoptions)(PetscOptionItems*,KSP); 25 PetscErrorCode (*publishoptions)(KSP); 26 PetscErrorCode (*computeextremesingularvalues)(KSP,PetscReal*,PetscReal*); 27 PetscErrorCode (*computeeigenvalues)(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt *); 28 PetscErrorCode (*computeritz)(KSP,PetscBool,PetscBool,PetscInt*,Vec[],PetscReal*,PetscReal*); 29 PetscErrorCode (*destroy)(KSP); 30 PetscErrorCode (*view)(KSP,PetscViewer); 31 PetscErrorCode (*reset)(KSP); 32 PetscErrorCode (*load)(KSP,PetscViewer); 33 }; 34 35 typedef struct _KSPGuessOps *KSPGuessOps; 36 37 struct _KSPGuessOps { 38 PetscErrorCode (*formguess)(KSPGuess,Vec,Vec); /* Form initial guess */ 39 PetscErrorCode (*update)(KSPGuess,Vec,Vec); /* Update database */ 40 PetscErrorCode (*setfromoptions)(KSPGuess); 41 PetscErrorCode (*setup)(KSPGuess); 42 PetscErrorCode (*destroy)(KSPGuess); 43 PetscErrorCode (*view)(KSPGuess,PetscViewer); 44 PetscErrorCode (*reset)(KSPGuess); 45 }; 46 47 /* 48 Defines the KSPGuess data structure. 49 */ 50 struct _p_KSPGuess { 51 PETSCHEADER(struct _KSPGuessOps); 52 KSP ksp; /* the parent KSP */ 53 Mat A; /* the current linear operator */ 54 PetscObjectState omatstate; /* previous linear operator state */ 55 void *data; /* pointer to the specific implementation */ 56 }; 57 58 PETSC_EXTERN PetscErrorCode KSPGuessCreate_Fischer(KSPGuess); 59 PETSC_EXTERN PetscErrorCode KSPGuessCreate_POD(KSPGuess); 60 61 /* 62 Maximum number of monitors you can run with a single KSP 63 */ 64 #define MAXKSPMONITORS 5 65 typedef enum {KSP_SETUP_NEW, KSP_SETUP_NEWMATRIX, KSP_SETUP_NEWRHS} KSPSetUpStage; 66 67 /* 68 Defines the KSP data structure. 69 */ 70 struct _p_KSP { 71 PETSCHEADER(struct _KSPOps); 72 DM dm; 73 PetscBool dmAuto; /* DM was created automatically by KSP */ 74 PetscBool dmActive; /* KSP should use DM for computing operators */ 75 /*------------------------- User parameters--------------------------*/ 76 PetscInt max_it; /* maximum number of iterations */ 77 KSPGuess guess; 78 PetscBool guess_zero, /* flag for whether initial guess is 0 */ 79 calc_sings, /* calculate extreme Singular Values */ 80 calc_ritz, /* calculate (harmonic) Ritz pairs */ 81 guess_knoll; /* use initial guess of PCApply(ksp->B,b */ 82 PCSide pc_side; /* flag for left, right, or symmetric preconditioning */ 83 PetscInt normsupporttable[KSP_NORM_MAX][PC_SIDE_MAX]; /* Table of supported norms and pc_side, see KSPSetSupportedNorm() */ 84 PetscReal rtol, /* relative tolerance */ 85 abstol, /* absolute tolerance */ 86 ttol, /* (not set by user) */ 87 divtol; /* divergence tolerance */ 88 PetscReal rnorm0; /* initial residual norm (used for divergence testing) */ 89 PetscReal rnorm; /* current residual norm */ 90 KSPConvergedReason reason; 91 PetscBool errorifnotconverged; /* create an error if the KSPSolve() does not converge */ 92 93 Vec vec_sol,vec_rhs; /* pointer to where user has stashed 94 the solution and rhs, these are 95 never touched by the code, only 96 passed back to the user */ 97 PetscReal *res_hist; /* If !0 stores residual at iterations*/ 98 PetscReal *res_hist_alloc; /* If !0 means user did not provide buffer, needs deallocation */ 99 PetscInt res_hist_len; /* current size of residual history array */ 100 PetscInt res_hist_max; /* actual amount of data in residual_history */ 101 PetscBool res_hist_reset; /* reset history to size zero for each new solve */ 102 103 PetscInt chknorm; /* only compute/check norm if iterations is great than this */ 104 PetscBool lagnorm; /* Lag the residual norm calculation so that it is computed as part of the 105 MPI_Allreduce() for computing the inner products for the next iteration. */ 106 /* --------User (or default) routines (most return -1 on error) --------*/ 107 PetscErrorCode (*monitor[MAXKSPMONITORS])(KSP,PetscInt,PetscReal,void*); /* returns control to user after */ 108 PetscErrorCode (*monitordestroy[MAXKSPMONITORS])(void**); /* */ 109 void *monitorcontext[MAXKSPMONITORS]; /* residual calculation, allows user */ 110 PetscInt numbermonitors; /* to, for instance, print residual norm, etc. */ 111 112 PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*); 113 PetscErrorCode (*convergeddestroy)(void*); 114 void *cnvP; 115 116 void *user; /* optional user-defined context */ 117 118 PC pc; 119 120 void *data; /* holder for misc stuff associated 121 with a particular iterative solver */ 122 123 PetscBool view, viewPre, viewReason, viewMat, viewPMat, viewRhs, viewSol, viewMatExp, viewEV, viewSV, viewEVExp, viewFinalRes, viewPOpExp; 124 PetscViewer viewer, viewerPre, viewerReason, viewerMat, viewerPMat, viewerRhs, viewerSol, viewerMatExp, viewerEV, viewerSV, viewerEVExp, viewerFinalRes, viewerPOpExp; 125 PetscViewerFormat format, formatPre, formatReason, formatMat, formatPMat, formatRhs, formatSol, formatMatExp, formatEV, formatSV, formatEVExp, formatFinalRes, formatPOpExp; 126 127 /* ----------------Default work-area management -------------------- */ 128 PetscInt nwork; 129 Vec *work; 130 131 KSPSetUpStage setupstage; 132 PetscBool setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */ 133 134 PetscInt its; /* number of iterations so far computed in THIS linear solve*/ 135 PetscInt totalits; /* number of iterations used by this KSP object since it was created */ 136 137 PetscBool transpose_solve; /* solve transpose system instead */ 138 139 KSPNormType normtype; /* type of norm used for convergence tests */ 140 141 PCSide pc_side_set; /* PC type set explicitly by user */ 142 KSPNormType normtype_set; /* Norm type set explicitly by user */ 143 144 /* Allow diagonally scaling the matrix before computing the preconditioner or using 145 the Krylov method. Note this is NOT just Jacobi preconditioning */ 146 147 PetscBool dscale; /* diagonal scale system; used with KSPSetDiagonalScale() */ 148 PetscBool dscalefix; /* unscale system after solve */ 149 PetscBool dscalefix2; /* system has been unscaled */ 150 Vec diagonal; /* 1/sqrt(diag of matrix) */ 151 Vec truediagonal; 152 153 PetscInt setfromoptionscalled; 154 PetscBool skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */ 155 156 PetscViewer eigviewer; /* Viewer where computed eigenvalues are displayed */ 157 158 PetscErrorCode (*presolve)(KSP,Vec,Vec,void*); 159 PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*); 160 void *prectx,*postctx; 161 }; 162 163 typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */ 164 PetscReal coef; 165 PetscReal bnrm; 166 } KSPDynTolCtx; 167 168 typedef struct { 169 PetscBool initialrtol; /* default relative residual decrease is computing from initial residual, not rhs */ 170 PetscBool mininitialrtol; /* default relative residual decrease is computing from min of initial residual and rhs */ 171 Vec work; 172 } KSPConvergedDefaultCtx; 173 174 PETSC_STATIC_INLINE PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm) 175 { 176 PetscErrorCode ierr; 177 178 PetscFunctionBegin; 179 ierr = PetscObjectSAWsTakeAccess((PetscObject)ksp);CHKERRQ(ierr); 180 if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) { 181 ksp->res_hist[ksp->res_hist_len++] = norm; 182 } 183 ierr = PetscObjectSAWsGrantAccess((PetscObject)ksp);CHKERRQ(ierr); 184 PetscFunctionReturn(0); 185 } 186 187 PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP,PetscBool,KSPNormType*,PCSide*); 188 189 PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP,PetscInt,const PetscReal*,const PetscReal*); 190 191 typedef struct _p_DMKSP *DMKSP; 192 typedef struct _DMKSPOps *DMKSPOps; 193 struct _DMKSPOps { 194 PetscErrorCode (*computeoperators)(KSP,Mat,Mat,void*); 195 PetscErrorCode (*computerhs)(KSP,Vec,void*); 196 PetscErrorCode (*computeinitialguess)(KSP,Vec,void*); 197 PetscErrorCode (*destroy)(DMKSP*); 198 PetscErrorCode (*duplicate)(DMKSP,DMKSP); 199 }; 200 201 struct _p_DMKSP { 202 PETSCHEADER(struct _DMKSPOps); 203 void *operatorsctx; 204 void *rhsctx; 205 void *initialguessctx; 206 void *data; 207 208 /* This is NOT reference counted. The DM on which this context was first created is cached here to implement one-way 209 * copy-on-write. When DMGetDMKSPWrite() sees a request using a different DM, it makes a copy. Thus, if a user 210 * only interacts directly with one level, e.g., using KSPSetComputeOperators(), then coarse levels are constructed by 211 * PCMG, then the user changes the routine with another call to KSPSetComputeOperators(), it automatically propagates 212 * to all the levels. If instead, they get out a specific level and set the routines on that level, subsequent changes 213 * to the original level will no longer propagate to that level. 214 */ 215 DM originaldm; 216 217 void (*fortran_func_pointers[3])(void); /* Store our own function pointers so they are associated with the DMKSP instead of the DM */ 218 }; 219 PETSC_EXTERN PetscErrorCode DMGetDMKSP(DM,DMKSP*); 220 PETSC_EXTERN PetscErrorCode DMGetDMKSPWrite(DM,DMKSP*); 221 PETSC_EXTERN PetscErrorCode DMCopyDMKSP(DM,DM); 222 223 /* 224 These allow the various Krylov methods to apply to either the linear system or its transpose. 225 */ 226 PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpace(KSP ksp,Vec y) 227 { 228 PetscErrorCode ierr; 229 PetscFunctionBegin; 230 if (ksp->pc_side == PC_LEFT) { 231 Mat A; 232 MatNullSpace nullsp; 233 ierr = PCGetOperators(ksp->pc,&A,NULL);CHKERRQ(ierr); 234 ierr = MatGetNullSpace(A,&nullsp);CHKERRQ(ierr); 235 if (nullsp) { 236 ierr = MatNullSpaceRemove(nullsp,y);CHKERRQ(ierr); 237 } 238 } 239 PetscFunctionReturn(0); 240 } 241 242 PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp,Vec y) 243 { 244 PetscErrorCode ierr; 245 PetscFunctionBegin; 246 if (ksp->pc_side == PC_LEFT) { 247 Mat A; 248 MatNullSpace nullsp; 249 ierr = PCGetOperators(ksp->pc,&A,NULL);CHKERRQ(ierr); 250 ierr = MatGetTransposeNullSpace(A,&nullsp);CHKERRQ(ierr); 251 if (nullsp) { 252 ierr = MatNullSpaceRemove(nullsp,y);CHKERRQ(ierr); 253 } 254 } 255 PetscFunctionReturn(0); 256 } 257 258 PETSC_STATIC_INLINE PetscErrorCode KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y) 259 { 260 PetscErrorCode ierr; 261 PetscFunctionBegin; 262 if (!ksp->transpose_solve) {ierr = MatMult(A,x,y);CHKERRQ(ierr);} 263 else {ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);} 264 PetscFunctionReturn(0); 265 } 266 267 PETSC_STATIC_INLINE PetscErrorCode KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y) 268 { 269 PetscErrorCode ierr; 270 PetscFunctionBegin; 271 if (!ksp->transpose_solve) {ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr);} 272 else {ierr = MatMult(A,x,y);CHKERRQ(ierr);} 273 PetscFunctionReturn(0); 274 } 275 276 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y) 277 { 278 PetscErrorCode ierr; 279 PetscFunctionBegin; 280 if (!ksp->transpose_solve) { 281 ierr = PCApply(ksp->pc,x,y);CHKERRQ(ierr); 282 ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr); 283 } else { 284 ierr = PCApplyTranspose(ksp->pc,x,y);CHKERRQ(ierr); 285 ierr = KSP_RemoveNullSpaceTranspose(ksp,y);CHKERRQ(ierr); 286 } 287 PetscFunctionReturn(0); 288 } 289 290 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y) 291 { 292 PetscErrorCode ierr; 293 PetscFunctionBegin; 294 if (!ksp->transpose_solve) { 295 ierr = PCApplyTranspose(ksp->pc,x,y);CHKERRQ(ierr); 296 ierr = KSP_RemoveNullSpaceTranspose(ksp,y);CHKERRQ(ierr); 297 } else { 298 ierr = PCApply(ksp->pc,x,y);CHKERRQ(ierr); 299 ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr); 300 } 301 PetscFunctionReturn(0); 302 } 303 304 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w) 305 { 306 PetscErrorCode ierr; 307 PetscFunctionBegin; 308 if (!ksp->transpose_solve) { 309 ierr = PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 310 ierr = KSP_RemoveNullSpace(ksp,y);CHKERRQ(ierr); 311 } else { 312 ierr = PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 313 ierr = KSP_RemoveNullSpaceTranspose(ksp,y);CHKERRQ(ierr); 314 } 315 PetscFunctionReturn(0); 316 } 317 318 PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w) 319 { 320 PetscErrorCode ierr; 321 PetscFunctionBegin; 322 if (!ksp->transpose_solve) { 323 ierr = PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 324 } else { 325 ierr = PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);CHKERRQ(ierr); 326 } 327 PetscFunctionReturn(0); 328 } 329 330 PETSC_EXTERN PetscLogEvent KSP_GMRESOrthogonalization; 331 PETSC_EXTERN PetscLogEvent KSP_SetUp; 332 PETSC_EXTERN PetscLogEvent KSP_Solve; 333 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_0; 334 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_1; 335 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_2; 336 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_3; 337 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_4; 338 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_S; 339 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_L; 340 PETSC_EXTERN PetscLogEvent KSP_Solve_FS_U; 341 342 PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*); 343 PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*); 344 345 /* 346 Either generate an error or mark as diverged when a scalar from an inner product is Nan or Inf 347 */ 348 #define KSPCheckDot(ksp,beta) \ 349 if (PetscIsInfOrNanScalar(beta)) { \ 350 if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf inner product");\ 351 else {\ 352 PetscErrorCode ierr;\ 353 PCFailedReason pcreason;\ 354 PetscInt sendbuf,pcreason_max; \ 355 ierr = PCGetSetUpFailedReason(ksp->pc,&pcreason);CHKERRQ(ierr);\ 356 sendbuf = (PetscInt)pcreason; \ 357 ierr = MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr); \ 358 if (pcreason_max) {\ 359 ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;\ 360 ierr = VecSetInf(ksp->vec_sol);CHKERRQ(ierr);\ 361 } else {\ 362 ksp->reason = KSP_DIVERGED_NANORINF;\ 363 }\ 364 PetscFunctionReturn(0);\ 365 }\ 366 } 367 368 /* 369 Either generate an error or mark as diverged when a real from a norm is Nan or Inf 370 */ 371 #define KSPCheckNorm(ksp,beta) \ 372 if (PetscIsInfOrNanReal(beta)) { \ 373 if (ksp->errorifnotconverged) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged due to Nan or Inf norm");\ 374 else {\ 375 PetscErrorCode ierr;\ 376 PCFailedReason pcreason;\ 377 PetscInt sendbuf,pcreason_max; \ 378 ierr = PCGetSetUpFailedReason(ksp->pc,&pcreason);CHKERRQ(ierr);\ 379 sendbuf = (PetscInt)pcreason; \ 380 ierr = MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));CHKERRQ(ierr); \ 381 if (pcreason_max) {\ 382 ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;\ 383 ierr = VecSetInf(ksp->vec_sol);CHKERRQ(ierr);\ 384 } else {\ 385 ksp->reason = KSP_DIVERGED_NANORINF;\ 386 }\ 387 PetscFunctionReturn(0);\ 388 }\ 389 } 390 391 #endif 392