Actual source code: minres.c
1: #include <petsc/private/kspimpl.h>
3: PetscBool QLPcite = PETSC_FALSE;
4: const char QLPCitation[] = "@article{choi2011minres,\n"
5: " title={MINRES-QLP: A Krylov subspace method for indefinite or singular symmetric systems},\n"
6: " author={Choi, Sou-Cheng T and Paige, Christopher C and Saunders, Michael A},\n"
7: " journal={SIAM Journal on Scientific Computing},\n"
8: " volume={33},\n"
9: " number={4},\n"
10: " pages={1810--1836},\n"
11: " year={2011},\n}\n";
13: typedef struct {
14: PetscReal haptol;
15: PetscBool qlp;
16: PetscReal maxxnorm;
17: PetscReal TranCond;
18: PetscBool monitor;
19: PetscViewer viewer;
20: PetscViewerFormat viewer_fmt;
21: } KSP_MINRES;
23: static PetscErrorCode KSPSetUp_MINRES(KSP ksp)
24: {
25: PetscFunctionBegin;
26: PetscCall(KSPSetWorkVecs(ksp, 9));
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: /* Convenience functions */
31: #define KSPMinresSwap3(V1, V2, V3) \
32: do { \
33: Vec T = V1; \
34: V1 = V2; \
35: V2 = V3; \
36: V3 = T; \
37: } while (0)
39: static inline PetscReal Norm3(PetscReal a, PetscReal b, PetscReal c)
40: {
41: return PetscSqrtReal(PetscSqr(a) + PetscSqr(b) + PetscSqr(c));
42: }
44: static inline void SymOrtho(PetscReal a, PetscReal b, PetscReal *c, PetscReal *s, PetscReal *r)
45: {
46: if (b == 0.0) {
47: if (a == 0.0) *c = 1.0;
48: else *c = PetscCopysignReal(1.0, a);
49: *s = 0.0;
50: *r = PetscAbsReal(a);
51: } else if (a == 0.0) {
52: *c = 0.0;
53: *s = PetscCopysignReal(1.0, b);
54: *r = PetscAbsReal(b);
55: } else if (PetscAbsReal(b) > PetscAbsReal(a)) {
56: PetscReal t = a / b;
58: *s = PetscCopysignReal(1.0, b) / PetscSqrtReal(1.0 + t * t);
59: *c = (*s) * t;
60: *r = b / (*s); // computationally better than d = a / c since |c| <= |s|
61: } else {
62: PetscReal t = b / a;
64: *c = PetscCopysignReal(1.0, a) / PetscSqrtReal(1.0 + t * t);
65: *s = (*c) * t;
66: *r = a / (*c); // computationally better than d = b / s since |s| <= |c|
67: }
68: }
70: /*
71: Code adapted from https://stanford.edu/group/SOL/software/minresqlp/minresqlp-matlab/CPS11.zip
72: CSP11/Algorithms/MINRESQLP/minresQLP.m
73: */
74: static PetscErrorCode KSPSolve_MINRES(KSP ksp)
75: {
76: Mat Amat;
77: Vec X, B, R1, R2, R3, V, W, WL, WL2, XL2, RN;
78: PetscReal alpha, beta, beta1, betan, betal;
79: PetscBool diagonalscale;
80: PetscReal zero = 0.0, dbar, dltan = 0.0, dlta, cs = -1.0, sn = 0.0, epln, eplnn = 0.0, gbar, dlta_QLP;
81: PetscReal gamal3 = 0.0, gamal2 = 0.0, gamal = 0.0, gama = 0.0, gama_tmp;
82: PetscReal taul2 = 0.0, taul = 0.0, tau = 0.0, phi;
83: PetscReal Axnorm, xnorm, xnorm_tmp, xl2norm = 0.0, pnorm, Anorm = 0.0, gmin = 0.0, gminl = 0.0, gminl2 = 0.0;
84: PetscReal Acond = 1.0, Acondl = 0.0, rnorml, rnorm, rootl, relAresl, relres, relresl, Arnorml, Anorml = 0.0, xnorml = 0.0;
85: PetscReal epsx, realmin = PETSC_REAL_MIN, eps = PETSC_MACHINE_EPSILON;
86: PetscReal veplnl2 = 0.0, veplnl = 0.0, vepln = 0.0, etal2 = 0.0, etal = 0.0, eta = 0.0;
87: PetscReal dlta_tmp, sr2 = 0.0, cr2 = -1.0, cr1 = -1.0, sr1 = 0.0;
88: PetscReal ul4 = 0.0, ul3 = 0.0, ul2 = 0.0, ul = 0.0, u = 0.0, ul_QLP = 0.0, u_QLP = 0.0;
89: PetscReal vepln_QLP = 0.0, gamal_QLP = 0.0, gama_QLP = 0.0, gamal_tmp, abs_gama;
90: PetscInt flag = -2, flag0 = -2, QLPiter = 0;
91: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
93: PetscFunctionBegin;
94: PetscCall(PetscCitationsRegister(QLPCitation, &QLPcite));
95: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
96: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
98: X = ksp->vec_sol;
99: B = ksp->vec_rhs;
100: R1 = ksp->work[0];
101: R2 = ksp->work[1];
102: R3 = ksp->work[2];
103: V = ksp->work[3];
104: W = ksp->work[4];
105: WL = ksp->work[5];
106: WL2 = ksp->work[6];
107: XL2 = ksp->work[7];
108: RN = ksp->work[8];
109: PetscCall(PCGetOperators(ksp->pc, &Amat, NULL));
111: ksp->its = 0;
112: ksp->rnorm = 0.0;
113: if (!ksp->guess_zero) {
114: PetscCall(KSP_MatMult(ksp, Amat, X, R2));
115: PetscCall(VecNorm(R2, NORM_2, &Axnorm));
116: PetscCall(VecNorm(X, NORM_2, &xnorm));
117: PetscCall(VecAYPX(R2, -1.0, B));
118: } else {
119: PetscCall(VecCopy(B, R2));
120: Axnorm = 0.0;
121: xnorm = 0.0;
122: }
123: if (ksp->converged_neg_curve) PetscCall(VecCopy(R2, RN));
124: PetscCall(KSP_PCApply(ksp, R2, R3));
125: PetscCall(VecDotRealPart(R3, R2, &beta1));
126: KSPCheckDot(ksp, beta1);
127: if (beta1 < 0.0) {
128: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g", (double)beta1);
129: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
132: beta1 = PetscSqrtReal(beta1);
134: rnorm = beta1;
135: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
136: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
137: PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
138: PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
139: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
141: relres = rnorm / beta1;
142: betan = beta1;
143: phi = beta1;
144: betan = beta1;
145: beta = 0.0;
146: do {
147: /* Lanczos */
148: ksp->its++;
149: betal = beta;
150: beta = betan;
151: PetscCall(VecAXPBY(V, 1.0 / beta, 0.0, R3));
152: PetscCall(KSP_MatMult(ksp, Amat, V, R3));
153: if (ksp->its > 1) PetscCall(VecAXPY(R3, -beta / betal, R1));
154: PetscCall(VecDotRealPart(R3, V, &alpha));
155: PetscCall(VecAXPY(R3, -alpha / beta, R2));
156: KSPMinresSwap3(R1, R2, R3);
158: PetscCall(KSP_PCApply(ksp, R2, R3));
159: PetscCall(VecDotRealPart(R3, R2, &betan));
160: KSPCheckDot(ksp, betan);
161: if (betan < 0.0) {
162: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite preconditioner %g", (double)betan);
163: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
166: betan = PetscSqrtReal(betan);
168: pnorm = Norm3(betal, alpha, betan);
170: // Apply previous left rotation Q_{k-1}
171: dbar = dltan;
172: epln = eplnn;
173: dlta = cs * dbar + sn * alpha;
174: gbar = sn * dbar - cs * alpha;
175: eplnn = sn * betan;
176: dltan = -cs * betan;
177: dlta_QLP = dlta;
179: // Stop if negative curvature is detected and return residual
180: // This is very experimental and maybe changed in the future
181: // based on https://arxiv.org/pdf/2208.07095.pdf
182: if (ksp->converged_neg_curve) {
183: if (cs * gbar >= 0.0) {
184: PetscCall(PetscInfo(ksp, "Detected negative curvature c_nm1 %g, gbar %g\n", (double)cs, (double)gbar));
185: ksp->reason = KSP_CONVERGED_NEG_CURVE;
186: PetscCall(VecCopy(RN, X));
187: break;
188: } else {
189: PetscCall(VecAXPBY(RN, -phi * cs, PetscSqr(sn), V));
190: }
191: }
193: // Compute the current left plane rotation Q_k
194: gamal3 = gamal2;
195: gamal2 = gamal;
196: gamal = gama;
197: SymOrtho(gbar, betan, &cs, &sn, &gama);
199: gama_tmp = gama;
200: taul2 = taul;
201: taul = tau;
202: tau = cs * phi;
203: Axnorm = Norm3(Axnorm, tau, zero);
204: phi = sn * phi;
206: //Apply the previous right plane rotation P{k-2,k}
207: if (ksp->its > 2) {
208: veplnl2 = veplnl;
209: etal2 = etal;
210: etal = eta;
211: dlta_tmp = sr2 * vepln - cr2 * dlta;
212: veplnl = cr2 * vepln + sr2 * dlta;
213: dlta = dlta_tmp;
214: eta = sr2 * gama;
215: gama = -cr2 * gama;
216: }
218: // Compute the current right plane rotation P{k-1,k}, P_12, P_23,...
219: if (ksp->its > 1) {
220: SymOrtho(gamal, dlta, &cr1, &sr1, &gamal);
221: vepln = sr1 * gama;
222: gama = -cr1 * gama;
223: }
225: // Update xnorm
226: xnorml = xnorm;
227: ul4 = ul3;
228: ul3 = ul2;
229: if (ksp->its > 2) ul2 = (taul2 - etal2 * ul4 - veplnl2 * ul3) / gamal2;
230: if (ksp->its > 1) ul = (taul - etal * ul3 - veplnl * ul2) / gamal;
231: xnorm_tmp = Norm3(xl2norm, ul2, ul);
232: if (PetscAbsReal(gama) > realmin && xnorm_tmp < minres->maxxnorm) {
233: u = (tau - eta * ul2 - vepln * ul) / gama;
234: if (Norm3(xnorm_tmp, u, zero) > minres->maxxnorm) {
235: u = 0;
236: flag = 6;
237: }
238: } else {
239: u = 0;
240: flag = 9;
241: }
242: xl2norm = Norm3(xl2norm, ul2, zero);
243: xnorm = Norm3(xl2norm, ul, u);
245: // Update w. Update x except if it will become too big
246: //if (Acond < minres->TranCond && flag != flag0 && QLPiter == 0) { // I believe they have a typo in the matlab code
247: if ((Acond < minres->TranCond || !minres->qlp) && flag == flag0 && QLPiter == 0) { // MINRES
248: KSPMinresSwap3(WL2, WL, W);
249: PetscCall(VecAXPBY(W, 1.0 / gama_tmp, 0.0, V));
250: if (ksp->its > 1) {
251: Vec T[] = {WL, WL2};
252: PetscScalar alphas[] = {-dlta_QLP / gama_tmp, -epln / gama_tmp};
253: PetscInt nv = (ksp->its == 2 ? 1 : 2);
255: PetscCall(VecMAXPY(W, nv, alphas, T));
256: }
257: if (xnorm < minres->maxxnorm) {
258: PetscCall(VecAXPY(X, tau, W));
259: } else {
260: flag = 6;
261: }
262: } else if (minres->qlp) { //MINRES-QLP updates
263: QLPiter = QLPiter + 1;
264: if (QLPiter == 1) {
265: // xl2 = x - wl*ul_QLP - w*u_QLP;
266: PetscScalar maxpys[] = {1.0, -ul_QLP, -u_QLP};
267: Vec maxpyv[] = {X, WL, W};
269: PetscCall(VecSet(XL2, 0.0));
270: // construct w_{k-3}, w_{k-2}, w_{k-1}
271: if (ksp->its > 1) {
272: if (ksp->its > 3) { // w_{k-3}
273: //wl2 = gamal3*wl2 + veplnl2*wl + etal*w;
274: PetscCall(VecAXPBYPCZ(WL2, veplnl2, etal, gamal3, WL, W));
275: }
276: if (ksp->its > 2) { // w_{k-2}
277: //wl = gamal_QLP*wl + vepln_QLP*w;
278: PetscCall(VecAXPBY(WL, vepln_QLP, gamal_QLP, W));
279: }
280: // w = gama_QLP*w;
281: PetscCall(VecScale(W, gama_QLP));
282: // xl2 = x - wl*ul_QLP - w*u_QLP;
283: PetscCall(VecMAXPY(XL2, 3, maxpys, maxpyv));
284: }
285: }
286: if (ksp->its == 1) {
287: //wl2 = wl; wl = v*sr1; w = -v*cr1;
288: PetscCall(VecCopy(WL, WL2));
289: PetscCall(VecAXPBY(WL, sr1, 0, V));
290: PetscCall(VecAXPBY(W, -cr1, 0, V));
291: } else if (ksp->its == 2) {
292: //wl2 = wl;
293: //wl = w*cr1 + v*sr1;
294: //w = w*sr1 - v*cr1;
295: PetscCall(VecCopy(WL, WL2));
296: PetscCall(VecAXPBYPCZ(WL, cr1, sr1, 0.0, W, V));
297: PetscCall(VecAXPBY(W, -cr1, sr1, V));
298: } else {
299: //wl2 = wl; wl = w; w = wl2*sr2 - v*cr2;
300: //wl2 = wl2*cr2 + v*sr2; v = wl *cr1 + w*sr1;
301: //w = wl *sr1 - w*cr1; wl = v;
302: PetscCall(VecCopy(WL, WL2));
303: PetscCall(VecCopy(W, WL));
304: PetscCall(VecAXPBYPCZ(W, sr2, -cr2, 0.0, WL2, V));
305: PetscCall(VecAXPBY(WL2, sr2, cr2, V));
306: PetscCall(VecAXPBYPCZ(V, cr1, sr1, 0.0, WL, W));
307: PetscCall(VecAXPBY(W, sr1, -cr1, WL));
308: PetscCall(VecCopy(V, WL));
309: }
311: //xl2 = xl2 + wl2*ul2;
312: PetscCall(VecAXPY(XL2, ul2, WL2));
313: //x = xl2 + wl *ul + w*u;
314: PetscCall(VecCopy(XL2, X));
315: PetscCall(VecAXPBYPCZ(X, ul, u, 1.0, WL, W));
316: }
317: // Compute the next right plane rotation P{k-1,k+1}
318: gamal_tmp = gamal;
319: SymOrtho(gamal, eplnn, &cr2, &sr2, &gamal);
321: //Store quantities for transferring from MINRES to MINRES-QLP
322: gamal_QLP = gamal_tmp;
323: vepln_QLP = vepln;
324: gama_QLP = gama;
325: ul_QLP = ul;
326: u_QLP = u;
328: // Estimate various norms
329: abs_gama = PetscAbsReal(gama);
330: Anorml = Anorm;
331: Anorm = PetscMax(PetscMax(Anorm, pnorm), PetscMax(gamal, abs_gama));
332: if (ksp->its == 1) {
333: gmin = gama;
334: gminl = gmin;
335: } else {
336: gminl2 = gminl;
337: gminl = gmin;
338: gmin = PetscMin(gminl2, PetscMin(gamal, abs_gama));
339: }
340: Acondl = Acond;
341: Acond = Anorm / gmin;
342: rnorml = rnorm;
343: relresl = relres;
344: if (flag != 9) rnorm = phi;
345: relres = rnorm / (Anorm * xnorm + beta1);
346: rootl = Norm3(gbar, dltan, zero);
347: Arnorml = rnorml * rootl;
348: relAresl = rootl / Anorm;
350: // See if any of the stopping criteria are satisfied.
351: epsx = Anorm * xnorm * eps;
352: if (flag == flag0 || flag == 9) {
353: //if (Acond >= Acondlim) flag = 7; // Huge Acond
354: if (epsx >= beta1) flag = 5; // x is an eigenvector
355: if (minres->qlp) { /* We use these indicators only if the QLP variant has been selected */
356: PetscReal t1 = 1.0 + relres;
357: PetscReal t2 = 1.0 + relAresl;
358: if (xnorm >= minres->maxxnorm) flag = 6; // xnorm exceeded its limit
359: if (t2 <= 1) flag = 4; // Accurate LS solution
360: if (t1 <= 1) flag = 3; // Accurate Ax=b solution
361: if (relAresl <= ksp->rtol) flag = 2; // Good enough LS solution
362: if (relres <= ksp->rtol) flag = 1; // Good enough Ax=b solution
363: }
364: }
366: if (flag == 2 || flag == 4 || flag == 6 || flag == 7) {
367: Acond = Acondl;
368: rnorm = rnorml;
369: relres = relresl;
370: }
372: if (minres->monitor) { /* Mimics matlab code with extra flag */
373: PetscCall(PetscViewerPushFormat(minres->viewer, minres->viewer_fmt));
374: if (ksp->its == 1) PetscCall(PetscViewerASCIIPrintf(minres->viewer, " flag rnorm Arnorm Compatible LS Anorm Acond xnorm\n"));
375: PetscCall(PetscViewerASCIIPrintf(minres->viewer, "%s %5d %2d %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\n", QLPiter == 1 ? "P" : " ", (int)ksp->its - 1, (int)flag, (double)rnorml, (double)Arnorml, (double)relresl, (double)relAresl, (double)Anorml, (double)Acondl, (double)xnorml));
376: PetscCall(PetscViewerPopFormat(minres->viewer));
377: }
379: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
380: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
381: PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
382: PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
383: if (!ksp->reason) {
384: switch (flag) {
385: case 1:
386: case 2:
387: case 5: /* XXX */
388: ksp->reason = KSP_CONVERGED_RTOL;
389: break;
390: case 3:
391: case 4:
392: ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN;
393: break;
394: case 6:
395: ksp->reason = KSP_CONVERGED_STEP_LENGTH;
396: break;
397: default:
398: break;
399: }
400: }
401: if (ksp->reason) break;
402: } while (ksp->its < ksp->max_it);
404: if (minres->monitor && flag != 2 && flag != 4 && flag != 6 && flag != 7) {
405: PetscCall(VecNorm(X, NORM_2, &xnorm));
406: PetscCall(KSP_MatMult(ksp, Amat, X, R1));
407: PetscCall(VecAYPX(R1, -1.0, B));
408: PetscCall(VecNorm(R1, NORM_2, &rnorml));
409: PetscCall(KSP_MatMult(ksp, Amat, R1, R2));
410: PetscCall(VecNorm(R2, NORM_2, &Arnorml));
411: relresl = rnorml / (Anorm * xnorm + beta1);
412: relAresl = rnorml > realmin ? Arnorml / (Anorm * rnorml) : 0.0;
413: PetscCall(PetscViewerPushFormat(minres->viewer, minres->viewer_fmt));
414: PetscCall(PetscViewerASCIIPrintf(minres->viewer, "%s %5d %2d %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\n", QLPiter == 1 ? "P" : " ", (int)ksp->its, (int)flag, (double)rnorml, (double)Arnorml, (double)relresl, (double)relAresl, (double)Anorml, (double)Acondl, (double)xnorml));
415: PetscCall(PetscViewerPopFormat(minres->viewer));
416: }
417: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: /* This was the original implementation provided by R. Scheichl */
422: static PetscErrorCode KSPSolve_MINRES_OLD(KSP ksp)
423: {
424: PetscInt i;
425: PetscScalar alpha, beta, betaold, eta, c = 1.0, ceta, cold = 1.0, coold, s = 0.0, sold = 0.0, soold;
426: PetscScalar rho0, rho1, rho2, rho3, dp = 0.0;
427: const PetscScalar none = -1.0;
428: PetscReal np;
429: Vec X, B, R, Z, U, V, W, UOLD, VOLD, WOLD, WOOLD;
430: Mat Amat;
431: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
432: PetscBool diagonalscale;
434: PetscFunctionBegin;
435: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
436: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
438: X = ksp->vec_sol;
439: B = ksp->vec_rhs;
440: R = ksp->work[0];
441: Z = ksp->work[1];
442: U = ksp->work[2];
443: V = ksp->work[3];
444: W = ksp->work[4];
445: UOLD = ksp->work[5];
446: VOLD = ksp->work[6];
447: WOLD = ksp->work[7];
448: WOOLD = ksp->work[8];
450: PetscCall(PCGetOperators(ksp->pc, &Amat, NULL));
452: ksp->its = 0;
454: if (!ksp->guess_zero) {
455: PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - A*x */
456: PetscCall(VecAYPX(R, -1.0, B));
457: } else {
458: PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
459: }
460: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */
461: PetscCall(VecNorm(Z, NORM_2, &np)); /* np <- ||z|| */
462: KSPCheckNorm(ksp, np);
463: PetscCall(VecDot(R, Z, &dp));
464: KSPCheckDot(ksp, dp);
466: if (PetscRealPart(dp) < minres->haptol && np > minres->haptol) {
467: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g tolerance %g", (double)PetscRealPart(dp), (double)minres->haptol);
468: PetscCall(PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
469: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
470: PetscFunctionReturn(PETSC_SUCCESS);
471: }
473: ksp->rnorm = 0.0;
474: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
475: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
476: PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
477: PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
478: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
480: dp = PetscAbsScalar(dp);
481: dp = PetscSqrtScalar(dp);
482: beta = dp; /* beta <- sqrt(r'*z) */
483: eta = beta;
484: PetscCall(VecAXPBY(V, 1.0 / beta, 0, R)); /* v <- r / beta */
485: PetscCall(VecAXPBY(U, 1.0 / beta, 0, Z)); /* u <- z / beta */
487: i = 0;
488: do {
489: ksp->its = i + 1;
491: /* Lanczos */
493: PetscCall(KSP_MatMult(ksp, Amat, U, R)); /* r <- A*u */
494: PetscCall(VecDot(U, R, &alpha)); /* alpha <- r'*u */
495: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */
497: if (ksp->its > 1) {
498: Vec T[2];
499: PetscScalar alphas[] = {-alpha, -beta};
500: /* r <- r - alpha v - beta v_old */
501: T[0] = V;
502: T[1] = VOLD;
503: PetscCall(VecMAXPY(R, 2, alphas, T));
504: /* z <- z - alpha u - beta u_old */
505: T[0] = U;
506: T[1] = UOLD;
507: PetscCall(VecMAXPY(Z, 2, alphas, T));
508: } else {
509: PetscCall(VecAXPY(R, -alpha, V)); /* r <- r - alpha v */
510: PetscCall(VecAXPY(Z, -alpha, U)); /* z <- z - alpha u */
511: }
513: betaold = beta;
515: PetscCall(VecDot(R, Z, &dp));
516: KSPCheckDot(ksp, dp);
517: dp = PetscAbsScalar(dp);
518: beta = PetscSqrtScalar(dp); /* beta <- sqrt(r'*z) */
520: /* QR factorisation */
522: coold = cold;
523: cold = c;
524: soold = sold;
525: sold = s;
527: rho0 = cold * alpha - coold * sold * betaold;
528: rho1 = PetscSqrtScalar(rho0 * rho0 + beta * beta);
529: rho2 = sold * alpha + coold * cold * betaold;
530: rho3 = soold * betaold;
532: /* Stop if negative curvature is detected */
533: if (ksp->converged_neg_curve && PetscRealPart(cold * rho0) <= 0.0) {
534: PetscCall(PetscInfo(ksp, "Detected negative curvature c_nm1=%g, gbar %g\n", (double)PetscRealPart(cold), -(double)PetscRealPart(rho0)));
535: ksp->reason = KSP_CONVERGED_NEG_CURVE;
536: break;
537: }
539: /* Givens rotation */
541: c = rho0 / rho1;
542: s = beta / rho1;
544: /* Update */
545: /* w_oold <- w_old */
546: /* w_old <- w */
547: KSPMinresSwap3(WOOLD, WOLD, W);
549: /* w <- (u - rho2 w_old - rho3 w_oold)/rho1 */
550: PetscCall(VecAXPBY(W, 1.0 / rho1, 0.0, U));
551: if (ksp->its > 1) {
552: Vec T[] = {WOLD, WOOLD};
553: PetscScalar alphas[] = {-rho2 / rho1, -rho3 / rho1};
554: PetscInt nv = (ksp->its == 2 ? 1 : 2);
556: PetscCall(VecMAXPY(W, nv, alphas, T));
557: }
559: ceta = c * eta;
560: PetscCall(VecAXPY(X, ceta, W)); /* x <- x + c eta w */
562: /*
563: when dp is really small we have either convergence or an indefinite operator so compute true
564: residual norm to check for convergence
565: */
566: if (PetscRealPart(dp) < minres->haptol) {
567: PetscCall(PetscInfo(ksp, "Possible indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
568: PetscCall(KSP_MatMult(ksp, Amat, X, VOLD));
569: PetscCall(VecAXPY(VOLD, none, B));
570: PetscCall(VecNorm(VOLD, NORM_2, &np));
571: KSPCheckNorm(ksp, np);
572: } else {
573: /* otherwise compute new residual norm via recurrence relation */
574: np *= PetscAbsScalar(s);
575: }
577: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
578: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
579: PetscCall(KSPMonitor(ksp, i + 1, ksp->rnorm));
580: PetscCall((*ksp->converged)(ksp, i + 1, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
581: if (ksp->reason) break;
583: if (PetscRealPart(dp) < minres->haptol) {
584: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g tolerance %g", (double)PetscRealPart(dp), (double)minres->haptol);
585: PetscCall(PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
586: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
587: break;
588: }
590: eta = -s * eta;
591: KSPMinresSwap3(VOLD, V, R);
592: KSPMinresSwap3(UOLD, U, Z);
593: PetscCall(VecScale(V, 1.0 / beta)); /* v <- r / beta */
594: PetscCall(VecScale(U, 1.0 / beta)); /* u <- z / beta */
596: i++;
597: } while (i < ksp->max_it);
598: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: static PetscErrorCode KSPDestroy_MINRES(KSP ksp)
603: {
604: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
606: PetscFunctionBegin;
607: PetscCall(PetscViewerDestroy(&minres->viewer));
608: PetscCall(PetscFree(ksp->data));
609: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetRadius_C", NULL));
610: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetUseQLP_C", NULL));
611: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESGetUseQLP_C", NULL));
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: static PetscErrorCode KSPMINRESSetUseQLP_MINRES(KSP ksp, PetscBool qlp)
616: {
617: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
619: PetscFunctionBegin;
620: minres->qlp = qlp;
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: static PetscErrorCode KSPMINRESSetRadius_MINRES(KSP ksp, PetscReal radius)
625: {
626: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
628: PetscFunctionBegin;
629: minres->maxxnorm = radius;
630: PetscFunctionReturn(PETSC_SUCCESS);
631: }
633: static PetscErrorCode KSPMINRESGetUseQLP_MINRES(KSP ksp, PetscBool *qlp)
634: {
635: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
637: PetscFunctionBegin;
638: *qlp = minres->qlp;
639: PetscFunctionReturn(PETSC_SUCCESS);
640: }
642: static PetscErrorCode KSPSetFromOptions_MINRES(KSP ksp, PetscOptionItems *PetscOptionsObject)
643: {
644: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
646: PetscFunctionBegin;
647: PetscOptionsHeadBegin(PetscOptionsObject, "KSP MINRES options");
648: { /* Allow comparing with the old code (to be removed in a few releases) */
649: PetscBool flg = PETSC_FALSE;
650: PetscCall(PetscOptionsBool("-ksp_minres_old", "Use old implementation (to be removed)", "None", flg, &flg, NULL));
651: if (flg) ksp->ops->solve = KSPSolve_MINRES_OLD;
652: else ksp->ops->solve = KSPSolve_MINRES;
653: }
654: PetscCall(PetscOptionsBool("-ksp_minres_qlp", "Solve with QLP variant", "KSPMINRESSetUseQLP", minres->qlp, &minres->qlp, NULL));
655: PetscCall(PetscOptionsReal("-ksp_minres_radius", "Maximum allowed norm of solution", "KSPMINRESSetRadius", minres->maxxnorm, &minres->maxxnorm, NULL));
656: PetscCall(PetscOptionsReal("-ksp_minres_trancond", "Threshold on condition number to dynamically switch to QLP", "None", minres->TranCond, &minres->TranCond, NULL));
657: PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp), PetscOptionsObject->options, PetscOptionsObject->prefix, "-ksp_minres_monitor", &minres->viewer, &minres->viewer_fmt, &minres->monitor));
658: PetscOptionsHeadEnd();
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: /*@
663: KSPMINRESSetUseQLP - Use the QLP variant of the algorithm.
665: Logically Collective
667: Input Parameters:
668: + ksp - the iterative context
669: - qlp - a Boolean indicating if the QLP variant should be used
671: Level: beginner
673: Note:
674: By default, the QLP variant is not used.
676: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESGetUseQLP()`
677: @*/
678: PetscErrorCode KSPMINRESSetUseQLP(KSP ksp, PetscBool qlp)
679: {
680: PetscFunctionBegin;
683: PetscTryMethod(ksp, "KSPMINRESSetUseQLP_C", (KSP, PetscBool), (ksp, qlp));
684: PetscFunctionReturn(PETSC_SUCCESS);
685: }
687: /*@
688: KSPMINRESSetRadius - Set the maximum solution norm allowed.
690: Logically Collective
692: Input Parameters:
693: + ksp - the iterative context
694: - radius - the value
696: Level: beginner
698: Options Database Key:
699: . -ksp_minres_radius <real> - maximum allowed solution norm
701: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESSetUseQLP()`
702: @*/
703: PetscErrorCode KSPMINRESSetRadius(KSP ksp, PetscReal radius)
704: {
705: PetscFunctionBegin;
708: PetscTryMethod(ksp, "KSPMINRESSetRadius_C", (KSP, PetscReal), (ksp, radius));
709: PetscFunctionReturn(PETSC_SUCCESS);
710: }
712: /*@
713: KSPMINRESGetUseQLP - Get the flag for the QLP variant.
715: Logically Collective
717: Input Parameter:
718: . ksp - the iterative context
720: Output Parameter:
721: . qlp - a Boolean indicating if the QLP variant is used
723: Level: beginner
725: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESSetUseQLP()`
726: @*/
727: PetscErrorCode KSPMINRESGetUseQLP(KSP ksp, PetscBool *qlp)
728: {
729: PetscFunctionBegin;
732: PetscUseMethod(ksp, "KSPMINRESGetUseQLP_C", (KSP, PetscBool *), (ksp, qlp));
733: PetscFunctionReturn(PETSC_SUCCESS);
734: }
736: /*MC
737: KSPMINRES - This code implements the MINRES (Minimum Residual) method and its QLP variant.
739: Options Database Keys:
740: + -ksp_minres_qlp <bool> - activates QLP code
741: . -ksp_minres_radius <real> - maximum allowed solution norm
742: . -ksp_minres_trancond <real> - threshold on condition number to dynamically switch to QLP iterations when QLP has been activated
743: - -ksp_minres_monitor - monitors convergence quantities
745: Level: beginner
747: Notes:
748: The operator and the preconditioner must be symmetric and the preconditioner must be positive definite for this method.
750: Supports only left preconditioning.
752: Reference:
753: + * - Paige & Saunders, Solution of sparse indefinite systems of linear equations, SIAM J. Numer. Anal. 12, 1975.
754: - * - S.-C. T. Choi, C. C. Paige and M. A. Saunders. MINRES-QLP: A Krylov subspace method for indefinite or singular symmetric systems, SIAM J. Sci. Comput. 33:4, 2011.
756: Original MINRES code contributed by: Robert Scheichl: maprs@maths.bath.ac.uk
757: QLP variant adapted from: https://stanford.edu/group/SOL/software/minresqlp/minresqlp-matlab/CPS11.zip
759: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`, `KSPCR`
760: M*/
761: PETSC_EXTERN PetscErrorCode KSPCreate_MINRES(KSP ksp)
762: {
763: KSP_MINRES *minres;
765: PetscFunctionBegin;
766: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
767: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
768: PetscCall(PetscNew(&minres));
770: /* this parameter is arbitrary and belongs to the old implementation; but e-50 didn't work for __float128 in one example */
771: #if defined(PETSC_USE_REAL___FLOAT128)
772: minres->haptol = 1.e-100;
773: #elif defined(PETSC_USE_REAL_SINGLE)
774: minres->haptol = 1.e-25;
775: #else
776: minres->haptol = 1.e-50;
777: #endif
778: /* those are set as 1.e7 in the matlab code -> use 1.0/sqrt(eps) to support single precision */
779: minres->maxxnorm = 1.0 / PETSC_SQRT_MACHINE_EPSILON;
780: minres->TranCond = 1.0 / PETSC_SQRT_MACHINE_EPSILON;
782: ksp->data = (void *)minres;
784: ksp->ops->setup = KSPSetUp_MINRES;
785: ksp->ops->solve = KSPSolve_MINRES;
786: ksp->ops->destroy = KSPDestroy_MINRES;
787: ksp->ops->setfromoptions = KSPSetFromOptions_MINRES;
788: ksp->ops->buildsolution = KSPBuildSolutionDefault;
789: ksp->ops->buildresidual = KSPBuildResidualDefault;
791: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetRadius_C", KSPMINRESSetRadius_MINRES));
792: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetUseQLP_C", KSPMINRESSetUseQLP_MINRES));
793: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESGetUseQLP_C", KSPMINRESGetUseQLP_MINRES));
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }