Actual source code: bcgsl.c
petsc-3.4.5 2014-06-29
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> /*I "petscksp.h" I*/
14: #include <../src/ksp/ksp/impls/bcgsl/bcgslimpl.h>
15: #include <petscblaslapack.h>
20: static PetscErrorCode KSPSolve_BCGSL(KSP ksp)
21: {
22: KSP_BCGSL *bcgsl = (KSP_BCGSL*) ksp->data;
23: PetscScalar alpha, beta, omega, sigma;
24: PetscScalar rho0, rho1;
25: PetscReal kappa0, kappaA, kappa1;
26: PetscReal ghat;
27: PetscReal zeta, zeta0, rnmax_computed, rnmax_true, nrm0;
28: PetscBool bUpdateX;
29: PetscInt maxit;
30: PetscInt h, i, j, k, vi, ell;
31: PetscBLASInt ldMZ,bierr;
32: PetscScalar utb;
33: PetscReal max_s, pinv_tol;
37: /* set up temporary vectors */
38: vi = 0;
39: ell = bcgsl->ell;
40: bcgsl->vB = ksp->work[vi]; vi++;
41: bcgsl->vRt = ksp->work[vi]; vi++;
42: bcgsl->vTm = ksp->work[vi]; vi++;
43: bcgsl->vvR = ksp->work+vi; vi += ell+1;
44: bcgsl->vvU = ksp->work+vi; vi += ell+1;
45: bcgsl->vXr = ksp->work[vi]; vi++;
46: PetscBLASIntCast(ell+1,&ldMZ);
48: /* Prime the iterative solver */
49: KSPInitialResidual(ksp, VX, VTM, VB, VVR[0], ksp->vec_rhs);
50: VecNorm(VVR[0], NORM_2, &zeta0);
51: rnmax_computed = zeta0;
52: rnmax_true = zeta0;
54: (*ksp->converged)(ksp, 0, zeta0, &ksp->reason, ksp->cnvP);
55: if (ksp->reason) {
56: PetscObjectAMSTakeAccess((PetscObject)ksp);
57: ksp->its = 0;
58: ksp->rnorm = zeta0;
59: PetscObjectAMSGrantAccess((PetscObject)ksp);
60: return(0);
61: }
63: VecSet(VVU[0],0.0);
64: alpha = 0.;
65: rho0 = omega = 1;
67: if (bcgsl->delta>0.0) {
68: VecCopy(VX, VXR);
69: VecSet(VX,0.0);
70: VecCopy(VVR[0], VB);
71: } else {
72: VecCopy(ksp->vec_rhs, VB);
73: }
75: /* Life goes on */
76: VecCopy(VVR[0], VRT);
77: zeta = zeta0;
79: KSPGetTolerances(ksp, NULL, NULL, NULL, &maxit);
81: for (k=0; k<maxit; k += bcgsl->ell) {
82: ksp->its = k;
83: ksp->rnorm = zeta;
85: KSPLogResidualHistory(ksp, zeta);
86: KSPMonitor(ksp, ksp->its, zeta);
88: (*ksp->converged)(ksp, k, zeta, &ksp->reason, ksp->cnvP);
89: if (ksp->reason < 0) return(0);
90: else if (ksp->reason) break;
92: /* BiCG part */
93: rho0 = -omega*rho0;
94: nrm0 = zeta;
95: for (j=0; j<bcgsl->ell; j++) {
96: /* rho1 <- r_j' * r_tilde */
97: VecDot(VVR[j], VRT, &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: if (sigma == 0.0) {
113: ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
114: return(0);
115: }
116: alpha = rho1/sigma;
118: /* x <- x + alpha*u_0 */
119: VecAXPY(VX, alpha, VVU[0]);
121: for (i=0; i<=j; i++) {
122: /* r_i <- r_i - alpha*u_{i+1} */
123: VecAXPY(VVR[i], -alpha, VVU[i+1]);
124: }
126: /* r_{j+1} <- inv(K)*A*r_j */
127: KSP_PCApplyBAorAB(ksp, VVR[j], VVR[j+1], VTM);
129: VecNorm(VVR[0], NORM_2, &nrm0);
130: if (bcgsl->delta>0.0) {
131: if (rnmax_computed<nrm0) rnmax_computed = nrm0;
132: if (rnmax_true<nrm0) rnmax_true = nrm0;
133: }
135: /* NEW: check for early exit */
136: (*ksp->converged)(ksp, k+j, nrm0, &ksp->reason, ksp->cnvP);
137: if (ksp->reason) {
138: PetscObjectAMSTakeAccess((PetscObject)ksp);
140: ksp->its = k+j;
141: ksp->rnorm = nrm0;
143: PetscObjectAMSGrantAccess((PetscObject)ksp);
144: if (ksp->reason < 0) return(0);
145: }
146: }
148: /* Polynomial part */
149: for (i = 0; i <= bcgsl->ell; ++i) {
150: VecMDot(VVR[i], i+1, VVR, &MZa[i*ldMZ]);
151: }
152: /* Symmetrize MZa */
153: for (i = 0; i <= bcgsl->ell; ++i) {
154: for (j = i+1; j <= bcgsl->ell; ++j) {
155: MZa[i*ldMZ+j] = MZa[j*ldMZ+i] = PetscConj(MZa[j*ldMZ+i]);
156: }
157: }
158: /* Copy MZa to MZb */
159: PetscMemcpy(MZb,MZa,ldMZ*ldMZ*sizeof(PetscScalar));
161: if (!bcgsl->bConvex || bcgsl->ell==1) {
162: PetscBLASInt ione = 1,bell;
163: PetscBLASIntCast(bcgsl->ell,&bell);
165: AY0c[0] = -1;
166: if (bcgsl->pinv) {
167: #if defined(PETSC_MISSING_LAPACK_GESVD)
168: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable.");
169: #else
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: #endif
176: if (bierr!=0) {
177: ksp->reason = KSP_DIVERGED_BREAKDOWN;
178: return(0);
179: }
180: /* Apply pseudo-inverse */
181: max_s = bcgsl->s[0];
182: for (i=1; i<bell; i++) {
183: if (bcgsl->s[i] > max_s) {
184: max_s = bcgsl->s[i];
185: }
186: }
187: /* tolerance is hardwired to bell*max(s)*PETSC_MACHINE_EPSILON */
188: pinv_tol = bell*max_s*PETSC_MACHINE_EPSILON;
189: PetscMemzero(&AY0c[1],bell*sizeof(PetscScalar));
190: for (i=0; i<bell; i++) {
191: if (bcgsl->s[i] >= pinv_tol) {
192: utb=0.;
193: for (j=0; j<bell; j++) {
194: utb += MZb[1+j]*bcgsl->u[i*bell+j];
195: }
197: for (j=0; j<bell; j++) {
198: AY0c[1+j] += utb/bcgsl->s[i]*bcgsl->v[j*bell+i];
199: }
200: }
201: }
202: } else {
203: #if defined(PETSC_MISSING_LAPACK_POTRF)
204: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
205: #else
206: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("Lower", &bell, &MZa[1+ldMZ], &ldMZ, &bierr));
207: #endif
208: if (bierr!=0) {
209: ksp->reason = KSP_DIVERGED_BREAKDOWN;
210: return(0);
211: }
212: PetscMemcpy(&AY0c[1],&MZb[1],bcgsl->ell*sizeof(PetscScalar));
213: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &bell, &ione, &MZa[1+ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
214: }
215: } else {
216: PetscBLASInt ione = 1;
217: PetscScalar aone = 1.0, azero = 0.0;
218: PetscBLASInt neqs;
219: PetscBLASIntCast(bcgsl->ell-1,&neqs);
221: #if defined(PETSC_MISSING_LAPACK_POTRF)
222: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable.");
223: #else
224: PetscStackCallBLAS("LAPACKpotrf",LAPACKpotrf_("Lower", &neqs, &MZa[1+ldMZ], &ldMZ, &bierr));
225: #endif
226: if (bierr!=0) {
227: ksp->reason = KSP_DIVERGED_BREAKDOWN;
228: return(0);
229: }
230: PetscMemcpy(&AY0c[1],&MZb[1],(bcgsl->ell-1)*sizeof(PetscScalar));
231: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1+ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
232: AY0c[0] = -1;
233: AY0c[bcgsl->ell] = 0.;
235: PetscMemcpy(&AYlc[1],&MZb[1+ldMZ*(bcgsl->ell)],(bcgsl->ell-1)*sizeof(PetscScalar));
236: PetscStackCallBLAS("LAPACKpotrs",LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1+ldMZ], &ldMZ, &AYlc[1], &ldMZ, &bierr));
238: AYlc[0] = 0.;
239: AYlc[bcgsl->ell] = -1;
241: PetscStackCallBLAS("BLASgemv",BLASgemv_("NoTr", &ldMZ, &ldMZ, &aone, MZb, &ldMZ, AY0c, &ione, &azero, AYtc, &ione));
243: kappa0 = PetscRealPart(BLASdot_(&ldMZ, AY0c, &ione, AYtc, &ione));
245: /* round-off can cause negative kappa's */
246: if (kappa0<0) kappa0 = -kappa0;
247: kappa0 = PetscSqrtReal(kappa0);
249: kappaA = PetscRealPart(BLASdot_(&ldMZ, AYlc, &ione, AYtc, &ione));
251: PetscStackCallBLAS("BLASgemv",BLASgemv_("noTr", &ldMZ, &ldMZ, &aone, MZb, &ldMZ, AYlc, &ione, &azero, AYtc, &ione));
253: kappa1 = PetscRealPart(BLASdot_(&ldMZ, AYlc, &ione, AYtc, &ione));
255: if (kappa1<0) kappa1 = -kappa1;
256: kappa1 = PetscSqrtReal(kappa1);
258: if (kappa0!=0.0 && kappa1!=0.0) {
259: if (kappaA<0.7*kappa0*kappa1) {
260: ghat = (kappaA<0.0) ? -0.7*kappa0/kappa1 : 0.7*kappa0/kappa1;
261: } else {
262: ghat = kappaA/(kappa1*kappa1);
263: }
264: for (i=0; i<=bcgsl->ell; i++) {
265: AY0c[i] = AY0c[i] - ghat* AYlc[i];
266: }
267: }
268: }
270: omega = AY0c[bcgsl->ell];
271: for (h=bcgsl->ell; h>0 && omega==0.0; h--) omega = AY0c[h];
272: if (omega==0.0) {
273: ksp->reason = KSP_DIVERGED_BREAKDOWN;
274: return(0);
275: }
278: VecMAXPY(VX, bcgsl->ell,AY0c+1, VVR);
279: for (i=1; i<=bcgsl->ell; i++) AY0c[i] *= -1.0;
280: VecMAXPY(VVU[0], bcgsl->ell,AY0c+1, VVU+1);
281: VecMAXPY(VVR[0], bcgsl->ell,AY0c+1, VVR+1);
282: for (i=1; i<=bcgsl->ell; i++) AY0c[i] *= -1.0;
283: VecNorm(VVR[0], NORM_2, &zeta);
285: /* Accurate Update */
286: if (bcgsl->delta>0.0) {
287: if (rnmax_computed<zeta) rnmax_computed = zeta;
288: if (rnmax_true<zeta) rnmax_true = zeta;
290: bUpdateX = (PetscBool) (zeta<bcgsl->delta*zeta0 && zeta0<=rnmax_computed);
291: if ((zeta<bcgsl->delta*rnmax_true && zeta0<=rnmax_true) || bUpdateX) {
292: /* r0 <- b-inv(K)*A*X */
293: KSP_PCApplyBAorAB(ksp, VX, VVR[0], VTM);
294: VecAYPX(VVR[0], -1.0, VB);
295: rnmax_true = zeta;
297: if (bUpdateX) {
298: VecAXPY(VXR,1.0,VX);
299: VecSet(VX,0.0);
300: VecCopy(VVR[0], VB);
301: rnmax_computed = zeta;
302: }
303: }
304: }
305: }
306: if (bcgsl->delta>0.0) {
307: VecAXPY(VX,1.0,VXR);
308: }
310: (*ksp->converged)(ksp, k, zeta, &ksp->reason, ksp->cnvP);
311: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
312: return(0);
313: }
317: /*@
318: KSPBCGSLSetXRes - Sets the parameter governing when
319: exact residuals will be used instead of computed residuals.
321: Logically Collective on KSP
323: Input Parameters:
324: + ksp - iterative context obtained from KSPCreate
325: - delta - computed residuals are used alone when delta is not positive
327: Options Database Keys:
329: . -ksp_bcgsl_xres delta
331: Level: intermediate
333: .keywords: KSP, BiCGStab(L), set, exact residuals
335: .seealso: KSPBCGSLSetEll(), KSPBCGSLSetPol()
336: @*/
337: PetscErrorCode KSPBCGSLSetXRes(KSP ksp, PetscReal delta)
338: {
339: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
344: if (ksp->setupstage) {
345: if ((delta<=0 && bcgsl->delta>0) || (delta>0 && bcgsl->delta<=0)) {
346: VecDestroyVecs(ksp->nwork,&ksp->work);
347: PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
348: PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);
349: ksp->setupstage = KSP_SETUP_NEW;
350: }
351: }
352: bcgsl->delta = delta;
353: return(0);
354: }
358: /*@
359: KSPBCGSLSetUsePseudoinverse - Use pseudoinverse (via SVD) to solve polynomial part of update
361: Logically Collective on KSP
363: Input Parameters:
364: + ksp - iterative context obtained from KSPCreate
365: - use_pinv - set to PETSC_TRUE when using pseudoinverse
367: Options Database Keys:
369: + -ksp_bcgsl_pinv - use pseudoinverse
371: Level: intermediate
373: .keywords: KSP, BiCGStab(L), set, polynomial
375: .seealso: KSPBCGSLSetEll()
376: @*/
377: PetscErrorCode KSPBCGSLSetUsePseudoinverse(KSP ksp,PetscBool use_pinv)
378: {
379: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
382: bcgsl->pinv = use_pinv;
383: return(0);
384: }
388: /*@
389: KSPBCGSLSetPol - Sets the type of polynomial part will
390: be used in the BiCGSTab(L) solver.
392: Logically Collective on KSP
394: Input Parameters:
395: + ksp - iterative context obtained from KSPCreate
396: - uMROR - set to PETSC_TRUE when the polynomial is a convex combination of an MR and an OR step.
398: Options Database Keys:
400: + -ksp_bcgsl_cxpoly - use enhanced polynomial
401: . -ksp_bcgsl_mrpoly - use standard polynomial
403: Level: intermediate
405: .keywords: KSP, BiCGStab(L), set, polynomial
407: .seealso: @()
408: @*/
409: PetscErrorCode KSPBCGSLSetPol(KSP ksp, PetscBool uMROR)
410: {
411: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
417: if (!ksp->setupstage) bcgsl->bConvex = uMROR;
418: else if (bcgsl->bConvex != uMROR) {
419: /* free the data structures,
420: then create them again
421: */
422: VecDestroyVecs(ksp->nwork,&ksp->work);
423: PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
424: PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);
426: bcgsl->bConvex = uMROR;
427: ksp->setupstage = KSP_SETUP_NEW;
428: }
429: return(0);
430: }
434: /*@
435: KSPBCGSLSetEll - Sets the number of search directions in BiCGStab(L).
437: Logically Collective on KSP
439: Input Parameters:
440: + ksp - iterative context obtained from KSPCreate
441: - ell - number of search directions
443: Options Database Keys:
445: . -ksp_bcgsl_ell ell
447: Level: intermediate
449: Notes:
450: For large ell it is common for the polynomial update problem to become singular (due to happy breakdown for smallish
451: test problems, but also for larger problems). Consequently, by default, the system is solved by pseudoinverse, which
452: allows the iteration to complete successfully. See KSPBCGSLSetUsePseudoinverse() to switch to a conventional solve.
454: .keywords: KSP, BiCGStab(L), set, exact residuals,
456: .seealso: KSPBCGSLSetUsePseudoinverse()
457: @*/
458: PetscErrorCode KSPBCGSLSetEll(KSP ksp, PetscInt ell)
459: {
460: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
464: if (ell < 1) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "KSPBCGSLSetEll: second argument must be positive");
467: if (!ksp->setupstage) bcgsl->ell = ell;
468: else if (bcgsl->ell != ell) {
469: /* free the data structures, then create them again */
470: VecDestroyVecs(ksp->nwork,&ksp->work);
471: PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
472: PetscFree4(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v);
474: bcgsl->ell = ell;
475: ksp->setupstage = KSP_SETUP_NEW;
476: }
477: return(0);
478: }
482: PetscErrorCode KSPView_BCGSL(KSP ksp, PetscViewer viewer)
483: {
484: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
486: PetscBool isascii;
489: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
491: if (isascii) {
492: PetscViewerASCIIPrintf(viewer, " BCGSL: Ell = %D\n", bcgsl->ell);
493: PetscViewerASCIIPrintf(viewer, " BCGSL: Delta = %lg\n", bcgsl->delta);
494: }
495: return(0);
496: }
500: PetscErrorCode KSPSetFromOptions_BCGSL(KSP ksp)
501: {
502: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
504: PetscInt this_ell;
505: PetscReal delta;
506: PetscBool flga = PETSC_FALSE, flg;
509: /* PetscOptionsBegin/End are called in KSPSetFromOptions. They
510: don't need to be called here.
511: */
512: PetscOptionsHead("KSP BiCGStab(L) Options");
514: /* Set number of search directions */
515: PetscOptionsInt("-ksp_bcgsl_ell","Number of Krylov search directions","KSPBCGSLSetEll",bcgsl->ell,&this_ell,&flg);
516: if (flg) {
517: KSPBCGSLSetEll(ksp, this_ell);
518: }
520: /* Set polynomial type */
521: PetscOptionsBool("-ksp_bcgsl_cxpoly", "Polynomial part of BiCGStabL is MinRes + OR", "KSPBCGSLSetPol", flga,&flga,NULL);
522: if (flga) {
523: KSPBCGSLSetPol(ksp, PETSC_TRUE);
524: } else {
525: flg = PETSC_FALSE;
526: PetscOptionsBool("-ksp_bcgsl_mrpoly", "Polynomial part of BiCGStabL is MinRes", "KSPBCGSLSetPol", flg,&flg,NULL);
527: KSPBCGSLSetPol(ksp, PETSC_FALSE);
528: }
530: /* Will computed residual be refreshed? */
531: PetscOptionsReal("-ksp_bcgsl_xres", "Threshold used to decide when to refresh computed residuals", "KSPBCGSLSetXRes", bcgsl->delta, &delta, &flg);
532: if (flg) {
533: KSPBCGSLSetXRes(ksp, delta);
534: }
536: /* Use pseudoinverse? */
537: flg = bcgsl->pinv;
538: PetscOptionsBool("-ksp_bcgsl_pinv", "Polynomial correction via pseudoinverse", "KSPBCGSLSetUsePseudoinverse",flg,&flg,NULL);
539: KSPBCGSLSetUsePseudoinverse(ksp,flg);
540: PetscOptionsTail();
541: return(0);
542: }
546: PetscErrorCode KSPSetUp_BCGSL(KSP ksp)
547: {
548: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
549: PetscInt ell = bcgsl->ell,ldMZ = ell+1;
553: KSPSetWorkVecs(ksp, 6+2*ell);
554: PetscMalloc5(ldMZ,PetscScalar,&AY0c,ldMZ,PetscScalar,&AYlc,ldMZ,PetscScalar,&AYtc,ldMZ*ldMZ,PetscScalar,&MZa,ldMZ*ldMZ,PetscScalar,&MZb);
555: PetscBLASIntCast(5*ell,&bcgsl->lwork);
556: PetscMalloc5(bcgsl->lwork,PetscScalar,&bcgsl->work,ell,PetscReal,&bcgsl->s,ell*ell,PetscScalar,&bcgsl->u,ell*ell,PetscScalar,&bcgsl->v,5*ell,PetscReal,&bcgsl->realwork);
557: return(0);
558: }
562: PetscErrorCode KSPReset_BCGSL(KSP ksp)
563: {
564: KSP_BCGSL *bcgsl = (KSP_BCGSL*)ksp->data;
568: VecDestroyVecs(ksp->nwork,&ksp->work);
569: PetscFree5(AY0c,AYlc,AYtc,MZa,MZb);
570: PetscFree5(bcgsl->work,bcgsl->s,bcgsl->u,bcgsl->v,bcgsl->realwork);
571: return(0);
572: }
576: PetscErrorCode KSPDestroy_BCGSL(KSP ksp)
577: {
581: KSPReset_BCGSL(ksp);
582: KSPDestroyDefault(ksp);
583: return(0);
584: }
586: /*MC
587: KSPBCGSL - Implements a slight variant of the Enhanced
588: BiCGStab(L) algorithm in (3) and (2). The variation
589: concerns cases when either kappa0**2 or kappa1**2 is
590: negative due to round-off. Kappa0 has also been pulled
591: out of the denominator in the formula for ghat.
593: References:
594: 1. G.L.G. Sleijpen, H.A. van der Vorst, "An overview of
595: approaches for the stable computation of hybrid BiCG
596: methods", Applied Numerical Mathematics: Transactions
597: f IMACS, 19(3), pp 235-54, 1996.
598: 2. G.L.G. Sleijpen, H.A. van der Vorst, D.R. Fokkema,
599: "BiCGStab(L) and other hybrid Bi-CG methods",
600: Numerical Algorithms, 7, pp 75-109, 1994.
601: 3. D.R. Fokkema, "Enhanced implementation of BiCGStab(L)
602: for solving linear systems of equations", preprint
603: from www.citeseer.com.
605: Contributed by: Joel M. Malard, email jm.malard@pnl.gov
607: Options Database Keys:
608: + -ksp_bcgsl_ell <ell> Number of Krylov search directions, defaults to 2 -- KSPBCGSLSetEll()
609: . -ksp_bcgsl_cxpol - Use a convex function of the MinRes and OR polynomials after the BiCG step instead of default MinRes -- KSPBCGSLSetPol()
610: . -ksp_bcgsl_mrpoly - Use the default MinRes polynomial after the BiCG step -- KSPBCGSLSetPol()
611: . -ksp_bcgsl_xres <res> Threshold used to decide when to refresh computed residuals -- KSPBCGSLSetXRes()
612: - -ksp_bcgsl_pinv <true/false> - (de)activate use of pseudoinverse -- KSPBCGSLSetUsePseudoinverse()
614: Level: beginner
616: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPFGMRES, KSPBCGS, KSPSetPCSide(), KSPBCGSLSetEll(), KSPBCGSLSetXRes()
618: M*/
621: PETSC_EXTERN PetscErrorCode KSPCreate_BCGSL(KSP ksp)
622: {
624: KSP_BCGSL *bcgsl;
627: /* allocate BiCGStab(L) context */
628: PetscNewLog(ksp, KSP_BCGSL, &bcgsl);
629: ksp->data = (void*)bcgsl;
631: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
632: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,1);
634: ksp->ops->setup = KSPSetUp_BCGSL;
635: ksp->ops->solve = KSPSolve_BCGSL;
636: ksp->ops->reset = KSPReset_BCGSL;
637: ksp->ops->destroy = KSPDestroy_BCGSL;
638: ksp->ops->buildsolution = KSPBuildSolutionDefault;
639: ksp->ops->buildresidual = KSPBuildResidualDefault;
640: ksp->ops->setfromoptions = KSPSetFromOptions_BCGSL;
641: ksp->ops->view = KSPView_BCGSL;
643: /* Let the user redefine the number of directions vectors */
644: bcgsl->ell = 2;
646: /*Choose between a single MR step or an averaged MR/OR */
647: bcgsl->bConvex = PETSC_FALSE;
649: bcgsl->pinv = PETSC_TRUE; /* Use the reliable method by default */
651: /* Set the threshold for when exact residuals will be used */
652: bcgsl->delta = 0.0;
653: return(0);
654: }