Actual source code: pipelcg.c

  1: #include <petsc/private/kspimpl.h>
  2: #include <petsc/private/vecimpl.h>

  4: #define offset(j)       PetscMax(((j) - (2 * l)), 0)
  5: #define shift(i, j)     ((i)-offset((j)))
  6: #define G(i, j)         (plcg->G[((j) * (2 * l + 1)) + (shift((i), (j)))])
  7: #define G_noshift(i, j) (plcg->G[((j) * (2 * l + 1)) + (i)])
  8: #define alpha(i)        (plcg->alpha[(i)])
  9: #define gamma(i)        (plcg->gamma[(i)])
 10: #define delta(i)        (plcg->delta[(i)])
 11: #define sigma(i)        (plcg->sigma[(i)])
 12: #define req(i)          (plcg->req[(i)])

 14: typedef struct KSP_CG_PIPE_L_s KSP_CG_PIPE_L;
 15: struct KSP_CG_PIPE_L_s {
 16:   PetscInt     l; /* pipeline depth */
 17:   Vec         *Z; /* Z vectors (shifted base) */
 18:   Vec         *U; /* U vectors (unpreconditioned shifted base) */
 19:   Vec         *V; /* V vectors (original base) */
 20:   Vec         *Q; /* Q vectors (auxiliary bases) */
 21:   Vec          p; /* work vector */
 22:   PetscScalar *G; /* such that Z = VG (band matrix)*/
 23:   PetscScalar *gamma, *delta, *alpha;
 24:   PetscReal    lmin, lmax; /* min and max eigen values estimates to compute base shifts */
 25:   PetscReal   *sigma;      /* base shifts */
 26:   MPI_Request *req;        /* request array for asynchronous global collective */
 27:   PetscBool    show_rstrt; /* flag to show restart information in output (default: not shown) */
 28: };

 30: /*
 31:   KSPSetUp_PIPELCG - Sets up the workspace needed by the PIPELCG method.

 33:   This is called once, usually automatically by KSPSolve() or KSPSetUp()
 34:   but can be called directly by KSPSetUp()
 35: */
 36: static PetscErrorCode KSPSetUp_PIPELCG(KSP ksp)
 37: {
 38:   KSP_CG_PIPE_L *plcg = (KSP_CG_PIPE_L *)ksp->data;
 39:   PetscInt       l = plcg->l, max_it = ksp->max_it;
 40:   MPI_Comm       comm;

 42:   comm = PetscObjectComm((PetscObject)ksp);

 47:   KSPSetWorkVecs(ksp, 1); /* get work vectors needed by PIPELCG */
 48:   plcg->p = ksp->work[0];

 50:   VecDuplicateVecs(plcg->p, PetscMax(3, l + 1), &plcg->Z);
 51:   VecDuplicateVecs(plcg->p, 3, &plcg->U);
 52:   VecDuplicateVecs(plcg->p, 3, &plcg->V);
 53:   VecDuplicateVecs(plcg->p, 3 * (l - 1) + 1, &plcg->Q);
 54:   PetscCalloc1(2, &plcg->alpha);
 55:   PetscCalloc1(l, &plcg->sigma);

 57:   return 0;
 58: }

 60: static PetscErrorCode KSPReset_PIPELCG(KSP ksp)
 61: {
 62:   KSP_CG_PIPE_L *plcg = (KSP_CG_PIPE_L *)ksp->data;
 63:   PetscInt       l    = plcg->l;

 65:   PetscFree(plcg->sigma);
 66:   PetscFree(plcg->alpha);
 67:   VecDestroyVecs(PetscMax(3, l + 1), &plcg->Z);
 68:   VecDestroyVecs(3, &plcg->U);
 69:   VecDestroyVecs(3, &plcg->V);
 70:   VecDestroyVecs(3 * (l - 1) + 1, &plcg->Q);
 71:   return 0;
 72: }

 74: static PetscErrorCode KSPDestroy_PIPELCG(KSP ksp)
 75: {
 76:   KSPReset_PIPELCG(ksp);
 77:   KSPDestroyDefault(ksp);
 78:   return 0;
 79: }

 81: static PetscErrorCode KSPSetFromOptions_PIPELCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
 82: {
 83:   KSP_CG_PIPE_L *plcg = (KSP_CG_PIPE_L *)ksp->data;
 84:   PetscBool      flag = PETSC_FALSE;

 86:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP PIPELCG options");
 87:   PetscOptionsInt("-ksp_pipelcg_pipel", "Pipeline length", "", plcg->l, &plcg->l, &flag);
 88:   if (!flag) plcg->l = 1;
 89:   PetscOptionsReal("-ksp_pipelcg_lmin", "Estimate for smallest eigenvalue", "", plcg->lmin, &plcg->lmin, &flag);
 90:   if (!flag) plcg->lmin = 0.0;
 91:   PetscOptionsReal("-ksp_pipelcg_lmax", "Estimate for largest eigenvalue", "", plcg->lmax, &plcg->lmax, &flag);
 92:   if (!flag) plcg->lmax = 0.0;
 93:   PetscOptionsBool("-ksp_pipelcg_monitor", "Output information on restarts when they occur? (default: 0)", "", plcg->show_rstrt, &plcg->show_rstrt, &flag);
 94:   if (!flag) plcg->show_rstrt = PETSC_FALSE;
 95:   PetscOptionsHeadEnd();
 96:   return 0;
 97: }

 99: static PetscErrorCode MPIPetsc_Iallreduce(void *sendbuf, void *recvbuf, PetscMPIInt count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm, MPI_Request *request)
100: {
101: #if defined(PETSC_HAVE_MPI_NONBLOCKING_COLLECTIVES)
102:   MPI_Iallreduce(sendbuf, recvbuf, count, datatype, op, comm, request);
103: #else
104:   MPIU_Allreduce(sendbuf, recvbuf, count, datatype, op, comm);
105:   *request = MPI_REQUEST_NULL;
106: #endif
107:   return 0;
108: }

110: static PetscErrorCode KSPView_PIPELCG(KSP ksp, PetscViewer viewer)
111: {
112:   KSP_CG_PIPE_L *plcg   = (KSP_CG_PIPE_L *)ksp->data;
113:   PetscBool      iascii = PETSC_FALSE, isstring = PETSC_FALSE;

115:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
116:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
117:   if (iascii) {
118:     PetscViewerASCIIPrintf(viewer, "  Pipeline depth: %" PetscInt_FMT "\n", plcg->l);
119:     PetscViewerASCIIPrintf(viewer, "  Minimal eigenvalue estimate %g\n", (double)plcg->lmin);
120:     PetscViewerASCIIPrintf(viewer, "  Maximal eigenvalue estimate %g\n", (double)plcg->lmax);
121:   } else if (isstring) {
122:     PetscViewerStringSPrintf(viewer, "  Pipeline depth: %" PetscInt_FMT "\n", plcg->l);
123:     PetscViewerStringSPrintf(viewer, "  Minimal eigenvalue estimate %g\n", (double)plcg->lmin);
124:     PetscViewerStringSPrintf(viewer, "  Maximal eigenvalue estimate %g\n", (double)plcg->lmax);
125:   }
126:   return 0;
127: }

129: static PetscErrorCode KSPSolve_InnerLoop_PIPELCG(KSP ksp)
130: {
131:   KSP_CG_PIPE_L *plcg = (KSP_CG_PIPE_L *)ksp->data;
132:   Mat            A = NULL, Pmat = NULL;
133:   PetscInt       it = 0, max_it = ksp->max_it, l = plcg->l, i = 0, j = 0, k = 0;
134:   PetscInt       start = 0, middle = 0, end = 0;
135:   Vec           *Z = plcg->Z, *U = plcg->U, *V = plcg->V, *Q = plcg->Q;
136:   Vec            x = NULL, p = NULL, temp = NULL;
137:   PetscScalar    sum_dummy = 0.0, eta = 0.0, zeta = 0.0, lambda = 0.0;
138:   PetscReal      dp = 0.0, tmp = 0.0, beta = 0.0, invbeta2 = 0.0;
139:   MPI_Comm       comm;

141:   x = ksp->vec_sol;
142:   p = plcg->p;

144:   comm = PetscObjectComm((PetscObject)ksp);
145:   PCGetOperators(ksp->pc, &A, &Pmat);

147:   for (it = 0; it < max_it + l; ++it) {
148:     /* ----------------------------------- */
149:     /* Multiplication  z_{it+1} =  Az_{it} */
150:     /* ----------------------------------- */
151:     /* Shift the U vector pointers */
152:     temp = U[2];
153:     for (i = 2; i > 0; i--) U[i] = U[i - 1];
154:     U[0] = temp;
155:     if (it < l) {
156:       /* SpMV and Sigma-shift and Prec */
157:       MatMult(A, Z[l - it], U[0]);
158:       VecAXPY(U[0], -sigma(it), U[1]);
159:       KSP_PCApply(ksp, U[0], Z[l - it - 1]);
160:       if (it < l - 1) VecCopy(Z[l - it - 1], Q[3 * it]);
161:     } else {
162:       /* Shift the Z vector pointers */
163:       temp = Z[PetscMax(l, 2)];
164:       for (i = PetscMax(l, 2); i > 0; --i) Z[i] = Z[i - 1];
165:       Z[0] = temp;
166:       /* SpMV and Prec */
167:       MatMult(A, Z[1], U[0]);
168:       KSP_PCApply(ksp, U[0], Z[0]);
169:     }

171:     /* ----------------------------------- */
172:     /* Adjust the G matrix */
173:     /* ----------------------------------- */
174:     if (it >= l) {
175:       if (it == l) {
176:         /* MPI_Wait for G(0,0),scale V0 and Z and U and Q vectors with 1/beta */
177:         MPI_Wait(&req(0), MPI_STATUS_IGNORE);
178:         beta    = PetscSqrtReal(PetscRealPart(G(0, 0)));
179:         G(0, 0) = 1.0;
180:         VecAXPY(V[0], 1.0 / beta, p); /* this assumes V[0] to be zero initially */
181:         for (j = 0; j <= PetscMax(l, 2); ++j) VecScale(Z[j], 1.0 / beta);
182:         for (j = 0; j <= 2; ++j) VecScale(U[j], 1.0 / beta);
183:         for (j = 0; j < l - 1; ++j) VecScale(Q[3 * j], 1.0 / beta);
184:       }

186:       /* MPI_Wait until the dot products,started l iterations ago,are completed */
187:       MPI_Wait(&req(it - l + 1), MPI_STATUS_IGNORE);
188:       if (it >= 2 * l) {
189:         for (j = PetscMax(0, it - 3 * l + 1); j <= it - 2 * l; j++) { G(j, it - l + 1) = G(it - 2 * l + 1, j + l); /* exploit symmetry in G matrix */ }
190:       }

192:       if (it <= 2 * l - 1) {
193:         invbeta2 = 1.0 / (beta * beta);
194:         /* Scale columns 1 up to l of G with 1/beta^2 */
195:         for (j = PetscMax(it - 3 * l + 1, 0); j <= it - l + 1; ++j) G(j, it - l + 1) *= invbeta2;
196:       }

198:       for (j = PetscMax(it - 2 * l + 2, 0); j <= it - l; ++j) {
199:         sum_dummy = 0.0;
200:         for (k = PetscMax(it - 3 * l + 1, 0); k <= j - 1; ++k) sum_dummy = sum_dummy + G(k, j) * G(k, it - l + 1);
201:         G(j, it - l + 1) = (G(j, it - l + 1) - sum_dummy) / G(j, j);
202:       }

204:       sum_dummy = 0.0;
205:       for (k = PetscMax(it - 3 * l + 1, 0); k <= it - l; ++k) sum_dummy = sum_dummy + G(k, it - l + 1) * G(k, it - l + 1);

207:       tmp = PetscRealPart(G(it - l + 1, it - l + 1) - sum_dummy);
208:       /* Breakdown check */
209:       if (tmp < 0) {
210:         if (plcg->show_rstrt) PetscPrintf(comm, "Sqrt breakdown in iteration %" PetscInt_FMT ": sqrt argument is %e. Iteration was restarted.\n", ksp->its + 1, (double)tmp);
211:         /* End hanging dot-products in the pipeline before exiting for-loop */
212:         start = it - l + 2;
213:         end   = PetscMin(it + 1, max_it + 1); /* !warning! 'it' can actually be greater than 'max_it' */
214:         for (i = start; i < end; ++i) MPI_Wait(&req(i), MPI_STATUS_IGNORE);
215:         break;
216:       }
217:       G(it - l + 1, it - l + 1) = PetscSqrtReal(tmp);

219:       if (it < 2 * l) {
220:         if (it == l) {
221:           gamma(it - l) = (G(it - l, it - l + 1) + sigma(it - l) * G(it - l, it - l)) / G(it - l, it - l);
222:         } else {
223:           gamma(it - l) = (G(it - l, it - l + 1) + sigma(it - l) * G(it - l, it - l) - delta(it - l - 1) * G(it - l - 1, it - l)) / G(it - l, it - l);
224:         }
225:         delta(it - l) = G(it - l + 1, it - l + 1) / G(it - l, it - l);
226:       } else {
227:         gamma(it - l) = (G(it - l, it - l) * gamma(it - 2 * l) + G(it - l, it - l + 1) * delta(it - 2 * l) - G(it - l - 1, it - l) * delta(it - l - 1)) / G(it - l, it - l);
228:         delta(it - l) = (G(it - l + 1, it - l + 1) * delta(it - 2 * l)) / G(it - l, it - l);
229:       }

231:       /* -------------------------------------------------- */
232:       /* Recursively compute the next V, Q, Z and U vectors */
233:       /* -------------------------------------------------- */
234:       /* Shift the V vector pointers */
235:       temp = V[2];
236:       for (i = 2; i > 0; i--) V[i] = V[i - 1];
237:       V[0] = temp;

239:       /* Recurrence V vectors */
240:       if (l == 1) {
241:         VecCopy(Z[1], V[0]);
242:       } else {
243:         VecCopy(Q[0], V[0]);
244:       }
245:       if (it == l) {
246:         VecAXPY(V[0], sigma(0) - gamma(it - l), V[1]);
247:       } else {
248:         alpha(0) = sigma(0) - gamma(it - l);
249:         alpha(1) = -delta(it - l - 1);
250:         VecMAXPY(V[0], 2, &alpha(0), &V[1]);
251:       }
252:       VecScale(V[0], 1.0 / delta(it - l));

254:       /* Recurrence Q vectors */
255:       for (j = 0; j < l - 1; ++j) {
256:         /* Shift the Q vector pointers */
257:         temp = Q[3 * j + 2];
258:         for (i = 2; i > 0; i--) Q[3 * j + i] = Q[3 * j + i - 1];
259:         Q[3 * j] = temp;

261:         if (j < l - 2) {
262:           VecCopy(Q[3 * (j + 1)], Q[3 * j]);
263:         } else {
264:           VecCopy(Z[1], Q[3 * j]);
265:         }
266:         if (it == l) {
267:           VecAXPY(Q[3 * j], sigma(j + 1) - gamma(it - l), Q[3 * j + 1]);
268:         } else {
269:           alpha(0) = sigma(j + 1) - gamma(it - l);
270:           alpha(1) = -delta(it - l - 1);
271:           VecMAXPY(Q[3 * j], 2, &alpha(0), &Q[3 * j + 1]);
272:         }
273:         VecScale(Q[3 * j], 1.0 / delta(it - l));
274:       }

276:       /* Recurrence Z and U vectors */
277:       if (it == l) {
278:         VecAXPY(Z[0], -gamma(it - l), Z[1]);
279:         VecAXPY(U[0], -gamma(it - l), U[1]);
280:       } else {
281:         alpha(0) = -gamma(it - l);
282:         alpha(1) = -delta(it - l - 1);
283:         VecMAXPY(Z[0], 2, &alpha(0), &Z[1]);
284:         VecMAXPY(U[0], 2, &alpha(0), &U[1]);
285:       }
286:       VecScale(Z[0], 1.0 / delta(it - l));
287:       VecScale(U[0], 1.0 / delta(it - l));
288:     }

290:     /* ---------------------------------------- */
291:     /* Compute and communicate the dot products */
292:     /* ---------------------------------------- */
293:     if (it < l) {
294:       for (j = 0; j < it + 2; ++j) { (*U[0]->ops->dot_local)(U[0], Z[l - j], &G(j, it + 1)); /* dot-products (U[0],Z[j]) */ }
295:       MPIPetsc_Iallreduce(MPI_IN_PLACE, &G(0, it + 1), it + 2, MPIU_SCALAR, MPIU_SUM, comm, &req(it + 1));
296:     } else if ((it >= l) && (it < max_it)) {
297:       middle = it - l + 2;
298:       end    = it + 2;
299:       (*U[0]->ops->dot_local)(U[0], V[0], &G(it - l + 1, it + 1)); /* dot-product (U[0],V[0]) */
300:       for (j = middle; j < end; ++j) { (*U[0]->ops->dot_local)(U[0], plcg->Z[it + 1 - j], &G(j, it + 1)); /* dot-products (U[0],Z[j]) */ }
301:       MPIPetsc_Iallreduce(MPI_IN_PLACE, &G(it - l + 1, it + 1), l + 1, MPIU_SCALAR, MPIU_SUM, comm, &req(it + 1));
302:     }

304:     /* ----------------------------------------- */
305:     /* Compute solution vector and residual norm */
306:     /* ----------------------------------------- */
307:     if (it >= l) {
308:       if (it == l) {
309:         if (ksp->its != 0) ++ksp->its;
310:         eta  = gamma(0);
311:         zeta = beta;
312:         VecCopy(V[1], p);
313:         VecScale(p, 1.0 / eta);
314:         VecAXPY(x, zeta, p);
315:         dp = beta;
316:       } else if (it > l) {
317:         k = it - l;
318:         ++ksp->its;
319:         lambda = delta(k - 1) / eta;
320:         eta    = gamma(k) - lambda * delta(k - 1);
321:         zeta   = -lambda * zeta;
322:         VecScale(p, -delta(k - 1) / eta);
323:         VecAXPY(p, 1.0 / eta, V[1]);
324:         VecAXPY(x, zeta, p);
325:         dp = PetscAbsScalar(zeta);
326:       }
327:       ksp->rnorm = dp;
328:       KSPLogResidualHistory(ksp, dp);
329:       KSPMonitor(ksp, ksp->its, dp);
330:       (*ksp->converged)(ksp, ksp->its, dp, &ksp->reason, ksp->cnvP);

332:       if (ksp->its >= max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
333:       if (ksp->reason) {
334:         /* End hanging dot-products in the pipeline before exiting for-loop */
335:         start = it - l + 2;
336:         end   = PetscMin(it + 2, max_it + 1); /* !warning! 'it' can actually be greater than 'max_it' */
337:         for (i = start; i < end; ++i) MPI_Wait(&req(i), MPI_STATUS_IGNORE);
338:         break;
339:       }
340:     }
341:   } /* End inner for loop */
342:   return 0;
343: }

345: static PetscErrorCode KSPSolve_ReInitData_PIPELCG(KSP ksp)
346: {
347:   KSP_CG_PIPE_L *plcg = (KSP_CG_PIPE_L *)ksp->data;
348:   PetscInt       i = 0, j = 0, l = plcg->l, max_it = ksp->max_it;

350:   for (i = 0; i < PetscMax(3, l + 1); ++i) VecSet(plcg->Z[i], 0.0);
351:   for (i = 1; i < 3; ++i) VecSet(plcg->U[i], 0.0);
352:   for (i = 0; i < 3; ++i) VecSet(plcg->V[i], 0.0);
353:   for (i = 0; i < 3 * (l - 1) + 1; ++i) VecSet(plcg->Q[i], 0.0);
354:   for (j = 0; j < (max_it + 1); ++j) {
355:     gamma(j) = 0.0;
356:     delta(j) = 0.0;
357:     for (i = 0; i < (2 * l + 1); ++i) G_noshift(i, j) = 0.0;
358:   }
359:   return 0;
360: }

362: /*
363:   KSPSolve_PIPELCG - This routine actually applies the pipelined(l) conjugate gradient method
364: */
365: static PetscErrorCode KSPSolve_PIPELCG(KSP ksp)
366: {
367:   KSP_CG_PIPE_L *plcg = (KSP_CG_PIPE_L *)ksp->data;
368:   Mat            A = NULL, Pmat = NULL;
369:   Vec            b = NULL, x = NULL, p = NULL;
370:   PetscInt       max_it = ksp->max_it, l = plcg->l;
371:   PetscInt       i = 0, outer_it = 0, curr_guess_zero = 0;
372:   PetscReal      lmin = plcg->lmin, lmax = plcg->lmax;
373:   PetscBool      diagonalscale = PETSC_FALSE;
374:   MPI_Comm       comm;

376:   comm = PetscObjectComm((PetscObject)ksp);
377:   PCGetDiagonalScale(ksp->pc, &diagonalscale);

380:   x = ksp->vec_sol;
381:   b = ksp->vec_rhs;
382:   p = plcg->p;

384:   PetscCalloc1((max_it + 1) * (2 * l + 1), &plcg->G);
385:   PetscCalloc1(max_it + 1, &plcg->gamma);
386:   PetscCalloc1(max_it + 1, &plcg->delta);
387:   PetscCalloc1(max_it + 1, &plcg->req);

389:   PCGetOperators(ksp->pc, &A, &Pmat);

391:   for (i = 0; i < l; ++i) sigma(i) = (0.5 * (lmin + lmax) + (0.5 * (lmax - lmin) * PetscCosReal(PETSC_PI * (2.0 * i + 1.0) / (2.0 * l))));

393:   ksp->its        = 0;
394:   outer_it        = 0;
395:   curr_guess_zero = !!ksp->guess_zero;

397:   while (ksp->its < max_it) { /* OUTER LOOP (gmres-like restart to handle breakdowns) */
398:     /* RESTART LOOP */
399:     if (!curr_guess_zero) {
400:       KSP_MatMult(ksp, A, x, plcg->U[0]); /* u <- b - Ax */
401:       VecAYPX(plcg->U[0], -1.0, b);
402:     } else {
403:       VecCopy(b, plcg->U[0]); /* u <- b (x is 0) */
404:     }
405:     KSP_PCApply(ksp, plcg->U[0], p); /* p <- Bu */

407:     if (outer_it > 0) {
408:       /* Re-initialize Z,U,V,Q,gamma,delta,G after restart occurred */
409:       KSPSolve_ReInitData_PIPELCG(ksp);
410:     }

412:     (*plcg->U[0]->ops->dot_local)(plcg->U[0], p, &G(0, 0));
413:     MPIPetsc_Iallreduce(MPI_IN_PLACE, &G(0, 0), 1, MPIU_SCALAR, MPIU_SUM, comm, &req(0));
414:     VecCopy(p, plcg->Z[l]);

416:     KSPSolve_InnerLoop_PIPELCG(ksp);

418:     if (ksp->reason) break; /* convergence or divergence */
419:     ++outer_it;
420:     curr_guess_zero = 0;
421:   }

423:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
424:   PetscFree(plcg->G);
425:   PetscFree(plcg->gamma);
426:   PetscFree(plcg->delta);
427:   PetscFree(plcg->req);
428:   return 0;
429: }

431: /*MC
432:     KSPPIPELCG - Deep pipelined (length l) Conjugate Gradient method. This method has only a single non-blocking global
433:     reduction per iteration, compared to 2 blocking reductions for standard `KSPCG`. The reduction is overlapped by the
434:     matrix-vector product and preconditioner application of the next l iterations. The pipeline length l is a parameter
435:     of the method. [](sec_pipelineksp)

437:     Options Database Keys:
438: +   -ksp_pipelcg_pipel - pipelined length
439: .   -ksp_pipelcg_lmin - approximation to the smallest eigenvalue of the preconditioned operator (default: 0.0)
440: .   -ksp_pipelcg_lmax - approximation to the largest eigenvalue of the preconditioned operator (default: 0.0)
441: -   -ksp_pipelcg_monitor - output where/why the method restarts when a sqrt breakdown occurs

443:     Level: advanced

445:     Notes:
446:     MPI configuration may be necessary for reductions to make asynchronous progress, which is important for
447:     performance of pipelined methods. See [](doc_faq_pipelined)

449:     Contributed by:
450:     Siegfried Cools, University of Antwerp, Dept. Mathematics and Computer Science,
451:     funded by Flemish Research Foundation (FWO) grant number 12H4617N.

453:     Example usage:
454: .vb
455:     KSP tutorials ex2, no preconditioner, pipel = 2, lmin = 0.0, lmax = 8.0 :
456:         $mpiexec -n 14 ./ex2 -m 1000 -n 1000 -ksp_type pipelcg -pc_type none -ksp_norm_type natural
457:            -ksp_rtol 1e-10 -ksp_max_it 1000 -ksp_pipelcg_pipel 2 -ksp_pipelcg_lmin 0.0 -ksp_pipelcg_lmax 8.0 -log_view
458:     SNES tutorials ex48, bjacobi preconditioner, pipel = 3, lmin = 0.0, lmax = 2.0, show restart information :
459:         $mpiexec -n 14 ./ex48 -M 150 -P 100 -ksp_type pipelcg -pc_type bjacobi -ksp_rtol 1e-10 -ksp_pipelcg_pipel 3
460:            -ksp_pipelcg_lmin 0.0 -ksp_pipelcg_lmax 2.0 -ksp_pipelcg_monitor -log_view
461: .ve

463:     References:
464: +   * - J. Cornelis, S. Cools and W. Vanroose,
465:         "The Communication-Hiding Conjugate Gradient Method with Deep Pipelines"
466:         Submitted to SIAM Journal on Scientific Computing (SISC), 2018.
467: -   * - S. Cools, J. Cornelis and W. Vanroose,
468:         "Numerically Stable Recurrence Relations for the Communication Hiding Pipelined Conjugate Gradient Method"
469:         Submitted to IEEE Transactions on Parallel and Distributed Systems, 2019.

471: .seealso: [](chapter_ksp), [](sec_pipelineksp), [](doc_faq_pipelined), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSPCG`, `KSPPIPECG`, `KSPPIPECGRR`, `KSPPGMRES`,
472:           `KSPPIPEBCGS`, `KSPSetPCSide()`, `KSPGROPPCG`
473: M*/
474: PETSC_EXTERN PetscErrorCode KSPCreate_PIPELCG(KSP ksp)
475: {
476:   KSP_CG_PIPE_L *plcg = NULL;

478:   PetscNew(&plcg);
479:   ksp->data = (void *)plcg;

481:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);
482:   KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2);

484:   ksp->ops->setup          = KSPSetUp_PIPELCG;
485:   ksp->ops->solve          = KSPSolve_PIPELCG;
486:   ksp->ops->reset          = KSPReset_PIPELCG;
487:   ksp->ops->destroy        = KSPDestroy_PIPELCG;
488:   ksp->ops->view           = KSPView_PIPELCG;
489:   ksp->ops->setfromoptions = KSPSetFromOptions_PIPELCG;
490:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
491:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
492:   return 0;
493: }