Actual source code: bcgsl.c

petsc-3.13.3 2020-07-01
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:   KSPCheckNorm(ksp,zeta0);
 50:   rnmax_computed = zeta0;
 51:   rnmax_true     = zeta0;

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

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

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

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

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

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

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

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

 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);

142:         ksp->its   = k+j;
143:         ksp->rnorm = nrm0;

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

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

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

167:       AY0c[0] = -1;
168:       if (bcgsl->pinv) {
169: #  if defined(PETSC_USE_COMPLEX)
170:         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));
171: #  else
172:         PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&bell,&bell,&MZa[1+ldMZ],&ldMZ,bcgsl->s,bcgsl->u,&bell,bcgsl->v,&bell,bcgsl->work,&bcgsl->lwork,&bierr));
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:         PetscArrayzero(&AY0c[1],bell);
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:         PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("Lower", &bell, &MZa[1+ldMZ], &ldMZ, &bierr));
202:         if (bierr!=0) {
203:           ksp->reason = KSP_DIVERGED_BREAKDOWN;
204:           return(0);
205:         }
206:         PetscArraycpy(&AY0c[1],&MZb[1],bcgsl->ell);
207:         PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &bell, &ione, &MZa[1+ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
208:       }
209:     } else {
210:       PetscBLASInt ione = 1;
211:       PetscScalar  aone = 1.0, azero = 0.0;
212:       PetscBLASInt neqs;
213:       PetscBLASIntCast(bcgsl->ell-1,&neqs);

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

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

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

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

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

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

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

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

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

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

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

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


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->converged)(ksp, k, zeta, &ksp->reason, ksp->cnvP);
302:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
303:   return(0);
304: }

306: /*@
307:    KSPBCGSLSetXRes - Sets the parameter governing when
308:    exact residuals will be used instead of computed residuals.

310:    Logically Collective on ksp

312:    Input Parameters:
313: +  ksp - iterative context obtained from KSPCreate
314: -  delta - computed residuals are used alone when delta is not positive

316:    Options Database Keys:

318: .  -ksp_bcgsl_xres delta

320:    Level: intermediate

322: .seealso: KSPBCGSLSetEll(), KSPBCGSLSetPol(), KSP
323: @*/
324: PetscErrorCode  KSPBCGSLSetXRes(KSP ksp, PetscReal delta)
325: {
326:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;

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

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

346:    Logically Collective on ksp

348:    Input Parameters:
349: +  ksp - iterative context obtained from KSPCreate
350: -  use_pinv - set to PETSC_TRUE when using pseudoinverse

352:    Options Database Keys:

354: .  -ksp_bcgsl_pinv - use pseudoinverse

356:    Level: intermediate

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

365:   bcgsl->pinv = use_pinv;
366:   return(0);
367: }

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

373:    Logically Collective on ksp

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

379:    Options Database Keys:

381: +  -ksp_bcgsl_cxpoly - use enhanced polynomial
382: -  -ksp_bcgsl_mrpoly - use standard polynomial

384:    Level: intermediate

386: .seealso: KSP, KSPBCGSL, KSPCreate(), KSPSetType()
387: @*/
388: PetscErrorCode  KSPBCGSLSetPol(KSP ksp, PetscBool uMROR)
389: {
390:   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

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;

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

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

449:     bcgsl->ell      = ell;
450:     ksp->setupstage = KSP_SETUP_NEW;
451:   }
452:   return(0);
453: }

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

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

464:   if (isascii) {
465:     PetscViewerASCIIPrintf(viewer, "  Ell = %D\n", bcgsl->ell);
466:     PetscViewerASCIIPrintf(viewer, "  Delta = %lg\n", bcgsl->delta);
467:   }
468:   return(0);
469: }

471: PetscErrorCode KSPSetFromOptions_BCGSL(PetscOptionItems *PetscOptionsObject,KSP ksp)
472: {
473:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
475:   PetscInt       this_ell;
476:   PetscReal      delta;
477:   PetscBool      flga = PETSC_FALSE, flg;

480:   /* PetscOptionsBegin/End are called in KSPSetFromOptions. They
481:      don't need to be called here.
482:   */
483:   PetscOptionsHead(PetscOptionsObject,"KSP BiCGStab(L) Options");

485:   /* Set number of search directions */
486:   PetscOptionsInt("-ksp_bcgsl_ell","Number of Krylov search directions","KSPBCGSLSetEll",bcgsl->ell,&this_ell,&flg);
487:   if (flg) {
488:     KSPBCGSLSetEll(ksp, this_ell);
489:   }

491:   /* Set polynomial type */
492:   PetscOptionsBool("-ksp_bcgsl_cxpoly", "Polynomial part of BiCGStabL is MinRes + OR", "KSPBCGSLSetPol", flga,&flga,NULL);
493:   if (flga) {
494:     KSPBCGSLSetPol(ksp, PETSC_TRUE);
495:   } else {
496:     flg  = PETSC_FALSE;
497:     PetscOptionsBool("-ksp_bcgsl_mrpoly", "Polynomial part of BiCGStabL is MinRes", "KSPBCGSLSetPol", flg,&flg,NULL);
498:     KSPBCGSLSetPol(ksp, PETSC_FALSE);
499:   }

501:   /* Will computed residual be refreshed? */
502:   PetscOptionsReal("-ksp_bcgsl_xres", "Threshold used to decide when to refresh computed residuals", "KSPBCGSLSetXRes", bcgsl->delta, &delta, &flg);
503:   if (flg) {
504:     KSPBCGSLSetXRes(ksp, delta);
505:   }

507:   /* Use pseudoinverse? */
508:   flg  = bcgsl->pinv;
509:   PetscOptionsBool("-ksp_bcgsl_pinv", "Polynomial correction via pseudoinverse", "KSPBCGSLSetUsePseudoinverse",flg,&flg,NULL);
510:   KSPBCGSLSetUsePseudoinverse(ksp,flg);
511:   PetscOptionsTail();
512:   return(0);
513: }

515: PetscErrorCode KSPSetUp_BCGSL(KSP ksp)
516: {
517:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
518:   PetscInt       ell    = bcgsl->ell,ldMZ = ell+1;

522:   KSPSetWorkVecs(ksp, 6+2*ell);
523:   PetscMalloc5(ldMZ,&AY0c,ldMZ,&AYlc,ldMZ,&AYtc,ldMZ*ldMZ,&MZa,ldMZ*ldMZ,&MZb);
524:   PetscBLASIntCast(5*ell,&bcgsl->lwork);
525:   PetscMalloc5(bcgsl->lwork,&bcgsl->work,ell,&bcgsl->s,ell*ell,&bcgsl->u,ell*ell,&bcgsl->v,5*ell,&bcgsl->realwork);
526:   return(0);
527: }

529: PetscErrorCode KSPReset_BCGSL(KSP ksp)
530: {
531:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;

535:   VecDestroyVecs(ksp->nwork,&ksp->work);
536:   PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
537:   PetscFree5(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v,bcgsl->realwork);
538:   return(0);
539: }

541: PetscErrorCode KSPDestroy_BCGSL(KSP ksp)
542: {

546:   KSPReset_BCGSL(ksp);
547:   KSPDestroyDefault(ksp);
548:   return(0);
549: }

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

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

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

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

579:    Level: beginner

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

583: M*/
584: PETSC_EXTERN PetscErrorCode KSPCreate_BCGSL(KSP ksp)
585: {
587:   KSP_BCGSL      *bcgsl;

590:   /* allocate BiCGStab(L) context */
591:   PetscNewLog(ksp,&bcgsl);
592:   ksp->data = (void*)bcgsl;

594:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
595:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);

597:   ksp->ops->setup          = KSPSetUp_BCGSL;
598:   ksp->ops->solve          = KSPSolve_BCGSL;
599:   ksp->ops->reset          = KSPReset_BCGSL;
600:   ksp->ops->destroy        = KSPDestroy_BCGSL;
601:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
602:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
603:   ksp->ops->setfromoptions = KSPSetFromOptions_BCGSL;
604:   ksp->ops->view           = KSPView_BCGSL;

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

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

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

614:   /* Set the threshold for when exact residuals will be used */
615:   bcgsl->delta = 0.0;
616:   return(0);
617: }