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;
34: /* set up temporary vectors */
35: vi = 0;
36: ell = bcgsl->ell;
37: bcgsl->vB = ksp->work[vi]; vi++;
38: bcgsl->vRt = ksp->work[vi]; vi++;
39: bcgsl->vTm = ksp->work[vi]; vi++;
40: bcgsl->vvR = ksp->work+vi; vi += ell+1;
41: bcgsl->vvU = ksp->work+vi; vi += ell+1;
42: bcgsl->vXr = ksp->work[vi]; vi++;
43: PetscBLASIntCast(ell+1,&ldMZ);
45: /* Prime the iterative solver */
46: KSPInitialResidual(ksp, VX, VTM, VB, VVR[0], ksp->vec_rhs);
47: VecNorm(VVR[0], NORM_2, &zeta0);
48: KSPCheckNorm(ksp,zeta0);
49: rnmax_computed = zeta0;
50: rnmax_true = zeta0;
52: PetscObjectSAWsTakeAccess((PetscObject)ksp);
53: ksp->its = 0;
54: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta0;
55: else ksp->rnorm = 0.0;
56: PetscObjectSAWsGrantAccess((PetscObject)ksp);
57: (*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP);
58: if (ksp->reason) return(0);
60: VecSet(VVU[0],0.0);
61: alpha = 0.;
62: rho0 = omega = 1;
64: if (bcgsl->delta>0.0) {
65: VecCopy(VX, VXR);
66: VecSet(VX,0.0);
67: VecCopy(VVR[0], VB);
68: } else {
69: VecCopy(ksp->vec_rhs, VB);
70: }
72: /* Life goes on */
73: VecCopy(VVR[0], VRT);
74: zeta = zeta0;
76: KSPGetTolerances(ksp, NULL, NULL, NULL, &maxit);
78: for (k=0; k<maxit; k += bcgsl->ell) {
79: ksp->its = k;
80: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta;
81: else ksp->rnorm = 0.0;
83: KSPLogResidualHistory(ksp, ksp->rnorm);
84: KSPMonitor(ksp, ksp->its, ksp->rnorm);
86: (*ksp->converged)(ksp, k, ksp->rnorm, &ksp->reason, ksp->cnvP);
87: if (ksp->reason < 0) return(0);
88: if (ksp->reason) {
89: if (bcgsl->delta>0.0) {
90: VecAXPY(VX,1.0,VXR);
91: }
92: return(0);
93: }
95: /* BiCG part */
96: rho0 = -omega*rho0;
97: nrm0 = zeta;
98: for (j=0; j<bcgsl->ell; j++) {
99: /* rho1 <- r_j' * r_tilde */
100: VecDot(VVR[j], VRT, &rho1);
101: KSPCheckDot(ksp,rho1);
102: if (rho1 == 0.0) {
103: ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
104: return(0);
105: }
106: beta = alpha*(rho1/rho0);
107: rho0 = rho1;
108: for (i=0; i<=j; i++) {
109: /* u_i <- r_i - beta*u_i */
110: VecAYPX(VVU[i], -beta, VVR[i]);
111: }
112: /* u_{j+1} <- inv(K)*A*u_j */
113: KSP_PCApplyBAorAB(ksp, VVU[j], VVU[j+1], VTM);
115: VecDot(VVU[j+1], VRT, &sigma);
116: KSPCheckDot(ksp,sigma);
117: if (sigma == 0.0) {
118: ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
119: return(0);
120: }
121: alpha = rho1/sigma;
123: /* x <- x + alpha*u_0 */
124: VecAXPY(VX, alpha, VVU[0]);
126: for (i=0; i<=j; i++) {
127: /* r_i <- r_i - alpha*u_{i+1} */
128: VecAXPY(VVR[i], -alpha, VVU[i+1]);
129: }
131: /* r_{j+1} <- inv(K)*A*r_j */
132: KSP_PCApplyBAorAB(ksp, VVR[j], VVR[j+1], VTM);
134: VecNorm(VVR[0], NORM_2, &nrm0);
135: KSPCheckNorm(ksp,nrm0);
136: if (bcgsl->delta>0.0) {
137: if (rnmax_computed<nrm0) rnmax_computed = nrm0;
138: if (rnmax_true<nrm0) rnmax_true = nrm0;
139: }
141: /* NEW: check for early exit */
142: (*ksp->converged)(ksp, k+j, nrm0, &ksp->reason, ksp->cnvP);
143: if (ksp->reason) {
144: PetscObjectSAWsTakeAccess((PetscObject)ksp);
145: ksp->its = k+j;
146: ksp->rnorm = nrm0;
148: PetscObjectSAWsGrantAccess((PetscObject)ksp);
149: if (ksp->reason < 0) return(0);
150: }
151: }
153: /* Polynomial part */
154: for (i = 0; i <= bcgsl->ell; ++i) {
155: VecMDot(VVR[i], i+1, VVR, &MZa[i*ldMZ]);
156: }
157: /* Symmetrize MZa */
158: for (i = 0; i <= bcgsl->ell; ++i) {
159: for (j = i+1; j <= bcgsl->ell; ++j) {
160: MZa[i*ldMZ+j] = MZa[j*ldMZ+i] = PetscConj(MZa[j*ldMZ+i]);
161: }
162: }
163: /* Copy MZa to MZb */
164: PetscArraycpy(MZb,MZa,ldMZ*ldMZ);
166: if (!bcgsl->bConvex || bcgsl->ell==1) {
167: PetscBLASInt ione = 1,bell;
168: PetscBLASIntCast(bcgsl->ell,&bell);
170: AY0c[0] = -1;
171: if (bcgsl->pinv) {
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: if (bierr!=0) {
178: ksp->reason = KSP_DIVERGED_BREAKDOWN;
179: return(0);
180: }
181: /* Apply pseudo-inverse */
182: max_s = bcgsl->s[0];
183: for (i=1; i<bell; i++) {
184: if (bcgsl->s[i] > max_s) {
185: max_s = bcgsl->s[i];
186: }
187: }
188: /* tolerance is hardwired to bell*max(s)*PETSC_MACHINE_EPSILON */
189: pinv_tol = bell*max_s*PETSC_MACHINE_EPSILON;
190: PetscArrayzero(&AY0c[1],bell);
191: for (i=0; i<bell; i++) {
192: if (bcgsl->s[i] >= pinv_tol) {
193: utb=0.;
194: for (j=0; j<bell; j++) {
195: utb += MZb[1+j]*bcgsl->u[i*bell+j];
196: }
198: for (j=0; j<bell; j++) {
199: AY0c[1+j] += utb/bcgsl->s[i]*bcgsl->v[j*bell+i];
200: }
201: }
202: }
203: } else {
204: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("Lower", &bell, &MZa[1+ldMZ], &ldMZ, &bierr));
205: if (bierr!=0) {
206: ksp->reason = KSP_DIVERGED_BREAKDOWN;
207: return(0);
208: }
209: PetscArraycpy(&AY0c[1],&MZb[1],bcgsl->ell);
210: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &bell, &ione, &MZa[1+ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
211: }
212: } else {
213: PetscBLASInt ione = 1;
214: PetscScalar aone = 1.0, azero = 0.0;
215: PetscBLASInt neqs;
216: PetscBLASIntCast(bcgsl->ell-1,&neqs);
218: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("Lower", &neqs, &MZa[1+ldMZ], &ldMZ, &bierr));
219: if (bierr!=0) {
220: ksp->reason = KSP_DIVERGED_BREAKDOWN;
221: return(0);
222: }
223: PetscArraycpy(&AY0c[1],&MZb[1],bcgsl->ell-1);
224: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1+ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
225: AY0c[0] = -1;
226: AY0c[bcgsl->ell] = 0.;
228: PetscArraycpy(&AYlc[1],&MZb[1+ldMZ*(bcgsl->ell)],bcgsl->ell-1);
229: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1+ldMZ], &ldMZ, &AYlc[1], &ldMZ, &bierr));
231: AYlc[0] = 0.;
232: AYlc[bcgsl->ell] = -1;
234: PetscStackCallBLAS("BLASgemv",BLASgemv_("NoTr", &ldMZ, &ldMZ, &aone, MZb, &ldMZ, AY0c, &ione, &azero, AYtc, &ione));
236: kappa0 = PetscRealPart(BLASdot_(&ldMZ, AY0c, &ione, AYtc, &ione));
238: /* round-off can cause negative kappa's */
239: if (kappa0<0) kappa0 = -kappa0;
240: kappa0 = PetscSqrtReal(kappa0);
242: kappaA = PetscRealPart(BLASdot_(&ldMZ, AYlc, &ione, AYtc, &ione));
244: PetscStackCallBLAS("BLASgemv",BLASgemv_("noTr", &ldMZ, &ldMZ, &aone, MZb, &ldMZ, AYlc, &ione, &azero, AYtc, &ione));
246: kappa1 = PetscRealPart(BLASdot_(&ldMZ, AYlc, &ione, AYtc, &ione));
248: if (kappa1<0) kappa1 = -kappa1;
249: kappa1 = PetscSqrtReal(kappa1);
251: if (kappa0!=0.0 && kappa1!=0.0) {
252: if (kappaA<0.7*kappa0*kappa1) {
253: ghat = (kappaA<0.0) ? -0.7*kappa0/kappa1 : 0.7*kappa0/kappa1;
254: } else {
255: ghat = kappaA/(kappa1*kappa1);
256: }
257: for (i=0; i<=bcgsl->ell; i++) {
258: AY0c[i] = AY0c[i] - ghat* AYlc[i];
259: }
260: }
261: }
263: omega = AY0c[bcgsl->ell];
264: for (h=bcgsl->ell; h>0 && omega==0.0; h--) omega = AY0c[h];
265: if (omega==0.0) {
266: ksp->reason = KSP_DIVERGED_BREAKDOWN;
267: return(0);
268: }
270: VecMAXPY(VX, bcgsl->ell,AY0c+1, VVR);
271: for (i=1; i<=bcgsl->ell; i++) AY0c[i] *= -1.0;
272: VecMAXPY(VVU[0], bcgsl->ell,AY0c+1, VVU+1);
273: VecMAXPY(VVR[0], bcgsl->ell,AY0c+1, VVR+1);
274: for (i=1; i<=bcgsl->ell; i++) AY0c[i] *= -1.0;
275: VecNorm(VVR[0], NORM_2, &zeta);
276: KSPCheckNorm(ksp,zeta);
278: /* Accurate Update */
279: if (bcgsl->delta>0.0) {
280: if (rnmax_computed<zeta) rnmax_computed = zeta;
281: if (rnmax_true<zeta) rnmax_true = zeta;
283: bUpdateX = (PetscBool) (zeta<bcgsl->delta*zeta0 && zeta0<=rnmax_computed);
284: if ((zeta<bcgsl->delta*rnmax_true && zeta0<=rnmax_true) || bUpdateX) {
285: /* r0 <- b-inv(K)*A*X */
286: KSP_PCApplyBAorAB(ksp, VX, VVR[0], VTM);
287: VecAYPX(VVR[0], -1.0, VB);
288: rnmax_true = zeta;
290: if (bUpdateX) {
291: VecAXPY(VXR,1.0,VX);
292: VecSet(VX,0.0);
293: VecCopy(VVR[0], VB);
294: rnmax_computed = zeta;
295: }
296: }
297: }
298: }
299: if (bcgsl->delta>0.0) {
300: VecAXPY(VX,1.0,VXR);
301: }
303: ksp->its = k;
304: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta;
305: else ksp->rnorm = 0.0;
306: KSPMonitor(ksp, ksp->its, ksp->rnorm);
307: KSPLogResidualHistory(ksp, ksp->rnorm);
308: (*ksp->converged)(ksp, k, ksp->rnorm, &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: .seealso: KSPBCGSLSetEll(), KSPBCGSLSetPol(), KSP
330: @*/
331: PetscErrorCode KSPBCGSLSetXRes(KSP ksp, PetscReal delta)
332: {
333: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
338: if (ksp->setupstage) {
339: if ((delta<=0 && bcgsl->delta>0) || (delta>0 && bcgsl->delta<=0)) {
340: VecDestroyVecs(ksp->nwork,&ksp->work);
341: PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
342: PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);
343: ksp->setupstage = KSP_SETUP_NEW;
344: }
345: }
346: bcgsl->delta = delta;
347: return(0);
348: }
350: /*@
351: KSPBCGSLSetUsePseudoinverse - Use pseudoinverse (via SVD) to solve polynomial part of update
353: Logically Collective on ksp
355: Input Parameters:
356: + ksp - iterative context obtained from KSPCreate
357: - use_pinv - set to PETSC_TRUE when using pseudoinverse
359: Options Database Keys:
361: . -ksp_bcgsl_pinv - use pseudoinverse
363: Level: intermediate
365: .seealso: KSPBCGSLSetEll(), KSP
366: @*/
367: PetscErrorCode KSPBCGSLSetUsePseudoinverse(KSP ksp,PetscBool use_pinv)
368: {
369: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
372: bcgsl->pinv = use_pinv;
373: return(0);
374: }
376: /*@
377: KSPBCGSLSetPol - Sets the type of polynomial part will
378: be used in the BiCGSTab(L) solver.
380: Logically Collective on ksp
382: Input Parameters:
383: + ksp - iterative context obtained from KSPCreate
384: - uMROR - set to PETSC_TRUE when the polynomial is a convex combination of an MR and an OR step.
386: Options Database Keys:
388: + -ksp_bcgsl_cxpoly - use enhanced polynomial
389: - -ksp_bcgsl_mrpoly - use standard polynomial
391: Level: intermediate
393: .seealso: KSP, KSPBCGSL, KSPCreate(), KSPSetType()
394: @*/
395: PetscErrorCode KSPBCGSLSetPol(KSP ksp, PetscBool uMROR)
396: {
397: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
403: if (!ksp->setupstage) bcgsl->bConvex = uMROR;
404: else if (bcgsl->bConvex != uMROR) {
405: /* free the data structures,
406: then create them again
407: */
408: VecDestroyVecs(ksp->nwork,&ksp->work);
409: PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
410: PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);
412: bcgsl->bConvex = uMROR;
413: ksp->setupstage = KSP_SETUP_NEW;
414: }
415: return(0);
416: }
418: /*@
419: KSPBCGSLSetEll - Sets the number of search directions in BiCGStab(L).
421: Logically Collective on ksp
423: Input Parameters:
424: + ksp - iterative context obtained from KSPCreate
425: - ell - number of search directions
427: Options Database Keys:
429: . -ksp_bcgsl_ell ell
431: Level: intermediate
433: Notes:
434: For large ell it is common for the polynomial update problem to become singular (due to happy breakdown for smallish
435: test problems, but also for larger problems). Consequently, by default, the system is solved by pseudoinverse, which
436: allows the iteration to complete successfully. See KSPBCGSLSetUsePseudoinverse() to switch to a conventional solve.
438: .seealso: KSPBCGSLSetUsePseudoinverse(), KSP, KSPBCGSL
439: @*/
440: PetscErrorCode KSPBCGSLSetEll(KSP ksp, PetscInt ell)
441: {
442: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
446: if (ell < 1) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "KSPBCGSLSetEll: second argument must be positive");
449: if (!ksp->setupstage) bcgsl->ell = ell;
450: else if (bcgsl->ell != ell) {
451: /* free the data structures, then create them again */
452: VecDestroyVecs(ksp->nwork,&ksp->work);
453: PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
454: PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);
456: bcgsl->ell = ell;
457: ksp->setupstage = KSP_SETUP_NEW;
458: }
459: return(0);
460: }
462: PetscErrorCode KSPView_BCGSL(KSP ksp, PetscViewer viewer)
463: {
464: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
466: PetscBool isascii;
469: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
471: if (isascii) {
472: PetscViewerASCIIPrintf(viewer, " Ell = %D\n", bcgsl->ell);
473: PetscViewerASCIIPrintf(viewer, " Delta = %lg\n", bcgsl->delta);
474: }
475: return(0);
476: }
478: PetscErrorCode KSPSetFromOptions_BCGSL(PetscOptionItems *PetscOptionsObject,KSP ksp)
479: {
480: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
482: PetscInt this_ell;
483: PetscReal delta;
484: PetscBool flga = PETSC_FALSE, flg;
487: /* PetscOptionsBegin/End are called in KSPSetFromOptions. They
488: don't need to be called here.
489: */
490: PetscOptionsHead(PetscOptionsObject,"KSP BiCGStab(L) Options");
492: /* Set number of search directions */
493: PetscOptionsInt("-ksp_bcgsl_ell","Number of Krylov search directions","KSPBCGSLSetEll",bcgsl->ell,&this_ell,&flg);
494: if (flg) {
495: KSPBCGSLSetEll(ksp, this_ell);
496: }
498: /* Set polynomial type */
499: PetscOptionsBool("-ksp_bcgsl_cxpoly", "Polynomial part of BiCGStabL is MinRes + OR", "KSPBCGSLSetPol", flga,&flga,NULL);
500: if (flga) {
501: KSPBCGSLSetPol(ksp, PETSC_TRUE);
502: } else {
503: flg = PETSC_FALSE;
504: PetscOptionsBool("-ksp_bcgsl_mrpoly", "Polynomial part of BiCGStabL is MinRes", "KSPBCGSLSetPol", flg,&flg,NULL);
505: KSPBCGSLSetPol(ksp, PETSC_FALSE);
506: }
508: /* Will computed residual be refreshed? */
509: PetscOptionsReal("-ksp_bcgsl_xres", "Threshold used to decide when to refresh computed residuals", "KSPBCGSLSetXRes", bcgsl->delta, &delta, &flg);
510: if (flg) {
511: KSPBCGSLSetXRes(ksp, delta);
512: }
514: /* Use pseudoinverse? */
515: flg = bcgsl->pinv;
516: PetscOptionsBool("-ksp_bcgsl_pinv", "Polynomial correction via pseudoinverse", "KSPBCGSLSetUsePseudoinverse",flg,&flg,NULL);
517: KSPBCGSLSetUsePseudoinverse(ksp,flg);
518: PetscOptionsTail();
519: return(0);
520: }
522: PetscErrorCode KSPSetUp_BCGSL(KSP ksp)
523: {
524: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
525: PetscInt ell = bcgsl->ell,ldMZ = ell+1;
529: KSPSetWorkVecs(ksp, 6+2*ell);
530: PetscMalloc5(ldMZ,&AY0c,ldMZ,&AYlc,ldMZ,&AYtc,ldMZ*ldMZ,&MZa,ldMZ*ldMZ,&MZb);
531: PetscBLASIntCast(5*ell,&bcgsl->lwork);
532: PetscMalloc5(bcgsl->lwork,&bcgsl->work,ell,&bcgsl->s,ell*ell,&bcgsl->u,ell*ell,&bcgsl->v,5*ell,&bcgsl->realwork);
533: return(0);
534: }
536: PetscErrorCode KSPReset_BCGSL(KSP ksp)
537: {
538: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
542: VecDestroyVecs(ksp->nwork,&ksp->work);
543: PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
544: PetscFree5(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v,bcgsl->realwork);
545: return(0);
546: }
548: PetscErrorCode KSPDestroy_BCGSL(KSP ksp)
549: {
553: KSPReset_BCGSL(ksp);
554: KSPDestroyDefault(ksp);
555: return(0);
556: }
558: /*MC
559: KSPBCGSL - Implements a slight variant of the Enhanced
560: BiCGStab(L) algorithm in (3) and (2). The variation
561: concerns cases when either kappa0**2 or kappa1**2 is
562: negative due to round-off. Kappa0 has also been pulled
563: out of the denominator in the formula for ghat.
565: References:
566: + 1. - G.L.G. Sleijpen, H.A. van der Vorst, "An overview of
567: approaches for the stable computation of hybrid BiCG
568: methods", Applied Numerical Mathematics: Transactions
569: f IMACS, 19(3), 1996.
570: . 2. - G.L.G. Sleijpen, H.A. van der Vorst, D.R. Fokkema,
571: "BiCGStab(L) and other hybrid BiCG methods",
572: Numerical Algorithms, 7, 1994.
573: - 3. - D.R. Fokkema, "Enhanced implementation of BiCGStab(L)
574: for solving linear systems of equations", preprint
575: from www.citeseer.com.
577: Contributed by: Joel M. Malard, email jm.malard@pnl.gov
579: Options Database Keys:
580: + -ksp_bcgsl_ell <ell> Number of Krylov search directions, defaults to 2 -- KSPBCGSLSetEll()
581: . -ksp_bcgsl_cxpol - Use a convex function of the MinRes and OR polynomials after the BiCG step instead of default MinRes -- KSPBCGSLSetPol()
582: . -ksp_bcgsl_mrpoly - Use the default MinRes polynomial after the BiCG step -- KSPBCGSLSetPol()
583: . -ksp_bcgsl_xres <res> Threshold used to decide when to refresh computed residuals -- KSPBCGSLSetXRes()
584: - -ksp_bcgsl_pinv <true/false> - (de)activate use of pseudoinverse -- KSPBCGSLSetUsePseudoinverse()
586: Level: beginner
588: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPFGMRES, KSPBCGS, KSPSetPCSide(), KSPBCGSLSetEll(), KSPBCGSLSetXRes()
590: M*/
591: PETSC_EXTERN PetscErrorCode KSPCreate_BCGSL(KSP ksp)
592: {
594: KSP_BCGSL *bcgsl;
597: /* allocate BiCGStab(L) context */
598: PetscNewLog(ksp,&bcgsl);
599: ksp->data = (void*)bcgsl;
601: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
602: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
603: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
605: ksp->ops->setup = KSPSetUp_BCGSL;
606: ksp->ops->solve = KSPSolve_BCGSL;
607: ksp->ops->reset = KSPReset_BCGSL;
608: ksp->ops->destroy = KSPDestroy_BCGSL;
609: ksp->ops->buildsolution = KSPBuildSolutionDefault;
610: ksp->ops->buildresidual = KSPBuildResidualDefault;
611: ksp->ops->setfromoptions = KSPSetFromOptions_BCGSL;
612: ksp->ops->view = KSPView_BCGSL;
614: /* Let the user redefine the number of directions vectors */
615: bcgsl->ell = 2;
617: /*Choose between a single MR step or an averaged MR/OR */
618: bcgsl->bConvex = PETSC_FALSE;
620: bcgsl->pinv = PETSC_TRUE; /* Use the reliable method by default */
622: /* Set the threshold for when exact residuals will be used */
623: bcgsl->delta = 0.0;
624: return(0);
625: }