Actual source code: kspimpl.h
petsc-master 2021-01-18
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