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:  * This has not been completely cleaned up into PETSc style.
  8:  *
  9:  * All the BLAS and LAPACK calls below should be removed and replaced with
 10:  * loops and the macros for block solvers converted from LINPACK; there is no way
 11:  * calls to BLAS/LAPACK make sense for size 2, 3, 4, etc.
 12:  */
 13: #include <petsc/private/kspimpl.h>
 14: #include <../src/ksp/ksp/impls/bcgsl/bcgslimpl.h>
 15: #include <petscblaslapack.h>

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

 34:   /* set up temporary vectors */
 35:   vi         = 0;
 36:   ell        = bcgsl->ell;
 37:   bcgsl->vB  = ksp->work[vi]; vi++;
 38:   bcgsl->vRt = ksp->work[vi]; vi++;
 39:   bcgsl->vTm = ksp->work[vi]; vi++;
 40:   bcgsl->vvR = ksp->work+vi; vi += ell+1;
 41:   bcgsl->vvU = ksp->work+vi; vi += ell+1;
 42:   bcgsl->vXr = ksp->work[vi]; vi++;
 43:   PetscBLASIntCast(ell+1,&ldMZ);

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

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

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

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

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

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

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

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

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

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

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

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

126:       for (i=0; i<=j; i++) {
127:         /* r_i <- r_i - alpha*u_{i+1} */
128:         VecAXPY(VVR[i], -alpha, VVU[i+1]);
129:       }

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

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

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

148:         PetscObjectSAWsGrantAccess((PetscObject)ksp);
149:         if (ksp->reason < 0) return(0);
150:       }
151:     }

153:     /* Polynomial part */
154:     for (i = 0; i <= bcgsl->ell; ++i) {
155:       VecMDot(VVR[i], i+1, VVR, &MZa[i*ldMZ]);
156:     }
157:     /* Symmetrize MZa */
158:     for (i = 0; i <= bcgsl->ell; ++i) {
159:       for (j = i+1; j <= bcgsl->ell; ++j) {
160:         MZa[i*ldMZ+j] = MZa[j*ldMZ+i] = PetscConj(MZa[j*ldMZ+i]);
161:       }
162:     }
163:     /* Copy MZa to MZb */
164:     PetscArraycpy(MZb,MZa,ldMZ*ldMZ);

166:     if (!bcgsl->bConvex || bcgsl->ell==1) {
167:       PetscBLASInt ione = 1,bell;
168:       PetscBLASIntCast(bcgsl->ell,&bell);

170:       AY0c[0] = -1;
171:       if (bcgsl->pinv) {
172: #  if defined(PETSC_USE_COMPLEX)
173:         PetscStackCallBLAS("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));
174: #  else
175:         PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&bell,&bell,&MZa[1+ldMZ],&ldMZ,bcgsl->s,bcgsl->u,&bell,bcgsl->v,&bell,bcgsl->work,&bcgsl->lwork,&bierr));
176: #  endif
177:         if (bierr!=0) {
178:           ksp->reason = KSP_DIVERGED_BREAKDOWN;
179:           return(0);
180:         }
181:         /* Apply pseudo-inverse */
182:         max_s = bcgsl->s[0];
183:         for (i=1; i<bell; i++) {
184:           if (bcgsl->s[i] > max_s) {
185:             max_s = bcgsl->s[i];
186:           }
187:         }
188:         /* tolerance is hardwired to bell*max(s)*PETSC_MACHINE_EPSILON */
189:         pinv_tol = bell*max_s*PETSC_MACHINE_EPSILON;
190:         PetscArrayzero(&AY0c[1],bell);
191:         for (i=0; i<bell; i++) {
192:           if (bcgsl->s[i] >= pinv_tol) {
193:             utb=0.;
194:             for (j=0; j<bell; j++) {
195:               utb += MZb[1+j]*bcgsl->u[i*bell+j];
196:             }

198:             for (j=0; j<bell; j++) {
199:               AY0c[1+j] += utb/bcgsl->s[i]*bcgsl->v[j*bell+i];
200:             }
201:           }
202:         }
203:       } else {
204:         PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("Lower", &bell, &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);
210:         PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &bell, &ione, &MZa[1+ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
211:       }
212:     } else {
213:       PetscBLASInt ione = 1;
214:       PetscScalar  aone = 1.0, azero = 0.0;
215:       PetscBLASInt neqs;
216:       PetscBLASIntCast(bcgsl->ell-1,&neqs);

218:       PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("Lower", &neqs, &MZa[1+ldMZ], &ldMZ, &bierr));
219:       if (bierr!=0) {
220:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
221:         return(0);
222:       }
223:       PetscArraycpy(&AY0c[1],&MZb[1],bcgsl->ell-1);
224:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1+ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
225:       AY0c[0]          = -1;
226:       AY0c[bcgsl->ell] = 0.;

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

231:       AYlc[0]          = 0.;
232:       AYlc[bcgsl->ell] = -1;

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

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

238:       /* round-off can cause negative kappa's */
239:       if (kappa0<0) kappa0 = -kappa0;
240:       kappa0 = PetscSqrtReal(kappa0);

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

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

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

248:       if (kappa1<0) kappa1 = -kappa1;
249:       kappa1 = PetscSqrtReal(kappa1);

251:       if (kappa0!=0.0 && kappa1!=0.0) {
252:         if (kappaA<0.7*kappa0*kappa1) {
253:           ghat = (kappaA<0.0) ?  -0.7*kappa0/kappa1 : 0.7*kappa0/kappa1;
254:         } else {
255:           ghat = kappaA/(kappa1*kappa1);
256:         }
257:         for (i=0; i<=bcgsl->ell; i++) {
258:           AY0c[i] = AY0c[i] - ghat* AYlc[i];
259:         }
260:       }
261:     }

263:     omega = AY0c[bcgsl->ell];
264:     for (h=bcgsl->ell; h>0 && omega==0.0; h--) omega = AY0c[h];
265:     if (omega==0.0) {
266:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
267:       return(0);
268:     }

270:     VecMAXPY(VX, bcgsl->ell,AY0c+1, VVR);
271:     for (i=1; i<=bcgsl->ell; i++) AY0c[i] *= -1.0;
272:     VecMAXPY(VVU[0], bcgsl->ell,AY0c+1, VVU+1);
273:     VecMAXPY(VVR[0], bcgsl->ell,AY0c+1, VVR+1);
274:     for (i=1; i<=bcgsl->ell; i++) AY0c[i] *= -1.0;
275:     VecNorm(VVR[0], NORM_2, &zeta);
276:     KSPCheckNorm(ksp,zeta);

278:     /* Accurate Update */
279:     if (bcgsl->delta>0.0) {
280:       if (rnmax_computed<zeta) rnmax_computed = zeta;
281:       if (rnmax_true<zeta) rnmax_true = zeta;

283:       bUpdateX = (PetscBool) (zeta<bcgsl->delta*zeta0 && zeta0<=rnmax_computed);
284:       if ((zeta<bcgsl->delta*rnmax_true && zeta0<=rnmax_true) || bUpdateX) {
285:         /* r0 <- b-inv(K)*A*X */
286:         KSP_PCApplyBAorAB(ksp, VX, VVR[0], VTM);
287:         VecAYPX(VVR[0], -1.0, VB);
288:         rnmax_true = zeta;

290:         if (bUpdateX) {
291:           VecAXPY(VXR,1.0,VX);
292:           VecSet(VX,0.0);
293:           VecCopy(VVR[0], VB);
294:           rnmax_computed = zeta;
295:         }
296:       }
297:     }
298:   }
299:   if (bcgsl->delta>0.0) {
300:     VecAXPY(VX,1.0,VXR);
301:   }

303:   ksp->its = k;
304:   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta;
305:   else ksp->rnorm = 0.0;
306:   KSPMonitor(ksp, ksp->its, ksp->rnorm);
307:   KSPLogResidualHistory(ksp, ksp->rnorm);
308:   (*ksp->converged)(ksp, k, ksp->rnorm, &ksp->reason, ksp->cnvP);
309:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
310:   return(0);
311: }

313: /*@
314:    KSPBCGSLSetXRes - Sets the parameter governing when
315:    exact residuals will be used instead of computed residuals.

317:    Logically Collective on ksp

319:    Input Parameters:
320: +  ksp - iterative context obtained from KSPCreate
321: -  delta - computed residuals are used alone when delta is not positive

323:    Options Database Keys:

325: .  -ksp_bcgsl_xres delta

327:    Level: intermediate

329: .seealso: KSPBCGSLSetEll(), KSPBCGSLSetPol(), KSP
330: @*/
331: PetscErrorCode  KSPBCGSLSetXRes(KSP ksp, PetscReal delta)
332: {
333:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;

338:   if (ksp->setupstage) {
339:     if ((delta<=0 && bcgsl->delta>0) || (delta>0 && bcgsl->delta<=0)) {
340:       VecDestroyVecs(ksp->nwork,&ksp->work);
341:       PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
342:       PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);
343:       ksp->setupstage = KSP_SETUP_NEW;
344:     }
345:   }
346:   bcgsl->delta = delta;
347:   return(0);
348: }

350: /*@
351:    KSPBCGSLSetUsePseudoinverse - Use pseudoinverse (via SVD) to solve polynomial part of update

353:    Logically Collective on ksp

355:    Input Parameters:
356: +  ksp - iterative context obtained from KSPCreate
357: -  use_pinv - set to PETSC_TRUE when using pseudoinverse

359:    Options Database Keys:

361: .  -ksp_bcgsl_pinv - use pseudoinverse

363:    Level: intermediate

365: .seealso: KSPBCGSLSetEll(), KSP
366: @*/
367: PetscErrorCode KSPBCGSLSetUsePseudoinverse(KSP ksp,PetscBool use_pinv)
368: {
369:   KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;

372:   bcgsl->pinv = use_pinv;
373:   return(0);
374: }

376: /*@
377:    KSPBCGSLSetPol - Sets the type of polynomial part will
378:    be used in the BiCGSTab(L) solver.

380:    Logically Collective on ksp

382:    Input Parameters:
383: +  ksp - iterative context obtained from KSPCreate
384: -  uMROR - set to PETSC_TRUE when the polynomial is a convex combination of an MR and an OR step.

386:    Options Database Keys:

388: +  -ksp_bcgsl_cxpoly - use enhanced polynomial
389: -  -ksp_bcgsl_mrpoly - use standard polynomial

391:    Level: intermediate

393: .seealso: KSP, KSPBCGSL, KSPCreate(), KSPSetType()
394: @*/
395: PetscErrorCode  KSPBCGSLSetPol(KSP ksp, PetscBool uMROR)
396: {
397:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;


403:   if (!ksp->setupstage) bcgsl->bConvex = uMROR;
404:   else if (bcgsl->bConvex != uMROR) {
405:     /* free the data structures,
406:        then create them again
407:      */
408:     VecDestroyVecs(ksp->nwork,&ksp->work);
409:     PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
410:     PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);

412:     bcgsl->bConvex  = uMROR;
413:     ksp->setupstage = KSP_SETUP_NEW;
414:   }
415:   return(0);
416: }

418: /*@
419:    KSPBCGSLSetEll - Sets the number of search directions in BiCGStab(L).

421:    Logically Collective on ksp

423:    Input Parameters:
424: +  ksp - iterative context obtained from KSPCreate
425: -  ell - number of search directions

427:    Options Database Keys:

429: .  -ksp_bcgsl_ell ell

431:    Level: intermediate

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

438: .seealso: KSPBCGSLSetUsePseudoinverse(), KSP, KSPBCGSL
439: @*/
440: PetscErrorCode  KSPBCGSLSetEll(KSP ksp, PetscInt ell)
441: {
442:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;

446:   if (ell < 1) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "KSPBCGSLSetEll: second argument must be positive");

449:   if (!ksp->setupstage) bcgsl->ell = ell;
450:   else if (bcgsl->ell != ell) {
451:     /* free the data structures, then create them again */
452:     VecDestroyVecs(ksp->nwork,&ksp->work);
453:     PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
454:     PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);

456:     bcgsl->ell      = ell;
457:     ksp->setupstage = KSP_SETUP_NEW;
458:   }
459:   return(0);
460: }

462: PetscErrorCode KSPView_BCGSL(KSP ksp, PetscViewer viewer)
463: {
464:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
466:   PetscBool      isascii;

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

471:   if (isascii) {
472:     PetscViewerASCIIPrintf(viewer, "  Ell = %D\n", bcgsl->ell);
473:     PetscViewerASCIIPrintf(viewer, "  Delta = %lg\n", bcgsl->delta);
474:   }
475:   return(0);
476: }

478: PetscErrorCode KSPSetFromOptions_BCGSL(PetscOptionItems *PetscOptionsObject,KSP ksp)
479: {
480:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
482:   PetscInt       this_ell;
483:   PetscReal      delta;
484:   PetscBool      flga = PETSC_FALSE, flg;

487:   /* PetscOptionsBegin/End are called in KSPSetFromOptions. They
488:      don't need to be called here.
489:   */
490:   PetscOptionsHead(PetscOptionsObject,"KSP BiCGStab(L) Options");

492:   /* Set number of search directions */
493:   PetscOptionsInt("-ksp_bcgsl_ell","Number of Krylov search directions","KSPBCGSLSetEll",bcgsl->ell,&this_ell,&flg);
494:   if (flg) {
495:     KSPBCGSLSetEll(ksp, this_ell);
496:   }

498:   /* Set polynomial type */
499:   PetscOptionsBool("-ksp_bcgsl_cxpoly", "Polynomial part of BiCGStabL is MinRes + OR", "KSPBCGSLSetPol", flga,&flga,NULL);
500:   if (flga) {
501:     KSPBCGSLSetPol(ksp, PETSC_TRUE);
502:   } else {
503:     flg  = PETSC_FALSE;
504:     PetscOptionsBool("-ksp_bcgsl_mrpoly", "Polynomial part of BiCGStabL is MinRes", "KSPBCGSLSetPol", flg,&flg,NULL);
505:     KSPBCGSLSetPol(ksp, PETSC_FALSE);
506:   }

508:   /* Will computed residual be refreshed? */
509:   PetscOptionsReal("-ksp_bcgsl_xres", "Threshold used to decide when to refresh computed residuals", "KSPBCGSLSetXRes", bcgsl->delta, &delta, &flg);
510:   if (flg) {
511:     KSPBCGSLSetXRes(ksp, delta);
512:   }

514:   /* Use pseudoinverse? */
515:   flg  = bcgsl->pinv;
516:   PetscOptionsBool("-ksp_bcgsl_pinv", "Polynomial correction via pseudoinverse", "KSPBCGSLSetUsePseudoinverse",flg,&flg,NULL);
517:   KSPBCGSLSetUsePseudoinverse(ksp,flg);
518:   PetscOptionsTail();
519:   return(0);
520: }

522: PetscErrorCode KSPSetUp_BCGSL(KSP ksp)
523: {
524:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
525:   PetscInt       ell    = bcgsl->ell,ldMZ = ell+1;

529:   KSPSetWorkVecs(ksp, 6+2*ell);
530:   PetscMalloc5(ldMZ,&AY0c,ldMZ,&AYlc,ldMZ,&AYtc,ldMZ*ldMZ,&MZa,ldMZ*ldMZ,&MZb);
531:   PetscBLASIntCast(5*ell,&bcgsl->lwork);
532:   PetscMalloc5(bcgsl->lwork,&bcgsl->work,ell,&bcgsl->s,ell*ell,&bcgsl->u,ell*ell,&bcgsl->v,5*ell,&bcgsl->realwork);
533:   return(0);
534: }

536: PetscErrorCode KSPReset_BCGSL(KSP ksp)
537: {
538:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;

542:   VecDestroyVecs(ksp->nwork,&ksp->work);
543:   PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
544:   PetscFree5(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v,bcgsl->realwork);
545:   return(0);
546: }

548: PetscErrorCode KSPDestroy_BCGSL(KSP ksp)
549: {

553:   KSPReset_BCGSL(ksp);
554:   KSPDestroyDefault(ksp);
555:   return(0);
556: }

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

565:     References:
566: +     1. - G.L.G. Sleijpen, H.A. van der Vorst, "An overview of
567:          approaches for the stable computation of hybrid BiCG
568:          methods", Applied Numerical Mathematics: Transactions
569:          f IMACS, 19(3), 1996.
570: .     2. -  G.L.G. Sleijpen, H.A. van der Vorst, D.R. Fokkema,
571:          "BiCGStab(L) and other hybrid BiCG methods",
572:           Numerical Algorithms, 7, 1994.
573: -     3. -  D.R. Fokkema, "Enhanced implementation of BiCGStab(L)
574:          for solving linear systems of equations", preprint
575:          from www.citeseer.com.

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

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

586:    Level: beginner

588: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPFGMRES, KSPBCGS, KSPSetPCSide(), KSPBCGSLSetEll(), KSPBCGSLSetXRes()

590: M*/
591: PETSC_EXTERN PetscErrorCode KSPCreate_BCGSL(KSP ksp)
592: {
594:   KSP_BCGSL      *bcgsl;

597:   /* allocate BiCGStab(L) context */
598:   PetscNewLog(ksp,&bcgsl);
599:   ksp->data = (void*)bcgsl;

601:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
602:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
603:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);

605:   ksp->ops->setup          = KSPSetUp_BCGSL;
606:   ksp->ops->solve          = KSPSolve_BCGSL;
607:   ksp->ops->reset          = KSPReset_BCGSL;
608:   ksp->ops->destroy        = KSPDestroy_BCGSL;
609:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
610:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
611:   ksp->ops->setfromoptions = KSPSetFromOptions_BCGSL;
612:   ksp->ops->view           = KSPView_BCGSL;

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

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

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

622:   /* Set the threshold for when exact residuals will be used */
623:   bcgsl->delta = 0.0;
624:   return(0);
625: }