Actual source code: pipefcg.c
2: #include <../src/ksp/ksp/impls/fcg/pipefcg/pipefcgimpl.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 KSPPIPEFCG_DEFAULT_MMAX 15
19: #define KSPPIPEFCG_DEFAULT_NPREALLOC 5
20: #define KSPPIPEFCG_DEFAULT_VECB 5
21: #define KSPPIPEFCG_DEFAULT_TRUNCSTRAT KSP_FCD_TRUNC_TYPE_NOTAY
23: static PetscErrorCode KSPAllocateVectors_PIPEFCG(KSP ksp, PetscInt nvecsneeded, PetscInt chunksize)
24: {
25: PetscInt i;
26: KSP_PIPEFCG *pipefcg;
27: PetscInt nnewvecs, nvecsprev;
29: pipefcg = (KSP_PIPEFCG *)ksp->data;
31: /* Allocate enough new vectors to add chunksize new vectors, reach nvecsneedtotal, or to reach mmax+1, whichever is smallest */
32: if (pipefcg->nvecs < PetscMin(pipefcg->mmax + 1, nvecsneeded)) {
33: nvecsprev = pipefcg->nvecs;
34: nnewvecs = PetscMin(PetscMax(nvecsneeded - pipefcg->nvecs, chunksize), pipefcg->mmax + 1 - pipefcg->nvecs);
35: KSPCreateVecs(ksp, nnewvecs, &pipefcg->pQvecs[pipefcg->nchunks], 0, NULL);
36: KSPCreateVecs(ksp, nnewvecs, &pipefcg->pZETAvecs[pipefcg->nchunks], 0, NULL);
37: KSPCreateVecs(ksp, nnewvecs, &pipefcg->pPvecs[pipefcg->nchunks], 0, NULL);
38: KSPCreateVecs(ksp, nnewvecs, &pipefcg->pSvecs[pipefcg->nchunks], 0, NULL);
39: pipefcg->nvecs += nnewvecs;
40: for (i = 0; i < nnewvecs; ++i) {
41: pipefcg->Qvecs[nvecsprev + i] = pipefcg->pQvecs[pipefcg->nchunks][i];
42: pipefcg->ZETAvecs[nvecsprev + i] = pipefcg->pZETAvecs[pipefcg->nchunks][i];
43: pipefcg->Pvecs[nvecsprev + i] = pipefcg->pPvecs[pipefcg->nchunks][i];
44: pipefcg->Svecs[nvecsprev + i] = pipefcg->pSvecs[pipefcg->nchunks][i];
45: }
46: pipefcg->chunksizes[pipefcg->nchunks] = nnewvecs;
47: ++pipefcg->nchunks;
48: }
49: return 0;
50: }
52: static PetscErrorCode KSPSetUp_PIPEFCG(KSP ksp)
53: {
54: KSP_PIPEFCG *pipefcg;
55: const PetscInt nworkstd = 5;
57: pipefcg = (KSP_PIPEFCG *)ksp->data;
59: /* Allocate "standard" work vectors (not including the basis and transformed basis vectors) */
60: 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: PetscMalloc4(pipefcg->mmax + 1, &(pipefcg->Pvecs), pipefcg->mmax + 1, &(pipefcg->pPvecs), pipefcg->mmax + 1, &(pipefcg->Svecs), pipefcg->mmax + 1, &(pipefcg->pSvecs));
66: PetscMalloc4(pipefcg->mmax + 1, &(pipefcg->Qvecs), pipefcg->mmax + 1, &(pipefcg->pQvecs), pipefcg->mmax + 1, &(pipefcg->ZETAvecs), pipefcg->mmax + 1, &(pipefcg->pZETAvecs));
67: PetscMalloc4(pipefcg->mmax + 1, &(pipefcg->Pold), pipefcg->mmax + 1, &(pipefcg->Sold), pipefcg->mmax + 1, &(pipefcg->Qold), pipefcg->mmax + 1, &(pipefcg->ZETAold));
68: PetscMalloc1(pipefcg->mmax + 1, &(pipefcg->chunksizes));
69: 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) 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: KSPAllocateVectors_PIPEFCG(ksp, pipefcg->nprealloc, pipefcg->nprealloc);
76: return 0;
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: /* We have not checked these routines for use with complex numbers. The inner products
89: are likely not defined correctly for that case */
92: #define VecXDot(x, y, a) (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x, y, a) : VecTDot(x, y, a))
93: #define VecXDotBegin(x, y, a) (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecDotBegin(x, y, a) : VecTDotBegin(x, y, a))
94: #define VecXDotEnd(x, y, a) (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecDotEnd(x, y, a) : VecTDotEnd(x, y, a))
95: #define VecMXDot(x, n, y, a) (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecMDot(x, n, y, a) : VecMTDot(x, n, y, a))
96: #define VecMXDotBegin(x, n, y, a) (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecMDotBegin(x, n, y, a) : VecMTDotBegin(x, n, y, a))
97: #define VecMXDotEnd(x, n, y, a) (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecMDotEnd(x, n, y, a) : VecMTDotEnd(x, n, y, a))
99: pipefcg = (KSP_PIPEFCG *)ksp->data;
100: X = ksp->vec_sol;
101: B = ksp->vec_rhs;
102: R = ksp->work[0];
103: Z = ksp->work[1];
104: W = ksp->work[2];
105: M = ksp->work[3];
106: N = ksp->work[4];
108: redux = pipefcg->redux;
109: dots = pipefcg->dots;
110: etas = pipefcg->etas;
111: betas = dots; /* dots takes the result of all dot products of which the betas are a subset */
113: PCGetOperators(ksp->pc, &Amat, &Pmat);
115: /* Compute cycle initial residual */
116: KSP_MatMult(ksp, Amat, X, R);
117: VecAYPX(R, -1.0, B); /* r <- b - Ax */
118: KSP_PCApply(ksp, R, Z); /* z <- Br */
120: Pcurr = pipefcg->Pvecs[0];
121: Scurr = pipefcg->Svecs[0];
122: Qcurr = pipefcg->Qvecs[0];
123: ZETAcurr = pipefcg->ZETAvecs[0];
124: VecCopy(Z, Pcurr);
125: KSP_MatMult(ksp, Amat, Pcurr, Scurr); /* S = Ap */
126: VecCopy(Scurr, W); /* w = s = Az */
128: /* Initial state of pipelining intermediates */
129: redux[0] = R;
130: redux[1] = W;
131: VecMXDotBegin(Z, 2, redux, dots);
132: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)Z)); /* perform asynchronous reduction */
133: KSP_PCApply(ksp, W, M); /* m = B(w) */
134: KSP_MatMult(ksp, Amat, M, N); /* n = Am */
135: VecCopy(M, Qcurr); /* q = m */
136: VecCopy(N, ZETAcurr); /* zeta = n */
137: VecMXDotEnd(Z, 2, redux, dots);
138: gamma = dots[0];
139: delta = PetscRealPart(dots[1]);
140: etas[0] = delta;
141: alpha = gamma / delta;
143: i = 0;
144: do {
145: ksp->its++;
147: /* Update X, R, Z, W */
148: VecAXPY(X, +alpha, Pcurr); /* x <- x + alpha * pi */
149: VecAXPY(R, -alpha, Scurr); /* r <- r - alpha * si */
150: VecAXPY(Z, -alpha, Qcurr); /* z <- z - alpha * qi */
151: VecAXPY(W, -alpha, ZETAcurr); /* w <- w - alpha * zetai */
153: /* Compute norm for convergence check */
154: switch (ksp->normtype) {
155: case KSP_NORM_PRECONDITIONED:
156: VecNorm(Z, NORM_2, &dp); /* dp <- sqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */
157: break;
158: case KSP_NORM_UNPRECONDITIONED:
159: VecNorm(R, NORM_2, &dp); /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e) */
160: break;
161: case KSP_NORM_NATURAL:
162: dp = PetscSqrtReal(PetscAbsScalar(gamma)); /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e) */
163: break;
164: case KSP_NORM_NONE:
165: dp = 0.0;
166: break;
167: default:
168: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
169: }
171: /* Check for convergence */
172: ksp->rnorm = dp;
173: KSPLogResidualHistory(ksp, dp);
174: KSPMonitor(ksp, ksp->its, dp);
175: (*ksp->converged)(ksp, ksp->its, dp, &ksp->reason, ksp->cnvP);
176: if (ksp->reason) return 0;
178: /* Computations of current iteration done */
179: ++i;
181: /* If needbe, allocate a new chunk of vectors in P and C */
182: KSPAllocateVectors_PIPEFCG(ksp, i + 1, pipefcg->vecb);
184: /* Note that we wrap around and start clobbering old vectors */
185: idx = i % (pipefcg->mmax + 1);
186: Pcurr = pipefcg->Pvecs[idx];
187: Scurr = pipefcg->Svecs[idx];
188: Qcurr = pipefcg->Qvecs[idx];
189: ZETAcurr = pipefcg->ZETAvecs[idx];
190: eta = pipefcg->etas + idx;
192: /* number of old directions to orthogonalize against */
193: switch (pipefcg->truncstrat) {
194: case KSP_FCD_TRUNC_TYPE_STANDARD:
195: mi = pipefcg->mmax;
196: break;
197: case KSP_FCD_TRUNC_TYPE_NOTAY:
198: mi = ((i - 1) % pipefcg->mmax) + 1;
199: break;
200: default:
201: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized Truncation Strategy");
202: }
204: /* Pick old p,s,q,zeta in a way suitable for VecMDot */
205: VecCopy(Z, Pcurr);
206: for (k = PetscMax(0, i - mi), j = 0; k < i; ++j, ++k) {
207: kdx = k % (pipefcg->mmax + 1);
208: pipefcg->Pold[j] = pipefcg->Pvecs[kdx];
209: pipefcg->Sold[j] = pipefcg->Svecs[kdx];
210: pipefcg->Qold[j] = pipefcg->Qvecs[kdx];
211: pipefcg->ZETAold[j] = pipefcg->ZETAvecs[kdx];
212: redux[j] = pipefcg->Svecs[kdx];
213: }
214: redux[j] = R; /* If the above loop is not executed redux contains only R => all beta_k = 0, only gamma, delta != 0 */
215: redux[j + 1] = W;
217: VecMXDotBegin(Z, j + 2, redux, betas); /* Start split reductions for beta_k = (z,s_k), gamma = (z,r), delta = (z,w) */
218: PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)Z)); /* perform asynchronous reduction */
219: VecWAXPY(N, -1.0, R, W); /* m = u + B(w-r): (a) ntmp = w-r */
220: KSP_PCApply(ksp, N, M); /* m = u + B(w-r): (b) mtmp = B(ntmp) = B(w-r) */
221: VecAXPY(M, 1.0, Z); /* m = u + B(w-r): (c) m = z + mtmp */
222: KSP_MatMult(ksp, Amat, M, N); /* n = Am */
223: VecMXDotEnd(Z, j + 2, redux, betas); /* Finish split reductions */
224: gamma = betas[j];
225: delta = PetscRealPart(betas[j + 1]);
227: *eta = 0.;
228: for (k = PetscMax(0, i - mi), j = 0; k < i; ++j, ++k) {
229: kdx = k % (pipefcg->mmax + 1);
230: betas[j] /= -etas[kdx]; /* betak /= etak */
231: *eta -= ((PetscReal)(PetscAbsScalar(betas[j]) * PetscAbsScalar(betas[j]))) * etas[kdx];
232: /* etaitmp = -betaik^2 * etak */
233: }
234: *eta += delta; /* etai = delta -betaik^2 * etak */
235: if (*eta < 0.) {
236: pipefcg->norm_breakdown = PETSC_TRUE;
237: PetscInfo(ksp, "Restart due to square root breakdown at it = %" PetscInt_FMT "\n", ksp->its);
238: break;
239: } else {
240: alpha = gamma / (*eta); /* alpha = gamma/etai */
241: }
243: /* project out stored search directions using classical G-S */
244: VecCopy(Z, Pcurr);
245: VecCopy(W, Scurr);
246: VecCopy(M, Qcurr);
247: VecCopy(N, ZETAcurr);
248: VecMAXPY(Pcurr, j, betas, pipefcg->Pold); /* pi <- ui - sum_k beta_k p_k */
249: VecMAXPY(Scurr, j, betas, pipefcg->Sold); /* si <- wi - sum_k beta_k s_k */
250: VecMAXPY(Qcurr, j, betas, pipefcg->Qold); /* qi <- m - sum_k beta_k q_k */
251: VecMAXPY(ZETAcurr, j, betas, pipefcg->ZETAold); /* zetai <- n - sum_k beta_k zeta_k */
253: } while (ksp->its < ksp->max_it);
254: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
255: return 0;
256: }
258: static PetscErrorCode KSPSolve_PIPEFCG(KSP ksp)
259: {
260: KSP_PIPEFCG *pipefcg;
261: PetscScalar gamma;
262: PetscReal dp = 0.0;
263: Vec B, R, Z, X;
264: Mat Amat, Pmat;
266: #define VecXDot(x, y, a) (((pipefcg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x, y, a) : VecTDot(x, y, a))
268: PetscCitationsRegister(citation, &cited);
270: pipefcg = (KSP_PIPEFCG *)ksp->data;
271: X = ksp->vec_sol;
272: B = ksp->vec_rhs;
273: R = ksp->work[0];
274: Z = ksp->work[1];
276: PCGetOperators(ksp->pc, &Amat, &Pmat);
278: /* Compute initial residual needed for convergence check*/
279: ksp->its = 0;
280: if (!ksp->guess_zero) {
281: KSP_MatMult(ksp, Amat, X, R);
282: VecAYPX(R, -1.0, B); /* r <- b - Ax */
283: } else {
284: VecCopy(B, R); /* r <- b (x is 0) */
285: }
286: switch (ksp->normtype) {
287: case KSP_NORM_PRECONDITIONED:
288: KSP_PCApply(ksp, R, Z); /* z <- Br */
289: VecNorm(Z, NORM_2, &dp); /* dp <- dqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */
290: break;
291: case KSP_NORM_UNPRECONDITIONED:
292: VecNorm(R, NORM_2, &dp); /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e) */
293: break;
294: case KSP_NORM_NATURAL:
295: KSP_PCApply(ksp, R, Z); /* z <- Br */
296: VecXDot(Z, R, &gamma);
297: dp = PetscSqrtReal(PetscAbsScalar(gamma)); /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e) */
298: break;
299: case KSP_NORM_NONE:
300: dp = 0.0;
301: break;
302: default:
303: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
304: }
306: /* Initial Convergence Check */
307: KSPLogResidualHistory(ksp, dp);
308: KSPMonitor(ksp, 0, dp);
309: ksp->rnorm = dp;
310: (*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP);
311: if (ksp->reason) return 0;
313: do {
314: /* A cycle is broken only if a norm breakdown occurs. If not the entire solve happens in a single cycle.
315: This is coded this way to allow both truncation and truncation-restart strategy
316: (see KSPFCDGetNumOldDirections()) */
317: KSPSolve_PIPEFCG_cycle(ksp);
318: if (ksp->reason) return 0;
319: if (pipefcg->norm_breakdown) {
320: pipefcg->n_restarts++;
321: pipefcg->norm_breakdown = PETSC_FALSE;
322: }
323: } while (ksp->its < ksp->max_it);
325: if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
326: return 0;
327: }
329: static PetscErrorCode KSPDestroy_PIPEFCG(KSP ksp)
330: {
331: PetscInt i;
332: KSP_PIPEFCG *pipefcg;
334: pipefcg = (KSP_PIPEFCG *)ksp->data;
336: /* Destroy "standard" work vecs */
337: VecDestroyVecs(ksp->nwork, &ksp->work);
339: /* Destroy vectors of old directions and the arrays that manage pointers to them */
340: if (pipefcg->nvecs) {
341: for (i = 0; i < pipefcg->nchunks; ++i) {
342: VecDestroyVecs(pipefcg->chunksizes[i], &pipefcg->pPvecs[i]);
343: VecDestroyVecs(pipefcg->chunksizes[i], &pipefcg->pSvecs[i]);
344: VecDestroyVecs(pipefcg->chunksizes[i], &pipefcg->pQvecs[i]);
345: VecDestroyVecs(pipefcg->chunksizes[i], &pipefcg->pZETAvecs[i]);
346: }
347: }
348: PetscFree4(pipefcg->Pvecs, pipefcg->Svecs, pipefcg->pPvecs, pipefcg->pSvecs);
349: PetscFree4(pipefcg->Qvecs, pipefcg->ZETAvecs, pipefcg->pQvecs, pipefcg->pZETAvecs);
350: PetscFree4(pipefcg->Pold, pipefcg->Sold, pipefcg->Qold, pipefcg->ZETAold);
351: PetscFree(pipefcg->chunksizes);
352: PetscFree3(pipefcg->dots, pipefcg->etas, pipefcg->redux);
353: KSPDestroyDefault(ksp);
354: return 0;
355: }
357: static PetscErrorCode KSPView_PIPEFCG(KSP ksp, PetscViewer viewer)
358: {
359: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
360: PetscBool iascii, isstring;
361: const char *truncstr;
363: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
364: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
366: if (pipefcg->truncstrat == KSP_FCD_TRUNC_TYPE_STANDARD) {
367: truncstr = "Using standard truncation strategy";
368: } else if (pipefcg->truncstrat == KSP_FCD_TRUNC_TYPE_NOTAY) {
369: truncstr = "Using Notay's truncation strategy";
370: } else {
371: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Undefined FCD truncation strategy");
372: }
374: if (iascii) {
375: PetscViewerASCIIPrintf(viewer, " max previous directions = %" PetscInt_FMT "\n", pipefcg->mmax);
376: PetscViewerASCIIPrintf(viewer, " preallocated %" PetscInt_FMT " directions\n", PetscMin(pipefcg->nprealloc, pipefcg->mmax + 1));
377: PetscViewerASCIIPrintf(viewer, " %s\n", truncstr);
378: PetscViewerASCIIPrintf(viewer, " restarts performed = %" PetscInt_FMT " \n", pipefcg->n_restarts);
379: } else if (isstring) {
380: PetscViewerStringSPrintf(viewer, "max previous directions = %" PetscInt_FMT ", preallocated %" PetscInt_FMT " directions, %s truncation strategy", pipefcg->mmax, pipefcg->nprealloc, truncstr);
381: }
382: return 0;
383: }
385: /*@
386: KSPPIPEFCGSetMmax - set the maximum number of previous directions `KSPPIPEFCG` will store for orthogonalization
388: Logically Collective
390: Input Parameters:
391: + ksp - the Krylov space context
392: - mmax - the maximum number of previous directions to orthogonalize against
394: Options Database Key:
395: . -ksp_pipefcg_mmax <N> - maximum number of previous directions
397: Level: intermediate
399: Note:
400: mmax + 1 directions are stored (mmax previous ones along with the current one)
401: and whether all are used in each iteration also depends on the truncation strategy, see `KSPPIPEFCGSetTruncationType()`
403: .seealso: [](chapter_ksp), `KSPPIPEFCG`, `KSPPIPEFCGSetTruncationType()`, `KSPPIPEFCGSetNprealloc()`
404: @*/
405: PetscErrorCode KSPPIPEFCGSetMmax(KSP ksp, PetscInt mmax)
406: {
407: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
411: pipefcg->mmax = mmax;
412: return 0;
413: }
415: /*@
416: KSPPIPEFCGGetMmax - get the maximum number of previous directions `KSPPIPEFCG` will store
418: Not Collective
420: Input Parameter:
421: . ksp - the Krylov space context
423: Output Parameter:
424: . mmax - the maximum number of previous directions allowed for orthogonalization
426: Options Database Key:
427: . -ksp_pipefcg_mmax <N> - maximum number of previous directions
429: Level: intermediate
431: .seealso: [](chapter_ksp), `KSPPIPEFCG`, `KSPPIPEFCGGetTruncationType()`, `KSPPIPEFCGGetNprealloc()`, `KSPPIPEFCGSetMmax()`
432: @*/
433: PetscErrorCode KSPPIPEFCGGetMmax(KSP ksp, PetscInt *mmax)
434: {
435: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
438: *mmax = pipefcg->mmax;
439: return 0;
440: }
442: /*@
443: KSPPIPEFCGSetNprealloc - set the number of directions to preallocate with `KSPPIPEFCG`
445: Logically Collective
447: Input Parameters:
448: + ksp - the Krylov space context
449: - nprealloc - the number of vectors to preallocate
451: Options Database Key:
452: . -ksp_pipefcg_nprealloc <N> - the number of vectors to preallocate
454: Level: advanced
456: .seealso: [](chapter_ksp), `KSPPIPEFCG`, `KSPPIPEFCGSetTruncationType()`, `KSPPIPEFCGGetNprealloc()`, `KSPPIPEFCGSetMmax()`, `KSPPIPEFCGGetMmax()`
457: @*/
458: PetscErrorCode KSPPIPEFCGSetNprealloc(KSP ksp, PetscInt nprealloc)
459: {
460: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
464: pipefcg->nprealloc = nprealloc;
465: return 0;
466: }
468: /*@
469: KSPPIPEFCGGetNprealloc - get the number of directions to preallocate by `KSPPIPEFCG`
471: Not Collective
473: Input Parameter:
474: . ksp - the Krylov space context
476: Output Parameter:
477: . nprealloc - the number of directions preallocated
479: Level: advanced
481: .seealso: [](chapter_ksp), `KSPPIPEFCG`, `KSPPIPEFCGGetTruncationType()`, `KSPPIPEFCGSetNprealloc()`, `KSPPIPEFCGSetMmax()`, `KSPPIPEFCGGetMmax()`
482: @*/
483: PetscErrorCode KSPPIPEFCGGetNprealloc(KSP ksp, PetscInt *nprealloc)
484: {
485: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
488: *nprealloc = pipefcg->nprealloc;
489: return 0;
490: }
492: /*@
493: KSPPIPEFCGSetTruncationType - specify how many of its stored previous directions `KSPPIPEFCG` uses during orthoganalization
495: Logically Collective
497: Input Parameters:
498: + ksp - the Krylov space context
499: - truncstrat - the choice of strategy
500: .vb
501: KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
502: KSP_FCD_TRUNC_TYPE_NOTAY uses max(1,mod(i,mmax)) stored directions at iteration i=0,1,..
503: .ve
505: Options Database Key:
506: . -ksp_pipefcg_truncation_type <standard,notay> - which stored search directions to orthogonalize against
508: Level: intermediate
510: .seealso: [](chapter_ksp), `KSPPIPEFCG`, `KSPPIPEFCGGetTruncationType`, `KSPFCDTruncationType`
511: @*/
512: PetscErrorCode KSPPIPEFCGSetTruncationType(KSP ksp, KSPFCDTruncationType truncstrat)
513: {
514: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
518: pipefcg->truncstrat = truncstrat;
519: return 0;
520: }
522: /*@
523: KSPPIPEFCGGetTruncationType - get the truncation strategy employed by `KSPPIPEFCG`
525: Not Collective
527: Input Parameter:
528: . ksp - the Krylov space context
530: Output Parameter:
531: . truncstrat - the strategy type
533: Options Database Key:
534: . -ksp_pipefcg_truncation_type <standard,notay> - which stored basis vectors to orthogonalize against
536: Level: intermediate
538: .seealso: [](chapter_ksp), `KSPPIPEFCG`, `KSPPIPEFCGSetTruncationType`, `KSPFCDTruncationType`
539: @*/
540: PetscErrorCode KSPPIPEFCGGetTruncationType(KSP ksp, KSPFCDTruncationType *truncstrat)
541: {
542: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
545: *truncstrat = pipefcg->truncstrat;
546: return 0;
547: }
549: static PetscErrorCode KSPSetFromOptions_PIPEFCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
550: {
551: KSP_PIPEFCG *pipefcg = (KSP_PIPEFCG *)ksp->data;
552: PetscInt mmax, nprealloc;
553: PetscBool flg;
555: PetscOptionsHeadBegin(PetscOptionsObject, "KSP PIPEFCG options");
556: PetscOptionsInt("-ksp_pipefcg_mmax", "Number of search directions to storue", "KSPPIPEFCGSetMmax", pipefcg->mmax, &mmax, &flg);
557: if (flg) KSPPIPEFCGSetMmax(ksp, mmax);
558: PetscOptionsInt("-ksp_pipefcg_nprealloc", "Number of directions to preallocate", "KSPPIPEFCGSetNprealloc", pipefcg->nprealloc, &nprealloc, &flg);
559: if (flg) KSPPIPEFCGSetNprealloc(ksp, nprealloc);
560: PetscOptionsEnum("-ksp_pipefcg_truncation_type", "Truncation approach for directions", "KSPFCGSetTruncationType", KSPFCDTruncationTypes, (PetscEnum)pipefcg->truncstrat, (PetscEnum *)&pipefcg->truncstrat, NULL);
561: PetscOptionsHeadEnd();
562: return 0;
563: }
565: /*MC
567: KSPPIPEFCG - Implements a Pipelined, Flexible Conjugate Gradient method. [](sec_pipelineksp). [](sec_flexibleksp)
569: Options Database Keys:
570: + -ksp_pipefcg_mmax <N> - The number of previous search directions to store
571: . -ksp_pipefcg_nprealloc <N> - The number of previous search directions to preallocate
572: - -ksp_pipefcg_truncation_type <standard,notay> - which stored search directions to orthogonalize against
574: Notes:
575: Supports left preconditioning only.
577: 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.
578: Choosing preconditioned or unpreconditioned norms involve an extra blocking global reduction, thus removing any benefit from pipelining.
580: MPI configuration may be necessary for reductions to make asynchronous progress, which is important for performance of pipelined methods.
581: See [](doc_faq_pipelined)
583: Contributed by:
584: Patrick Sanan and Sascha M. Schnepp
586: Reference:
587: . * - P. Sanan, S.M. Schnepp, and D.A. May, "Pipelined, Flexible Krylov Subspace Methods", SIAM Journal on Scientific Computing 2016 38:5, C441-C470,
588: DOI: 10.1137/15M1049130
590: Level: intermediate
592: .seealso: [](chapter_ksp), [](doc_faq_pipelined), [](sec_pipelineksp), [](sec_flexibleksp), `KSPFCG`, `KSPPIPECG`, `KSPPIPECR`, `KSPGCR`, `KSPPIPEGCR`, `KSPFGMRES`, `KSPCG`, `KSPPIPEFCGSetMmax()`, `KSPPIPEFCGGetMmax()`, `KSPPIPEFCGSetNprealloc()`,
593: `KSPPIPEFCGGetNprealloc()`, `KSPPIPEFCGSetTruncationType()`, `KSPPIPEFCGGetTruncationType()`
594: M*/
595: PETSC_EXTERN PetscErrorCode KSPCreate_PIPEFCG(KSP ksp)
596: {
597: KSP_PIPEFCG *pipefcg;
599: PetscNew(&pipefcg);
600: #if !defined(PETSC_USE_COMPLEX)
601: pipefcg->type = KSP_CG_SYMMETRIC;
602: #else
603: pipefcg->type = KSP_CG_HERMITIAN;
604: #endif
605: pipefcg->mmax = KSPPIPEFCG_DEFAULT_MMAX;
606: pipefcg->nprealloc = KSPPIPEFCG_DEFAULT_NPREALLOC;
607: pipefcg->nvecs = 0;
608: pipefcg->vecb = KSPPIPEFCG_DEFAULT_VECB;
609: pipefcg->nchunks = 0;
610: pipefcg->truncstrat = KSPPIPEFCG_DEFAULT_TRUNCSTRAT;
611: pipefcg->n_restarts = 0;
613: ksp->data = (void *)pipefcg;
615: KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
616: KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 1);
617: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 1);
618: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);
620: ksp->ops->setup = KSPSetUp_PIPEFCG;
621: ksp->ops->solve = KSPSolve_PIPEFCG;
622: ksp->ops->destroy = KSPDestroy_PIPEFCG;
623: ksp->ops->view = KSPView_PIPEFCG;
624: ksp->ops->setfromoptions = KSPSetFromOptions_PIPEFCG;
625: ksp->ops->buildsolution = KSPBuildSolutionDefault;
626: ksp->ops->buildresidual = KSPBuildResidualDefault;
627: return 0;
628: }