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>

     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;

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

     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:     if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta;
     83:     else ksp->rnorm = 0.0;

     85:     KSPLogResidualHistory(ksp, ksp->rnorm);
     86:     KSPMonitor(ksp, ksp->its, ksp->rnorm);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    233:       AYlc[0]          = 0.;
    234:       AYlc[bcgsl->ell] = -1;

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

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

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

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

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

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

    250:       if (kappa1<0) kappa1 = -kappa1;
    251:       kappa1 = PetscSqrtReal(kappa1);

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

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

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

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

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

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

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

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

    320:    Logically Collective on ksp

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

    326:    Options Database Keys:

    328: .  -ksp_bcgsl_xres delta

    330:    Level: intermediate

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

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

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

    356:    Logically Collective on ksp

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

    362:    Options Database Keys:

    364: .  -ksp_bcgsl_pinv - use pseudoinverse

    366:    Level: intermediate

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

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

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

    383:    Logically Collective on ksp

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

    389:    Options Database Keys:

    391: +  -ksp_bcgsl_cxpoly - use enhanced polynomial
    392: -  -ksp_bcgsl_mrpoly - use standard polynomial

    394:    Level: intermediate

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

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

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

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

    424:    Logically Collective on ksp

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

    430:    Options Database Keys:

    432: .  -ksp_bcgsl_ell ell

    434:    Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    551: PetscErrorCode KSPDestroy_BCGSL(KSP ksp)
    552: {

    556:   KSPReset_BCGSL(ksp);
    557:   KSPDestroyDefault(ksp);
    558:   return(0);
    559: }

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

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

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

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

    589:    Level: beginner

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

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

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

    604:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
    605:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
    606:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);

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

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

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

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

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