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;

 32:   /* set up temporary vectors */
 33:   vi         = 0;
 34:   ell        = bcgsl->ell;
 35:   bcgsl->vB  = ksp->work[vi]; vi++;
 36:   bcgsl->vRt = ksp->work[vi]; vi++;
 37:   bcgsl->vTm = ksp->work[vi]; vi++;
 38:   bcgsl->vvR = ksp->work+vi; vi += ell+1;
 39:   bcgsl->vvU = ksp->work+vi; vi += ell+1;
 40:   bcgsl->vXr = ksp->work[vi]; 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) {
 88:         VecAXPY(VX,1.0,VXR);
 89:       }
 90:       return 0;
 91:     }

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

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

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

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

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

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

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

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

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

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

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

196:             for (j=0; j<bell; j++) {
197:               AY0c[1+j] += utb/bcgsl->s[i]*bcgsl->v[j*bell+i];
198:             }
199:           }
200:         }
201:       } else {
202:         PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("Lower", &bell, &MZa[1+ldMZ], &ldMZ, &bierr));
203:         if (bierr!=0) {
204:           ksp->reason = KSP_DIVERGED_BREAKDOWN;
205:           return 0;
206:         }
207:         PetscArraycpy(&AY0c[1],&MZb[1],bcgsl->ell);
208:         PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &bell, &ione, &MZa[1+ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
209:       }
210:     } else {
211:       PetscBLASInt ione = 1;
212:       PetscScalar  aone = 1.0, azero = 0.0;
213:       PetscBLASInt neqs;
214:       PetscBLASIntCast(bcgsl->ell-1,&neqs);

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

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

229:       AYlc[0]          = 0.;
230:       AYlc[bcgsl->ell] = -1;

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

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

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

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

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

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

246:       if (kappa1<0) kappa1 = -kappa1;
247:       kappa1 = PetscSqrtReal(kappa1);

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

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

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

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

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

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

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

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

315:    Logically Collective on ksp

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

321:    Options Database Keys:

323: .  -ksp_bcgsl_xres delta - Threshold used to decide when to refresh computed residuals

325:    Level: intermediate

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

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

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

349:    Logically Collective on ksp

351:    Input Parameters:
352: +  ksp - iterative context obtained from KSPCreate
353: -  use_pinv - set to PETSC_TRUE when using pseudoinverse

355:    Options Database Keys:

357: .  -ksp_bcgsl_pinv - use pseudoinverse

359:    Level: intermediate

361: .seealso: KSPBCGSLSetEll(), KSP
362: @*/
363: PetscErrorCode KSPBCGSLSetUsePseudoinverse(KSP ksp,PetscBool use_pinv)
364: {
365:   KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;

367:   bcgsl->pinv = use_pinv;
368:   return 0;
369: }

371: /*@
372:    KSPBCGSLSetPol - Sets the type of polynomial part will
373:    be used in the BiCGSTab(L) solver.

375:    Logically Collective on ksp

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

381:    Options Database Keys:

383: +  -ksp_bcgsl_cxpoly - use enhanced polynomial
384: -  -ksp_bcgsl_mrpoly - use standard polynomial

386:    Level: intermediate

388: .seealso: KSP, KSPBCGSL, KSPCreate(), KSPSetType()
389: @*/
390: PetscErrorCode  KSPBCGSLSetPol(KSP ksp, PetscBool uMROR)
391: {
392:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;


396:   if (!ksp->setupstage) bcgsl->bConvex = uMROR;
397:   else if (bcgsl->bConvex != uMROR) {
398:     /* free the data structures,
399:        then create them again
400:      */
401:     VecDestroyVecs(ksp->nwork,&ksp->work);
402:     PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
403:     PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);

405:     bcgsl->bConvex  = uMROR;
406:     ksp->setupstage = KSP_SETUP_NEW;
407:   }
408:   return 0;
409: }

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

414:    Logically Collective on ksp

416:    Input Parameters:
417: +  ksp - iterative context obtained from KSPCreate
418: -  ell - number of search directions

420:    Options Database Keys:

422: .  -ksp_bcgsl_ell ell - Number of Krylov search directions

424:    Level: intermediate

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

431: .seealso: KSPBCGSLSetUsePseudoinverse(), KSP, KSPBCGSL
432: @*/
433: PetscErrorCode  KSPBCGSLSetEll(KSP ksp, PetscInt ell)
434: {
435:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;


440:   if (!ksp->setupstage) bcgsl->ell = ell;
441:   else if (bcgsl->ell != ell) {
442:     /* free the data structures, then create them again */
443:     VecDestroyVecs(ksp->nwork,&ksp->work);
444:     PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
445:     PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);

447:     bcgsl->ell      = ell;
448:     ksp->setupstage = KSP_SETUP_NEW;
449:   }
450:   return 0;
451: }

453: PetscErrorCode KSPView_BCGSL(KSP ksp, PetscViewer viewer)
454: {
455:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
456:   PetscBool      isascii;

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

460:   if (isascii) {
461:     PetscViewerASCIIPrintf(viewer, "  Ell = %D\n", bcgsl->ell);
462:     PetscViewerASCIIPrintf(viewer, "  Delta = %lg\n", bcgsl->delta);
463:   }
464:   return 0;
465: }

467: PetscErrorCode KSPSetFromOptions_BCGSL(PetscOptionItems *PetscOptionsObject,KSP ksp)
468: {
469:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
470:   PetscInt       this_ell;
471:   PetscReal      delta;
472:   PetscBool      flga = PETSC_FALSE, flg;

474:   /* PetscOptionsBegin/End are called in KSPSetFromOptions. They
475:      don't need to be called here.
476:   */
477:   PetscOptionsHead(PetscOptionsObject,"KSP BiCGStab(L) Options");

479:   /* Set number of search directions */
480:   PetscOptionsInt("-ksp_bcgsl_ell","Number of Krylov search directions","KSPBCGSLSetEll",bcgsl->ell,&this_ell,&flg);
481:   if (flg) {
482:     KSPBCGSLSetEll(ksp, this_ell);
483:   }

485:   /* Set polynomial type */
486:   PetscOptionsBool("-ksp_bcgsl_cxpoly", "Polynomial part of BiCGStabL is MinRes + OR", "KSPBCGSLSetPol", flga,&flga,NULL);
487:   if (flga) {
488:     KSPBCGSLSetPol(ksp, PETSC_TRUE);
489:   } else {
490:     flg  = PETSC_FALSE;
491:     PetscOptionsBool("-ksp_bcgsl_mrpoly", "Polynomial part of BiCGStabL is MinRes", "KSPBCGSLSetPol", flg,&flg,NULL);
492:     KSPBCGSLSetPol(ksp, PETSC_FALSE);
493:   }

495:   /* Will computed residual be refreshed? */
496:   PetscOptionsReal("-ksp_bcgsl_xres", "Threshold used to decide when to refresh computed residuals", "KSPBCGSLSetXRes", bcgsl->delta, &delta, &flg);
497:   if (flg) {
498:     KSPBCGSLSetXRes(ksp, delta);
499:   }

501:   /* Use pseudoinverse? */
502:   flg  = bcgsl->pinv;
503:   PetscOptionsBool("-ksp_bcgsl_pinv", "Polynomial correction via pseudoinverse", "KSPBCGSLSetUsePseudoinverse",flg,&flg,NULL);
504:   KSPBCGSLSetUsePseudoinverse(ksp,flg);
505:   PetscOptionsTail();
506:   return 0;
507: }

509: PetscErrorCode KSPSetUp_BCGSL(KSP ksp)
510: {
511:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
512:   PetscInt       ell    = bcgsl->ell,ldMZ = ell+1;

514:   KSPSetWorkVecs(ksp, 6+2*ell);
515:   PetscMalloc5(ldMZ,&AY0c,ldMZ,&AYlc,ldMZ,&AYtc,ldMZ*ldMZ,&MZa,ldMZ*ldMZ,&MZb);
516:   PetscBLASIntCast(5*ell,&bcgsl->lwork);
517:   PetscMalloc5(bcgsl->lwork,&bcgsl->work,ell,&bcgsl->s,ell*ell,&bcgsl->u,ell*ell,&bcgsl->v,5*ell,&bcgsl->realwork);
518:   return 0;
519: }

521: PetscErrorCode KSPReset_BCGSL(KSP ksp)
522: {
523:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;

525:   VecDestroyVecs(ksp->nwork,&ksp->work);
526:   PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
527:   PetscFree5(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v,bcgsl->realwork);
528:   return 0;
529: }

531: PetscErrorCode KSPDestroy_BCGSL(KSP ksp)
532: {
533:   KSPReset_BCGSL(ksp);
534:   KSPDestroyDefault(ksp);
535:   return 0;
536: }

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

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:    Contributed by: Joel M. Malard, email jm.malard@pnl.gov

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

566:    Level: beginner

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

570: M*/
571: PETSC_EXTERN PetscErrorCode KSPCreate_BCGSL(KSP ksp)
572: {
573:   KSP_BCGSL      *bcgsl;

575:   /* allocate BiCGStab(L) context */
576:   PetscNewLog(ksp,&bcgsl);
577:   ksp->data = (void*)bcgsl;

579:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
580:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
581:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);

583:   ksp->ops->setup          = KSPSetUp_BCGSL;
584:   ksp->ops->solve          = KSPSolve_BCGSL;
585:   ksp->ops->reset          = KSPReset_BCGSL;
586:   ksp->ops->destroy        = KSPDestroy_BCGSL;
587:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
588:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
589:   ksp->ops->setfromoptions = KSPSetFromOptions_BCGSL;
590:   ksp->ops->view           = KSPView_BCGSL;

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

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

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

600:   /* Set the threshold for when exact residuals will be used */
601:   bcgsl->delta = 0.0;
602:   return 0;
603: }