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;
23: PetscScalar utb;
24: PetscReal max_s, pinv_tol;
26: PetscFunctionBegin;
27: PetscCheck(ksp->nwork == 6 + 2 * bcgsl->ell, PetscObjectComm((PetscObject)ksp), PETSC_ERR_COR, "Unexpected number of work vectors %" PetscInt_FMT " != %" PetscInt_FMT, ksp->nwork, 6 + 2 * bcgsl->ell);
28: vi = 0;
29: ell = bcgsl->ell;
30: bcgsl->vB = ksp->work[vi];
31: vi++;
32: bcgsl->vRt = ksp->work[vi];
33: vi++;
34: bcgsl->vTm = ksp->work[vi];
35: vi++;
36: bcgsl->vvR = ksp->work + vi;
37: vi += ell + 1;
38: bcgsl->vvU = ksp->work + vi;
39: vi += ell + 1;
40: bcgsl->vXr = ksp->work[vi];
41: vi++;
42: PetscCall(PetscBLASIntCast(ell + 1, &ldMZ));
44: /* Prime the iterative solver */
45: PetscCall(KSPInitialResidual(ksp, VX, VTM, VB, VVR[0], ksp->vec_rhs));
46: PetscCall(VecNorm(VVR[0], NORM_2, &zeta0));
47: KSPCheckNorm(ksp, zeta0);
48: rnmax_computed = zeta0;
49: rnmax_true = zeta0;
51: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
52: ksp->its = 0;
53: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta0;
54: else ksp->rnorm = 0.0;
55: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
56: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
57: PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
58: PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP));
59: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
61: PetscCall(VecSet(VVU[0], 0.0));
62: alpha = 0.;
63: rho0 = omega = 1;
65: if (bcgsl->delta > 0.0) {
66: PetscCall(VecCopy(VX, VXR));
67: PetscCall(VecSet(VX, 0.0));
68: PetscCall(VecCopy(VVR[0], VB));
69: } else {
70: PetscCall(VecCopy(ksp->vec_rhs, VB));
71: }
73: /* Life goes on */
74: PetscCall(VecCopy(VVR[0], VRT));
75: zeta = zeta0;
77: PetscCall(KSPGetTolerances(ksp, NULL, NULL, NULL, &maxit));
79: for (k = 0; k < maxit; k += bcgsl->ell) {
80: ksp->its = k;
81: if (k > 0) {
82: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta;
83: else ksp->rnorm = 0.0;
85: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
86: PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
88: PetscCall((*ksp->converged)(ksp, k, ksp->rnorm, &ksp->reason, ksp->cnvP));
89: if (ksp->reason < 0) PetscFunctionReturn(PETSC_SUCCESS);
90: if (ksp->reason) {
91: if (bcgsl->delta > 0.0) PetscCall(VecAXPY(VX, 1.0, VXR));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
94: }
96: /* BiCG part */
97: rho0 = -omega * rho0;
98: nrm0 = zeta;
99: for (j = 0; j < bcgsl->ell; j++) {
100: /* rho1 <- r_j' * r_tilde */
101: PetscCall(VecDot(VVR[j], VRT, &rho1));
102: KSPCheckDot(ksp, rho1);
103: if (rho1 == 0.0) {
104: ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
107: beta = alpha * (rho1 / rho0);
108: rho0 = rho1;
109: for (i = 0; i <= j; i++) {
110: /* u_i <- r_i - beta*u_i */
111: PetscCall(VecAYPX(VVU[i], -beta, VVR[i]));
112: }
113: /* u_{j+1} <- inv(K)*A*u_j */
114: PetscCall(KSP_PCApplyBAorAB(ksp, VVU[j], VVU[j + 1], VTM));
116: PetscCall(VecDot(VVU[j + 1], VRT, &sigma));
117: KSPCheckDot(ksp, sigma);
118: if (sigma == 0.0) {
119: ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
122: alpha = rho1 / sigma;
124: /* x <- x + alpha*u_0 */
125: PetscCall(VecAXPY(VX, alpha, VVU[0]));
127: for (i = 0; i <= j; i++) {
128: /* r_i <- r_i - alpha*u_{i+1} */
129: PetscCall(VecAXPY(VVR[i], -alpha, VVU[i + 1]));
130: }
132: /* r_{j+1} <- inv(K)*A*r_j */
133: PetscCall(KSP_PCApplyBAorAB(ksp, VVR[j], VVR[j + 1], VTM));
135: PetscCall(VecNorm(VVR[0], NORM_2, &nrm0));
136: KSPCheckNorm(ksp, nrm0);
137: if (bcgsl->delta > 0.0) {
138: if (rnmax_computed < nrm0) rnmax_computed = nrm0;
139: if (rnmax_true < nrm0) rnmax_true = nrm0;
140: }
142: /* NEW: check for early exit */
143: PetscCall((*ksp->converged)(ksp, k + j, nrm0, &ksp->reason, ksp->cnvP));
144: if (ksp->reason) {
145: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
146: ksp->its = k + j;
147: ksp->rnorm = nrm0;
149: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
150: if (ksp->reason < 0) PetscFunctionReturn(PETSC_SUCCESS);
151: }
152: }
154: /* Polynomial part */
155: for (i = 0; i <= bcgsl->ell; ++i) PetscCall(VecMDot(VVR[i], i + 1, VVR, &MZa[i * ldMZ]));
156: /* Symmetrize MZa */
157: for (i = 0; i <= bcgsl->ell; ++i) {
158: for (j = i + 1; j <= bcgsl->ell; ++j) MZa[i * ldMZ + j] = MZa[j * ldMZ + i] = PetscConj(MZa[j * ldMZ + i]);
159: }
160: /* Copy MZa to MZb */
161: PetscCall(PetscArraycpy(MZb, MZa, ldMZ * ldMZ));
163: if (!bcgsl->bConvex || bcgsl->ell == 1) {
164: PetscBLASInt ione = 1, bell, info;
165: PetscCall(PetscBLASIntCast(bcgsl->ell, &bell));
167: AY0c[0] = -1;
168: if (bcgsl->pinv) {
169: #if PetscDefined(USE_COMPLEX)
170: 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, &info));
171: #else
172: PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("A", "A", &bell, &bell, &MZa[1 + ldMZ], &ldMZ, bcgsl->s, bcgsl->u, &bell, bcgsl->v, &bell, bcgsl->work, &bcgsl->lwork, &info));
173: #endif
174: if (info != 0) {
175: ksp->reason = KSP_DIVERGED_BREAKDOWN;
176: PetscFunctionReturn(PETSC_SUCCESS);
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) max_s = bcgsl->s[i];
182: }
183: /* tolerance is hardwired to bell*max(s)*PETSC_MACHINE_EPSILON */
184: pinv_tol = bell * max_s * PETSC_MACHINE_EPSILON;
185: PetscCall(PetscArrayzero(&AY0c[1], bell));
186: for (i = 0; i < bell; i++) {
187: if (bcgsl->s[i] >= pinv_tol) {
188: utb = 0.;
189: for (j = 0; j < bell; j++) utb += MZb[1 + j] * bcgsl->u[i * bell + j];
191: for (j = 0; j < bell; j++) AY0c[1 + j] += utb / bcgsl->s[i] * bcgsl->v[j * bell + i];
192: }
193: }
194: } else {
195: PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("Lower", &bell, &MZa[1 + ldMZ], &ldMZ, &info));
196: if (info != 0) {
197: ksp->reason = KSP_DIVERGED_BREAKDOWN;
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
200: PetscCall(PetscArraycpy(&AY0c[1], &MZb[1], bcgsl->ell));
201: PetscCallLAPACKInfo("LAPACKpotrs", LAPACKpotrs_("Lower", &bell, &ione, &MZa[1 + ldMZ], &ldMZ, &AY0c[1], &ldMZ, &info));
202: }
203: } else {
204: PetscBLASInt ione = 1, info;
205: PetscScalar aone = 1.0, azero = 0.0;
206: PetscBLASInt neqs;
207: PetscCall(PetscBLASIntCast(bcgsl->ell - 1, &neqs));
209: PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("Lower", &neqs, &MZa[1 + ldMZ], &ldMZ, &info));
210: if (info != 0) {
211: ksp->reason = KSP_DIVERGED_BREAKDOWN;
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
214: PetscCall(PetscArraycpy(&AY0c[1], &MZb[1], bcgsl->ell - 1));
215: PetscCallLAPACKInfo("LAPACKpotrs", LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1 + ldMZ], &ldMZ, &AY0c[1], &ldMZ, &info));
216: AY0c[0] = -1;
217: AY0c[bcgsl->ell] = 0.;
219: PetscCall(PetscArraycpy(&AYlc[1], &MZb[1 + ldMZ * (bcgsl->ell)], bcgsl->ell - 1));
220: PetscCallLAPACKInfo("LAPACKpotrs", LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1 + ldMZ], &ldMZ, &AYlc[1], &ldMZ, &info));
222: AYlc[0] = 0.;
223: AYlc[bcgsl->ell] = -1;
225: PetscCallBLAS("BLASgemv", BLASgemv_("NoTr", &ldMZ, &ldMZ, &aone, MZb, &ldMZ, AY0c, &ione, &azero, AYtc, &ione));
227: kappa0 = PetscRealPart(BLASdot_(&ldMZ, AY0c, &ione, AYtc, &ione));
229: /* round-off can cause negative kappa's */
230: if (kappa0 < 0) kappa0 = -kappa0;
231: kappa0 = PetscSqrtReal(kappa0);
233: kappaA = PetscRealPart(BLASdot_(&ldMZ, AYlc, &ione, AYtc, &ione));
235: PetscCallBLAS("BLASgemv", BLASgemv_("noTr", &ldMZ, &ldMZ, &aone, MZb, &ldMZ, AYlc, &ione, &azero, AYtc, &ione));
237: kappa1 = PetscRealPart(BLASdot_(&ldMZ, AYlc, &ione, AYtc, &ione));
239: if (kappa1 < 0) kappa1 = -kappa1;
240: kappa1 = PetscSqrtReal(kappa1);
242: if (kappa0 != 0.0 && kappa1 != 0.0) {
243: if (kappaA < 0.7 * kappa0 * kappa1) {
244: ghat = (kappaA < 0.0) ? -0.7 * kappa0 / kappa1 : 0.7 * kappa0 / kappa1;
245: } else {
246: ghat = kappaA / (kappa1 * kappa1);
247: }
248: for (i = 0; i <= bcgsl->ell; i++) AY0c[i] = AY0c[i] - ghat * AYlc[i];
249: }
250: }
252: omega = AY0c[bcgsl->ell];
253: for (h = bcgsl->ell; h > 0 && omega == 0.0; h--) omega = AY0c[h];
254: if (omega == 0.0) {
255: ksp->reason = KSP_DIVERGED_BREAKDOWN;
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: PetscCall(VecMAXPY(VX, bcgsl->ell, AY0c + 1, VVR));
260: for (i = 1; i <= bcgsl->ell; i++) AY0c[i] *= -1.0;
261: PetscCall(VecMAXPY(VVU[0], bcgsl->ell, AY0c + 1, VVU + 1));
262: PetscCall(VecMAXPY(VVR[0], bcgsl->ell, AY0c + 1, VVR + 1));
263: for (i = 1; i <= bcgsl->ell; i++) AY0c[i] *= -1.0;
264: PetscCall(VecNorm(VVR[0], NORM_2, &zeta));
265: KSPCheckNorm(ksp, zeta);
267: /* Accurate Update */
268: if (bcgsl->delta > 0.0) {
269: if (rnmax_computed < zeta) rnmax_computed = zeta;
270: if (rnmax_true < zeta) rnmax_true = zeta;
272: bUpdateX = (PetscBool)(zeta < bcgsl->delta * zeta0 && zeta0 <= rnmax_computed);
273: if ((zeta < bcgsl->delta * rnmax_true && zeta0 <= rnmax_true) || bUpdateX) {
274: /* r0 <- b-inv(K)*A*X */
275: PetscCall(KSP_PCApplyBAorAB(ksp, VX, VVR[0], VTM));
276: PetscCall(VecAYPX(VVR[0], -1.0, VB));
277: rnmax_true = zeta;
279: if (bUpdateX) {
280: PetscCall(VecAXPY(VXR, 1.0, VX));
281: PetscCall(VecSet(VX, 0.0));
282: PetscCall(VecCopy(VVR[0], VB));
283: rnmax_computed = zeta;
284: }
285: }
286: }
287: }
288: if (bcgsl->delta > 0.0) PetscCall(VecAXPY(VX, 1.0, VXR));
290: ksp->its = k;
291: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta;
292: else ksp->rnorm = 0.0;
293: PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
294: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
295: PetscCall((*ksp->converged)(ksp, k, ksp->rnorm, &ksp->reason, ksp->cnvP));
296: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: static PetscErrorCode KSPReset_BCGSL_Private(KSP ksp)
301: {
302: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
304: PetscFunctionBegin;
305: PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));
306: PetscCall(PetscFree5(AY0c, AYlc, AYtc, MZa, MZb));
307: PetscCall(PetscFree5(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v, bcgsl->realwork));
309: ksp->setupstage = KSP_SETUP_NEW;
310: ksp->nwork = 0;
311: PetscFunctionReturn(PETSC_SUCCESS);
312: }
314: /*@
315: KSPBCGSLSetXRes - Sets the parameter governing when
316: exact residuals will be used instead of computed residuals for `KSPCBGSL`.
318: Logically Collective
320: Input Parameters:
321: + ksp - iterative context of type `KSPBCGSL`
322: - delta - computed residuals are used alone when delta is not positive
324: Options Database Key:
325: . -ksp_bcgsl_xres delta - Threshold used to decide when to refresh computed residuals
327: Level: intermediate
329: .seealso: [](ch_ksp), `KSPBCGSLSetEll()`, `KSPBCGSLSetPol()`, `KSP`, `KSPCBGSL`, `KSPBCGSLSetUsePseudoinverse()`
330: @*/
331: PetscErrorCode KSPBCGSLSetXRes(KSP ksp, PetscReal delta)
332: {
333: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
335: PetscFunctionBegin;
337: if (ksp->setupstage && ((delta <= 0 && bcgsl->delta > 0) || (delta > 0 && bcgsl->delta <= 0))) PetscCall(KSPReset_BCGSL_Private(ksp));
338: bcgsl->delta = delta;
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: /*@
343: KSPBCGSLSetUsePseudoinverse - Use pseudoinverse (via SVD) to solve polynomial part of the update in `KSPCBGSL` solver
345: Logically Collective
347: Input Parameters:
348: + ksp - iterative context of type `KSPCBGSL`
349: - use_pinv - set to `PETSC_TRUE` when using pseudoinverse
351: Options Database Key:
352: . -ksp_bcgsl_pinv (true|false) - use pseudoinverse
354: Level: intermediate
356: .seealso: [](ch_ksp), `KSPBCGSLSetEll()`, `KSP`, `KSPCBGSL`, `KSPBCGSLSetPol()`, `KSPBCGSLSetXRes()`
357: @*/
358: PetscErrorCode KSPBCGSLSetUsePseudoinverse(KSP ksp, PetscBool use_pinv)
359: {
360: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
362: PetscFunctionBegin;
363: bcgsl->pinv = use_pinv;
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: /*@
368: KSPBCGSLSetPol - Sets the type of polynomial part that will
369: be used in the `KSPCBGSL` `KSPSolve()`
371: Logically Collective
373: Input Parameters:
374: + ksp - iterative context of type `KSPCBGSL`
375: - uMROR - set to `PETSC_TRUE` when the polynomial is a convex combination of an MR and an OR step.
377: Options Database Keys:
378: + -ksp_bcgsl_cxpoly - use enhanced polynomial
379: - -ksp_bcgsl_mrpoly - use standard polynomial
381: Level: intermediate
383: .seealso: [](ch_ksp), `KSP`, `KSPBCGSL`, `KSPCreate()`, `KSPSetType()`, `KSPCBGSL`, `KSPBCGSLSetUsePseudoinverse()`, `KSPBCGSLSetEll()`, `KSPBCGSLSetXRes()`
384: @*/
385: PetscErrorCode KSPBCGSLSetPol(KSP ksp, PetscBool uMROR)
386: {
387: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
389: PetscFunctionBegin;
391: if (ksp->setupstage && bcgsl->bConvex != uMROR) PetscCall(KSPReset_BCGSL_Private(ksp));
392: bcgsl->bConvex = uMROR;
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /*@
397: KSPBCGSLSetEll - Sets the number of search directions to use in the `KSPBCGSL` Krylov solver
399: Logically Collective
401: Input Parameters:
402: + ksp - iterative context, `KSP`, of type `KSPBCGSL`
403: - ell - number of search directions to use
405: Options Database Key:
406: . -ksp_bcgsl_ell ell - Number of Krylov search directions
408: Level: intermediate
410: Notes:
411: For large `ell` it is common for the polynomial update problem to become singular (due to happy breakdown for smallish
412: test problems, but also for larger problems). Consequently, by default, the system is solved by using the pseudoinverse, which
413: allows the iteration to complete successfully. See `KSPBCGSLSetUsePseudoinverse()` to switch to a conventional solve.
415: .seealso: [](ch_ksp), `KSPBCGSLSetUsePseudoinverse()`, `KSP`, `KSPBCGSL`, `KSPBCGSLSetPol()`, `KSPBCGSLSetXRes()`
416: @*/
417: PetscErrorCode KSPBCGSLSetEll(KSP ksp, PetscInt ell)
418: {
419: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
421: PetscFunctionBegin;
423: PetscCheck(ell > 0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "KSPBCGSLSetEll(): second argument must be positive");
424: if (ksp->setupstage && bcgsl->ell != ell) PetscCall(KSPReset_BCGSL_Private(ksp));
425: bcgsl->ell = ell;
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: static PetscErrorCode KSPView_BCGSL(KSP ksp, PetscViewer viewer)
430: {
431: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
432: PetscBool isascii;
434: PetscFunctionBegin;
435: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
437: if (isascii) {
438: PetscCall(PetscViewerASCIIPrintf(viewer, " Ell = %" PetscInt_FMT "\n", bcgsl->ell));
439: PetscCall(PetscViewerASCIIPrintf(viewer, " Delta = %g\n", (double)bcgsl->delta));
440: }
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: static PetscErrorCode KSPSetFromOptions_BCGSL(KSP ksp, PetscOptionItems PetscOptionsObject)
445: {
446: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
447: PetscInt this_ell;
448: PetscReal delta;
449: PetscBool flga = PETSC_FALSE, flg;
451: PetscFunctionBegin;
452: PetscOptionsHeadBegin(PetscOptionsObject, "KSPBCGSL Options");
454: /* Set number of search directions */
455: PetscCall(PetscOptionsInt("-ksp_bcgsl_ell", "Number of Krylov search directions", "KSPBCGSLSetEll", bcgsl->ell, &this_ell, &flg));
456: if (flg) PetscCall(KSPBCGSLSetEll(ksp, this_ell));
458: /* Set polynomial type */
459: PetscCall(PetscOptionsBool("-ksp_bcgsl_cxpoly", "Polynomial part of BiCGStabL is MinRes + OR", "KSPBCGSLSetPol", flga, &flga, NULL));
460: if (flga) PetscCall(KSPBCGSLSetPol(ksp, PETSC_TRUE));
461: else {
462: flg = PETSC_FALSE;
463: PetscCall(PetscOptionsBool("-ksp_bcgsl_mrpoly", "Polynomial part of BiCGStabL is MinRes", "KSPBCGSLSetPol", flg, &flg, NULL));
464: PetscCall(KSPBCGSLSetPol(ksp, PETSC_FALSE));
465: }
467: /* Will computed residual be refreshed? */
468: PetscCall(PetscOptionsReal("-ksp_bcgsl_xres", "Threshold used to decide when to refresh computed residuals", "KSPBCGSLSetXRes", bcgsl->delta, &delta, &flg));
469: if (flg) PetscCall(KSPBCGSLSetXRes(ksp, delta));
471: /* Use pseudoinverse? */
472: flg = bcgsl->pinv;
473: PetscCall(PetscOptionsBool("-ksp_bcgsl_pinv", "Polynomial correction via pseudoinverse", "KSPBCGSLSetUsePseudoinverse", flg, &flg, NULL));
474: PetscCall(KSPBCGSLSetUsePseudoinverse(ksp, flg));
475: PetscOptionsHeadEnd();
476: PetscFunctionReturn(PETSC_SUCCESS);
477: }
479: static PetscErrorCode KSPSetUp_BCGSL(KSP ksp)
480: {
481: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
482: PetscInt ell = bcgsl->ell, ldMZ = ell + 1;
484: PetscFunctionBegin;
485: PetscCall(KSPSetWorkVecs(ksp, 6 + 2 * ell));
486: PetscCall(PetscMalloc5(ldMZ, &AY0c, ldMZ, &AYlc, ldMZ, &AYtc, ldMZ * ldMZ, &MZa, ldMZ * ldMZ, &MZb));
487: PetscCall(PetscBLASIntCast(5 * ell, &bcgsl->lwork));
488: PetscCall(PetscMalloc5(bcgsl->lwork, &bcgsl->work, ell, &bcgsl->s, ell * ell, &bcgsl->u, ell * ell, &bcgsl->v, 5 * ell, &bcgsl->realwork));
489: PetscFunctionReturn(PETSC_SUCCESS);
490: }
492: static PetscErrorCode KSPReset_BCGSL(KSP ksp)
493: {
494: PetscFunctionBegin;
495: PetscCall(KSPReset_BCGSL_Private(ksp));
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: /*MC
500: KSPBCGSL - Implements a slight variant of the Enhanced BiCGStab(L) algorithm in {cite}`fokkema1996enhanced`
501: and {cite}`sleijpen1994bicgstab`, see also {cite}`sleijpen1995overview`. The variation
502: concerns cases when either kappa0**2 or kappa1**2 is
503: negative due to round-off. Kappa0 has also been pulled
504: out of the denominator in the formula for ghat.
506: Options Database Keys:
507: + -ksp_bcgsl_ell ell - Number of Krylov search directions to use, defaults to 2, cf. `KSPBCGSLSetEll()`
508: . -ksp_bcgsl_cxpol (true|false) - Use a convex function of the MinRes and OR polynomials after the BiCG step instead of default MinRes, cf. `KSPBCGSLSetPol()`
509: . -ksp_bcgsl_mrpoly (true|false) - Use the default MinRes polynomial after the BiCG step, cf. `KSPBCGSLSetPol()`
510: . -ksp_bcgsl_xres res - Threshold used to decide when to refresh computed residuals, cf. `KSPBCGSLSetXRes()`
511: - -ksp_bcgsl_pinv (true|false) - (de)activate use of pseudoinverse, cf. `KSPBCGSLSetUsePseudoinverse()`
513: Level: intermediate
515: Note:
516: The "sub-iterations" of this solver are not reported by `-ksp_monitor` or recorded in `KSPSetResidualHistory()` since the solution is not directly computed for
517: these sub-iterations.
519: Contributed by:
520: Joel M. Malard, email jm.malard@pnl.gov
522: Developer Notes:
523: This has not been completely cleaned up into PETSc style.
525: 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.
527: .seealso: [](ch_ksp), `KSPFBCGS`, `KSPFBCGSR`, `KSPBCGS`, `KSPPIPEBCGS`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPSetPCSide()`,
528: `KSPBCGSLSetEll()`, `KSPBCGSLSetXRes()`, `KSPBCGSLSetUsePseudoinverse()`, `KSPBCGSLSetPol()`
529: M*/
530: PETSC_EXTERN PetscErrorCode KSPCreate_BCGSL(KSP ksp)
531: {
532: KSP_BCGSL *bcgsl;
534: PetscFunctionBegin;
535: /* allocate BiCGStab(L) context */
536: PetscCall(PetscNew(&bcgsl));
537: ksp->data = (void *)bcgsl;
539: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
540: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
541: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
543: ksp->ops->setup = KSPSetUp_BCGSL;
544: ksp->ops->solve = KSPSolve_BCGSL;
545: ksp->ops->reset = KSPReset_BCGSL;
546: ksp->ops->destroy = KSPDestroyDefault;
547: ksp->ops->buildsolution = KSPBuildSolutionDefault;
548: ksp->ops->buildresidual = KSPBuildResidualDefault;
549: ksp->ops->setfromoptions = KSPSetFromOptions_BCGSL;
550: ksp->ops->view = KSPView_BCGSL;
552: /* Let the user redefine the number of directions vectors */
553: bcgsl->ell = 2;
555: /*Choose between a single MR step or an averaged MR/OR */
556: bcgsl->bConvex = PETSC_FALSE;
558: bcgsl->pinv = PETSC_TRUE; /* Use the reliable method by default */
560: /* Set the threshold for when exact residuals will be used */
561: bcgsl->delta = 0.0;
562: PetscFunctionReturn(PETSC_SUCCESS);
563: }