Actual source code: bcgsl.c

petsc-3.8.3 2017-12-09
Report Typos and Errors
  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>


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

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

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

 52:   (*ksp->converged)(ksp, 0, zeta0, &ksp->reason, ksp->cnvP);
 53:   if (ksp->reason) {
 54:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
 55:     ksp->its   = 0;
 56:     ksp->rnorm = zeta0;
 57:     PetscObjectSAWsGrantAccess((PetscObject)ksp);
 58:     return(0);
 59:   }

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

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

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

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

 79:   for (k=0; k<maxit; k += bcgsl->ell) {
 80:     ksp->its   = k;
 81:     ksp->rnorm = zeta;

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

 86:     (*ksp->converged)(ksp, k, zeta, &ksp->reason, ksp->cnvP);
 87:     if (ksp->reason < 0) return(0);
 88:     else if (ksp->reason) break;

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

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

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

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

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

127:       VecNorm(VVR[0], NORM_2, &nrm0);
128:       if (bcgsl->delta>0.0) {
129:         if (rnmax_computed<nrm0) rnmax_computed = nrm0;
130:         if (rnmax_true<nrm0) rnmax_true = nrm0;
131:       }

133:       /* NEW: check for early exit */
134:       (*ksp->converged)(ksp, k+j, nrm0, &ksp->reason, ksp->cnvP);
135:       if (ksp->reason) {
136:         PetscObjectSAWsTakeAccess((PetscObject)ksp);

138:         ksp->its   = k+j;
139:         ksp->rnorm = nrm0;

141:         PetscObjectSAWsGrantAccess((PetscObject)ksp);
142:         if (ksp->reason < 0) return(0);
143:       }
144:     }

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

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

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

195:             for (j=0; j<bell; j++) {
196:               AY0c[1+j] += utb/bcgsl->s[i]*bcgsl->v[j*bell+i];
197:             }
198:           }
199:         }
200:       } else {
201: #if defined(PETSC_MISSING_LAPACK_POTRF)
202:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
203: #else
204:         PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("Lower", &bell, &MZa[1+ldMZ], &ldMZ, &bierr));
205: #endif
206:         if (bierr!=0) {
207:           ksp->reason = KSP_DIVERGED_BREAKDOWN;
208:           return(0);
209:         }
210:         PetscMemcpy(&AY0c[1],&MZb[1],bcgsl->ell*sizeof(PetscScalar));
211:         PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &bell, &ione, &MZa[1+ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
212:       }
213:     } else {
214:       PetscBLASInt ione = 1;
215:       PetscScalar  aone = 1.0, azero = 0.0;
216:       PetscBLASInt neqs;
217:       PetscBLASIntCast(bcgsl->ell-1,&neqs);

219: #if defined(PETSC_MISSING_LAPACK_POTRF)
220:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
221: #else
222:       PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("Lower", &neqs, &MZa[1+ldMZ], &ldMZ, &bierr));
223: #endif
224:       if (bierr!=0) {
225:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
226:         return(0);
227:       }
228:       PetscMemcpy(&AY0c[1],&MZb[1],(bcgsl->ell-1)*sizeof(PetscScalar));
229:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1+ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
230:       AY0c[0]          = -1;
231:       AY0c[bcgsl->ell] = 0.;

233:       PetscMemcpy(&AYlc[1],&MZb[1+ldMZ*(bcgsl->ell)],(bcgsl->ell-1)*sizeof(PetscScalar));
234:       PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1+ldMZ], &ldMZ, &AYlc[1], &ldMZ, &bierr));

236:       AYlc[0]          = 0.;
237:       AYlc[bcgsl->ell] = -1;

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

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

243:       /* round-off can cause negative kappa's */
244:       if (kappa0<0) kappa0 = -kappa0;
245:       kappa0 = PetscSqrtReal(kappa0);

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

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

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

253:       if (kappa1<0) kappa1 = -kappa1;
254:       kappa1 = PetscSqrtReal(kappa1);

256:       if (kappa0!=0.0 && kappa1!=0.0) {
257:         if (kappaA<0.7*kappa0*kappa1) {
258:           ghat = (kappaA<0.0) ?  -0.7*kappa0/kappa1 : 0.7*kappa0/kappa1;
259:         } else {
260:           ghat = kappaA/(kappa1*kappa1);
261:         }
262:         for (i=0; i<=bcgsl->ell; i++) {
263:           AY0c[i] = AY0c[i] - ghat* AYlc[i];
264:         }
265:       }
266:     }

268:     omega = AY0c[bcgsl->ell];
269:     for (h=bcgsl->ell; h>0 && omega==0.0; h--) omega = AY0c[h];
270:     if (omega==0.0) {
271:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
272:       return(0);
273:     }


276:     VecMAXPY(VX, bcgsl->ell,AY0c+1, VVR);
277:     for (i=1; i<=bcgsl->ell; i++) AY0c[i] *= -1.0;
278:     VecMAXPY(VVU[0], bcgsl->ell,AY0c+1, VVU+1);
279:     VecMAXPY(VVR[0], bcgsl->ell,AY0c+1, VVR+1);
280:     for (i=1; i<=bcgsl->ell; i++) AY0c[i] *= -1.0;
281:     VecNorm(VVR[0], NORM_2, &zeta);

283:     /* Accurate Update */
284:     if (bcgsl->delta>0.0) {
285:       if (rnmax_computed<zeta) rnmax_computed = zeta;
286:       if (rnmax_true<zeta) rnmax_true = zeta;

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

295:         if (bUpdateX) {
296:           VecAXPY(VXR,1.0,VX);
297:           VecSet(VX,0.0);
298:           VecCopy(VVR[0], VB);
299:           rnmax_computed = zeta;
300:         }
301:       }
302:     }
303:   }
304:   if (bcgsl->delta>0.0) {
305:     VecAXPY(VX,1.0,VXR);
306:   }

308:   (*ksp->converged)(ksp, k, zeta, &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: .keywords: KSP, BiCGStab(L), set, exact residuals

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

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

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

355:    Logically Collective on KSP

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

361:    Options Database Keys:

363: +  -ksp_bcgsl_pinv - use pseudoinverse

365:    Level: intermediate

367: .keywords: KSP, BiCGStab(L), set, polynomial

369: .seealso: KSPBCGSLSetEll()
370: @*/
371: PetscErrorCode KSPBCGSLSetUsePseudoinverse(KSP ksp,PetscBool use_pinv)
372: {
373:   KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;

376:   bcgsl->pinv = use_pinv;
377:   return(0);
378: }

380: /*@
381:    KSPBCGSLSetPol - Sets the type of polynomial part will
382:    be used in the BiCGSTab(L) solver.

384:    Logically Collective on KSP

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

390:    Options Database Keys:

392: +  -ksp_bcgsl_cxpoly - use enhanced polynomial
393: .  -ksp_bcgsl_mrpoly - use standard polynomial

395:    Level: intermediate

397: .keywords: KSP, BiCGStab(L), set, polynomial

399: .seealso: @()
400: @*/
401: PetscErrorCode  KSPBCGSLSetPol(KSP ksp, PetscBool uMROR)
402: {
403:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;


409:   if (!ksp->setupstage) bcgsl->bConvex = uMROR;
410:   else if (bcgsl->bConvex != uMROR) {
411:     /* free the data structures,
412:        then create them again
413:      */
414:     VecDestroyVecs(ksp->nwork,&ksp->work);
415:     PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
416:     PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);

418:     bcgsl->bConvex  = uMROR;
419:     ksp->setupstage = KSP_SETUP_NEW;
420:   }
421:   return(0);
422: }

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

427:    Logically Collective on KSP

429:    Input Parameters:
430: +  ksp - iterative context obtained from KSPCreate
431: -  ell - number of search directions

433:    Options Database Keys:

435: .  -ksp_bcgsl_ell ell

437:    Level: intermediate

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

444: .keywords: KSP, BiCGStab(L), set, exact residuals,

446: .seealso: KSPBCGSLSetUsePseudoinverse()
447: @*/
448: PetscErrorCode  KSPBCGSLSetEll(KSP ksp, PetscInt ell)
449: {
450:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;

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

457:   if (!ksp->setupstage) bcgsl->ell = ell;
458:   else if (bcgsl->ell != ell) {
459:     /* free the data structures, then create them again */
460:     VecDestroyVecs(ksp->nwork,&ksp->work);
461:     PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
462:     PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);

464:     bcgsl->ell      = ell;
465:     ksp->setupstage = KSP_SETUP_NEW;
466:   }
467:   return(0);
468: }

470: PetscErrorCode KSPView_BCGSL(KSP ksp, PetscViewer viewer)
471: {
472:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
474:   PetscBool      isascii;

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

479:   if (isascii) {
480:     PetscViewerASCIIPrintf(viewer, "  Ell = %D\n", bcgsl->ell);
481:     PetscViewerASCIIPrintf(viewer, "  Delta = %lg\n", bcgsl->delta);
482:   }
483:   return(0);
484: }

486: PetscErrorCode KSPSetFromOptions_BCGSL(PetscOptionItems *PetscOptionsObject,KSP ksp)
487: {
488:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
490:   PetscInt       this_ell;
491:   PetscReal      delta;
492:   PetscBool      flga = PETSC_FALSE, flg;

495:   /* PetscOptionsBegin/End are called in KSPSetFromOptions. They
496:      don't need to be called here.
497:   */
498:   PetscOptionsHead(PetscOptionsObject,"KSP BiCGStab(L) Options");

500:   /* Set number of search directions */
501:   PetscOptionsInt("-ksp_bcgsl_ell","Number of Krylov search directions","KSPBCGSLSetEll",bcgsl->ell,&this_ell,&flg);
502:   if (flg) {
503:     KSPBCGSLSetEll(ksp, this_ell);
504:   }

506:   /* Set polynomial type */
507:   PetscOptionsBool("-ksp_bcgsl_cxpoly", "Polynomial part of BiCGStabL is MinRes + OR", "KSPBCGSLSetPol", flga,&flga,NULL);
508:   if (flga) {
509:     KSPBCGSLSetPol(ksp, PETSC_TRUE);
510:   } else {
511:     flg  = PETSC_FALSE;
512:     PetscOptionsBool("-ksp_bcgsl_mrpoly", "Polynomial part of BiCGStabL is MinRes", "KSPBCGSLSetPol", flg,&flg,NULL);
513:     KSPBCGSLSetPol(ksp, PETSC_FALSE);
514:   }

516:   /* Will computed residual be refreshed? */
517:   PetscOptionsReal("-ksp_bcgsl_xres", "Threshold used to decide when to refresh computed residuals", "KSPBCGSLSetXRes", bcgsl->delta, &delta, &flg);
518:   if (flg) {
519:     KSPBCGSLSetXRes(ksp, delta);
520:   }

522:   /* Use pseudoinverse? */
523:   flg  = bcgsl->pinv;
524:   PetscOptionsBool("-ksp_bcgsl_pinv", "Polynomial correction via pseudoinverse", "KSPBCGSLSetUsePseudoinverse",flg,&flg,NULL);
525:   KSPBCGSLSetUsePseudoinverse(ksp,flg);
526:   PetscOptionsTail();
527:   return(0);
528: }

530: PetscErrorCode KSPSetUp_BCGSL(KSP ksp)
531: {
532:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
533:   PetscInt       ell    = bcgsl->ell,ldMZ = ell+1;

537:   KSPSetWorkVecs(ksp, 6+2*ell);
538:   PetscMalloc5(ldMZ,&AY0c,ldMZ,&AYlc,ldMZ,&AYtc,ldMZ*ldMZ,&MZa,ldMZ*ldMZ,&MZb);
539:   PetscBLASIntCast(5*ell,&bcgsl->lwork);
540:   PetscMalloc5(bcgsl->lwork,&bcgsl->work,ell,&bcgsl->s,ell*ell,&bcgsl->u,ell*ell,&bcgsl->v,5*ell,&bcgsl->realwork);
541:   return(0);
542: }

544: PetscErrorCode KSPReset_BCGSL(KSP ksp)
545: {
546:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;

550:   VecDestroyVecs(ksp->nwork,&ksp->work);
551:   PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
552:   PetscFree5(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v,bcgsl->realwork);
553:   return(0);
554: }

556: PetscErrorCode KSPDestroy_BCGSL(KSP ksp)
557: {

561:   KSPReset_BCGSL(ksp);
562:   KSPDestroyDefault(ksp);
563:   return(0);
564: }

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

573:     References:
574: +     1. - G.L.G. Sleijpen, H.A. van der Vorst, "An overview of
575:          approaches for the stable computation of hybrid BiCG
576:          methods", Applied Numerical Mathematics: Transactions
577:          f IMACS, 19(3), 1996.
578: .     2. -  G.L.G. Sleijpen, H.A. van der Vorst, D.R. Fokkema,
579:          "BiCGStab(L) and other hybrid BiCG methods",
580:           Numerical Algorithms, 7, 1994.
581: -     3. -  D.R. Fokkema, "Enhanced implementation of BiCGStab(L)
582:          for solving linear systems of equations", preprint
583:          from www.citeseer.com.

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

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

594:    Level: beginner

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

598: M*/
599: PETSC_EXTERN PetscErrorCode KSPCreate_BCGSL(KSP ksp)
600: {
602:   KSP_BCGSL      *bcgsl;

605:   /* allocate BiCGStab(L) context */
606:   PetscNewLog(ksp,&bcgsl);
607:   ksp->data = (void*)bcgsl;

609:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
610:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);

612:   ksp->ops->setup          = KSPSetUp_BCGSL;
613:   ksp->ops->solve          = KSPSolve_BCGSL;
614:   ksp->ops->reset          = KSPReset_BCGSL;
615:   ksp->ops->destroy        = KSPDestroy_BCGSL;
616:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
617:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
618:   ksp->ops->setfromoptions = KSPSetFromOptions_BCGSL;
619:   ksp->ops->view           = KSPView_BCGSL;

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

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

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

629:   /* Set the threshold for when exact residuals will be used */
630:   bcgsl->delta = 0.0;
631:   return(0);
632: }