Actual source code: kspimpl.h

petsc-master 2021-01-18
Report Typos and Errors

  2: #ifndef _KSPIMPL_H
  3: #define _KSPIMPL_H

  5: #include <petscksp.h>
  6: #include <petscds.h>
  7: #include <petsc/private/petscimpl.h>

  9: PETSC_EXTERN PetscBool KSPRegisterAllCalled;
 10: PETSC_EXTERN PetscErrorCode KSPRegisterAll(void);
 11: PETSC_EXTERN PetscErrorCode KSPGuessRegisterAll(void);
 12: PETSC_EXTERN PetscErrorCode KSPMatRegisterAll(void);

 14: typedef struct _KSPOps *KSPOps;

 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: };

 37: typedef struct _KSPGuessOps *KSPGuessOps;

 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: };

 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: };

 60: PETSC_EXTERN PetscErrorCode KSPGuessCreate_Fischer(KSPGuess);
 61: PETSC_EXTERN PetscErrorCode KSPGuessCreate_POD(KSPGuess);

 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;

 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 */

 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 */

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. */

119:   PetscErrorCode (*converged)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
120:   PetscErrorCode (*convergeddestroy)(void*);
121:   void       *cnvP;

123:   void       *user;             /* optional user-defined context */

125:   PC         pc;

127:   void       *data;                      /* holder for misc stuff associated
128:                                    with a particular iterative solver */

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;

134:   /* ----------------Default work-area management -------------------- */
135:   PetscInt       nwork;
136:   Vec            *work;

138:   KSPSetUpStage  setupstage;
139:   PetscBool      setupnewmatrix; /* true if we need to call ksp->ops->setup with KSP_SETUP_NEWMATRIX */

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 */

144:   PetscBool      transpose_solve;    /* solve transpose system instead */

146:   KSPNormType    normtype;          /* type of norm used for convergence tests */

148:   PCSide         pc_side_set;   /* PC type set explicitly by user */
149:   KSPNormType    normtype_set;  /* Norm type set explicitly by user */

151:   /*   Allow diagonally scaling the matrix before computing the preconditioner or using
152:        the Krylov method. Note this is NOT just Jacobi preconditioning */

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;

160:   PetscInt     setfromoptionscalled;
161:   PetscBool    skippcsetfromoptions; /* if set then KSPSetFromOptions() does not call PCSetFromOptions() */

163:   PetscViewer  eigviewer;   /* Viewer where computed eigenvalues are displayed */

165:   PetscErrorCode (*presolve)(KSP,Vec,Vec,void*);
166:   PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*);
167:   void           *prectx,*postctx;
168: };

170: typedef struct { /* dummy data structure used in KSPMonitorDynamicTolerance() */
171:   PetscReal coef;
172:   PetscReal bnrm;
173: } KSPDynTolCtx;

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;

182: PETSC_STATIC_INLINE PetscErrorCode KSPLogResidualHistory(KSP ksp,PetscReal norm)
183: {

187:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
188:   if (ksp->res_hist && ksp->res_hist_max > ksp->res_hist_len) {
189:     ksp->res_hist[ksp->res_hist_len++] = norm;
190:   }
191:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
192:   return(0);
193: }

195: PETSC_STATIC_INLINE PetscErrorCode KSPLogErrorHistory(KSP ksp)
196: {
197:   DM             dm;

201:   PetscObjectSAWsTakeAccess((PetscObject) ksp);
202:   KSPGetDM(ksp, &dm);
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;

211:     KSPBuildSolution(ksp, NULL, &u);
212:     /* TODO Was needed to correct for Newton solution, but I just need to set a solution */
213:     //VecScale(u, -1.0);
214:     /* TODO Case when I have a solution */
215:     if (0) {
216:       DMGetDS(dm, &ds);
217:       PetscDSGetNumFields(ds, &Nf);
218:       if (Nf > 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle number of fields %D > 1 right now", Nf);
219:       PetscDSGetExactSolution(ds, 0, &exactSol, &exactCtx);
220:       DMComputeL2FieldDiff(dm, 0.0, &exactSol, &exactCtx, u, &error);
221:     } else {
222:       /* The null solution A 0 = 0 */
223:       VecNorm(u, NORM_2, &error);
224:     }
225:     ksp->err_hist[ksp->err_hist_len++] = error;
226:   }
227:   PetscObjectSAWsGrantAccess((PetscObject) ksp);
228:   return(0);
229: }

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: }

240: PETSC_STATIC_INLINE PetscErrorCode KSPSetNoisy_Private(Vec v)
241: {
242:   PetscScalar   *a;
243:   PetscInt       n, istart, i;

246:   VecGetOwnershipRange(v, &istart, NULL);
247:   VecGetLocalSize(v, &n);
248:   VecGetArrayWrite(v, &a);
249:   for (i = 0; i < n; ++i) a[i] = KSPNoisyHash_Private(i+istart);
250:   VecRestoreArrayWrite(v, &a);
251:   return(0);
252: }

254: PETSC_INTERN PetscErrorCode KSPSetUpNorms_Private(KSP,PetscBool,KSPNormType*,PCSide*);

256: PETSC_INTERN PetscErrorCode KSPPlotEigenContours_Private(KSP,PetscInt,const PetscReal*,const PetscReal*);

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: };

268: struct _p_DMKSP {
269:   PETSCHEADER(struct _DMKSPOps);
270:   void *operatorsctx;
271:   void *rhsctx;
272:   void *initialguessctx;
273:   void *data;

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;

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);

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: {
297:   if (ksp->pc_side == PC_LEFT) {
298:     Mat          A;
299:     MatNullSpace nullsp;
300:     PCGetOperators(ksp->pc,&A,NULL);
301:     MatGetNullSpace(A,&nullsp);
302:     if (nullsp) {
303:       MatNullSpaceRemove(nullsp,y);
304:     }
305:   }
306:   return(0);
307: }

309: PETSC_STATIC_INLINE PetscErrorCode KSP_RemoveNullSpaceTranspose(KSP ksp,Vec y)
310: {
313:   if (ksp->pc_side == PC_LEFT) {
314:     Mat          A;
315:     MatNullSpace nullsp;
316:     PCGetOperators(ksp->pc,&A,NULL);
317:     MatGetTransposeNullSpace(A,&nullsp);
318:     if (nullsp) {
319:       MatNullSpaceRemove(nullsp,y);
320:     }
321:   }
322:   return(0);
323: }

325: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMult(KSP ksp,Mat A,Vec x,Vec y)
326: {
329:   if (!ksp->transpose_solve) {MatMult(A,x,y);}
330:   else                       {MatMultTranspose(A,x,y);}
331:   return(0);
332: }

334: PETSC_STATIC_INLINE PetscErrorCode KSP_MatMultTranspose(KSP ksp,Mat A,Vec x,Vec y)
335: {
338:   if (!ksp->transpose_solve) {MatMultTranspose(A,x,y);}
339:   else                       {MatMult(A,x,y);}
340:   return(0);
341: }

343: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApply(KSP ksp,Vec x,Vec y)
344: {
347:   if (!ksp->transpose_solve) {
348:     PCApply(ksp->pc,x,y);
349:     KSP_RemoveNullSpace(ksp,y);
350:   } else {
351:     PCApplyTranspose(ksp->pc,x,y);
352:     KSP_RemoveNullSpaceTranspose(ksp,y);
353:   }
354:   return(0);
355: }

357: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyTranspose(KSP ksp,Vec x,Vec y)
358: {
361:   if (!ksp->transpose_solve) {
362:     PCApplyTranspose(ksp->pc,x,y);
363:     KSP_RemoveNullSpaceTranspose(ksp,y);
364:   } else {
365:     PCApply(ksp->pc,x,y);
366:     KSP_RemoveNullSpace(ksp,y);
367:   }
368:   return(0);
369: }

371: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorAB(KSP ksp,Vec x,Vec y,Vec w)
372: {
375:   if (!ksp->transpose_solve) {
376:     PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
377:     KSP_RemoveNullSpace(ksp,y);
378:   } else {
379:     PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
380:     KSP_RemoveNullSpaceTranspose(ksp,y);
381:   }
382:   return(0);
383: }

385: PETSC_STATIC_INLINE PetscErrorCode KSP_PCApplyBAorABTranspose(KSP ksp,Vec x,Vec y,Vec w)
386: {
389:   if (!ksp->transpose_solve) {
390:     PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w);
391:   } else {
392:     PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w);
393:   }
394:   return(0);
395: }

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;

411: PETSC_INTERN PetscErrorCode MatGetSchurComplement_Basic(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatSchurComplementAinvType,MatReuse,Mat*);
412: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC,PetscBool*);

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

418:    Collective on ksp

420:    Input Parameter:
421: .  ksp - the linear solver (KSP) context.

423:    Output Parameter:
424: .  beta - the result of the inner product

426:    Level: developer

428:    Developer Note:
429:    this is used to manage returning from KSP solvers whose preconditioners have failed in some way

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 {\
438:       PCFailedReason pcreason;\
439:       PetscInt       sendbuf,recvbuf; \
440:       PCGetFailedReasonRank(ksp->pc,&pcreason);\
441:       sendbuf = (PetscInt)pcreason; \
442:       MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));CHKERRMPI(ierr);\
443:       if (recvbuf) {                                                           \
444:         PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf); \
445:         ksp->reason = KSP_DIVERGED_PC_FAILED;\
446:         VecSetInf(ksp->vec_sol);\
447:       } else {\
448:         ksp->reason = KSP_DIVERGED_NANORINF;\
449:       }\
450:       return(0);\
451:     }\
452:   } } while (0)

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

458:    Collective on ksp

460:    Input Parameter:
461: .  ksp - the linear solver (KSP) context.

463:    Output Parameter:
464: .  beta - the result of the norm

466:    Level: developer

468:    Developer Note:
469:    this is used to manage returning from KSP solvers whose preconditioners have failed in some way

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 {\
478:       PCFailedReason pcreason;\
479:       PetscInt       sendbuf,recvbuf; \
480:       PCGetFailedReasonRank(ksp->pc,&pcreason);\
481:       sendbuf = (PetscInt)pcreason; \
482:       MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));CHKERRMPI(ierr);\
483:       if (recvbuf) {                                                           \
484:         PCSetFailedReason(ksp->pc,(PCFailedReason)recvbuf); \
485:         ksp->reason = KSP_DIVERGED_PC_FAILED;                         \
486:         VecSetInf(ksp->vec_sol);\
487:       } else {\
488:         PCSetFailedReason(ksp->pc,PC_NOERROR); \
489:         ksp->reason = KSP_DIVERGED_NANORINF;\
490:       }\
491:       return(0);\
492:     }\
493:   } } while (0)

495: #endif