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: #include <petsc/private/kspimpl.h>
8: #include <../src/ksp/ksp/impls/bcgsl/bcgslimpl.h>
9: #include <petscblaslapack.h>
11: static PetscErrorCode KSPSolve_BCGSL(KSP ksp)
12: {
13: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
14: PetscScalar alpha, beta, omega, sigma;
15: PetscScalar rho0, rho1;
16: PetscReal kappa0, kappaA, kappa1;
17: PetscReal ghat;
18: PetscReal zeta, zeta0, rnmax_computed, rnmax_true, nrm0;
19: PetscBool bUpdateX;
20: PetscInt maxit;
21: PetscInt h, i, j, k, vi, ell;
22: PetscBLASInt ldMZ, bierr;
23: PetscScalar utb;
24: PetscReal max_s, pinv_tol;
26: /* set up temporary vectors */
27: vi = 0;
28: ell = bcgsl->ell;
29: bcgsl->vB = ksp->work[vi];
30: vi++;
31: bcgsl->vRt = ksp->work[vi];
32: vi++;
33: bcgsl->vTm = ksp->work[vi];
34: vi++;
35: bcgsl->vvR = ksp->work + vi;
36: vi += ell + 1;
37: bcgsl->vvU = ksp->work + vi;
38: vi += ell + 1;
39: bcgsl->vXr = ksp->work[vi];
40: 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) VecAXPY(VX, 1.0, VXR);
88: return 0;
89: }
91: /* BiCG part */
92: rho0 = -omega * rho0;
93: nrm0 = zeta;
94: for (j = 0; j < bcgsl->ell; j++) {
95: /* rho1 <- r_j' * r_tilde */
96: VecDot(VVR[j], VRT, &rho1);
97: KSPCheckDot(ksp, 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: KSPCheckDot(ksp, sigma);
113: if (sigma == 0.0) {
114: ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
115: return 0;
116: }
117: alpha = rho1 / sigma;
119: /* x <- x + alpha*u_0 */
120: VecAXPY(VX, alpha, VVU[0]);
122: for (i = 0; i <= j; i++) {
123: /* r_i <- r_i - alpha*u_{i+1} */
124: VecAXPY(VVR[i], -alpha, VVU[i + 1]);
125: }
127: /* r_{j+1} <- inv(K)*A*r_j */
128: KSP_PCApplyBAorAB(ksp, VVR[j], VVR[j + 1], VTM);
130: VecNorm(VVR[0], NORM_2, &nrm0);
131: KSPCheckNorm(ksp, nrm0);
132: if (bcgsl->delta > 0.0) {
133: if (rnmax_computed < nrm0) rnmax_computed = nrm0;
134: if (rnmax_true < nrm0) rnmax_true = nrm0;
135: }
137: /* NEW: check for early exit */
138: (*ksp->converged)(ksp, k + j, nrm0, &ksp->reason, ksp->cnvP);
139: if (ksp->reason) {
140: PetscObjectSAWsTakeAccess((PetscObject)ksp);
141: ksp->its = k + j;
142: ksp->rnorm = nrm0;
144: PetscObjectSAWsGrantAccess((PetscObject)ksp);
145: if (ksp->reason < 0) return 0;
146: }
147: }
149: /* Polynomial part */
150: for (i = 0; i <= bcgsl->ell; ++i) VecMDot(VVR[i], i + 1, VVR, &MZa[i * ldMZ]);
151: /* Symmetrize MZa */
152: for (i = 0; i <= bcgsl->ell; ++i) {
153: for (j = i + 1; j <= bcgsl->ell; ++j) MZa[i * ldMZ + j] = MZa[j * ldMZ + i] = PetscConj(MZa[j * ldMZ + i]);
154: }
155: /* Copy MZa to MZb */
156: PetscArraycpy(MZb, MZa, ldMZ * ldMZ);
158: if (!bcgsl->bConvex || bcgsl->ell == 1) {
159: PetscBLASInt ione = 1, bell;
160: PetscBLASIntCast(bcgsl->ell, &bell);
162: AY0c[0] = -1;
163: if (bcgsl->pinv) {
164: #if defined(PETSC_USE_COMPLEX)
165: PetscCallBLAS("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));
166: #else
167: PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("A", "A", &bell, &bell, &MZa[1 + ldMZ], &ldMZ, bcgsl->s, bcgsl->u, &bell, bcgsl->v, &bell, bcgsl->work, &bcgsl->lwork, &bierr));
168: #endif
169: if (bierr != 0) {
170: ksp->reason = KSP_DIVERGED_BREAKDOWN;
171: return 0;
172: }
173: /* Apply pseudo-inverse */
174: max_s = bcgsl->s[0];
175: for (i = 1; i < bell; i++) {
176: if (bcgsl->s[i] > max_s) max_s = bcgsl->s[i];
177: }
178: /* tolerance is hardwired to bell*max(s)*PETSC_MACHINE_EPSILON */
179: pinv_tol = bell * max_s * PETSC_MACHINE_EPSILON;
180: PetscArrayzero(&AY0c[1], bell);
181: for (i = 0; i < bell; i++) {
182: if (bcgsl->s[i] >= pinv_tol) {
183: utb = 0.;
184: for (j = 0; j < bell; j++) utb += MZb[1 + j] * bcgsl->u[i * bell + j];
186: for (j = 0; j < bell; j++) AY0c[1 + j] += utb / bcgsl->s[i] * bcgsl->v[j * bell + i];
187: }
188: }
189: } else {
190: PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("Lower", &bell, &MZa[1 + ldMZ], &ldMZ, &bierr));
191: if (bierr != 0) {
192: ksp->reason = KSP_DIVERGED_BREAKDOWN;
193: return 0;
194: }
195: PetscArraycpy(&AY0c[1], &MZb[1], bcgsl->ell);
196: PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("Lower", &bell, &ione, &MZa[1 + ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
197: }
198: } else {
199: PetscBLASInt ione = 1;
200: PetscScalar aone = 1.0, azero = 0.0;
201: PetscBLASInt neqs;
202: PetscBLASIntCast(bcgsl->ell - 1, &neqs);
204: PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("Lower", &neqs, &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 - 1);
210: PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1 + ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
211: AY0c[0] = -1;
212: AY0c[bcgsl->ell] = 0.;
214: PetscArraycpy(&AYlc[1], &MZb[1 + ldMZ * (bcgsl->ell)], bcgsl->ell - 1);
215: PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1 + ldMZ], &ldMZ, &AYlc[1], &ldMZ, &bierr));
217: AYlc[0] = 0.;
218: AYlc[bcgsl->ell] = -1;
220: PetscCallBLAS("BLASgemv", BLASgemv_("NoTr", &ldMZ, &ldMZ, &aone, MZb, &ldMZ, AY0c, &ione, &azero, AYtc, &ione));
222: kappa0 = PetscRealPart(BLASdot_(&ldMZ, AY0c, &ione, AYtc, &ione));
224: /* round-off can cause negative kappa's */
225: if (kappa0 < 0) kappa0 = -kappa0;
226: kappa0 = PetscSqrtReal(kappa0);
228: kappaA = PetscRealPart(BLASdot_(&ldMZ, AYlc, &ione, AYtc, &ione));
230: PetscCallBLAS("BLASgemv", BLASgemv_("noTr", &ldMZ, &ldMZ, &aone, MZb, &ldMZ, AYlc, &ione, &azero, AYtc, &ione));
232: kappa1 = PetscRealPart(BLASdot_(&ldMZ, AYlc, &ione, AYtc, &ione));
234: if (kappa1 < 0) kappa1 = -kappa1;
235: kappa1 = PetscSqrtReal(kappa1);
237: if (kappa0 != 0.0 && kappa1 != 0.0) {
238: if (kappaA < 0.7 * kappa0 * kappa1) {
239: ghat = (kappaA < 0.0) ? -0.7 * kappa0 / kappa1 : 0.7 * kappa0 / kappa1;
240: } else {
241: ghat = kappaA / (kappa1 * kappa1);
242: }
243: for (i = 0; i <= bcgsl->ell; i++) AY0c[i] = AY0c[i] - ghat * AYlc[i];
244: }
245: }
247: omega = AY0c[bcgsl->ell];
248: for (h = bcgsl->ell; h > 0 && omega == 0.0; h--) omega = AY0c[h];
249: if (omega == 0.0) {
250: ksp->reason = KSP_DIVERGED_BREAKDOWN;
251: return 0;
252: }
254: VecMAXPY(VX, bcgsl->ell, AY0c + 1, VVR);
255: for (i = 1; i <= bcgsl->ell; i++) AY0c[i] *= -1.0;
256: VecMAXPY(VVU[0], bcgsl->ell, AY0c + 1, VVU + 1);
257: VecMAXPY(VVR[0], bcgsl->ell, AY0c + 1, VVR + 1);
258: for (i = 1; i <= bcgsl->ell; i++) AY0c[i] *= -1.0;
259: VecNorm(VVR[0], NORM_2, &zeta);
260: KSPCheckNorm(ksp, zeta);
262: /* Accurate Update */
263: if (bcgsl->delta > 0.0) {
264: if (rnmax_computed < zeta) rnmax_computed = zeta;
265: if (rnmax_true < zeta) rnmax_true = zeta;
267: bUpdateX = (PetscBool)(zeta < bcgsl->delta * zeta0 && zeta0 <= rnmax_computed);
268: if ((zeta < bcgsl->delta * rnmax_true && zeta0 <= rnmax_true) || bUpdateX) {
269: /* r0 <- b-inv(K)*A*X */
270: KSP_PCApplyBAorAB(ksp, VX, VVR[0], VTM);
271: VecAYPX(VVR[0], -1.0, VB);
272: rnmax_true = zeta;
274: if (bUpdateX) {
275: VecAXPY(VXR, 1.0, VX);
276: VecSet(VX, 0.0);
277: VecCopy(VVR[0], VB);
278: rnmax_computed = zeta;
279: }
280: }
281: }
282: }
283: if (bcgsl->delta > 0.0) VecAXPY(VX, 1.0, VXR);
285: ksp->its = k;
286: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta;
287: else ksp->rnorm = 0.0;
288: KSPMonitor(ksp, ksp->its, ksp->rnorm);
289: KSPLogResidualHistory(ksp, ksp->rnorm);
290: (*ksp->converged)(ksp, k, ksp->rnorm, &ksp->reason, ksp->cnvP);
291: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
292: return 0;
293: }
295: /*@
296: KSPBCGSLSetXRes - Sets the parameter governing when
297: exact residuals will be used instead of computed residuals.
299: Logically Collective
301: Input Parameters:
302: + ksp - iterative context of type `KSPBCGSL`
303: - delta - computed residuals are used alone when delta is not positive
305: Options Database Key:
306: . -ksp_bcgsl_xres delta - Threshold used to decide when to refresh computed residuals
308: Level: intermediate
310: .seealso: [](chapter_ksp), `KSPBCGSLSetEll()`, `KSPBCGSLSetPol()`, `KSP`, `KSPCBGSL`, `KSPBCGSLSetUsePseudoinverse()`
311: @*/
312: PetscErrorCode KSPBCGSLSetXRes(KSP ksp, PetscReal delta)
313: {
314: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
317: if (ksp->setupstage) {
318: if ((delta <= 0 && bcgsl->delta > 0) || (delta > 0 && bcgsl->delta <= 0)) {
319: VecDestroyVecs(ksp->nwork, &ksp->work);
320: PetscFree5(AY0c, AYlc, AYtc, MZa, MZb);
321: PetscFree4(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v);
322: ksp->setupstage = KSP_SETUP_NEW;
323: }
324: }
325: bcgsl->delta = delta;
326: return 0;
327: }
329: /*@
330: KSPBCGSLSetUsePseudoinverse - Use pseudoinverse (via SVD) to solve polynomial part of update in `KSPCBGSL` solver
332: Logically Collective
334: Input Parameters:
335: + ksp - iterative context of type `KSPCBGSL`
336: - use_pinv - set to PETSC_TRUE when using pseudoinverse
338: Options Database Key:
339: . -ksp_bcgsl_pinv - use pseudoinverse
341: Level: intermediate
343: .seealso: [](chapter_ksp), `KSPBCGSLSetEll()`, `KSP`, `KSPCBGSL`, `KSPBCGSLSetPol()`, `KSPBCGSLSetXRes()`
344: @*/
345: PetscErrorCode KSPBCGSLSetUsePseudoinverse(KSP ksp, PetscBool use_pinv)
346: {
347: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
349: bcgsl->pinv = use_pinv;
350: return 0;
351: }
353: /*@
354: KSPBCGSLSetPol - Sets the type of polynomial part that will
355: be used in the `KSPCBGSL` solver.
357: Logically Collective
359: Input Parameters:
360: + ksp - iterative context of type `KSPCBGSL`
361: - uMROR - set to PETSC_TRUE when the polynomial is a convex combination of an MR and an OR step.
363: Options Database Keys:
364: + -ksp_bcgsl_cxpoly - use enhanced polynomial
365: - -ksp_bcgsl_mrpoly - use standard polynomial
367: Level: intermediate
369: .seealso: [](chapter_ksp), `KSP`, `KSPBCGSL`, `KSPCreate()`, `KSPSetType()`, `KSPCBGSL`, `KSPBCGSLSetUsePseudoinverse()`, `KSPBCGSLSetEll()`, `KSPBCGSLSetXRes()`
370: @*/
371: PetscErrorCode KSPBCGSLSetPol(KSP ksp, PetscBool uMROR)
372: {
373: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
377: if (!ksp->setupstage) bcgsl->bConvex = uMROR;
378: else if (bcgsl->bConvex != uMROR) {
379: /* free the data structures,
380: then create them again
381: */
382: VecDestroyVecs(ksp->nwork, &ksp->work);
383: PetscFree5(AY0c, AYlc, AYtc, MZa, MZb);
384: PetscFree4(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v);
386: bcgsl->bConvex = uMROR;
387: ksp->setupstage = KSP_SETUP_NEW;
388: }
389: return 0;
390: }
392: /*@
393: KSPBCGSLSetEll - Sets the number of search directions in `KSPCBGSL` solver
395: Logically Collective
397: Input Parameters:
398: + ksp - iterative context of type `KSPCBGSL`
399: - ell - number of search directions
401: Options Database Keys:
402: . -ksp_bcgsl_ell ell - Number of Krylov search directions
404: Level: intermediate
406: Notes:
407: For large ell it is common for the polynomial update problem to become singular (due to happy breakdown for smallish
408: test problems, but also for larger problems). Consequently, by default, the system is solved by pseudoinverse, which
409: allows the iteration to complete successfully. See `KSPBCGSLSetUsePseudoinverse()` to switch to a conventional solve.
411: .seealso: [](chapter_ksp), `KSPBCGSLSetUsePseudoinverse()`, `KSP`, `KSPBCGSL`, `KSPBCGSLSetPol()`, `KSPBCGSLSetXRes()`
412: @*/
413: PetscErrorCode KSPBCGSLSetEll(KSP ksp, PetscInt ell)
414: {
415: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
420: if (!ksp->setupstage) bcgsl->ell = ell;
421: else if (bcgsl->ell != ell) {
422: /* free the data structures, then create them again */
423: VecDestroyVecs(ksp->nwork, &ksp->work);
424: PetscFree5(AY0c, AYlc, AYtc, MZa, MZb);
425: PetscFree4(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v);
427: bcgsl->ell = ell;
428: ksp->setupstage = KSP_SETUP_NEW;
429: }
430: return 0;
431: }
433: PetscErrorCode KSPView_BCGSL(KSP ksp, PetscViewer viewer)
434: {
435: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
436: PetscBool isascii;
438: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
440: if (isascii) {
441: PetscViewerASCIIPrintf(viewer, " Ell = %" PetscInt_FMT "\n", bcgsl->ell);
442: PetscViewerASCIIPrintf(viewer, " Delta = %g\n", (double)bcgsl->delta);
443: }
444: return 0;
445: }
447: PetscErrorCode KSPSetFromOptions_BCGSL(KSP ksp, PetscOptionItems *PetscOptionsObject)
448: {
449: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
450: PetscInt this_ell;
451: PetscReal delta;
452: PetscBool flga = PETSC_FALSE, flg;
454: PetscOptionsHeadBegin(PetscOptionsObject, "KSPBCGSL Options");
456: /* Set number of search directions */
457: PetscOptionsInt("-ksp_bcgsl_ell", "Number of Krylov search directions", "KSPBCGSLSetEll", bcgsl->ell, &this_ell, &flg);
458: if (flg) KSPBCGSLSetEll(ksp, this_ell);
460: /* Set polynomial type */
461: PetscOptionsBool("-ksp_bcgsl_cxpoly", "Polynomial part of BiCGStabL is MinRes + OR", "KSPBCGSLSetPol", flga, &flga, NULL);
462: if (flga) {
463: KSPBCGSLSetPol(ksp, PETSC_TRUE);
464: } else {
465: flg = PETSC_FALSE;
466: PetscOptionsBool("-ksp_bcgsl_mrpoly", "Polynomial part of BiCGStabL is MinRes", "KSPBCGSLSetPol", flg, &flg, NULL);
467: KSPBCGSLSetPol(ksp, PETSC_FALSE);
468: }
470: /* Will computed residual be refreshed? */
471: PetscOptionsReal("-ksp_bcgsl_xres", "Threshold used to decide when to refresh computed residuals", "KSPBCGSLSetXRes", bcgsl->delta, &delta, &flg);
472: if (flg) KSPBCGSLSetXRes(ksp, delta);
474: /* Use pseudoinverse? */
475: flg = bcgsl->pinv;
476: PetscOptionsBool("-ksp_bcgsl_pinv", "Polynomial correction via pseudoinverse", "KSPBCGSLSetUsePseudoinverse", flg, &flg, NULL);
477: KSPBCGSLSetUsePseudoinverse(ksp, flg);
478: PetscOptionsHeadEnd();
479: return 0;
480: }
482: PetscErrorCode KSPSetUp_BCGSL(KSP ksp)
483: {
484: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
485: PetscInt ell = bcgsl->ell, ldMZ = ell + 1;
487: KSPSetWorkVecs(ksp, 6 + 2 * ell);
488: PetscMalloc5(ldMZ, &AY0c, ldMZ, &AYlc, ldMZ, &AYtc, ldMZ * ldMZ, &MZa, ldMZ * ldMZ, &MZb);
489: PetscBLASIntCast(5 * ell, &bcgsl->lwork);
490: PetscMalloc5(bcgsl->lwork, &bcgsl->work, ell, &bcgsl->s, ell * ell, &bcgsl->u, ell * ell, &bcgsl->v, 5 * ell, &bcgsl->realwork);
491: return 0;
492: }
494: PetscErrorCode KSPReset_BCGSL(KSP ksp)
495: {
496: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
498: VecDestroyVecs(ksp->nwork, &ksp->work);
499: PetscFree5(AY0c, AYlc, AYtc, MZa, MZb);
500: PetscFree5(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v, bcgsl->realwork);
501: return 0;
502: }
504: PetscErrorCode KSPDestroy_BCGSL(KSP ksp)
505: {
506: KSPReset_BCGSL(ksp);
507: KSPDestroyDefault(ksp);
508: return 0;
509: }
511: /*MC
512: KSPBCGSL - Implements a slight variant of the Enhanced
513: BiCGStab(L) algorithm in (3) and (2). The variation
514: concerns cases when either kappa0**2 or kappa1**2 is
515: negative due to round-off. Kappa0 has also been pulled
516: out of the denominator in the formula for ghat.
518: Options Database Keys:
519: + -ksp_bcgsl_ell <ell> - Number of Krylov search directions, defaults to 2, cf. `KSPBCGSLSetEll()`
520: . -ksp_bcgsl_cxpol - Use a convex function of the MinRes and OR polynomials after the BiCG step instead of default MinRes, cf. `KSPBCGSLSetPol()`
521: . -ksp_bcgsl_mrpoly - Use the default MinRes polynomial after the BiCG step, cf. `KSPBCGSLSetPol()`
522: . -ksp_bcgsl_xres <res> - Threshold used to decide when to refresh computed residuals, cf. `KSPBCGSLSetXRes()`
523: - -ksp_bcgsl_pinv <true/false> - (de)activate use of pseudoinverse, cf. `KSPBCGSLSetUsePseudoinverse()`
525: Level: intermediate
527: References:
528: + * - G.L.G. Sleijpen, H.A. van der Vorst, "An overview of
529: approaches for the stable computation of hybrid BiCG
530: methods", Applied Numerical Mathematics: Transactions
531: f IMACS, 19(3), 1996.
532: . * - G.L.G. Sleijpen, H.A. van der Vorst, D.R. Fokkema,
533: "BiCGStab(L) and other hybrid BiCG methods",
534: Numerical Algorithms, 7, 1994.
535: - * - D.R. Fokkema, "Enhanced implementation of BiCGStab(L)
536: for solving linear systems of equations", preprint
537: from www.citeseer.com.
539: Contributed by:
540: Joel M. Malard, email jm.malard@pnl.gov
542: Developer Notes:
543: This has not been completely cleaned up into PETSc style.
545: All the BLAS and LAPACK calls in the source should be removed and replaced with loops and the macros for block solvers converted from LINPACK.
547: .seealso: [](chapter_ksp), `KSPFBCGS`, `KSPFBCGSR`, `KSPBCGS`, `KSPPIPEBCGS`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPBCGS`, `KSPSetPCSide()`,
548: `KSPBCGSLSetEll()`, `KSPBCGSLSetXRes()`, `KSPBCGSLSetUsePseudoinverse()`, `KSPBCGSLSetPol()`
549: M*/
550: PETSC_EXTERN PetscErrorCode KSPCreate_BCGSL(KSP ksp)
551: {
552: KSP_BCGSL *bcgsl;
554: /* allocate BiCGStab(L) context */
555: PetscNew(&bcgsl);
556: ksp->data = (void *)bcgsl;
558: KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3);
559: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2);
560: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1);
562: ksp->ops->setup = KSPSetUp_BCGSL;
563: ksp->ops->solve = KSPSolve_BCGSL;
564: ksp->ops->reset = KSPReset_BCGSL;
565: ksp->ops->destroy = KSPDestroy_BCGSL;
566: ksp->ops->buildsolution = KSPBuildSolutionDefault;
567: ksp->ops->buildresidual = KSPBuildResidualDefault;
568: ksp->ops->setfromoptions = KSPSetFromOptions_BCGSL;
569: ksp->ops->view = KSPView_BCGSL;
571: /* Let the user redefine the number of directions vectors */
572: bcgsl->ell = 2;
574: /*Choose between a single MR step or an averaged MR/OR */
575: bcgsl->bConvex = PETSC_FALSE;
577: bcgsl->pinv = PETSC_TRUE; /* Use the reliable method by default */
579: /* Set the threshold for when exact residuals will be used */
580: bcgsl->delta = 0.0;
581: return 0;
582: }