Actual source code: bcgsl.c

petsc-master 2019-06-17
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:     PetscMemcpy(MZb,MZa,ldMZ*ldMZ*sizeof(PetscScalar));

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_MISSING_LAPACK_GESVD)
170:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable.");
171: #else
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: #endif
178:         if (bierr!=0) {
179:           ksp->reason = KSP_DIVERGED_BREAKDOWN;
180:           return(0);
181:         }
182:         /* Apply pseudo-inverse */
183:         max_s = bcgsl->s[0];
184:         for (i=1; i<bell; i++) {
185:           if (bcgsl->s[i] > max_s) {
186:             max_s = bcgsl->s[i];
187:           }
188:         }
189:         /* tolerance is hardwired to bell*max(s)*PETSC_MACHINE_EPSILON */
190:         pinv_tol = bell*max_s*PETSC_MACHINE_EPSILON;
191:         PetscMemzero(&AY0c[1],bell*sizeof(PetscScalar));
192:         for (i=0; i<bell; i++) {
193:           if (bcgsl->s[i] >= pinv_tol) {
194:             utb=0.;
195:             for (j=0; j<bell; j++) {
196:               utb += MZb[1+j]*bcgsl->u[i*bell+j];
197:             }

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

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

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

240:       AYlc[0]          = 0.;
241:       AYlc[bcgsl->ell] = -1;

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

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

247:       /* round-off can cause negative kappa's */
248:       if (kappa0<0) kappa0 = -kappa0;
249:       kappa0 = PetscSqrtReal(kappa0);

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

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

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

257:       if (kappa1<0) kappa1 = -kappa1;
258:       kappa1 = PetscSqrtReal(kappa1);

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

272:     omega = AY0c[bcgsl->ell];
273:     for (h=bcgsl->ell; h>0 && omega==0.0; h--) omega = AY0c[h];
274:     if (omega==0.0) {
275:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
276:       return(0);
277:     }


280:     VecMAXPY(VX, bcgsl->ell,AY0c+1, VVR);
281:     for (i=1; i<=bcgsl->ell; i++) AY0c[i] *= -1.0;
282:     VecMAXPY(VVU[0], bcgsl->ell,AY0c+1, VVU+1);
283:     VecMAXPY(VVR[0], bcgsl->ell,AY0c+1, VVR+1);
284:     for (i=1; i<=bcgsl->ell; i++) AY0c[i] *= -1.0;
285:     VecNorm(VVR[0], NORM_2, &zeta);
286:     KSPCheckNorm(ksp,zeta);

288:     /* Accurate Update */
289:     if (bcgsl->delta>0.0) {
290:       if (rnmax_computed<zeta) rnmax_computed = zeta;
291:       if (rnmax_true<zeta) rnmax_true = zeta;

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

300:         if (bUpdateX) {
301:           VecAXPY(VXR,1.0,VX);
302:           VecSet(VX,0.0);
303:           VecCopy(VVR[0], VB);
304:           rnmax_computed = zeta;
305:         }
306:       }
307:     }
308:   }
309:   if (bcgsl->delta>0.0) {
310:     VecAXPY(VX,1.0,VXR);
311:   }

313:   (*ksp->converged)(ksp, k, zeta, &ksp->reason, ksp->cnvP);
314:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
315:   return(0);
316: }

318: /*@
319:    KSPBCGSLSetXRes - Sets the parameter governing when
320:    exact residuals will be used instead of computed residuals.

322:    Logically Collective on ksp

324:    Input Parameters:
325: +  ksp - iterative context obtained from KSPCreate
326: -  delta - computed residuals are used alone when delta is not positive

328:    Options Database Keys:

330: .  -ksp_bcgsl_xres delta

332:    Level: intermediate

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

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

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

358:    Logically Collective on ksp

360:    Input Parameters:
361: +  ksp - iterative context obtained from KSPCreate
362: -  use_pinv - set to PETSC_TRUE when using pseudoinverse

364:    Options Database Keys:

366: +  -ksp_bcgsl_pinv - use pseudoinverse

368:    Level: intermediate

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

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

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

385:    Logically Collective on ksp

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

391:    Options Database Keys:

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

396:    Level: intermediate

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


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

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

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

426:    Logically Collective on ksp

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

432:    Options Database Keys:

434: .  -ksp_bcgsl_ell ell

436:    Level: intermediate

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

443: .seealso: KSPBCGSLSetUsePseudoinverse(), KSP, KSPBCGSL
444: @*/
445: PetscErrorCode  KSPBCGSLSetEll(KSP ksp, PetscInt ell)
446: {
447:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;

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

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

461:     bcgsl->ell      = ell;
462:     ksp->setupstage = KSP_SETUP_NEW;
463:   }
464:   return(0);
465: }

467: PetscErrorCode KSPView_BCGSL(KSP ksp, PetscViewer viewer)
468: {
469:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
471:   PetscBool      isascii;

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

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

483: PetscErrorCode KSPSetFromOptions_BCGSL(PetscOptionItems *PetscOptionsObject,KSP ksp)
484: {
485:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
487:   PetscInt       this_ell;
488:   PetscReal      delta;
489:   PetscBool      flga = PETSC_FALSE, flg;

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

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

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

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

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

527: PetscErrorCode KSPSetUp_BCGSL(KSP ksp)
528: {
529:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;
530:   PetscInt       ell    = bcgsl->ell,ldMZ = ell+1;

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

541: PetscErrorCode KSPReset_BCGSL(KSP ksp)
542: {
543:   KSP_BCGSL      *bcgsl = (KSP_BCGSL*)ksp->data;

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

553: PetscErrorCode KSPDestroy_BCGSL(KSP ksp)
554: {

558:   KSPReset_BCGSL(ksp);
559:   KSPDestroyDefault(ksp);
560:   return(0);
561: }

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

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

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

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

591:    Level: beginner

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

595: M*/
596: PETSC_EXTERN PetscErrorCode KSPCreate_BCGSL(KSP ksp)
597: {
599:   KSP_BCGSL      *bcgsl;

602:   /* allocate BiCGStab(L) context */
603:   PetscNewLog(ksp,&bcgsl);
604:   ksp->data = (void*)bcgsl;

606:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
607:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);

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

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

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

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

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