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