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

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

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

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

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

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

 74:   KSPGetTolerances(ksp, NULL, NULL, NULL, &maxit);

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

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

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

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

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

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

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

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

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

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

144:         PetscObjectSAWsGrantAccess((PetscObject)ksp);
145:         if (ksp->reason < 0) return 0;
146:       }
147:     }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

299:    Logically Collective

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

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

308:    Level: intermediate

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

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

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

332:    Logically Collective

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

338:    Options Database Key:
339: .  -ksp_bcgsl_pinv - use pseudoinverse

341:    Level: intermediate

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

349:   bcgsl->pinv = use_pinv;
350:   return 0;
351: }

353: /*@
354:    KSPBCGSLSetPol - Sets the type of polynomial part that will
355:    be used in the  `KSPCBGSL` solver.

357:    Logically Collective

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

363:    Options Database Keys:
364: +  -ksp_bcgsl_cxpoly - use enhanced polynomial
365: -  -ksp_bcgsl_mrpoly - use standard polynomial

367:    Level: intermediate

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


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

386:     bcgsl->bConvex  = uMROR;
387:     ksp->setupstage = KSP_SETUP_NEW;
388:   }
389:   return 0;
390: }

392: /*@
393:    KSPBCGSLSetEll - Sets the number of search directions in `KSPCBGSL` solver

395:    Logically Collective

397:    Input Parameters:
398: +  ksp - iterative context of type `KSPCBGSL`
399: -  ell - number of search directions

401:    Options Database Keys:
402: .  -ksp_bcgsl_ell ell - Number of Krylov search directions

404:    Level: intermediate

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

411: .seealso: [](chapter_ksp), `KSPBCGSLSetUsePseudoinverse()`, `KSP`, `KSPBCGSL`, `KSPBCGSLSetPol()`, `KSPBCGSLSetXRes()`
412: @*/
413: PetscErrorCode KSPBCGSLSetEll(KSP ksp, PetscInt ell)
414: {
415:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;


420:   if (!ksp->setupstage) bcgsl->ell = ell;
421:   else if (bcgsl->ell != ell) {
422:     /* free the data structures, then create them again */
423:     VecDestroyVecs(ksp->nwork, &ksp->work);
424:     PetscFree5(AY0c, AYlc, AYtc, MZa, MZb);
425:     PetscFree4(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v);

427:     bcgsl->ell      = ell;
428:     ksp->setupstage = KSP_SETUP_NEW;
429:   }
430:   return 0;
431: }

433: PetscErrorCode KSPView_BCGSL(KSP ksp, PetscViewer viewer)
434: {
435:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
436:   PetscBool  isascii;

438:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);

440:   if (isascii) {
441:     PetscViewerASCIIPrintf(viewer, "  Ell = %" PetscInt_FMT "\n", bcgsl->ell);
442:     PetscViewerASCIIPrintf(viewer, "  Delta = %g\n", (double)bcgsl->delta);
443:   }
444:   return 0;
445: }

447: PetscErrorCode KSPSetFromOptions_BCGSL(KSP ksp, PetscOptionItems *PetscOptionsObject)
448: {
449:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
450:   PetscInt   this_ell;
451:   PetscReal  delta;
452:   PetscBool  flga = PETSC_FALSE, flg;

454:   PetscOptionsHeadBegin(PetscOptionsObject, "KSPBCGSL Options");

456:   /* Set number of search directions */
457:   PetscOptionsInt("-ksp_bcgsl_ell", "Number of Krylov search directions", "KSPBCGSLSetEll", bcgsl->ell, &this_ell, &flg);
458:   if (flg) KSPBCGSLSetEll(ksp, this_ell);

460:   /* Set polynomial type */
461:   PetscOptionsBool("-ksp_bcgsl_cxpoly", "Polynomial part of BiCGStabL is MinRes + OR", "KSPBCGSLSetPol", flga, &flga, NULL);
462:   if (flga) {
463:     KSPBCGSLSetPol(ksp, PETSC_TRUE);
464:   } else {
465:     flg = PETSC_FALSE;
466:     PetscOptionsBool("-ksp_bcgsl_mrpoly", "Polynomial part of BiCGStabL is MinRes", "KSPBCGSLSetPol", flg, &flg, NULL);
467:     KSPBCGSLSetPol(ksp, PETSC_FALSE);
468:   }

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

474:   /* Use pseudoinverse? */
475:   flg = bcgsl->pinv;
476:   PetscOptionsBool("-ksp_bcgsl_pinv", "Polynomial correction via pseudoinverse", "KSPBCGSLSetUsePseudoinverse", flg, &flg, NULL);
477:   KSPBCGSLSetUsePseudoinverse(ksp, flg);
478:   PetscOptionsHeadEnd();
479:   return 0;
480: }

482: PetscErrorCode KSPSetUp_BCGSL(KSP ksp)
483: {
484:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
485:   PetscInt   ell = bcgsl->ell, ldMZ = ell + 1;

487:   KSPSetWorkVecs(ksp, 6 + 2 * ell);
488:   PetscMalloc5(ldMZ, &AY0c, ldMZ, &AYlc, ldMZ, &AYtc, ldMZ * ldMZ, &MZa, ldMZ * ldMZ, &MZb);
489:   PetscBLASIntCast(5 * ell, &bcgsl->lwork);
490:   PetscMalloc5(bcgsl->lwork, &bcgsl->work, ell, &bcgsl->s, ell * ell, &bcgsl->u, ell * ell, &bcgsl->v, 5 * ell, &bcgsl->realwork);
491:   return 0;
492: }

494: PetscErrorCode KSPReset_BCGSL(KSP ksp)
495: {
496:   KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;

498:   VecDestroyVecs(ksp->nwork, &ksp->work);
499:   PetscFree5(AY0c, AYlc, AYtc, MZa, MZb);
500:   PetscFree5(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v, bcgsl->realwork);
501:   return 0;
502: }

504: PetscErrorCode KSPDestroy_BCGSL(KSP ksp)
505: {
506:   KSPReset_BCGSL(ksp);
507:   KSPDestroyDefault(ksp);
508:   return 0;
509: }

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

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

525:    Level: intermediate

527:     References:
528: +   * - G.L.G. Sleijpen, H.A. van der Vorst, "An overview of
529:          approaches for the stable computation of hybrid BiCG
530:          methods", Applied Numerical Mathematics: Transactions
531:          f IMACS, 19(3), 1996.
532: .   * - G.L.G. Sleijpen, H.A. van der Vorst, D.R. Fokkema,
533:          "BiCGStab(L) and other hybrid BiCG methods",
534:           Numerical Algorithms, 7, 1994.
535: -   * - D.R. Fokkema, "Enhanced implementation of BiCGStab(L)
536:          for solving linear systems of equations", preprint
537:          from www.citeseer.com.

539:    Contributed by:
540:    Joel M. Malard, email jm.malard@pnl.gov

542:    Developer Notes:
543:     This has not been completely cleaned up into PETSc style.

545:     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.

547: .seealso: [](chapter_ksp), `KSPFBCGS`, `KSPFBCGSR`, `KSPBCGS`, `KSPPIPEBCGS`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPBCGS`, `KSPSetPCSide()`,
548:           `KSPBCGSLSetEll()`, `KSPBCGSLSetXRes()`, `KSPBCGSLSetUsePseudoinverse()`, `KSPBCGSLSetPol()`
549: M*/
550: PETSC_EXTERN PetscErrorCode KSPCreate_BCGSL(KSP ksp)
551: {
552:   KSP_BCGSL *bcgsl;

554:   /* allocate BiCGStab(L) context */
555:   PetscNew(&bcgsl);
556:   ksp->data = (void *)bcgsl;

558:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3);
559:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2);
560:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1);

562:   ksp->ops->setup          = KSPSetUp_BCGSL;
563:   ksp->ops->solve          = KSPSolve_BCGSL;
564:   ksp->ops->reset          = KSPReset_BCGSL;
565:   ksp->ops->destroy        = KSPDestroy_BCGSL;
566:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
567:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
568:   ksp->ops->setfromoptions = KSPSetFromOptions_BCGSL;
569:   ksp->ops->view           = KSPView_BCGSL;

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

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

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

579:   /* Set the threshold for when exact residuals will be used */
580:   bcgsl->delta = 0.0;
581:   return 0;
582: }