Actual source code: pipebcgs.c


  2: /*
  3:     This file implements pipelined BiCGStab (pipe-BiCGStab).
  4: */
  5: #include <../src/ksp/ksp/impls/bcgs/bcgsimpl.h>

  7: static PetscErrorCode KSPSetUp_PIPEBCGS(KSP ksp)
  8: {
  9:   KSPSetWorkVecs(ksp, 15);
 10:   return 0;
 11: }

 13: /* Only need a few hacks from KSPSolve_BCGS */
 14: #include <petsc/private/pcimpl.h>
 15: static PetscErrorCode KSPSolve_PIPEBCGS(KSP ksp)
 16: {
 17:   PetscInt    i;
 18:   PetscScalar rho, rhoold, alpha, beta, omega = 0.0, d1, d2, d3;
 19:   Vec         X, B, S, R, RP, Y, Q, P2, Q2, R2, S2, W, Z, W2, Z2, T, V;
 20:   PetscReal   dp   = 0.0;
 21:   KSP_BCGS   *bcgs = (KSP_BCGS *)ksp->data;
 22:   PC          pc;

 24:   X  = ksp->vec_sol;
 25:   B  = ksp->vec_rhs;
 26:   R  = ksp->work[0];
 27:   RP = ksp->work[1];
 28:   S  = ksp->work[2];
 29:   Y  = ksp->work[3];
 30:   Q  = ksp->work[4];
 31:   Q2 = ksp->work[5];
 32:   P2 = ksp->work[6];
 33:   R2 = ksp->work[7];
 34:   S2 = ksp->work[8];
 35:   W  = ksp->work[9];
 36:   Z  = ksp->work[10];
 37:   W2 = ksp->work[11];
 38:   Z2 = ksp->work[12];
 39:   T  = ksp->work[13];
 40:   V  = ksp->work[14];

 42:   if (!ksp->guess_zero) {
 43:     if (!bcgs->guess) VecDuplicate(X, &bcgs->guess);
 44:     VecCopy(X, bcgs->guess);
 45:   } else {
 46:     VecSet(X, 0.0);
 47:   }

 49:   /* Compute initial residual */
 50:   KSPGetPC(ksp, &pc);
 51:   PCSetUp(pc);
 52:   if (!ksp->guess_zero) {
 53:     KSP_MatMult(ksp, pc->mat, X, Q2);
 54:     VecCopy(B, R);
 55:     VecAXPY(R, -1.0, Q2);
 56:   } else {
 57:     VecCopy(B, R);
 58:   }

 60:   /* Test for nothing to do */
 61:   if (ksp->normtype != KSP_NORM_NONE) {
 62:     VecNorm(R, NORM_2, &dp);
 63:   } else dp = 0.0;
 64:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
 65:   ksp->its   = 0;
 66:   ksp->rnorm = dp;
 67:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
 68:   KSPLogResidualHistory(ksp, dp);
 69:   KSPMonitor(ksp, 0, dp);
 70:   (*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP);
 71:   if (ksp->reason) return 0;

 73:   /* Initialize */
 74:   VecCopy(R, RP); /* rp <- r */

 76:   VecDotBegin(R, RP, &rho); /* rho <- (r,rp) */
 77:   PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R));
 78:   KSP_PCApply(ksp, R, R2);          /* r2 <- K r */
 79:   KSP_MatMult(ksp, pc->mat, R2, W); /* w <- A r2 */
 80:   VecDotEnd(R, RP, &rho);

 82:   VecDotBegin(W, RP, &d2); /* d2 <- (w,rp) */
 83:   PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)W));
 84:   KSP_PCApply(ksp, W, W2);          /* w2 <- K w */
 85:   KSP_MatMult(ksp, pc->mat, W2, T); /* t <- A w2 */
 86:   VecDotEnd(W, RP, &d2);

 88:   alpha = rho / d2;
 89:   beta  = 0.0;

 91:   /* Main loop */
 92:   i = 0;
 93:   do {
 94:     if (i == 0) {
 95:       VecCopy(R2, P2); /* p2 <- r2 */
 96:       VecCopy(W, S);   /* s  <- w  */
 97:       VecCopy(W2, S2); /* s2 <- w2 */
 98:       VecCopy(T, Z);   /* z  <- t  */
 99:     } else {
100:       VecAXPBYPCZ(P2, 1.0, -beta * omega, beta, R2, S2); /* p2 <- beta * p2 + r2 - beta * omega * s2 */
101:       VecAXPBYPCZ(S, 1.0, -beta * omega, beta, W, Z);    /* s  <- beta * s  + w  - beta * omega * z  */
102:       VecAXPBYPCZ(S2, 1.0, -beta * omega, beta, W2, Z2); /* s2 <- beta * s2 + w2 - beta * omega * z2 */
103:       VecAXPBYPCZ(Z, 1.0, -beta * omega, beta, T, V);    /* z  <- beta * z  + t  - beta * omega * v  */
104:     }
105:     VecWAXPY(Q, -alpha, S, R);    /* q  <- r  - alpha s  */
106:     VecWAXPY(Q2, -alpha, S2, R2); /* q2 <- r2 - alpha s2 */
107:     VecWAXPY(Y, -alpha, Z, W);    /* y  <- w  - alpha z  */

109:     VecDotBegin(Q, Y, &d1); /* d1 <- (q,y) */
110:     VecDotBegin(Y, Y, &d2); /* d2 <- (y,y) */

112:     PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)Q));
113:     KSP_PCApply(ksp, Z, Z2);          /* z2 <- K z */
114:     KSP_MatMult(ksp, pc->mat, Z2, V); /* v <- A z2 */

116:     VecDotEnd(Q, Y, &d1);
117:     VecDotEnd(Y, Y, &d2);

119:     if (d2 == 0.0) {
120:       /* y is 0. if q is 0, then alpha s == r, and hence alpha p may be our solution. Give it a try? */
121:       VecDot(Q, Q, &d1);
122:       if (d1 != 0.0) {
123:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
124:         break;
125:       }
126:       VecAXPY(X, alpha, P2); /* x <- x + alpha p2 */
127:       PetscObjectSAWsTakeAccess((PetscObject)ksp);
128:       ksp->its++;
129:       ksp->rnorm  = 0.0;
130:       ksp->reason = KSP_CONVERGED_RTOL;
131:       PetscObjectSAWsGrantAccess((PetscObject)ksp);
132:       KSPLogResidualHistory(ksp, dp);
133:       KSPMonitor(ksp, i + 1, 0.0);
134:       break;
135:     }
136:     omega = d1 / d2;                                      /* omega <- (y'q) / (y'y) */
137:     VecAXPBYPCZ(X, alpha, omega, 1.0, P2, Q2); /* x <- alpha * p2 + omega * q2 + x */
138:     VecWAXPY(R, -omega, Y, Q);                 /* r <- q - omega y */
139:     VecWAXPY(R2, -alpha, Z2, W2);              /* r2 <- w2 - alpha z2 */
140:     VecAYPX(R2, -omega, Q2);                   /* r2 <- q2 - omega r2 */
141:     VecWAXPY(W, -alpha, V, T);                 /* w <- t - alpha v */
142:     VecAYPX(W, -omega, Y);                     /* w <- y - omega w */
143:     rhoold = rho;

145:     if (ksp->normtype != KSP_NORM_NONE && ksp->chknorm < i + 2) { VecNormBegin(R, NORM_2, &dp); /* dp <- norm(r) */ }
146:     VecDotBegin(R, RP, &rho); /* rho <- (r,rp) */
147:     VecDotBegin(S, RP, &d1);  /* d1 <- (s,rp) */
148:     VecDotBegin(W, RP, &d2);  /* d2 <- (w,rp) */
149:     VecDotBegin(Z, RP, &d3);  /* d3 <- (z,rp) */

151:     PetscCommSplitReductionBegin(PetscObjectComm((PetscObject)R));
152:     KSP_PCApply(ksp, W, W2);          /* w2 <- K w */
153:     KSP_MatMult(ksp, pc->mat, W2, T); /* t <- A w2 */

155:     if (ksp->normtype != KSP_NORM_NONE && ksp->chknorm < i + 2) VecNormEnd(R, NORM_2, &dp);
156:     VecDotEnd(R, RP, &rho);
157:     VecDotEnd(S, RP, &d1);
158:     VecDotEnd(W, RP, &d2);
159:     VecDotEnd(Z, RP, &d3);


163:     beta  = (rho / rhoold) * (alpha / omega);
164:     alpha = rho / (d2 + beta * d1 - beta * omega * d3); /* alpha <- rho / (d2 + beta * d1 - beta * omega * d3) */

166:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
167:     ksp->its++;

169:     /* Residual replacement step  */
170:     if (i > 0 && i % 100 == 0 && i < 1001) {
171:       KSP_MatMult(ksp, pc->mat, X, R);
172:       VecAYPX(R, -1.0, B);              /* r  <- b - Ax */
173:       KSP_PCApply(ksp, R, R2);          /* r2 <- K r */
174:       KSP_MatMult(ksp, pc->mat, R2, W); /* w  <- A r2 */
175:       KSP_PCApply(ksp, W, W2);          /* w2 <- K w */
176:       KSP_MatMult(ksp, pc->mat, W2, T); /* t  <- A w2 */
177:       KSP_MatMult(ksp, pc->mat, P2, S); /* s  <- A p2 */
178:       KSP_PCApply(ksp, S, S2);          /* s2 <- K s */
179:       KSP_MatMult(ksp, pc->mat, S2, Z); /* z  <- A s2 */
180:       KSP_PCApply(ksp, Z, Z2);          /* z2 <- K z */
181:       KSP_MatMult(ksp, pc->mat, Z2, V); /* v  <- A z2 */
182:     }

184:     ksp->rnorm = dp;
185:     PetscObjectSAWsGrantAccess((PetscObject)ksp);
186:     KSPLogResidualHistory(ksp, dp);
187:     KSPMonitor(ksp, i + 1, dp);
188:     (*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP);
189:     if (ksp->reason) break;
190:     if (rho == 0.0) {
191:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
192:       break;
193:     }
194:     i++;
195:   } while (i < ksp->max_it);

197:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
198:   return 0;
199: }

201: /*MC
202:     KSPPIPEBCGS - Implements a pipelined BiCGStab method. [](sec_pipelineksp)

204:     Level: intermediate

206:     Notes:
207:     This method has only two non-blocking reductions per iteration, compared to 3 blocking for standard FBCGS.  The
208:     non-blocking reductions are overlapped by matrix-vector products and preconditioner applications.

210:     Periodic residual replacement may be used to increase robustness and maximal attainable accuracy.

212:     Like `KSPFBCGS`, the `KSPPIPEBCGS` implementation only allows for right preconditioning.

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

217:     Contributed by:
218:     Siegfried Cools, Universiteit Antwerpen,
219:     EXA2CT European Project on EXascale Algorithms and Advanced Computational Techniques, 2016.

221:     Reference:
222: .   * - S. Cools and W. Vanroose,
223:     "The communication-hiding pipelined BiCGStab method for the parallel solution of large unsymmetric linear systems",
224:     Parallel Computing, 65:1-20, 2017.

226: .seealso: [](chapter_ksp), `KSPFBCGS`, `KSPFBCGSR`, `KSPBCGS`, `KSPBCGSL`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBICG`, `KSPFBCGS`, `KSPFBCGSL`, `KSPSetPCSide()`,
227:            [](sec_pipelineksp), [](doc_faq_pipelined)
228: M*/
229: PETSC_EXTERN PetscErrorCode KSPCreate_PIPEBCGS(KSP ksp)
230: {
231:   KSP_BCGS *bcgs;

233:   PetscNew(&bcgs);

235:   ksp->data                = bcgs;
236:   ksp->ops->setup          = KSPSetUp_PIPEBCGS;
237:   ksp->ops->solve          = KSPSolve_PIPEBCGS;
238:   ksp->ops->destroy        = KSPDestroy_BCGS;
239:   ksp->ops->reset          = KSPReset_BCGS;
240:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
241:   ksp->ops->setfromoptions = KSPSetFromOptions_BCGS;

243:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2);
244:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1);
245:   return 0;
246: }