Actual source code: pipefcg.c
1: #include <../src/ksp/ksp/impls/fcg/pipefcg/pipefcgimpl.h>
3: static PetscBool cited = PETSC_FALSE;
4: static const char citation[] = "@article{SSM2016,\n"
5: " author = {P. Sanan and S.M. Schnepp and D.A. May},\n"
6: " title = {Pipelined, Flexible Krylov Subspace Methods},\n"
7: " journal = {SIAM Journal on Scientific Computing},\n"
8: " volume = {38},\n"
9: " number = {5},\n"
10: " pages = {C441-C470},\n"
11: " year = {2016},\n"
12: " doi = {10.1137/15M1049130},\n"
13: " URL = {http://dx.doi.org/10.1137/15M1049130},\n"
14: " eprint = {http://dx.doi.org/10.1137/15M1049130}\n"
15: "}\n";
17: #define KSPPIPEFCG_DEFAULT_MMAX 15
18: #define KSPPIPEFCG_DEFAULT_NPREALLOC 5
19: #define KSPPIPEFCG_DEFAULT_VECB 5
20: #define KSPPIPEFCG_DEFAULT_TRUNCSTRAT KSP_FCD_TRUNC_TYPE_NOTAY
22: static PetscErrorCode KSPAllocateVectors_PIPEFCG(KSP ksp, PetscInt nvecsneeded, PetscInt chunksize)
23: {
24: KSP_PIPEFCG *pipefcg;
25: PetscInt nnewvecs, nvecsprev;
27: PetscFunctionBegin;
28: pipefcg = (KSP_PIPEFCG *)ksp->data;
30: /* Allocate enough new vectors to add chunksize new vectors, reach nvecsneedtotal, or to reach mmax+1, whichever is smallest */
31: if (pipefcg->nvecs < PetscMin(pipefcg->mmax + 1, nvecsneeded)) {
32: nvecsprev = pipefcg->nvecs;
33: nnewvecs = PetscMin(PetscMax(nvecsneeded - pipefcg->nvecs, chunksize), pipefcg->mmax + 1 - pipefcg->nvecs);
34: PetscCall(KSPCreateVecs(ksp, nnewvecs, &pipefcg->pQvecs[pipefcg->nchunks], 0, NULL));
35: PetscCall(KSPCreateVecs(ksp, nnewvecs, &pipefcg->pZETAvecs[pipefcg->nchunks], 0, NULL));
36: PetscCall(KSPCreateVecs(ksp, nnewvecs, &pipefcg->pPvecs[pipefcg->nchunks], 0, NULL));
37: PetscCall(KSPCreateVecs(ksp, nnewvecs, &pipefcg->pSvecs[pipefcg->nchunks], 0, NULL));
38: pipefcg->nvecs += nnewvecs;
39: for (PetscInt i = 0; i < nnewvecs; ++i) {
40: pipefcg->Qvecs[nvecsprev + i] = pipefcg->pQvecs[pipefcg->nchunks][i];
41: pipefcg->ZETAvecs[nvecsprev + i] = pipefcg->pZETAvecs[pipefcg->nchunks][i];
42: pipefcg->Pvecs[nvecsprev + i] = pipefcg->pPvecs[pipefcg->nchunks][i];
43: pipefcg->Svecs[nvecsprev + i] = pipefcg->pSvecs[pipefcg->nchunks][i];
44: }
45: pipefcg->chunksizes[pipefcg->nchunks] = nnewvecs;
46: ++pipefcg->nchunks;
47: }
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: static PetscErrorCode KSPSetUp_PIPEFCG(KSP ksp)
52: {
53: KSP_PIPEFCG *pipefcg;
54: const PetscInt nworkstd = 5;
56: PetscFunctionBegin;
57: pipefcg = (KSP_PIPEFCG *)ksp->data;
59: /* Allocate "standard" work vectors (not including the basis and transformed basis vectors) */
60: PetscCall(KSPSetWorkVecs(ksp, nworkstd));
62: /* Allocated space for pointers to additional work vectors
63: note that mmax is the number of previous directions, so we add 1 for the current direction,
64: and an extra 1 for the prealloc (which might be empty) */
65: PetscCall(PetscMalloc4(pipefcg->mmax + 1, &pipefcg->Pvecs, pipefcg->mmax + 1, &pipefcg->pPvecs, pipefcg->mmax + 1, &pipefcg->Svecs, pipefcg->mmax + 1, &pipefcg->pSvecs));
66: PetscCall(PetscMalloc4(pipefcg->mmax + 1, &pipefcg->Qvecs, pipefcg->mmax + 1, &pipefcg->pQvecs, pipefcg->mmax + 1, &pipefcg->ZETAvecs, pipefcg->mmax + 1, &pipefcg->pZETAvecs));
67: PetscCall(PetscMalloc4(pipefcg->mmax + 1, &pipefcg->Pold, pipefcg->mmax + 1, &pipefcg->Sold, pipefcg->mmax + 1, &pipefcg->Qold, pipefcg->mmax + 1, &pipefcg->ZETAold));
68: PetscCall(PetscMalloc1(pipefcg->mmax + 1, &pipefcg->chunksizes));
69: PetscCall(PetscMalloc3(pipefcg->mmax + 2, &pipefcg->dots, pipefcg->mmax + 1, &pipefcg->etas, pipefcg->mmax + 2, &pipefcg->redux));
71: /* If the requested number of preallocated vectors is greater than mmax reduce nprealloc */
72: if (pipefcg->nprealloc > pipefcg->mmax + 1) PetscCall(PetscInfo(NULL, "Requested nprealloc=%" PetscInt_FMT " is greater than m_max+1=%" PetscInt_FMT ". Resetting nprealloc = m_max+1.\n", pipefcg->nprealloc, pipefcg->mmax + 1));
74: /* Preallocate additional work vectors */
75: PetscCall(KSPAllocateVectors_PIPEFCG(ksp, pipefcg->nprealloc, pipefcg->nprealloc));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: static PetscErrorCode KSPSolve_PIPEFCG_cycle(KSP ksp)
80: {
81: PetscInt i, j, k, idx, kdx, mi;
82: KSP_PIPEFCG *pipefcg;
83: PetscScalar alpha = 0.0, gamma, *betas, *dots;
84: PetscReal dp = 0.0, delta, *eta, *etas;
85: Vec B, R, Z, X, Qcurr, W, ZETAcurr, M, N, Pcurr, Scurr, *redux;
86: Mat Amat, Pmat;
88: PetscFunctionBegin;
89: /* We have not checked these routines for use with complex numbers. The inner products
90: are likely not defined correctly for that case */
91: PetscCheck(!PetscDefined(USE_COMPLEX) || PetscDefined(SKIP_COMPLEX), PETSC_COMM_WORLD, PETSC_ERR_SUP, "PIPEFGMRES has not been implemented for use with complex scalars");
93: #define VecXDot(x, y, a) (pipefcg->type == KSP_CG_HERMITIAN ? VecDot(x, y, a) : VecTDot(x, y, a))
94: #define VecXDotBegin(x, y, a) (pipefcg->type == KSP_CG_HERMITIAN ? VecDotBegin(x, y, a) : VecTDotBegin(x, y, a))
95: #define VecXDotEnd(x, y, a) (pipefcg->type == KSP_CG_HERMITIAN ? VecDotEnd(x, y, a) : VecTDotEnd(x, y, a))
96: #define VecMXDot(x, n, y, a) (pipefcg->type == KSP_CG_HERMITIAN ? VecMDot(x, n, y, a) : VecMTDot(x, n, y, a))
97: #define VecMXDotBegin(x, n, y, a) (pipefcg->type == KSP_CG_HERMITIAN ? VecMDotBegin(x, n, y, a) : VecMTDotBegin(x, n, y, a))
98: #define VecMXDotEnd(x, n, y, a) (pipefcg->type == KSP_CG_HERMITIAN ? VecMDotEnd(x, n, y, a) : VecMTDotEnd(x, n, y, a))
100: pipefcg = (KSP_PIPEFCG *)ksp->data;
101: X = ksp->vec_sol;
102: B = ksp->vec_rhs;
103: R = ksp->work[0];
104: Z = ksp->work[1];
105: W = ksp->work[2];
106: M = ksp->work[3];
107: N = ksp->work[4];
109: redux = pipefcg->redux;
110: dots = pipefcg->dots;
111: etas = pipefcg->etas;
112: betas = dots; /* dots takes the result of all dot products of which the betas are a subset */
114: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
116: /* Compute cycle initial residual */
117: PetscCall(KSP_MatMult(ksp, Amat, X, R));
118: PetscCall(VecAYPX(R, -1.0, B)); /* r <- b - Ax */
119: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
121: Pcurr = pipefcg->Pvecs[0];
122: Scurr = pipefcg->Svecs[0];
123: Qcurr = pipefcg->Qvecs[0];
124: ZETAcurr = pipefcg->ZETAvecs[0];
125: PetscCall(VecCopy(Z, Pcurr));
126: PetscCall(KSP_MatMult(ksp, Amat, Pcurr, Scurr)); /* S = Ap */
127: PetscCall(VecCopy(Scurr, W)); /* w = s = Az */
129: /* Initial state of pipelining intermediates */
130: redux[0] = R;
131: redux[1] = W;
132: PetscCall(VecMXDotBegin(Z, 2, redux, dots));
133: PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)Z))); /* perform asynchronous reduction */
134: PetscCall(KSP_PCApply(ksp, W, M)); /* m = B(w) */
135: PetscCall(KSP_MatMult(ksp, Amat, M, N)); /* n = Am */
136: PetscCall(VecCopy(M, Qcurr)); /* q = m */
137: PetscCall(VecCopy(N, ZETAcurr)); /* zeta = n */
138: PetscCall(VecMXDotEnd(Z, 2, redux, dots));
139: gamma = dots[0];
140: delta = PetscRealPart(dots[1]);
141: etas[0] = delta;
142: alpha = gamma / delta;
144: i = 0;
145: do {
146: ksp->its++;
148: /* Update X, R, Z, W */
149: PetscCall(VecAXPY(X, +alpha, Pcurr)); /* x <- x + alpha * pi */
150: PetscCall(VecAXPY(R, -alpha, Scurr)); /* r <- r - alpha * si */
151: PetscCall(VecAXPY(Z, -alpha, Qcurr)); /* z <- z - alpha * qi */
152: PetscCall(VecAXPY(W, -alpha, ZETAcurr)); /* w <- w - alpha * zetai */
154: /* Compute norm for convergence check */
155: switch (ksp->normtype) {
156: case KSP_NORM_PRECONDITIONED:
157: PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- sqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */
158: break;
159: case KSP_NORM_UNPRECONDITIONED:
160: PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e) */
161: break;
162: case KSP_NORM_NATURAL:
163: dp = PetscSqrtReal(PetscAbsScalar(gamma)); /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e) */
164: break;
165: case KSP_NORM_NONE:
166: dp = 0.0;
167: break;
168: default:
169: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
170: }
172: /* Check for convergence */
173: ksp->rnorm = dp;
174: PetscCall(KSPLogResidualHistory(ksp, dp));
175: PetscCall(KSPMonitor(ksp, ksp->its, dp));
176: PetscCall((*ksp->converged)(ksp, ksp->its, dp, &ksp->reason, ksp->cnvP));
177: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
179: /* Computations of current iteration done */
180: ++i;
182: /* If needbe, allocate a new chunk of vectors in P and C */
183: PetscCall(KSPAllocateVectors_PIPEFCG(ksp, i + 1, pipefcg->vecb));
185: /* Note that we wrap around and start clobbering old vectors */
186: idx = i % (pipefcg->mmax + 1);
187: Pcurr = pipefcg->Pvecs[idx];
188: Scurr = pipefcg->Svecs[idx];
189: Qcurr = pipefcg->Qvecs[idx];
190: ZETAcurr = pipefcg->ZETAvecs[idx];
191: eta = pipefcg->etas + idx;
193: /* number of old directions to orthogonalize against */
194: switch (pipefcg->truncstrat) {
195: case KSP_FCD_TRUNC_TYPE_STANDARD:
196: mi = pipefcg->mmax;
197: break;
198: case KSP_FCD_TRUNC_TYPE_NOTAY:
199: mi = ((i - 1) % pipefcg->mmax) + 1;
200: break;
201: default:
202: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized Truncation Strategy");
203: }
205: /* Pick old p,s,q,zeta in a way suitable for VecMDot */
206: PetscCall(VecCopy(Z, Pcurr));
207: for (k = PetscMax(0, i - mi), j = 0; k < i; ++j, ++k) {
208: kdx = k % (pipefcg->mmax + 1);
209: pipefcg->Pold[j] = pipefcg->Pvecs[kdx];
210: pipefcg->Sold[j] = pipefcg->Svecs[kdx];
211: pipefcg->Qold[j] = pipefcg->Qvecs[kdx];
212: pipefcg->ZETAold[j] = pipefcg->ZETAvecs[kdx];
213: redux[j] = pipefcg->Svecs[kdx];
214: }
215: redux[j] = R; /* If the above loop is not executed redux contains only R => all beta_k = 0, only gamma, delta != 0 */
216: redux[j + 1] = W;
218: PetscCall(VecMXDotBegin(Z, j + 2, redux, betas)); /* Start split reductions for beta_k = (z,s_k), gamma = (z,r), delta = (z,w) */
219: PetscCall(PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)Z))); /* perform asynchronous reduction */
220: PetscCall(VecWAXPY(N, -1.0, R, W)); /* m = u + B(w-r): (a) ntmp = w-r */
221: PetscCall(KSP_PCApply(ksp, N, M)); /* m = u + B(w-r): (b) mtmp = B(ntmp) = B(w-r) */
222: PetscCall(VecAXPY(M, 1.0, Z)); /* m = u + B(w-r): (c) m = z + mtmp */
223: PetscCall(KSP_MatMult(ksp, Amat, M, N)); /* n = Am */
224: PetscCall(VecMXDotEnd(Z, j + 2, redux, betas)); /* Finish split reductions */
225: gamma = betas[j];
226: delta = PetscRealPart(betas[j + 1]);
228: *eta = 0.;
229: for (k = PetscMax(0, i - mi), j = 0; k < i; ++j, ++k) {
230: kdx = k % (pipefcg->mmax + 1);
231: betas[j] /= -etas[kdx]; /* betak /= etak */
232: *eta -= PetscAbsScalar(betas[j]) * PetscAbsScalar(betas[j]) * etas[kdx];
233: /* etaitmp = -betaik^2 * etak */
234: }
235: *eta += delta; /* etai = delta -betaik^2 * etak */
236: if (*eta < 0.) {
237: pipefcg->norm_breakdown = PETSC_TRUE;
238: PetscCall(PetscInfo(ksp, "Restart due to square root breakdown at it = %" PetscInt_FMT "\n", ksp->its));
239: break;
240: } else {
241: alpha = gamma / (*eta); /* alpha = gamma/etai */
242: }
244: /* project out stored search directions using classical G-S */
245: PetscCall(VecCopy(Z, Pcurr));
246: PetscCall(VecCopy(W, Scurr));
247: PetscCall(VecCopy(M, Qcurr));
248: PetscCall(VecCopy(N, ZETAcurr));
249: PetscCall(VecMAXPY(Pcurr, j, betas, pipefcg->Pold)); /* pi <- ui - sum_k beta_k p_k */
250: PetscCall(VecMAXPY(Scurr, j, betas, pipefcg->Sold)); /* si <- wi - sum_k beta_k s_k */
251: PetscCall(VecMAXPY(Qcurr, j, betas, pipefcg->Qold)); /* qi <- m - sum_k beta_k q_k */
252: PetscCall(VecMAXPY(ZETAcurr, j, betas, pipefcg->ZETAold)); /* zetai <- n - sum_k beta_k zeta_k */
254: } while (ksp->its < ksp->max_it);
255: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: static PetscErrorCode KSPSolve_PIPEFCG(KSP ksp)
260: {
261: KSP_PIPEFCG *pipefcg;
262: PetscScalar gamma;
263: PetscReal dp = 0.0;
264: Vec B, R, Z, X;
265: Mat Amat, Pmat;
267: #define VecXDot(x, y, a) (pipefcg->type == KSP_CG_HERMITIAN ? VecDot(x, y, a) : VecTDot(x, y, a))
269: PetscFunctionBegin;
270: PetscCall(PetscCitationsRegister(citation, &cited));
272: pipefcg = (KSP_PIPEFCG *)ksp->data;
273: X = ksp->vec_sol;
274: B = ksp->vec_rhs;
275: R = ksp->work[0];
276: Z = ksp->work[1];
278: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
280: /* Compute initial residual needed for convergence check*/
281: ksp->its = 0;
282: if (!ksp->guess_zero) {
283: PetscCall(KSP_MatMult(ksp, Amat, X, R));
284: PetscCall(VecAYPX(R, -1.0, B)); /* r <- b - Ax */
285: } else {
286: PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
287: }
288: switch (ksp->normtype) {
289: case KSP_NORM_PRECONDITIONED:
290: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
291: PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- dqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */
292: break;
293: case KSP_NORM_UNPRECONDITIONED:
294: PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e) */
295: break;
296: case KSP_NORM_NATURAL:
297: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
298: PetscCall(VecXDot(Z, R, &gamma));
299: dp = PetscSqrtReal(PetscAbsScalar(gamma)); /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e) */
300: break;
301: case KSP_NORM_NONE:
302: dp = 0.0;
303: break;
304: default:
305: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
306: }
308: /* Initial Convergence Check */
309: PetscCall(KSPLogResidualHistory(ksp, dp));
310: PetscCall(KSPMonitor(ksp, 0, dp));
311: ksp->rnorm = dp;
312: PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP));
313: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
315: do {
316: /* A cycle is broken only if a norm breakdown occurs. If not the entire solve happens in a single cycle.
317: This is coded this way to allow both truncation and truncation-restart strategy
318: (see KSPFCDGetNumOldDirections()) */
319: PetscCall(KSPSolve_PIPEFCG_cycle(ksp));
320: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
321: if (pipefcg->norm_breakdown) {
322: pipefcg->n_restarts++;
323: pipefcg->norm_breakdown = PETSC_FALSE;
324: }
325: } while (ksp->its < ksp->max_it);
327: if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: static PetscErrorCode KSPDestroy_PIPEFCG(KSP ksp)
332: {
333: KSP_PIPEFCG *pipefcg;
335: PetscFunctionBegin;
336: pipefcg = (KSP_PIPEFCG *)ksp->data;
338: /* Destroy "standard" work vecs */
339: PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));
341: /* Destroy vectors of old directions and the arrays that manage pointers to them */
342: if (pipefcg->nvecs) {
343: for (PetscInt i = 0; i < pipefcg->nchunks; ++i) {
344: PetscCall(VecDestroyVecs(pipefcg->chunksizes[i], &pipefcg->pPvecs[i]));
345: PetscCall(VecDestroyVecs(pipefcg->chunksizes[i], &pipefcg->pSvecs[i]));
346: PetscCall(VecDestroyVecs(pipefcg->chunksizes[i], &pipefcg->pQvecs[i]));
347: PetscCall(VecDestroyVecs(pipefcg->chunksizes[i], &pipefcg->pZETAvecs[i]));
348: }
349: }
350: PetscCall(PetscFree4(pipefcg->Pvecs, pipefcg->Svecs, pipefcg->pPvecs, pipefcg->pSvecs));
351: PetscCall(PetscFree4(pipefcg->Qvecs, pipefcg->ZETAvecs, pipefcg->pQvecs, pipefcg->pZETAvecs));
352: PetscCall(PetscFree4(pipefcg->Pold, pipefcg->Sold, pipefcg->Qold, pipefcg->ZETAold));
353: PetscCall(PetscFree(pipefcg->chunksizes));
354: PetscCall(PetscFree3(pipefcg->dots, pipefcg->etas, pipefcg->redux));
355: PetscCall(KSPDestroyDefault(ksp));
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: static PetscErrorCode KSPView_PIPEFCG(KSP ksp, PetscViewer viewer)
360: {
361: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
362: PetscBool isascii, isstring;
363: const char *truncstr;
365: PetscFunctionBegin;
366: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
367: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
369: if (pipefcg->truncstrat == KSP_FCD_TRUNC_TYPE_STANDARD) {
370: truncstr = "Using standard truncation strategy";
371: } else if (pipefcg->truncstrat == KSP_FCD_TRUNC_TYPE_NOTAY) {
372: truncstr = "Using Notay's truncation strategy";
373: } else {
374: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Undefined FCD truncation strategy");
375: }
377: if (isascii) {
378: PetscCall(PetscViewerASCIIPrintf(viewer, " max previous directions = %" PetscInt_FMT "\n", pipefcg->mmax));
379: PetscCall(PetscViewerASCIIPrintf(viewer, " preallocated %" PetscInt_FMT " directions\n", PetscMin(pipefcg->nprealloc, pipefcg->mmax + 1)));
380: PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", truncstr));
381: PetscCall(PetscViewerASCIIPrintf(viewer, " restarts performed = %" PetscInt_FMT " \n", pipefcg->n_restarts));
382: } else if (isstring) {
383: PetscCall(PetscViewerStringSPrintf(viewer, "max previous directions = %" PetscInt_FMT ", preallocated %" PetscInt_FMT " directions, %s truncation strategy", pipefcg->mmax, pipefcg->nprealloc, truncstr));
384: }
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: /*@
389: KSPPIPEFCGSetMmax - set the maximum number of previous directions `KSPPIPEFCG` will store for orthogonalization
391: Logically Collective
393: Input Parameters:
394: + ksp - the Krylov space context
395: - mmax - the maximum number of previous directions to orthogonalize against
397: Options Database Key:
398: . -ksp_pipefcg_mmax N - maximum number of previous directions
400: Level: intermediate
402: Note:
403: `mmax` + 1 directions are stored (`mmax` previous ones along with the current one)
404: and whether all are used in each iteration also depends on the truncation strategy, see `KSPPIPEFCGSetTruncationType()`
406: .seealso: [](ch_ksp), `KSPPIPEFCG`, `KSPPIPEFCGSetTruncationType()`, `KSPPIPEFCGSetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
407: @*/
408: PetscErrorCode KSPPIPEFCGSetMmax(KSP ksp, PetscInt mmax)
409: {
410: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
412: PetscFunctionBegin;
415: pipefcg->mmax = mmax;
416: PetscFunctionReturn(PETSC_SUCCESS);
417: }
419: /*@
420: KSPPIPEFCGGetMmax - get the maximum number of previous directions `KSPPIPEFCG` will store
422: Not Collective
424: Input Parameter:
425: . ksp - the Krylov space context
427: Output Parameter:
428: . mmax - the maximum number of previous directions allowed for orthogonalization
430: Level: intermediate
432: .seealso: [](ch_ksp), `KSPPIPEFCG`, `KSPPIPEFCGGetTruncationType()`, `KSPPIPEFCGGetNprealloc()`, `KSPPIPEFCGSetMmax()`, `KSPFCGGetMmax()`, `KSPFCGSetMmax()`
433: @*/
434: PetscErrorCode KSPPIPEFCGGetMmax(KSP ksp, PetscInt *mmax)
435: {
436: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
438: PetscFunctionBegin;
440: *mmax = pipefcg->mmax;
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: /*@
445: KSPPIPEFCGSetNprealloc - set the number of directions to preallocate with `KSPPIPEFCG`
447: Logically Collective
449: Input Parameters:
450: + ksp - the Krylov space context
451: - nprealloc - the number of vectors to preallocate
453: Options Database Key:
454: . -ksp_pipefcg_nprealloc N - the number of vectors to preallocate
456: Level: advanced
458: .seealso: [](ch_ksp), `KSPPIPEFCG`, `KSPPIPEFCGSetTruncationType()`, `KSPPIPEFCGGetNprealloc()`, `KSPPIPEFCGSetMmax()`, `KSPPIPEFCGGetMmax()`
459: @*/
460: PetscErrorCode KSPPIPEFCGSetNprealloc(KSP ksp, PetscInt nprealloc)
461: {
462: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
464: PetscFunctionBegin;
467: pipefcg->nprealloc = nprealloc;
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: /*@
472: KSPPIPEFCGGetNprealloc - get the number of directions to preallocate by `KSPPIPEFCG`
474: Not Collective
476: Input Parameter:
477: . ksp - the Krylov space context
479: Output Parameter:
480: . nprealloc - the number of directions preallocated
482: Level: advanced
484: .seealso: [](ch_ksp), `KSPPIPEFCG`, `KSPPIPEFCGGetTruncationType()`, `KSPPIPEFCGSetNprealloc()`, `KSPPIPEFCGSetMmax()`, `KSPPIPEFCGGetMmax()`
485: @*/
486: PetscErrorCode KSPPIPEFCGGetNprealloc(KSP ksp, PetscInt *nprealloc)
487: {
488: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
490: PetscFunctionBegin;
492: *nprealloc = pipefcg->nprealloc;
493: PetscFunctionReturn(PETSC_SUCCESS);
494: }
496: /*@
497: KSPPIPEFCGSetTruncationType - specify how many of its stored previous directions `KSPPIPEFCG` uses during orthogonalization
499: Logically Collective
501: Input Parameters:
502: + ksp - the Krylov space context
503: - truncstrat - the choice of strategy
504: .vb
505: KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to `mmax`) stored directions
506: KSP_FCD_TRUNC_TYPE_NOTAY uses `max(1,mod(i,mmax))` stored directions at iteration i = 0, 1, ...
507: .ve
509: Options Database Key:
510: . -ksp_pipefcg_truncation_type (standard|notay) - which stored search directions to orthogonalize against
512: Level: intermediate
514: .seealso: [](ch_ksp), `KSPPIPEFCG`, `KSPPIPEFCGGetTruncationType`, `KSPFCDTruncationType`, `KSP_FCD_TRUNC_TYPE_STANDARD`, `KSP_FCD_TRUNC_TYPE_NOTAY`
515: @*/
516: PetscErrorCode KSPPIPEFCGSetTruncationType(KSP ksp, KSPFCDTruncationType truncstrat)
517: {
518: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
520: PetscFunctionBegin;
523: pipefcg->truncstrat = truncstrat;
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: /*@
528: KSPPIPEFCGGetTruncationType - get the truncation strategy employed by `KSPPIPEFCG`
530: Not Collective
532: Input Parameter:
533: . ksp - the Krylov space context
535: Output Parameter:
536: . truncstrat - the strategy type
538: Level: intermediate
540: .seealso: [](ch_ksp), `KSPPIPEFCG`, `KSPPIPEFCGSetTruncationType()`, `KSPFCDTruncationType`, `KSP_FCD_TRUNC_TYPE_STANDARD`, `KSP_FCD_TRUNC_TYPE_NOTAY`
541: @*/
542: PetscErrorCode KSPPIPEFCGGetTruncationType(KSP ksp, KSPFCDTruncationType *truncstrat)
543: {
544: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
546: PetscFunctionBegin;
548: *truncstrat = pipefcg->truncstrat;
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: static PetscErrorCode KSPSetFromOptions_PIPEFCG(KSP ksp, PetscOptionItems PetscOptionsObject)
553: {
554: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
555: PetscInt mmax, nprealloc;
556: PetscBool flg;
558: PetscFunctionBegin;
559: PetscOptionsHeadBegin(PetscOptionsObject, "KSP PIPEFCG options");
560: PetscCall(PetscOptionsInt("-ksp_pipefcg_mmax", "Number of search directions to storue", "KSPPIPEFCGSetMmax", pipefcg->mmax, &mmax, &flg));
561: if (flg) PetscCall(KSPPIPEFCGSetMmax(ksp, mmax));
562: PetscCall(PetscOptionsInt("-ksp_pipefcg_nprealloc", "Number of directions to preallocate", "KSPPIPEFCGSetNprealloc", pipefcg->nprealloc, &nprealloc, &flg));
563: if (flg) PetscCall(KSPPIPEFCGSetNprealloc(ksp, nprealloc));
564: PetscCall(PetscOptionsEnum("-ksp_pipefcg_truncation_type", "Truncation approach for directions", "KSPFCGSetTruncationType", KSPFCDTruncationTypes, (PetscEnum)pipefcg->truncstrat, (PetscEnum *)&pipefcg->truncstrat, NULL));
565: PetscOptionsHeadEnd();
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: /*MC
570: KSPPIPEFCG - Implements a Pipelined, Flexible Conjugate Gradient method {cite}`sananschneppmay2016`. [](sec_pipelineksp). [](sec_flexibleksp)
572: Options Database Keys:
573: + -ksp_pipefcg_mmax N - The number of previous search directions to store
574: . -ksp_pipefcg_nprealloc N - The number of previous search directions to preallocate
575: - -ksp_pipefcg_truncation_type (standard|notay) - which stored search directions to orthogonalize against
577: Level: intermediate
579: Notes:
580: Pipelined version of `KSPFCG` that overlaps communication (global reductions) with computation (preconditioner application and matrix-vector products) to reduce the impact of communication latency on parallel performance.
582: Supports left preconditioning only.
584: The natural "norm" for this method is $(u,Au)$, where $u$ is the preconditioned residual. As with standard `KSPCG`, this norm is available at no additional computational cost.
585: Choosing preconditioned or unpreconditioned norms involve an extra blocking global reduction, thus removing any benefit from pipelining.
587: MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
588: See [](doc_faq_pipelined)
590: Contributed by:
591: Patrick Sanan and Sascha M. Schnepp
593: .seealso: [](ch_ksp), [](doc_faq_pipelined), [](sec_pipelineksp), [](sec_flexibleksp), `KSPFCG`, `KSPPIPECG`, `KSPPIPECR`, `KSPGCR`, `KSPPIPEGCR`, `KSPFGMRES`,
594: `KSPCG`, `KSPPIPEFCGSetMmax()`, `KSPPIPEFCGGetMmax()`, `KSPPIPEFCGSetNprealloc()`,
595: `KSPPIPEFCGGetNprealloc()`, `KSPPIPEFCGSetTruncationType()`, `KSPPIPEFCGGetTruncationType()`
596: M*/
597: PETSC_EXTERN PetscErrorCode KSPCreate_PIPEFCG(KSP ksp)
598: {
599: KSP_PIPEFCG *pipefcg;
601: PetscFunctionBegin;
602: PetscCall(PetscNew(&pipefcg));
603: pipefcg->type = !PetscDefined(USE_COMPLEX) ? KSP_CG_SYMMETRIC : KSP_CG_HERMITIAN;
604: pipefcg->mmax = KSPPIPEFCG_DEFAULT_MMAX;
605: pipefcg->nprealloc = KSPPIPEFCG_DEFAULT_NPREALLOC;
606: pipefcg->nvecs = 0;
607: pipefcg->vecb = KSPPIPEFCG_DEFAULT_VECB;
608: pipefcg->nchunks = 0;
609: pipefcg->truncstrat = KSPPIPEFCG_DEFAULT_TRUNCSTRAT;
610: pipefcg->n_restarts = 0;
612: ksp->data = (void *)pipefcg;
614: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
615: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 1));
616: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 1));
617: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
619: ksp->ops->setup = KSPSetUp_PIPEFCG;
620: ksp->ops->solve = KSPSolve_PIPEFCG;
621: ksp->ops->destroy = KSPDestroy_PIPEFCG;
622: ksp->ops->view = KSPView_PIPEFCG;
623: ksp->ops->setfromoptions = KSPSetFromOptions_PIPEFCG;
624: ksp->ops->buildsolution = KSPBuildSolutionDefault;
625: ksp->ops->buildresidual = KSPBuildResidualDefault;
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }