Actual source code: pipefgmres.c
2: #include <../src/ksp/ksp/impls/gmres/pipefgmres/pipefgmresimpl.h>
4: static PetscBool cited = PETSC_FALSE;
5: static const char citation[] = "@article{SSM2016,\n"
6: " author = {P. Sanan and S.M. Schnepp and D.A. May},\n"
7: " title = {Pipelined, Flexible Krylov Subspace Methods},\n"
8: " journal = {SIAM Journal on Scientific Computing},\n"
9: " volume = {38},\n"
10: " number = {5},\n"
11: " pages = {C441-C470},\n"
12: " year = {2016},\n"
13: " doi = {10.1137/15M1049130},\n"
14: " URL = {http://dx.doi.org/10.1137/15M1049130},\n"
15: " eprint = {http://dx.doi.org/10.1137/15M1049130}\n"
16: "}\n";
18: #define PIPEFGMRES_DELTA_DIRECTIONS 10
19: #define PIPEFGMRES_DEFAULT_MAXK 30
21: static PetscErrorCode KSPPIPEFGMRESGetNewVectors(KSP, PetscInt);
22: static PetscErrorCode KSPPIPEFGMRESUpdateHessenberg(KSP, PetscInt, PetscBool *, PetscReal *);
23: static PetscErrorCode KSPPIPEFGMRESBuildSoln(PetscScalar *, Vec, Vec, KSP, PetscInt);
24: extern PetscErrorCode KSPReset_PIPEFGMRES(KSP);
26: /*
28: KSPSetUp_PIPEFGMRES - Sets up the workspace needed by pipefgmres.
30: This is called once, usually automatically by KSPSolve() or KSPSetUp(),
31: but can be called directly by KSPSetUp().
33: */
34: static PetscErrorCode KSPSetUp_PIPEFGMRES(KSP ksp)
35: {
36: PetscInt k;
37: KSP_PIPEFGMRES *pipefgmres = (KSP_PIPEFGMRES *)ksp->data;
38: const PetscInt max_k = pipefgmres->max_k;
40: KSPSetUp_GMRES(ksp);
42: PetscMalloc1((VEC_OFFSET + max_k), &pipefgmres->prevecs);
43: PetscMalloc1((VEC_OFFSET + max_k), &pipefgmres->prevecs_user_work);
45: KSPCreateVecs(ksp, pipefgmres->vv_allocated, &pipefgmres->prevecs_user_work[0], 0, NULL);
46: for (k = 0; k < pipefgmres->vv_allocated; k++) pipefgmres->prevecs[k] = pipefgmres->prevecs_user_work[0][k];
48: PetscMalloc1((VEC_OFFSET + max_k), &pipefgmres->zvecs);
49: PetscMalloc1((VEC_OFFSET + max_k), &pipefgmres->zvecs_user_work);
51: PetscMalloc1((VEC_OFFSET + max_k), &pipefgmres->redux);
53: KSPCreateVecs(ksp, pipefgmres->vv_allocated, &pipefgmres->zvecs_user_work[0], 0, NULL);
54: for (k = 0; k < pipefgmres->vv_allocated; k++) pipefgmres->zvecs[k] = pipefgmres->zvecs_user_work[0][k];
56: return 0;
57: }
59: /*
61: KSPPIPEFGMRESCycle - Run pipefgmres, possibly with restart. Return residual
62: history if requested.
64: input parameters:
65: . pipefgmres - structure containing parameters and work areas
67: output parameters:
68: . itcount - number of iterations used. If null, ignored.
69: . converged - 0 if not converged
71: Notes:
72: On entry, the value in vector VEC_VV(0) should be
73: the initial residual.
75: */
76: static PetscErrorCode KSPPIPEFGMRESCycle(PetscInt *itcount, KSP ksp)
77: {
78: KSP_PIPEFGMRES *pipefgmres = (KSP_PIPEFGMRES *)(ksp->data);
79: PetscReal res_norm;
80: PetscReal hapbnd, tt;
81: PetscScalar *hh, *hes, *lhh, shift = pipefgmres->shift;
82: PetscBool hapend = PETSC_FALSE; /* indicates happy breakdown ending */
83: PetscInt loc_it; /* local count of # of dir. in Krylov space */
84: PetscInt max_k = pipefgmres->max_k; /* max # of directions Krylov space */
85: PetscInt i, j, k;
86: Mat Amat, Pmat;
87: Vec Q, W; /* Pipelining vectors */
88: Vec *redux = pipefgmres->redux; /* workspace for single reduction */
90: if (itcount) *itcount = 0;
92: /* Assign simpler names to these vectors, allocated as pipelining workspace */
93: Q = VEC_Q;
94: W = VEC_W;
96: /* Allocate memory for orthogonalization work (freed in the GMRES Destroy routine)*/
97: /* Note that we add an extra value here to allow for a single reduction */
98: if (!pipefgmres->orthogwork) PetscMalloc1(pipefgmres->max_k + 2, &pipefgmres->orthogwork);
99: lhh = pipefgmres->orthogwork;
101: /* Number of pseudo iterations since last restart is the number
102: of prestart directions */
103: loc_it = 0;
105: /* note: (pipefgmres->it) is always set one less than (loc_it) It is used in
106: KSPBUILDSolution_PIPEFGMRES, where it is passed to KSPPIPEFGMRESBuildSoln.
107: Note that when KSPPIPEFGMRESBuildSoln is called from this function,
108: (loc_it -1) is passed, so the two are equivalent */
109: pipefgmres->it = (loc_it - 1);
111: /* initial residual is in VEC_VV(0) - compute its norm*/
112: VecNorm(VEC_VV(0), NORM_2, &res_norm);
114: /* first entry in right-hand-side of hessenberg system is just
115: the initial residual norm */
116: *RS(0) = res_norm;
118: PetscObjectSAWsTakeAccess((PetscObject)ksp);
119: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res_norm;
120: else ksp->rnorm = 0;
121: PetscObjectSAWsGrantAccess((PetscObject)ksp);
122: KSPLogResidualHistory(ksp, ksp->rnorm);
123: KSPMonitor(ksp, ksp->its, ksp->rnorm);
125: /* check for the convergence - maybe the current guess is good enough */
126: (*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP);
127: if (ksp->reason) {
128: if (itcount) *itcount = 0;
129: return 0;
130: }
132: /* scale VEC_VV (the initial residual) */
133: VecScale(VEC_VV(0), 1.0 / res_norm);
135: /* Fill the pipeline */
136: KSP_PCApply(ksp, VEC_VV(loc_it), PREVEC(loc_it));
137: PCGetOperators(ksp->pc, &Amat, &Pmat);
138: KSP_MatMult(ksp, Amat, PREVEC(loc_it), ZVEC(loc_it));
139: VecAXPY(ZVEC(loc_it), -shift, VEC_VV(loc_it)); /* Note shift */
141: /* MAIN ITERATION LOOP BEGINNING*/
142: /* keep iterating until we have converged OR generated the max number
143: of directions OR reached the max number of iterations for the method */
144: while (!ksp->reason && loc_it < max_k && ksp->its < ksp->max_it) {
145: if (loc_it) {
146: KSPLogResidualHistory(ksp, res_norm);
147: KSPMonitor(ksp, ksp->its, res_norm);
148: }
149: pipefgmres->it = (loc_it - 1);
151: /* see if more space is needed for work vectors */
152: if (pipefgmres->vv_allocated <= loc_it + VEC_OFFSET + 1) {
153: KSPPIPEFGMRESGetNewVectors(ksp, loc_it + 1);
154: /* (loc_it+1) is passed in as number of the first vector that should
155: be allocated */
156: }
158: /* Note that these inner products are with "Z" now, so
159: in particular, lhh[loc_it] is the 'barred' or 'shifted' value,
160: not the value from the equivalent FGMRES run (even in exact arithmetic)
161: That is, the H we need for the Arnoldi relation is different from the
162: coefficients we use in the orthogonalization process,because of the shift */
164: /* Do some local twiddling to allow for a single reduction */
165: for (i = 0; i < loc_it + 1; i++) redux[i] = VEC_VV(i);
166: redux[loc_it + 1] = ZVEC(loc_it);
168: /* note the extra dot product which ends up in lh[loc_it+1], which computes ||z||^2 */
169: VecMDotBegin(ZVEC(loc_it), loc_it + 2, redux, lhh);
171: /* Start the split reduction (This actually calls the MPI_Iallreduce, otherwise, the reduction is simply delayed until the "end" call)*/
172: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)ZVEC(loc_it)));
174: /* The work to be overlapped with the inner products follows.
175: This is application of the preconditioner and the operator
176: to compute intermediate quantites which will be combined (locally)
177: with the results of the inner products.
178: */
179: KSP_PCApply(ksp, ZVEC(loc_it), Q);
180: PCGetOperators(ksp->pc, &Amat, &Pmat);
181: KSP_MatMult(ksp, Amat, Q, W);
183: /* Compute inner products of the new direction with previous directions,
184: and the norm of the to-be-orthogonalized direction "Z".
185: This information is enough to build the required entries
186: of H. The inner product with VEC_VV(it_loc) is
187: *different* than in the standard FGMRES and need to be dealt with specially.
188: That is, for standard FGMRES the orthogonalization coefficients are the same
189: as the coefficients used in the Arnoldi relation to reconstruct, but here this
190: is not true (albeit only for the one entry of H which we "unshift" below. */
192: /* Finish the dot product, retrieving the extra entry */
193: VecMDotEnd(ZVEC(loc_it), loc_it + 2, redux, lhh);
194: tt = PetscRealPart(lhh[loc_it + 1]);
196: /* Hessenberg entries, and entries for (naive) classical Graham-Schmidt
197: Note that the Hessenberg entries require a shift, as these are for the
198: relation AU = VH, which is wrt unshifted basis vectors */
199: hh = HH(0, loc_it);
200: hes = HES(0, loc_it);
201: for (j = 0; j < loc_it; j++) {
202: hh[j] = lhh[j];
203: hes[j] = lhh[j];
204: }
205: hh[loc_it] = lhh[loc_it] + shift;
206: hes[loc_it] = lhh[loc_it] + shift;
208: /* we delay applying the shift here */
209: for (j = 0; j <= loc_it; j++) { lhh[j] = -lhh[j]; /* flip sign */ }
211: /* Compute the norm of the un-normalized new direction using the rearranged formula
212: Note that these are shifted ("barred") quantities */
213: for (k = 0; k <= loc_it; k++) tt -= ((PetscReal)(PetscAbsScalar(lhh[k]) * PetscAbsScalar(lhh[k])));
214: /* On AVX512 this is accumulating roundoff errors for eg: tt=-2.22045e-16 */
215: if ((tt < 0.0) && tt > -PETSC_SMALL) tt = 0.0;
216: if (tt < 0.0) {
217: /* If we detect square root breakdown in the norm, we must restart the algorithm.
218: Here this means we simply break the current loop and reconstruct the solution
219: using the basis we have computed thus far. Note that by breaking immediately,
220: we do not update the iteration count, so computation done in this iteration
221: should be disregarded.
222: */
223: PetscInfo(ksp, "Restart due to square root breakdown at it = %" PetscInt_FMT ", tt=%g\n", ksp->its, (double)tt);
224: break;
225: } else {
226: tt = PetscSqrtReal(tt);
227: }
229: /* new entry in hessenburg is the 2-norm of our new direction */
230: hh[loc_it + 1] = tt;
231: hes[loc_it + 1] = tt;
233: /* The recurred computation for the new direction
234: The division by tt is delayed to the happy breakdown check later
235: Note placement BEFORE the unshift
236: */
237: VecCopy(ZVEC(loc_it), VEC_VV(loc_it + 1));
238: VecMAXPY(VEC_VV(loc_it + 1), loc_it + 1, lhh, &VEC_VV(0));
239: /* (VEC_VV(loc_it+1) is not normalized yet) */
241: /* The recurred computation for the preconditioned vector (u) */
242: VecCopy(Q, PREVEC(loc_it + 1));
243: VecMAXPY(PREVEC(loc_it + 1), loc_it + 1, lhh, &PREVEC(0));
244: VecScale(PREVEC(loc_it + 1), 1.0 / tt);
246: /* Unshift an entry in the GS coefficients ("removing the bar") */
247: lhh[loc_it] -= shift;
249: /* The recurred computation for z (Au)
250: Note placement AFTER the "unshift" */
251: VecCopy(W, ZVEC(loc_it + 1));
252: VecMAXPY(ZVEC(loc_it + 1), loc_it + 1, lhh, &ZVEC(0));
253: VecScale(ZVEC(loc_it + 1), 1.0 / tt);
255: /* Happy Breakdown Check */
256: hapbnd = PetscAbsScalar((tt) / *RS(loc_it));
257: /* RS(loc_it) contains the res_norm from the last iteration */
258: hapbnd = PetscMin(pipefgmres->haptol, hapbnd);
259: if (tt > hapbnd) {
260: /* scale new direction by its norm */
261: VecScale(VEC_VV(loc_it + 1), 1.0 / tt);
262: } else {
263: /* This happens when the solution is exactly reached. */
264: /* So there is no new direction... */
265: VecSet(VEC_TEMP, 0.0); /* set VEC_TEMP to 0 */
266: hapend = PETSC_TRUE;
267: }
268: /* note that for pipefgmres we could get HES(loc_it+1, loc_it) = 0 and the
269: current solution would not be exact if HES was singular. Note that
270: HH non-singular implies that HES is not singular, and HES is guaranteed
271: to be nonsingular when PREVECS are linearly independent and A is
272: nonsingular (in GMRES, the nonsingularity of A implies the nonsingularity
273: of HES). So we should really add a check to verify that HES is nonsingular.*/
275: /* Note that to be thorough, in debug mode, one could call a LAPACK routine
276: here to check that the hessenberg matrix is indeed non-singular (since
277: FGMRES does not guarantee this) */
279: /* Now apply rotations to new col of hessenberg (and right side of system),
280: calculate new rotation, and get new residual norm at the same time*/
281: KSPPIPEFGMRESUpdateHessenberg(ksp, loc_it, &hapend, &res_norm);
282: if (ksp->reason) break;
284: loc_it++;
285: pipefgmres->it = (loc_it - 1); /* Add this here in case it has converged */
287: PetscObjectSAWsTakeAccess((PetscObject)ksp);
288: ksp->its++;
289: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = res_norm;
290: else ksp->rnorm = 0;
291: PetscObjectSAWsGrantAccess((PetscObject)ksp);
293: (*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP);
295: /* Catch error in happy breakdown and signal convergence and break from loop */
296: if (hapend) {
297: if (!ksp->reason) {
299: ksp->reason = KSP_DIVERGED_BREAKDOWN;
300: break;
301: }
302: }
303: }
304: /* END OF ITERATION LOOP */
305: KSPLogResidualHistory(ksp, ksp->rnorm);
307: /*
308: Monitor if we know that we will not return for a restart */
309: if (loc_it && (ksp->reason || ksp->its >= ksp->max_it)) KSPMonitor(ksp, ksp->its, ksp->rnorm);
311: if (itcount) *itcount = loc_it;
313: /*
314: Down here we have to solve for the "best" coefficients of the Krylov
315: columns, add the solution values together, and possibly unwind the
316: preconditioning from the solution
317: */
319: /* Form the solution (or the solution so far) */
320: /* Note: must pass in (loc_it-1) for iteration count so that KSPPIPEGMRESIIBuildSoln
321: properly navigates */
323: KSPPIPEFGMRESBuildSoln(RS(0), ksp->vec_sol, ksp->vec_sol, ksp, loc_it - 1);
325: return 0;
326: }
328: /*
329: KSPSolve_PIPEFGMRES - This routine applies the PIPEFGMRES method.
331: Input Parameter:
332: . ksp - the Krylov space object that was set to use pipefgmres
334: Output Parameter:
335: . outits - number of iterations used
337: */
338: static PetscErrorCode KSPSolve_PIPEFGMRES(KSP ksp)
339: {
340: PetscInt its, itcount;
341: KSP_PIPEFGMRES *pipefgmres = (KSP_PIPEFGMRES *)ksp->data;
342: PetscBool guess_zero = ksp->guess_zero;
344: /* We have not checked these routines for use with complex numbers. The inner products
345: are likely not defined correctly for that case */
348: PetscCitationsRegister(citation, &cited);
351: PetscObjectSAWsTakeAccess((PetscObject)ksp);
352: ksp->its = 0;
353: PetscObjectSAWsGrantAccess((PetscObject)ksp);
355: itcount = 0;
356: ksp->reason = KSP_CONVERGED_ITERATING;
357: while (!ksp->reason) {
358: KSPInitialResidual(ksp, ksp->vec_sol, VEC_TEMP, VEC_TEMP_MATOP, VEC_VV(0), ksp->vec_rhs);
359: KSPPIPEFGMRESCycle(&its, ksp);
360: itcount += its;
361: if (itcount >= ksp->max_it) {
362: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
363: break;
364: }
365: ksp->guess_zero = PETSC_FALSE; /* every future call to KSPInitialResidual() will have nonzero guess */
366: }
367: ksp->guess_zero = guess_zero; /* restore if user provided nonzero initial guess */
368: return 0;
369: }
371: static PetscErrorCode KSPDestroy_PIPEFGMRES(KSP ksp)
372: {
373: KSPReset_PIPEFGMRES(ksp);
374: KSPDestroy_GMRES(ksp);
375: return 0;
376: }
378: /*
379: KSPPIPEFGMRESBuildSoln - create the solution from the starting vector and the
380: current iterates.
382: Input parameters:
383: nrs - work area of size it + 1.
384: vguess - index of initial guess
385: vdest - index of result. Note that vguess may == vdest (replace
386: guess with the solution).
387: it - HH upper triangular part is a block of size (it+1) x (it+1)
389: This is an internal routine that knows about the PIPEFGMRES internals.
390: */
391: static PetscErrorCode KSPPIPEFGMRESBuildSoln(PetscScalar *nrs, Vec vguess, Vec vdest, KSP ksp, PetscInt it)
392: {
393: PetscScalar tt;
394: PetscInt k, j;
395: KSP_PIPEFGMRES *pipefgmres = (KSP_PIPEFGMRES *)(ksp->data);
397: /* Solve for solution vector that minimizes the residual */
399: if (it < 0) { /* no pipefgmres steps have been performed */
400: VecCopy(vguess, vdest); /* VecCopy() is smart, exits immediately if vguess == vdest */
401: return 0;
402: }
404: /* solve the upper triangular system - RS is the right side and HH is
405: the upper triangular matrix - put soln in nrs */
406: if (*HH(it, it) != 0.0) nrs[it] = *RS(it) / *HH(it, it);
407: else nrs[it] = 0.0;
409: for (k = it - 1; k >= 0; k--) {
410: tt = *RS(k);
411: for (j = k + 1; j <= it; j++) tt -= *HH(k, j) * nrs[j];
412: nrs[k] = tt / *HH(k, k);
413: }
415: /* Accumulate the correction to the solution of the preconditioned problem in VEC_TEMP */
416: VecZeroEntries(VEC_TEMP);
417: VecMAXPY(VEC_TEMP, it + 1, nrs, &PREVEC(0));
419: /* add solution to previous solution */
420: if (vdest == vguess) {
421: VecAXPY(vdest, 1.0, VEC_TEMP);
422: } else {
423: VecWAXPY(vdest, 1.0, VEC_TEMP, vguess);
424: }
425: return 0;
426: }
428: /*
430: KSPPIPEFGMRESUpdateHessenberg - Do the scalar work for the orthogonalization.
431: Return new residual.
433: input parameters:
435: . ksp - Krylov space object
436: . it - plane rotations are applied to the (it+1)th column of the
437: modified hessenberg (i.e. HH(:,it))
438: . hapend - PETSC_FALSE not happy breakdown ending.
440: output parameters:
441: . res - the new residual
443: */
444: /*
445: . it - column of the Hessenberg that is complete, PIPEFGMRES is actually computing two columns ahead of this
446: */
447: static PetscErrorCode KSPPIPEFGMRESUpdateHessenberg(KSP ksp, PetscInt it, PetscBool *hapend, PetscReal *res)
448: {
449: PetscScalar *hh, *cc, *ss, *rs;
450: PetscInt j;
451: PetscReal hapbnd;
452: KSP_PIPEFGMRES *pipefgmres = (KSP_PIPEFGMRES *)(ksp->data);
454: hh = HH(0, it); /* pointer to beginning of column to update */
455: cc = CC(0); /* beginning of cosine rotations */
456: ss = SS(0); /* beginning of sine rotations */
457: rs = RS(0); /* right hand side of least squares system */
459: /* The Hessenberg matrix is now correct through column it, save that form for possible spectral analysis */
460: for (j = 0; j <= it + 1; j++) *HES(j, it) = hh[j];
462: /* check for the happy breakdown */
463: hapbnd = PetscMin(PetscAbsScalar(hh[it + 1] / rs[it]), pipefgmres->haptol);
464: if (PetscAbsScalar(hh[it + 1]) < hapbnd) {
465: PetscInfo(ksp, "Detected happy breakdown, current hapbnd = %14.12e H(%" PetscInt_FMT ",%" PetscInt_FMT ") = %14.12e\n", (double)hapbnd, it + 1, it, (double)PetscAbsScalar(*HH(it + 1, it)));
466: *hapend = PETSC_TRUE;
467: }
469: /* Apply all the previously computed plane rotations to the new column
470: of the Hessenberg matrix */
471: /* Note: this uses the rotation [conj(c) s ; -s c], c= cos(theta), s= sin(theta),
472: and some refs have [c s ; -conj(s) c] (don't be confused!) */
474: for (j = 0; j < it; j++) {
475: PetscScalar hhj = hh[j];
476: hh[j] = PetscConj(cc[j]) * hhj + ss[j] * hh[j + 1];
477: hh[j + 1] = -ss[j] * hhj + cc[j] * hh[j + 1];
478: }
480: /*
481: compute the new plane rotation, and apply it to:
482: 1) the right-hand-side of the Hessenberg system (RS)
483: note: it affects RS(it) and RS(it+1)
484: 2) the new column of the Hessenberg matrix
485: note: it affects HH(it,it) which is currently pointed to
486: by hh and HH(it+1, it) (*(hh+1))
487: thus obtaining the updated value of the residual...
488: */
490: /* compute new plane rotation */
492: if (!*hapend) {
493: PetscReal delta = PetscSqrtReal(PetscSqr(PetscAbsScalar(hh[it])) + PetscSqr(PetscAbsScalar(hh[it + 1])));
494: if (delta == 0.0) {
495: ksp->reason = KSP_DIVERGED_NULL;
496: return 0;
497: }
499: cc[it] = hh[it] / delta; /* new cosine value */
500: ss[it] = hh[it + 1] / delta; /* new sine value */
502: hh[it] = PetscConj(cc[it]) * hh[it] + ss[it] * hh[it + 1];
503: rs[it + 1] = -ss[it] * rs[it];
504: rs[it] = PetscConj(cc[it]) * rs[it];
505: *res = PetscAbsScalar(rs[it + 1]);
506: } else { /* happy breakdown: HH(it+1, it) = 0, therefore we don't need to apply
507: another rotation matrix (so RH doesn't change). The new residual is
508: always the new sine term times the residual from last time (RS(it)),
509: but now the new sine rotation would be zero...so the residual should
510: be zero...so we will multiply "zero" by the last residual. This might
511: not be exactly what we want to do here -could just return "zero". */
513: *res = 0.0;
514: }
515: return 0;
516: }
518: /*
519: KSPBuildSolution_PIPEFGMRES
521: Input Parameter:
522: . ksp - the Krylov space object
523: . ptr-
525: Output Parameter:
526: . result - the solution
528: Note: this calls KSPPIPEFGMRESBuildSoln - the same function that KSPPIPEFGMRESCycle
529: calls directly.
531: */
532: PetscErrorCode KSPBuildSolution_PIPEFGMRES(KSP ksp, Vec ptr, Vec *result)
533: {
534: KSP_PIPEFGMRES *pipefgmres = (KSP_PIPEFGMRES *)ksp->data;
536: if (!ptr) {
537: if (!pipefgmres->sol_temp) { VecDuplicate(ksp->vec_sol, &pipefgmres->sol_temp); }
538: ptr = pipefgmres->sol_temp;
539: }
540: if (!pipefgmres->nrs) {
541: /* allocate the work area */
542: PetscMalloc1(pipefgmres->max_k, &pipefgmres->nrs);
543: }
545: KSPPIPEFGMRESBuildSoln(pipefgmres->nrs, ksp->vec_sol, ptr, ksp, pipefgmres->it);
546: if (result) *result = ptr;
547: return 0;
548: }
550: PetscErrorCode KSPSetFromOptions_PIPEFGMRES(KSP ksp, PetscOptionItems *PetscOptionsObject)
551: {
552: KSP_PIPEFGMRES *pipefgmres = (KSP_PIPEFGMRES *)ksp->data;
553: PetscBool flg;
554: PetscScalar shift;
556: KSPSetFromOptions_GMRES(ksp, PetscOptionsObject);
557: PetscOptionsHeadBegin(PetscOptionsObject, "KSP pipelined FGMRES Options");
558: PetscOptionsScalar("-ksp_pipefgmres_shift", "shift parameter", "KSPPIPEFGMRESSetShift", pipefgmres->shift, &shift, &flg);
559: if (flg) KSPPIPEFGMRESSetShift(ksp, shift);
560: PetscOptionsHeadEnd();
561: return 0;
562: }
564: PetscErrorCode KSPView_PIPEFGMRES(KSP ksp, PetscViewer viewer)
565: {
566: KSP_PIPEFGMRES *pipefgmres = (KSP_PIPEFGMRES *)ksp->data;
567: PetscBool iascii, isstring;
569: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
570: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
572: if (iascii) {
573: PetscViewerASCIIPrintf(viewer, " restart=%" PetscInt_FMT "\n", pipefgmres->max_k);
574: PetscViewerASCIIPrintf(viewer, " happy breakdown tolerance %g\n", (double)pipefgmres->haptol);
575: #if defined(PETSC_USE_COMPLEX)
576: PetscViewerASCIIPrintf(viewer, " shift=%g+%gi\n", (double)PetscRealPart(pipefgmres->shift), (double)PetscImaginaryPart(pipefgmres->shift));
577: #else
578: PetscViewerASCIIPrintf(viewer, " shift=%g\n", (double)pipefgmres->shift);
579: #endif
580: } else if (isstring) {
581: PetscViewerStringSPrintf(viewer, "restart %" PetscInt_FMT, pipefgmres->max_k);
582: #if defined(PETSC_USE_COMPLEX)
583: PetscViewerStringSPrintf(viewer, " shift=%g+%gi\n", (double)PetscRealPart(pipefgmres->shift), (double)PetscImaginaryPart(pipefgmres->shift));
584: #else
585: PetscViewerStringSPrintf(viewer, " shift=%g\n", (double)pipefgmres->shift);
586: #endif
587: }
588: return 0;
589: }
591: PetscErrorCode KSPReset_PIPEFGMRES(KSP ksp)
592: {
593: KSP_PIPEFGMRES *pipefgmres = (KSP_PIPEFGMRES *)ksp->data;
594: PetscInt i;
596: PetscFree(pipefgmres->prevecs);
597: PetscFree(pipefgmres->zvecs);
598: for (i = 0; i < pipefgmres->nwork_alloc; i++) {
599: VecDestroyVecs(pipefgmres->mwork_alloc[i], &pipefgmres->prevecs_user_work[i]);
600: VecDestroyVecs(pipefgmres->mwork_alloc[i], &pipefgmres->zvecs_user_work[i]);
601: }
602: PetscFree(pipefgmres->prevecs_user_work);
603: PetscFree(pipefgmres->zvecs_user_work);
604: PetscFree(pipefgmres->redux);
605: KSPReset_GMRES(ksp);
606: return 0;
607: }
609: /*MC
610: KSPPIPEFGMRES - Implements the Pipelined (1-stage) Flexible Generalized Minimal Residual method. [](sec_pipelineksp). [](sec_flexibleksp)
612: Options Database Keys:
613: + -ksp_gmres_restart <restart> - the number of Krylov directions to orthogonalize against
614: . -ksp_gmres_haptol <tol> - sets the tolerance for "happy ending" (exact convergence)
615: . -ksp_gmres_preallocate - preallocate all the Krylov search directions initially (otherwise groups of
616: . -ksp_pipefgmres_shift - the shift to use (defaults to 1. See KSPPIPEFGMRESSetShift()
617: vectors are allocated as needed)
618: - -ksp_gmres_krylov_monitor - plot the Krylov space generated
620: Level: intermediate
622: Notes:
623: This variant is not "explicitly normalized" like `KSPPGMRES`, and requires a shift parameter.
625: A heuristic for choosing the shift parameter is the largest eigenvalue of the preconditioned operator.
627: Only right preconditioning is supported (but this preconditioner may be nonlinear/variable/inexact, as with `KSPFGMRES`).
629: MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
630: See [](doc_faq_pipelined)
632: Developer Note:
633: This class is subclassed off of `KSPGMRES`.
635: Contributed by:
636: P. Sanan and S.M. Schnepp
638: Reference:
639: P. Sanan, S.M. Schnepp, and D.A. May,
640: "Pipelined, Flexible Krylov Subspace Methods,"
641: SIAM Journal on Scientific Computing 2016 38:5, C441-C470,
642: DOI: 10.1137/15M1049130
644: .seealso: [](chapter_ksp), [](doc_faq_pipelined), [](sec_pipelineksp), [](sec_flexibleksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPLGMRES`, `KSPPIPECG`, `KSPPIPECR`, `KSPPGMRES`, `KSPFGMRES`
645: `KSPGMRESSetRestart()`, `KSPGMRESSetHapTol()`, `KSPGMRESSetPreAllocateVectors()`, `KSPGMRESMonitorKrylov()`, `KSPPIPEFGMRESSetShift()`
646: M*/
648: PETSC_EXTERN PetscErrorCode KSPCreate_PIPEFGMRES(KSP ksp)
649: {
650: KSP_PIPEFGMRES *pipefgmres;
652: PetscNew(&pipefgmres);
654: ksp->data = (void *)pipefgmres;
655: ksp->ops->buildsolution = KSPBuildSolution_PIPEFGMRES;
656: ksp->ops->setup = KSPSetUp_PIPEFGMRES;
657: ksp->ops->solve = KSPSolve_PIPEFGMRES;
658: ksp->ops->reset = KSPReset_PIPEFGMRES;
659: ksp->ops->destroy = KSPDestroy_PIPEFGMRES;
660: ksp->ops->view = KSPView_PIPEFGMRES;
661: ksp->ops->setfromoptions = KSPSetFromOptions_PIPEFGMRES;
662: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_GMRES;
663: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_GMRES;
665: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 3);
666: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1);
668: PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetPreAllocateVectors_C", KSPGMRESSetPreAllocateVectors_GMRES);
669: PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESSetRestart_C", KSPGMRESSetRestart_GMRES);
670: PetscObjectComposeFunction((PetscObject)ksp, "KSPGMRESGetRestart_C", KSPGMRESGetRestart_GMRES);
672: pipefgmres->nextra_vecs = 1;
673: pipefgmres->haptol = 1.0e-30;
674: pipefgmres->q_preallocate = 0;
675: pipefgmres->delta_allocate = PIPEFGMRES_DELTA_DIRECTIONS;
676: pipefgmres->orthog = NULL;
677: pipefgmres->nrs = NULL;
678: pipefgmres->sol_temp = NULL;
679: pipefgmres->max_k = PIPEFGMRES_DEFAULT_MAXK;
680: pipefgmres->Rsvd = NULL;
681: pipefgmres->orthogwork = NULL;
682: pipefgmres->cgstype = KSP_GMRES_CGS_REFINE_NEVER;
683: pipefgmres->shift = 1.0;
684: return 0;
685: }
687: static PetscErrorCode KSPPIPEFGMRESGetNewVectors(KSP ksp, PetscInt it)
688: {
689: KSP_PIPEFGMRES *pipefgmres = (KSP_PIPEFGMRES *)ksp->data;
690: PetscInt nwork = pipefgmres->nwork_alloc; /* number of work vector chunks allocated */
691: PetscInt nalloc; /* number to allocate */
692: PetscInt k;
694: nalloc = pipefgmres->delta_allocate; /* number of vectors to allocate
695: in a single chunk */
697: /* Adjust the number to allocate to make sure that we don't exceed the
698: number of available slots (pipefgmres->vecs_allocated)*/
699: if (it + VEC_OFFSET + nalloc >= pipefgmres->vecs_allocated) nalloc = pipefgmres->vecs_allocated - it - VEC_OFFSET;
700: if (!nalloc) return 0;
702: pipefgmres->vv_allocated += nalloc; /* vv_allocated is the number of vectors allocated */
704: /* work vectors */
705: KSPCreateVecs(ksp, nalloc, &pipefgmres->user_work[nwork], 0, NULL);
706: for (k = 0; k < nalloc; k++) pipefgmres->vecs[it + VEC_OFFSET + k] = pipefgmres->user_work[nwork][k];
707: /* specify size of chunk allocated */
708: pipefgmres->mwork_alloc[nwork] = nalloc;
710: /* preconditioned vectors (note we don't use VEC_OFFSET) */
711: KSPCreateVecs(ksp, nalloc, &pipefgmres->prevecs_user_work[nwork], 0, NULL);
712: for (k = 0; k < nalloc; k++) pipefgmres->prevecs[it + k] = pipefgmres->prevecs_user_work[nwork][k];
714: KSPCreateVecs(ksp, nalloc, &pipefgmres->zvecs_user_work[nwork], 0, NULL);
715: for (k = 0; k < nalloc; k++) pipefgmres->zvecs[it + k] = pipefgmres->zvecs_user_work[nwork][k];
717: /* increment the number of work vector chunks */
718: pipefgmres->nwork_alloc++;
719: return 0;
720: }
722: /*@
723: KSPPIPEFGMRESSetShift - Set the shift parameter for the flexible, pipelined `KSPPIPEFGMRES` solver.
725: Logically Collective
727: Input Parameters:
728: + ksp - the Krylov space context
729: - shift - the shift
731: Options Database Key:
732: . -ksp_pipefgmres_shift <shift> - set the shift parameter
734: Level: intermediate
736: Note:
737: A heuristic is to set this to be comparable to the largest eigenvalue of the preconditioned operator.
738: This can be achieved with PETSc itself by using a few iterations of a Krylov method.
739: See `KSPComputeEigenvalues()` (and note the caveats there).
741: .seealso: [](chapter_ksp), `KSPPIPEFGMRES`, `KSPComputeEigenvalues()`
742: @*/
743: PetscErrorCode KSPPIPEFGMRESSetShift(KSP ksp, PetscScalar shift)
744: {
745: KSP_PIPEFGMRES *pipefgmres = (KSP_PIPEFGMRES *)ksp->data;
749: pipefgmres->shift = shift;
750: return 0;
751: }