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: }