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: PetscFunctionBegin;
27: /* set up temporary vectors */
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((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP));
57: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
59: PetscCall(VecSet(VVU[0], 0.0));
60: alpha = 0.;
61: rho0 = omega = 1;
63: if (bcgsl->delta > 0.0) {
64: PetscCall(VecCopy(VX, VXR));
65: PetscCall(VecSet(VX, 0.0));
66: PetscCall(VecCopy(VVR[0], VB));
67: } else {
68: PetscCall(VecCopy(ksp->vec_rhs, VB));
69: }
71: /* Life goes on */
72: PetscCall(VecCopy(VVR[0], VRT));
73: zeta = zeta0;
75: PetscCall(KSPGetTolerances(ksp, NULL, NULL, NULL, &maxit));
77: for (k = 0; k < maxit; k += bcgsl->ell) {
78: ksp->its = k;
79: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta;
80: else ksp->rnorm = 0.0;
82: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
83: PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
85: PetscCall((*ksp->converged)(ksp, k, ksp->rnorm, &ksp->reason, ksp->cnvP));
86: if (ksp->reason < 0) PetscFunctionReturn(PETSC_SUCCESS);
87: if (ksp->reason) {
88: if (bcgsl->delta > 0.0) PetscCall(VecAXPY(VX, 1.0, VXR));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
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: PetscCall(VecDot(VVR[j], VRT, &rho1));
98: KSPCheckDot(ksp, rho1);
99: if (rho1 == 0.0) {
100: ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
103: beta = alpha * (rho1 / rho0);
104: rho0 = rho1;
105: for (i = 0; i <= j; i++) {
106: /* u_i <- r_i - beta*u_i */
107: PetscCall(VecAYPX(VVU[i], -beta, VVR[i]));
108: }
109: /* u_{j+1} <- inv(K)*A*u_j */
110: PetscCall(KSP_PCApplyBAorAB(ksp, VVU[j], VVU[j + 1], VTM));
112: PetscCall(VecDot(VVU[j + 1], VRT, &sigma));
113: KSPCheckDot(ksp, sigma);
114: if (sigma == 0.0) {
115: ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
118: alpha = rho1 / sigma;
120: /* x <- x + alpha*u_0 */
121: PetscCall(VecAXPY(VX, alpha, VVU[0]));
123: for (i = 0; i <= j; i++) {
124: /* r_i <- r_i - alpha*u_{i+1} */
125: PetscCall(VecAXPY(VVR[i], -alpha, VVU[i + 1]));
126: }
128: /* r_{j+1} <- inv(K)*A*r_j */
129: PetscCall(KSP_PCApplyBAorAB(ksp, VVR[j], VVR[j + 1], VTM));
131: PetscCall(VecNorm(VVR[0], NORM_2, &nrm0));
132: KSPCheckNorm(ksp, nrm0);
133: if (bcgsl->delta > 0.0) {
134: if (rnmax_computed < nrm0) rnmax_computed = nrm0;
135: if (rnmax_true < nrm0) rnmax_true = nrm0;
136: }
138: /* NEW: check for early exit */
139: PetscCall((*ksp->converged)(ksp, k + j, nrm0, &ksp->reason, ksp->cnvP));
140: if (ksp->reason) {
141: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
142: ksp->its = k + j;
143: ksp->rnorm = nrm0;
145: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
146: if (ksp->reason < 0) PetscFunctionReturn(PETSC_SUCCESS);
147: }
148: }
150: /* Polynomial part */
151: for (i = 0; i <= bcgsl->ell; ++i) PetscCall(VecMDot(VVR[i], i + 1, VVR, &MZa[i * ldMZ]));
152: /* Symmetrize MZa */
153: for (i = 0; i <= bcgsl->ell; ++i) {
154: for (j = i + 1; j <= bcgsl->ell; ++j) MZa[i * ldMZ + j] = MZa[j * ldMZ + i] = PetscConj(MZa[j * ldMZ + i]);
155: }
156: /* Copy MZa to MZb */
157: PetscCall(PetscArraycpy(MZb, MZa, ldMZ * ldMZ));
159: if (!bcgsl->bConvex || bcgsl->ell == 1) {
160: PetscBLASInt ione = 1, bell;
161: PetscCall(PetscBLASIntCast(bcgsl->ell, &bell));
163: AY0c[0] = -1;
164: if (bcgsl->pinv) {
165: #if defined(PETSC_USE_COMPLEX)
166: 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));
167: #else
168: PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("A", "A", &bell, &bell, &MZa[1 + ldMZ], &ldMZ, bcgsl->s, bcgsl->u, &bell, bcgsl->v, &bell, bcgsl->work, &bcgsl->lwork, &bierr));
169: #endif
170: if (bierr != 0) {
171: ksp->reason = KSP_DIVERGED_BREAKDOWN;
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
174: /* Apply pseudo-inverse */
175: max_s = bcgsl->s[0];
176: for (i = 1; i < bell; i++) {
177: if (bcgsl->s[i] > max_s) max_s = bcgsl->s[i];
178: }
179: /* tolerance is hardwired to bell*max(s)*PETSC_MACHINE_EPSILON */
180: pinv_tol = bell * max_s * PETSC_MACHINE_EPSILON;
181: PetscCall(PetscArrayzero(&AY0c[1], bell));
182: for (i = 0; i < bell; i++) {
183: if (bcgsl->s[i] >= pinv_tol) {
184: utb = 0.;
185: for (j = 0; j < bell; j++) utb += MZb[1 + j] * bcgsl->u[i * bell + j];
187: for (j = 0; j < bell; j++) AY0c[1 + j] += utb / bcgsl->s[i] * bcgsl->v[j * bell + i];
188: }
189: }
190: } else {
191: PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("Lower", &bell, &MZa[1 + ldMZ], &ldMZ, &bierr));
192: if (bierr != 0) {
193: ksp->reason = KSP_DIVERGED_BREAKDOWN;
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
196: PetscCall(PetscArraycpy(&AY0c[1], &MZb[1], bcgsl->ell));
197: PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("Lower", &bell, &ione, &MZa[1 + ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
198: }
199: } else {
200: PetscBLASInt ione = 1;
201: PetscScalar aone = 1.0, azero = 0.0;
202: PetscBLASInt neqs;
203: PetscCall(PetscBLASIntCast(bcgsl->ell - 1, &neqs));
205: PetscCallBLAS("LAPACKpotrf", LAPACKpotrf_("Lower", &neqs, &MZa[1 + ldMZ], &ldMZ, &bierr));
206: if (bierr != 0) {
207: ksp->reason = KSP_DIVERGED_BREAKDOWN;
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
210: PetscCall(PetscArraycpy(&AY0c[1], &MZb[1], bcgsl->ell - 1));
211: PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1 + ldMZ], &ldMZ, &AY0c[1], &ldMZ, &bierr));
212: AY0c[0] = -1;
213: AY0c[bcgsl->ell] = 0.;
215: PetscCall(PetscArraycpy(&AYlc[1], &MZb[1 + ldMZ * (bcgsl->ell)], bcgsl->ell - 1));
216: PetscCallBLAS("LAPACKpotrs", LAPACKpotrs_("Lower", &neqs, &ione, &MZa[1 + ldMZ], &ldMZ, &AYlc[1], &ldMZ, &bierr));
218: AYlc[0] = 0.;
219: AYlc[bcgsl->ell] = -1;
221: PetscCallBLAS("BLASgemv", BLASgemv_("NoTr", &ldMZ, &ldMZ, &aone, MZb, &ldMZ, AY0c, &ione, &azero, AYtc, &ione));
223: kappa0 = PetscRealPart(BLASdot_(&ldMZ, AY0c, &ione, AYtc, &ione));
225: /* round-off can cause negative kappa's */
226: if (kappa0 < 0) kappa0 = -kappa0;
227: kappa0 = PetscSqrtReal(kappa0);
229: kappaA = PetscRealPart(BLASdot_(&ldMZ, AYlc, &ione, AYtc, &ione));
231: PetscCallBLAS("BLASgemv", BLASgemv_("noTr", &ldMZ, &ldMZ, &aone, MZb, &ldMZ, AYlc, &ione, &azero, AYtc, &ione));
233: kappa1 = PetscRealPart(BLASdot_(&ldMZ, AYlc, &ione, AYtc, &ione));
235: if (kappa1 < 0) kappa1 = -kappa1;
236: kappa1 = PetscSqrtReal(kappa1);
238: if (kappa0 != 0.0 && kappa1 != 0.0) {
239: if (kappaA < 0.7 * kappa0 * kappa1) {
240: ghat = (kappaA < 0.0) ? -0.7 * kappa0 / kappa1 : 0.7 * kappa0 / kappa1;
241: } else {
242: ghat = kappaA / (kappa1 * kappa1);
243: }
244: for (i = 0; i <= bcgsl->ell; i++) AY0c[i] = AY0c[i] - ghat * AYlc[i];
245: }
246: }
248: omega = AY0c[bcgsl->ell];
249: for (h = bcgsl->ell; h > 0 && omega == 0.0; h--) omega = AY0c[h];
250: if (omega == 0.0) {
251: ksp->reason = KSP_DIVERGED_BREAKDOWN;
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: PetscCall(VecMAXPY(VX, bcgsl->ell, AY0c + 1, VVR));
256: for (i = 1; i <= bcgsl->ell; i++) AY0c[i] *= -1.0;
257: PetscCall(VecMAXPY(VVU[0], bcgsl->ell, AY0c + 1, VVU + 1));
258: PetscCall(VecMAXPY(VVR[0], bcgsl->ell, AY0c + 1, VVR + 1));
259: for (i = 1; i <= bcgsl->ell; i++) AY0c[i] *= -1.0;
260: PetscCall(VecNorm(VVR[0], NORM_2, &zeta));
261: KSPCheckNorm(ksp, zeta);
263: /* Accurate Update */
264: if (bcgsl->delta > 0.0) {
265: if (rnmax_computed < zeta) rnmax_computed = zeta;
266: if (rnmax_true < zeta) rnmax_true = zeta;
268: bUpdateX = (PetscBool)(zeta < bcgsl->delta * zeta0 && zeta0 <= rnmax_computed);
269: if ((zeta < bcgsl->delta * rnmax_true && zeta0 <= rnmax_true) || bUpdateX) {
270: /* r0 <- b-inv(K)*A*X */
271: PetscCall(KSP_PCApplyBAorAB(ksp, VX, VVR[0], VTM));
272: PetscCall(VecAYPX(VVR[0], -1.0, VB));
273: rnmax_true = zeta;
275: if (bUpdateX) {
276: PetscCall(VecAXPY(VXR, 1.0, VX));
277: PetscCall(VecSet(VX, 0.0));
278: PetscCall(VecCopy(VVR[0], VB));
279: rnmax_computed = zeta;
280: }
281: }
282: }
283: }
284: if (bcgsl->delta > 0.0) PetscCall(VecAXPY(VX, 1.0, VXR));
286: ksp->its = k;
287: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = zeta;
288: else ksp->rnorm = 0.0;
289: PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
290: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
291: PetscCall((*ksp->converged)(ksp, k, ksp->rnorm, &ksp->reason, ksp->cnvP));
292: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: /*@
297: KSPBCGSLSetXRes - Sets the parameter governing when
298: exact residuals will be used instead of computed residuals for `KSPCBGSL`.
300: Logically Collective
302: Input Parameters:
303: + ksp - iterative context of type `KSPBCGSL`
304: - delta - computed residuals are used alone when delta is not positive
306: Options Database Key:
307: . -ksp_bcgsl_xres delta - Threshold used to decide when to refresh computed residuals
309: Level: intermediate
311: .seealso: [](ch_ksp), `KSPBCGSLSetEll()`, `KSPBCGSLSetPol()`, `KSP`, `KSPCBGSL`, `KSPBCGSLSetUsePseudoinverse()`
312: @*/
313: PetscErrorCode KSPBCGSLSetXRes(KSP ksp, PetscReal delta)
314: {
315: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
317: PetscFunctionBegin;
319: if (ksp->setupstage) {
320: if ((delta <= 0 && bcgsl->delta > 0) || (delta > 0 && bcgsl->delta <= 0)) {
321: PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));
322: PetscCall(PetscFree5(AY0c, AYlc, AYtc, MZa, MZb));
323: PetscCall(PetscFree4(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v));
324: ksp->setupstage = KSP_SETUP_NEW;
325: }
326: }
327: bcgsl->delta = delta;
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: /*@
332: KSPBCGSLSetUsePseudoinverse - Use pseudoinverse (via SVD) to solve polynomial part of the update in `KSPCBGSL` solver
334: Logically Collective
336: Input Parameters:
337: + ksp - iterative context of type `KSPCBGSL`
338: - use_pinv - set to `PETSC_TRUE` when using pseudoinverse
340: Options Database Key:
341: . -ksp_bcgsl_pinv <true,false> - use pseudoinverse
343: Level: intermediate
345: .seealso: [](ch_ksp), `KSPBCGSLSetEll()`, `KSP`, `KSPCBGSL`, `KSPBCGSLSetPol()`, `KSPBCGSLSetXRes()`
346: @*/
347: PetscErrorCode KSPBCGSLSetUsePseudoinverse(KSP ksp, PetscBool use_pinv)
348: {
349: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
351: PetscFunctionBegin;
352: bcgsl->pinv = use_pinv;
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /*@
357: KSPBCGSLSetPol - Sets the type of polynomial part that will
358: be used in the `KSPCBGSL` `KSPSolve()`
360: Logically Collective
362: Input Parameters:
363: + ksp - iterative context of type `KSPCBGSL`
364: - uMROR - set to `PETSC_TRUE` when the polynomial is a convex combination of an MR and an OR step.
366: Options Database Keys:
367: + -ksp_bcgsl_cxpoly - use enhanced polynomial
368: - -ksp_bcgsl_mrpoly - use standard polynomial
370: Level: intermediate
372: .seealso: [](ch_ksp), `KSP`, `KSPBCGSL`, `KSPCreate()`, `KSPSetType()`, `KSPCBGSL`, `KSPBCGSLSetUsePseudoinverse()`, `KSPBCGSLSetEll()`, `KSPBCGSLSetXRes()`
373: @*/
374: PetscErrorCode KSPBCGSLSetPol(KSP ksp, PetscBool uMROR)
375: {
376: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
378: PetscFunctionBegin;
381: if (!ksp->setupstage) bcgsl->bConvex = uMROR;
382: else if (bcgsl->bConvex != uMROR) {
383: /* free the data structures,
384: then create them again
385: */
386: PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));
387: PetscCall(PetscFree5(AY0c, AYlc, AYtc, MZa, MZb));
388: PetscCall(PetscFree4(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v));
390: bcgsl->bConvex = uMROR;
391: ksp->setupstage = KSP_SETUP_NEW;
392: }
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /*@
397: KSPBCGSLSetEll - Sets the number of search directions to use in the `KSPCBGSL` solver
399: Logically Collective
401: Input Parameters:
402: + ksp - iterative context of type `KSPCBGSL`
403: - ell - number of search directions
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 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;
422: PetscCheck(ell >= 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "KSPBCGSLSetEll: second argument must be positive");
425: if (!ksp->setupstage) bcgsl->ell = ell;
426: else if (bcgsl->ell != ell) {
427: /* free the data structures, then create them again */
428: PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));
429: PetscCall(PetscFree5(AY0c, AYlc, AYtc, MZa, MZb));
430: PetscCall(PetscFree4(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v));
432: bcgsl->ell = ell;
433: ksp->setupstage = KSP_SETUP_NEW;
434: }
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: static PetscErrorCode KSPView_BCGSL(KSP ksp, PetscViewer viewer)
439: {
440: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
441: PetscBool isascii;
443: PetscFunctionBegin;
444: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
446: if (isascii) {
447: PetscCall(PetscViewerASCIIPrintf(viewer, " Ell = %" PetscInt_FMT "\n", bcgsl->ell));
448: PetscCall(PetscViewerASCIIPrintf(viewer, " Delta = %g\n", (double)bcgsl->delta));
449: }
450: PetscFunctionReturn(PETSC_SUCCESS);
451: }
453: static PetscErrorCode KSPSetFromOptions_BCGSL(KSP ksp, PetscOptionItems *PetscOptionsObject)
454: {
455: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
456: PetscInt this_ell;
457: PetscReal delta;
458: PetscBool flga = PETSC_FALSE, flg;
460: PetscFunctionBegin;
461: PetscOptionsHeadBegin(PetscOptionsObject, "KSPBCGSL Options");
463: /* Set number of search directions */
464: PetscCall(PetscOptionsInt("-ksp_bcgsl_ell", "Number of Krylov search directions", "KSPBCGSLSetEll", bcgsl->ell, &this_ell, &flg));
465: if (flg) PetscCall(KSPBCGSLSetEll(ksp, this_ell));
467: /* Set polynomial type */
468: PetscCall(PetscOptionsBool("-ksp_bcgsl_cxpoly", "Polynomial part of BiCGStabL is MinRes + OR", "KSPBCGSLSetPol", flga, &flga, NULL));
469: if (flga) {
470: PetscCall(KSPBCGSLSetPol(ksp, PETSC_TRUE));
471: } else {
472: flg = PETSC_FALSE;
473: PetscCall(PetscOptionsBool("-ksp_bcgsl_mrpoly", "Polynomial part of BiCGStabL is MinRes", "KSPBCGSLSetPol", flg, &flg, NULL));
474: PetscCall(KSPBCGSLSetPol(ksp, PETSC_FALSE));
475: }
477: /* Will computed residual be refreshed? */
478: PetscCall(PetscOptionsReal("-ksp_bcgsl_xres", "Threshold used to decide when to refresh computed residuals", "KSPBCGSLSetXRes", bcgsl->delta, &delta, &flg));
479: if (flg) PetscCall(KSPBCGSLSetXRes(ksp, delta));
481: /* Use pseudoinverse? */
482: flg = bcgsl->pinv;
483: PetscCall(PetscOptionsBool("-ksp_bcgsl_pinv", "Polynomial correction via pseudoinverse", "KSPBCGSLSetUsePseudoinverse", flg, &flg, NULL));
484: PetscCall(KSPBCGSLSetUsePseudoinverse(ksp, flg));
485: PetscOptionsHeadEnd();
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: static PetscErrorCode KSPSetUp_BCGSL(KSP ksp)
490: {
491: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
492: PetscInt ell = bcgsl->ell, ldMZ = ell + 1;
494: PetscFunctionBegin;
495: PetscCall(KSPSetWorkVecs(ksp, 6 + 2 * ell));
496: PetscCall(PetscMalloc5(ldMZ, &AY0c, ldMZ, &AYlc, ldMZ, &AYtc, ldMZ * ldMZ, &MZa, ldMZ * ldMZ, &MZb));
497: PetscCall(PetscBLASIntCast(5 * ell, &bcgsl->lwork));
498: PetscCall(PetscMalloc5(bcgsl->lwork, &bcgsl->work, ell, &bcgsl->s, ell * ell, &bcgsl->u, ell * ell, &bcgsl->v, 5 * ell, &bcgsl->realwork));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: static PetscErrorCode KSPReset_BCGSL(KSP ksp)
503: {
504: KSP_BCGSL *bcgsl = (KSP_BCGSL *)ksp->data;
506: PetscFunctionBegin;
507: PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));
508: PetscCall(PetscFree5(AY0c, AYlc, AYtc, MZa, MZb));
509: PetscCall(PetscFree5(bcgsl->work, bcgsl->s, bcgsl->u, bcgsl->v, bcgsl->realwork));
510: PetscFunctionReturn(PETSC_SUCCESS);
511: }
513: static PetscErrorCode KSPDestroy_BCGSL(KSP ksp)
514: {
515: PetscFunctionBegin;
516: PetscCall(KSPReset_BCGSL(ksp));
517: PetscCall(KSPDestroyDefault(ksp));
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
521: /*MC
522: KSPBCGSL - Implements a slight variant of the Enhanced
523: BiCGStab(L) algorithm in (3) and (2). The variation
524: concerns cases when either kappa0**2 or kappa1**2 is
525: negative due to round-off. Kappa0 has also been pulled
526: out of the denominator in the formula for ghat.
528: Options Database Keys:
529: + -ksp_bcgsl_ell <ell> - Number of Krylov search directions, defaults to 2, cf. `KSPBCGSLSetEll()`
530: . -ksp_bcgsl_cxpol - Use a convex function of the MinRes and OR polynomials after the BiCG step instead of default MinRes, cf. `KSPBCGSLSetPol()`
531: . -ksp_bcgsl_mrpoly - Use the default MinRes polynomial after the BiCG step, cf. `KSPBCGSLSetPol()`
532: . -ksp_bcgsl_xres <res> - Threshold used to decide when to refresh computed residuals, cf. `KSPBCGSLSetXRes()`
533: - -ksp_bcgsl_pinv <true/false> - (de)activate use of pseudoinverse, cf. `KSPBCGSLSetUsePseudoinverse()`
535: Level: intermediate
537: Contributed by:
538: Joel M. Malard, email jm.malard@pnl.gov
540: Developer Notes:
541: This has not been completely cleaned up into PETSc style.
543: 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.
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: .seealso: [](ch_ksp), `KSPFBCGS`, `KSPFBCGSR`, `KSPBCGS`, `KSPPIPEBCGS`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPFGMRES`, `KSPBCGS`, `KSPSetPCSide()`,
558: `KSPBCGSLSetEll()`, `KSPBCGSLSetXRes()`, `KSPBCGSLSetUsePseudoinverse()`, `KSPBCGSLSetPol()`
559: M*/
560: PETSC_EXTERN PetscErrorCode KSPCreate_BCGSL(KSP ksp)
561: {
562: KSP_BCGSL *bcgsl;
564: PetscFunctionBegin;
565: /* allocate BiCGStab(L) context */
566: PetscCall(PetscNew(&bcgsl));
567: ksp->data = (void *)bcgsl;
569: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
570: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
571: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
573: ksp->ops->setup = KSPSetUp_BCGSL;
574: ksp->ops->solve = KSPSolve_BCGSL;
575: ksp->ops->reset = KSPReset_BCGSL;
576: ksp->ops->destroy = KSPDestroy_BCGSL;
577: ksp->ops->buildsolution = KSPBuildSolutionDefault;
578: ksp->ops->buildresidual = KSPBuildResidualDefault;
579: ksp->ops->setfromoptions = KSPSetFromOptions_BCGSL;
580: ksp->ops->view = KSPView_BCGSL;
582: /* Let the user redefine the number of directions vectors */
583: bcgsl->ell = 2;
585: /*Choose between a single MR step or an averaged MR/OR */
586: bcgsl->bConvex = PETSC_FALSE;
588: bcgsl->pinv = PETSC_TRUE; /* Use the reliable method by default */
590: /* Set the threshold for when exact residuals will be used */
591: bcgsl->delta = 0.0;
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }