Actual source code: bcgsl.c

  1: /*
  2:    Implementation of BiCGstab(L) the paper by D.R. Fokkema,
  3:    "Enhanced implementation of BiCGStab(L) for solving linear systems
  4:    of equations". This uses tricky delayed updating ideas to prevent
  5:    round-off buildup.
  6: */
  7: #include <petsc/private/kspimpl.h>
  8: #include <../src/ksp/ksp/impls/bcgsl/bcgslimpl.h>
  9: #include <petscblaslapack.h>

 11: static PetscErrorCode KSPSolve_BCGSL(KSP ksp)
 12: {
 13:   KSP_BCGSL   *bcgsl = (KSP_BCGSL *)ksp->data;
 14:   PetscScalar  alpha, beta, omega, sigma;
 15:   PetscScalar  rho0, rho1;
 16:   PetscReal    kappa0, kappaA, kappa1;
 17:   PetscReal    ghat;
 18:   PetscReal    zeta, zeta0, rnmax_computed, rnmax_true, nrm0;
 19:   PetscBool    bUpdateX;
 20:   PetscInt     maxit;
 21:   PetscInt     h, i, j, k, vi, ell;
 22:   PetscBLASInt ldMZ, bierr;
 23:   PetscScalar  utb;
 24:   PetscReal    max_s, pinv_tol;

 26:   PetscFunctionBegin;
 27:   /* set up temporary vectors */
 28:   vi        = 0;
 29:   ell       = bcgsl->ell;
 30:   bcgsl->vB = ksp->work[vi];
 31:   vi++;
 32:   bcgsl->vRt = ksp->work[vi];
 33:   vi++;
 34:   bcgsl->vTm = ksp->work[vi];
 35:   vi++;
 36:   bcgsl->vvR = ksp->work + vi;
 37:   vi += ell + 1;
 38:   bcgsl->vvU = ksp->work + vi;
 39:   vi += ell + 1;
 40:   bcgsl->vXr = ksp->work[vi];
 41:   vi++;
 42:   PetscCall(PetscBLASIntCast(ell + 1, &ldMZ));

 44:   /* Prime the iterative solver */
 45:   PetscCall(KSPInitialResidual(ksp, VX, VTM, VB, VVR[0], ksp->vec_rhs));
 46:   PetscCall(VecNorm(VVR[0], NORM_2, &zeta0));
 47:   KSPCheckNorm(ksp, zeta0);
 48:   rnmax_computed = zeta0;
 49:   rnmax_true     = zeta0;

 51:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
 52:   ksp->its = 0;
 53:   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta0;
 54:   else ksp->rnorm = 0.0;
 55:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
 56:   PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP));
 57:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

 59:   PetscCall(VecSet(VVU[0], 0.0));
 60:   alpha = 0.;
 61:   rho0 = omega = 1;

 63:   if (bcgsl->delta > 0.0) {
 64:     PetscCall(VecCopy(VX, VXR));
 65:     PetscCall(VecSet(VX, 0.0));
 66:     PetscCall(VecCopy(VVR[0], VB));
 67:   } else {
 68:     PetscCall(VecCopy(ksp->vec_rhs, VB));
 69:   }

 71:   /* Life goes on */
 72:   PetscCall(VecCopy(VVR[0], VRT));
 73:   zeta = zeta0;

 75:   PetscCall(KSPGetTolerances(ksp, NULL, NULL, NULL, &maxit));

 77:   for (k = 0; k < maxit; k += bcgsl->ell) {
 78:     ksp->its = k;
 79:     if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta;
 80:     else ksp->rnorm = 0.0;

 82:     PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
 83:     PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));

 85:     PetscCall((*ksp->converged)(ksp, k, ksp->rnorm, &ksp->reason, ksp->cnvP));
 86:     if (ksp->reason < 0) PetscFunctionReturn(PETSC_SUCCESS);
 87:     if (ksp->reason) {
 88:       if (bcgsl->delta > 0.0) PetscCall(VecAXPY(VX, 1.0, VXR));
 89:       PetscFunctionReturn(PETSC_SUCCESS);
 90:     }

 92:     /* BiCG part */
 93:     rho0 = -omega * rho0;
 94:     nrm0 = zeta;
 95:     for (j = 0; j < bcgsl->ell; j++) {
 96:       /* rho1 <- r_j' * r_tilde */
 97:       PetscCall(VecDot(VVR[j], VRT, &rho1));
 98:       KSPCheckDot(ksp, rho1);
 99:       if (rho1 == 0.0) {
100:         ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
101:         PetscFunctionReturn(PETSC_SUCCESS);
102:       }
103:       beta = alpha * (rho1 / rho0);
104:       rho0 = rho1;
105:       for (i = 0; i <= j; i++) {
106:         /* u_i <- r_i - beta*u_i */
107:         PetscCall(VecAYPX(VVU[i], -beta, VVR[i]));
108:       }
109:       /* u_{j+1} <- inv(K)*A*u_j */
110:       PetscCall(KSP_PCApplyBAorAB(ksp, VVU[j], VVU[j + 1], VTM));

112:       PetscCall(VecDot(VVU[j + 1], VRT, &sigma));
113:       KSPCheckDot(ksp, sigma);
114:       if (sigma == 0.0) {
115:         ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
116:         PetscFunctionReturn(PETSC_SUCCESS);
117:       }
118:       alpha = rho1 / sigma;

120:       /* x <- x + alpha*u_0 */
121:       PetscCall(VecAXPY(VX, alpha, VVU[0]));

123:       for (i = 0; i <= j; i++) {
124:         /* r_i <- r_i - alpha*u_{i+1} */
125:         PetscCall(VecAXPY(VVR[i], -alpha, VVU[i + 1]));
126:       }

128:       /* r_{j+1} <- inv(K)*A*r_j */
129:       PetscCall(KSP_PCApplyBAorAB(ksp, VVR[j], VVR[j + 1], VTM));

131:       PetscCall(VecNorm(VVR[0], NORM_2, &nrm0));
132:       KSPCheckNorm(ksp, nrm0);
133:       if (bcgsl->delta > 0.0) {
134:         if (rnmax_computed < nrm0) rnmax_computed = nrm0;
135:         if (rnmax_true < nrm0) rnmax_true = nrm0;
136:       }

138:       /* NEW: check for early exit */
139:       PetscCall((*ksp->converged)(ksp, k + j, nrm0, &ksp->reason, ksp->cnvP));
140:       if (ksp->reason) {
141:         PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
142:         ksp->its   = k + j;
143:         ksp->rnorm = nrm0;

145:         PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
146:         if (ksp->reason < 0) PetscFunctionReturn(PETSC_SUCCESS);
147:       }
148:     }

150:     /* Polynomial part */
151:     for (i = 0; i <= bcgsl->ell; ++i) PetscCall(VecMDot(VVR[i], i + 1, VVR, &MZa[i * ldMZ]));
152:     /* Symmetrize MZa */
153:     for (i = 0; i <= bcgsl->ell; ++i) {
154:       for (j = i + 1; j <= bcgsl->ell; ++j) MZa[i * ldMZ + j] = MZa[j * ldMZ + i] = PetscConj(MZa[j * ldMZ + i]);
155:     }
156:     /* Copy MZa to MZb */
157:     PetscCall(PetscArraycpy(MZb, MZa, ldMZ * ldMZ));

159:     if (!bcgsl->bConvex || bcgsl->ell == 1) {
160:       PetscBLASInt ione = 1, bell;
161:       PetscCall(PetscBLASIntCast(bcgsl->ell, &bell));

163:       AY0c[0] = -1;
164:       if (bcgsl->pinv) {
165: #if defined(PETSC_USE_COMPLEX)
166:         PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("A", "A", &bell, &bell, &MZa[1 + ldMZ], &ldMZ, bcgsl->s, bcgsl->u, &bell, bcgsl->v, &bell, bcgsl->work, &bcgsl->lwork, bcgsl->realwork, &bierr));
167: #else
168:         PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("A", "A", &bell, &bell, &MZa[1 + ldMZ], &ldMZ, bcgsl->s, bcgsl->u, &bell, bcgsl->v, &bell, bcgsl->work, &bcgsl->lwork, &bierr));
169: #endif
170:         if (bierr != 0) {
171:           ksp->reason = KSP_DIVERGED_BREAKDOWN;
172:           PetscFunctionReturn(PETSC_SUCCESS);
173:         }
174:         /* Apply pseudo-inverse */
175:         max_s = bcgsl->s[0];
176:         for (i = 1; i < bell; i++) {
177:           if (bcgsl->s[i] > max_s) max_s = bcgsl->s[i];
178:         }
179:         /* tolerance is hardwired to bell*max(s)*PETSC_MACHINE_EPSILON */
180:         pinv_tol = bell * max_s * PETSC_MACHINE_EPSILON;
181:         PetscCall(PetscArrayzero(&AY0c[1], bell));
182:         for (i = 0; i < bell; i++) {
183:           if (bcgsl->s[i] >= pinv_tol) {
184:             utb = 0.;
185:             for (j = 0; j < bell; j++) utb += MZb[1 + j] * bcgsl->u[i * bell + j];

187:             for (j = 0; j < bell; j++) AY0c[1 + j] += utb / bcgsl->s[i] * bcgsl->v[j * bell + i];
188:           }
189:         }
190:       } else {
191:         PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("Lower", &bell, &MZa[1 + ldMZ], &ldMZ, &bierr));
192:         if (bierr != 0) {
193:           ksp->reason = KSP_DIVERGED_BREAKDOWN;
194:           PetscFunctionReturn(PETSC_SUCCESS);
195:         }
196:         PetscCall(PetscArraycpy(&AY0c[1], &MZb[1], bcgsl->ell));
197:         PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("Lower", &bell, &ione, &MZa[1 + ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
198:       }
199:     } else {
200:       PetscBLASInt ione = 1;
201:       PetscScalar  aone = 1.0, azero = 0.0;
202:       PetscBLASInt neqs;
203:       PetscCall(PetscBLASIntCast(bcgsl->ell - 1, &neqs));

205:       PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("Lower", &neqs, &MZa[1 + ldMZ], &ldMZ, &bierr));
206:       if (bierr != 0) {
207:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
208:         PetscFunctionReturn(PETSC_SUCCESS);
209:       }
210:       PetscCall(PetscArraycpy(&AY0c[1], &MZb[1], bcgsl->ell - 1));
211:       PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1 + ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
212:       AY0c[0]          = -1;
213:       AY0c[bcgsl->ell] = 0.;

215:       PetscCall(PetscArraycpy(&AYlc[1], &MZb[1 + ldMZ * (bcgsl->ell)], bcgsl->ell - 1));
216:       PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1 + ldMZ], &ldMZ, &AYlc[1], &ldMZ, &bierr));

218:       AYlc[0]          = 0.;
219:       AYlc[bcgsl->ell] = -1;

221:       PetscCallBLAS("BLASgemv", BLASgemv_("NoTr", &ldMZ, &ldMZ, &aone, MZb, &ldMZ, AY0c, &ione, &azero, AYtc, &ione));

223:       kappa0 = PetscRealPart(BLASdot_(&ldMZ, AY0c, &ione, AYtc, &ione));

225:       /* round-off can cause negative kappa's */
226:       if (kappa0 < 0) kappa0 = -kappa0;
227:       kappa0 = PetscSqrtReal(kappa0);

229:       kappaA = PetscRealPart(BLASdot_(&ldMZ, AYlc, &ione, AYtc, &ione));

231:       PetscCallBLAS("BLASgemv", BLASgemv_("noTr", &ldMZ, &ldMZ, &aone, MZb, &ldMZ, AYlc, &ione, &azero, AYtc, &ione));

233:       kappa1 = PetscRealPart(BLASdot_(&ldMZ, AYlc, &ione, AYtc, &ione));

235:       if (kappa1 < 0) kappa1 = -kappa1;
236:       kappa1 = PetscSqrtReal(kappa1);

238:       if (kappa0 != 0.0 && kappa1 != 0.0) {
239:         if (kappaA < 0.7 * kappa0 * kappa1) {
240:           ghat = (kappaA < 0.0) ? -0.7 * kappa0 / kappa1 : 0.7 * kappa0 / kappa1;
241:         } else {
242:           ghat = kappaA / (kappa1 * kappa1);
243:         }
244:         for (i = 0; i <= bcgsl->ell; i++) AY0c[i] = AY0c[i] - ghat * AYlc[i];
245:       }
246:     }

248:     omega = AY0c[bcgsl->ell];
249:     for (h = bcgsl->ell; h > 0 && omega == 0.0; h--) omega = AY0c[h];
250:     if (omega == 0.0) {
251:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
252:       PetscFunctionReturn(PETSC_SUCCESS);
253:     }

255:     PetscCall(VecMAXPY(VX, bcgsl->ell, AY0c + 1, VVR));
256:     for (i = 1; i <= bcgsl->ell; i++) AY0c[i] *= -1.0;
257:     PetscCall(VecMAXPY(VVU[0], bcgsl->ell, AY0c + 1, VVU + 1));
258:     PetscCall(VecMAXPY(VVR[0], bcgsl->ell, AY0c + 1, VVR + 1));
259:     for (i = 1; i <= bcgsl->ell; i++) AY0c[i] *= -1.0;
260:     PetscCall(VecNorm(VVR[0], NORM_2, &zeta));
261:     KSPCheckNorm(ksp, zeta);

263:     /* Accurate Update */
264:     if (bcgsl->delta > 0.0) {
265:       if (rnmax_computed < zeta) rnmax_computed = zeta;
266:       if (rnmax_true < zeta) rnmax_true = zeta;

268:       bUpdateX = (PetscBool)(zeta < bcgsl->delta * zeta0 && zeta0 <= rnmax_computed);
269:       if ((zeta < bcgsl->delta * rnmax_true && zeta0 <= rnmax_true) || bUpdateX) {
270:         /* r0 <- b-inv(K)*A*X */
271:         PetscCall(KSP_PCApplyBAorAB(ksp, VX, VVR[0], VTM));
272:         PetscCall(VecAYPX(VVR[0], -1.0, VB));
273:         rnmax_true = zeta;

275:         if (bUpdateX) {
276:           PetscCall(VecAXPY(VXR, 1.0, VX));
277:           PetscCall(VecSet(VX, 0.0));
278:           PetscCall(VecCopy(VVR[0], VB));
279:           rnmax_computed = zeta;
280:         }
281:       }
282:     }
283:   }
284:   if (bcgsl->delta > 0.0) PetscCall(VecAXPY(VX, 1.0, VXR));

286:   ksp->its = k;
287:   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta;
288:   else ksp->rnorm = 0.0;
289:   PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
290:   PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
291:   PetscCall((*ksp->converged)(ksp, k, ksp->rnorm, &ksp->reason, ksp->cnvP));
292:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
293:   PetscFunctionReturn(PETSC_SUCCESS);
294: }

296: /*@
297:   KSPBCGSLSetXRes - Sets the parameter governing when
298:   exact residuals will be used instead of computed residuals for `KSPCBGSL`.

300:   Logically Collective

302:   Input Parameters:
303: + ksp   - iterative context of type `KSPBCGSL`
304: - delta - computed residuals are used alone when delta is not positive

306:   Options Database Key:
307: . -ksp_bcgsl_xres delta - Threshold used to decide when to refresh computed residuals

309:   Level: intermediate

311: .seealso: [](ch_ksp), `KSPBCGSLSetEll()`, `KSPBCGSLSetPol()`, `KSP`, `KSPCBGSL`, `KSPBCGSLSetUsePseudoinverse()`
312: @*/
313: PetscErrorCode KSPBCGSLSetXRes(KSP ksp, PetscReal delta)
314: {
315:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;

317:   PetscFunctionBegin;
319:   if (ksp->setupstage) {
320:     if ((delta <= 0 && bcgsl->delta > 0) || (delta > 0 && bcgsl->delta <= 0)) {
321:       PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));
322:       PetscCall(PetscFree5(AY0c, AYlc, AYtc, MZa, MZb));
323:       PetscCall(PetscFree4(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v));
324:       ksp->setupstage = KSP_SETUP_NEW;
325:     }
326:   }
327:   bcgsl->delta = delta;
328:   PetscFunctionReturn(PETSC_SUCCESS);
329: }

331: /*@
332:   KSPBCGSLSetUsePseudoinverse - Use pseudoinverse (via SVD) to solve polynomial part of the update in `KSPCBGSL` solver

334:   Logically Collective

336:   Input Parameters:
337: + ksp      - iterative context of type `KSPCBGSL`
338: - use_pinv - set to `PETSC_TRUE` when using pseudoinverse

340:   Options Database Key:
341: . -ksp_bcgsl_pinv <true,false> - use pseudoinverse

343:   Level: intermediate

345: .seealso: [](ch_ksp), `KSPBCGSLSetEll()`, `KSP`, `KSPCBGSL`, `KSPBCGSLSetPol()`, `KSPBCGSLSetXRes()`
346: @*/
347: PetscErrorCode KSPBCGSLSetUsePseudoinverse(KSP ksp, PetscBool use_pinv)
348: {
349:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;

351:   PetscFunctionBegin;
352:   bcgsl->pinv = use_pinv;
353:   PetscFunctionReturn(PETSC_SUCCESS);
354: }

356: /*@
357:   KSPBCGSLSetPol - Sets the type of polynomial part that will
358:   be used in the `KSPCBGSL` `KSPSolve()`

360:   Logically Collective

362:   Input Parameters:
363: + ksp   - iterative context of type `KSPCBGSL`
364: - uMROR - set to `PETSC_TRUE` when the polynomial is a convex combination of an MR and an OR step.

366:   Options Database Keys:
367: + -ksp_bcgsl_cxpoly - use enhanced polynomial
368: - -ksp_bcgsl_mrpoly - use standard polynomial

370:   Level: intermediate

372: .seealso: [](ch_ksp), `KSP`, `KSPBCGSL`, `KSPCreate()`, `KSPSetType()`, `KSPCBGSL`, `KSPBCGSLSetUsePseudoinverse()`, `KSPBCGSLSetEll()`, `KSPBCGSLSetXRes()`
373: @*/
374: PetscErrorCode KSPBCGSLSetPol(KSP ksp, PetscBool uMROR)
375: {
376:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;

378:   PetscFunctionBegin;

381:   if (!ksp->setupstage) bcgsl->bConvex = uMROR;
382:   else if (bcgsl->bConvex != uMROR) {
383:     /* free the data structures,
384:        then create them again
385:      */
386:     PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));
387:     PetscCall(PetscFree5(AY0c, AYlc, AYtc, MZa, MZb));
388:     PetscCall(PetscFree4(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v));

390:     bcgsl->bConvex  = uMROR;
391:     ksp->setupstage = KSP_SETUP_NEW;
392:   }
393:   PetscFunctionReturn(PETSC_SUCCESS);
394: }

396: /*@
397:   KSPBCGSLSetEll - Sets the number of search directions to use in the `KSPCBGSL` solver

399:   Logically Collective

401:   Input Parameters:
402: + ksp - iterative context of type `KSPCBGSL`
403: - ell - number of search directions

405:   Options Database Key:
406: . -ksp_bcgsl_ell ell - Number of Krylov search directions

408:   Level: intermediate

410:   Notes:
411:   For large `ell` it is common for the polynomial update problem to become singular (due to happy breakdown for smallish
412:   test problems, but also for larger problems). Consequently, by default, the system is solved by pseudoinverse, which
413:   allows the iteration to complete successfully. See `KSPBCGSLSetUsePseudoinverse()` to switch to a conventional solve.

415: .seealso: [](ch_ksp), `KSPBCGSLSetUsePseudoinverse()`, `KSP`, `KSPBCGSL`, `KSPBCGSLSetPol()`, `KSPBCGSLSetXRes()`
416: @*/
417: PetscErrorCode KSPBCGSLSetEll(KSP ksp, PetscInt ell)
418: {
419:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;

421:   PetscFunctionBegin;
422:   PetscCheck(ell >= 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "KSPBCGSLSetEll: second argument must be positive");

425:   if (!ksp->setupstage) bcgsl->ell = ell;
426:   else if (bcgsl->ell != ell) {
427:     /* free the data structures, then create them again */
428:     PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));
429:     PetscCall(PetscFree5(AY0c, AYlc, AYtc, MZa, MZb));
430:     PetscCall(PetscFree4(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v));

432:     bcgsl->ell      = ell;
433:     ksp->setupstage = KSP_SETUP_NEW;
434:   }
435:   PetscFunctionReturn(PETSC_SUCCESS);
436: }

438: static PetscErrorCode KSPView_BCGSL(KSP ksp, PetscViewer viewer)
439: {
440:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
441:   PetscBool  isascii;

443:   PetscFunctionBegin;
444:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));

446:   if (isascii) {
447:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Ell = %" PetscInt_FMT "\n", bcgsl->ell));
448:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Delta = %g\n", (double)bcgsl->delta));
449:   }
450:   PetscFunctionReturn(PETSC_SUCCESS);
451: }

453: static PetscErrorCode KSPSetFromOptions_BCGSL(KSP ksp, PetscOptionItems *PetscOptionsObject)
454: {
455:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
456:   PetscInt   this_ell;
457:   PetscReal  delta;
458:   PetscBool  flga = PETSC_FALSE, flg;

460:   PetscFunctionBegin;
461:   PetscOptionsHeadBegin(PetscOptionsObject, "KSPBCGSL Options");

463:   /* Set number of search directions */
464:   PetscCall(PetscOptionsInt("-ksp_bcgsl_ell", "Number of Krylov search directions", "KSPBCGSLSetEll", bcgsl->ell, &this_ell, &flg));
465:   if (flg) PetscCall(KSPBCGSLSetEll(ksp, this_ell));

467:   /* Set polynomial type */
468:   PetscCall(PetscOptionsBool("-ksp_bcgsl_cxpoly", "Polynomial part of BiCGStabL is MinRes + OR", "KSPBCGSLSetPol", flga, &flga, NULL));
469:   if (flga) {
470:     PetscCall(KSPBCGSLSetPol(ksp, PETSC_TRUE));
471:   } else {
472:     flg = PETSC_FALSE;
473:     PetscCall(PetscOptionsBool("-ksp_bcgsl_mrpoly", "Polynomial part of BiCGStabL is MinRes", "KSPBCGSLSetPol", flg, &flg, NULL));
474:     PetscCall(KSPBCGSLSetPol(ksp, PETSC_FALSE));
475:   }

477:   /* Will computed residual be refreshed? */
478:   PetscCall(PetscOptionsReal("-ksp_bcgsl_xres", "Threshold used to decide when to refresh computed residuals", "KSPBCGSLSetXRes", bcgsl->delta, &delta, &flg));
479:   if (flg) PetscCall(KSPBCGSLSetXRes(ksp, delta));

481:   /* Use pseudoinverse? */
482:   flg = bcgsl->pinv;
483:   PetscCall(PetscOptionsBool("-ksp_bcgsl_pinv", "Polynomial correction via pseudoinverse", "KSPBCGSLSetUsePseudoinverse", flg, &flg, NULL));
484:   PetscCall(KSPBCGSLSetUsePseudoinverse(ksp, flg));
485:   PetscOptionsHeadEnd();
486:   PetscFunctionReturn(PETSC_SUCCESS);
487: }

489: static PetscErrorCode KSPSetUp_BCGSL(KSP ksp)
490: {
491:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
492:   PetscInt   ell = bcgsl->ell, ldMZ = ell + 1;

494:   PetscFunctionBegin;
495:   PetscCall(KSPSetWorkVecs(ksp, 6 + 2 * ell));
496:   PetscCall(PetscMalloc5(ldMZ, &AY0c, ldMZ, &AYlc, ldMZ, &AYtc, ldMZ * ldMZ, &MZa, ldMZ * ldMZ, &MZb));
497:   PetscCall(PetscBLASIntCast(5 * ell, &bcgsl->lwork));
498:   PetscCall(PetscMalloc5(bcgsl->lwork, &bcgsl->work, ell, &bcgsl->s, ell * ell, &bcgsl->u, ell * ell, &bcgsl->v, 5 * ell, &bcgsl->realwork));
499:   PetscFunctionReturn(PETSC_SUCCESS);
500: }

502: static PetscErrorCode KSPReset_BCGSL(KSP ksp)
503: {
504:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;

506:   PetscFunctionBegin;
507:   PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));
508:   PetscCall(PetscFree5(AY0c, AYlc, AYtc, MZa, MZb));
509:   PetscCall(PetscFree5(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v, bcgsl->realwork));
510:   PetscFunctionReturn(PETSC_SUCCESS);
511: }

513: static PetscErrorCode KSPDestroy_BCGSL(KSP ksp)
514: {
515:   PetscFunctionBegin;
516:   PetscCall(KSPReset_BCGSL(ksp));
517:   PetscCall(KSPDestroyDefault(ksp));
518:   PetscFunctionReturn(PETSC_SUCCESS);
519: }

521: /*MC
522:     KSPBCGSL - Implements a slight variant of the Enhanced
523:                 BiCGStab(L) algorithm in (3) and (2).  The variation
524:                 concerns cases when either kappa0**2 or kappa1**2 is
525:                 negative due to round-off. Kappa0 has also been pulled
526:                 out of the denominator in the formula for ghat.

528:    Options Database Keys:
529: +  -ksp_bcgsl_ell <ell> - Number of Krylov search directions, defaults to 2, cf. `KSPBCGSLSetEll()`
530: .  -ksp_bcgsl_cxpol - Use a convex function of the MinRes and OR polynomials after the BiCG step instead of default MinRes, cf. `KSPBCGSLSetPol()`
531: .  -ksp_bcgsl_mrpoly - Use the default MinRes polynomial after the BiCG step, cf. `KSPBCGSLSetPol()`
532: .  -ksp_bcgsl_xres <res> - Threshold used to decide when to refresh computed residuals, cf. `KSPBCGSLSetXRes()`
533: -  -ksp_bcgsl_pinv <true/false> - (de)activate use of pseudoinverse, cf. `KSPBCGSLSetUsePseudoinverse()`

535:    Level: intermediate

537:    Contributed by:
538:    Joel M. Malard, email jm.malard@pnl.gov

540:    Developer Notes:
541:    This has not been completely cleaned up into PETSc style.

543:    All the BLAS and LAPACK calls in the source should be removed and replaced with loops and the macros for block solvers converted from LINPACK.

545:    References:
546: +  * - G.L.G. Sleijpen, H.A. van der Vorst, "An overview of
547:          approaches for the stable computation of hybrid BiCG
548:          methods", Applied Numerical Mathematics: Transactions
549:          f IMACS, 19(3), 1996.
550: .  * - G.L.G. Sleijpen, H.A. van der Vorst, D.R. Fokkema,
551:          "BiCGStab(L) and other hybrid BiCG methods",
552:           Numerical Algorithms, 7, 1994.
553: -  * - D.R. Fokkema, "Enhanced implementation of BiCGStab(L)
554:          for solving linear systems of equations", preprint
555:          from www.citeseer.com.

557: .seealso: [](ch_ksp), `KSPFBCGS`, `KSPFBCGSR`, `KSPBCGS`, `KSPPIPEBCGS`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPBCGS`, `KSPSetPCSide()`,
558:           `KSPBCGSLSetEll()`, `KSPBCGSLSetXRes()`, `KSPBCGSLSetUsePseudoinverse()`, `KSPBCGSLSetPol()`
559: M*/
560: PETSC_EXTERN PetscErrorCode KSPCreate_BCGSL(KSP ksp)
561: {
562:   KSP_BCGSL *bcgsl;

564:   PetscFunctionBegin;
565:   /* allocate BiCGStab(L) context */
566:   PetscCall(PetscNew(&bcgsl));
567:   ksp->data = (void *)bcgsl;

569:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
570:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
571:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));

573:   ksp->ops->setup          = KSPSetUp_BCGSL;
574:   ksp->ops->solve          = KSPSolve_BCGSL;
575:   ksp->ops->reset          = KSPReset_BCGSL;
576:   ksp->ops->destroy        = KSPDestroy_BCGSL;
577:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
578:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
579:   ksp->ops->setfromoptions = KSPSetFromOptions_BCGSL;
580:   ksp->ops->view           = KSPView_BCGSL;

582:   /* Let the user redefine the number of directions vectors */
583:   bcgsl->ell = 2;

585:   /*Choose between a single MR step or an averaged MR/OR */
586:   bcgsl->bConvex = PETSC_FALSE;

588:   bcgsl->pinv = PETSC_TRUE; /* Use the reliable method by default */

590:   /* Set the threshold for when exact residuals will be used */
591:   bcgsl->delta = 0.0;
592:   PetscFunctionReturn(PETSC_SUCCESS);
593: }