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: }