Actual source code: bcgsl.c
petsc-3.10.5 2019-03-28
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: BiCGStab(L), set, exact residuals
331: .seealso: KSPBCGSLSetEll(), KSPBCGSLSetPol(), KSP
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: BiCGStab(L), set, polynomial
369: .seealso: KSPBCGSLSetEll(), KSP
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: BiCGStab(L), set, polynomial
399: .seealso: KSP, KSPBCGSL, KSPCreate(), KSPSetType()
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: BiCGStab(L), set, exact residuals,
446: .seealso: KSPBCGSLSetUsePseudoinverse(), KSP, KSPBCGSL
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: }