Actual source code: minres.c
1: #include <petsc/private/kspimpl.h>
2: #include <petscblaslapack.h>
3: PETSC_INTERN PetscErrorCode KSPComputeExtremeSingularValues_MINRES(KSP, PetscReal *, PetscReal *);
4: PETSC_INTERN PetscErrorCode KSPComputeEigenvalues_MINRES(KSP, PetscInt, PetscReal *, PetscReal *, PetscInt *);
6: PetscBool QLPcite = PETSC_FALSE;
7: const char QLPCitation[] = "@article{choi2011minres,\n"
8: " title={MINRES-QLP: A Krylov subspace method for indefinite or singular symmetric systems},\n"
9: " author={Choi, Sou-Cheng T and Paige, Christopher C and Saunders, Michael A},\n"
10: " journal={SIAM Journal on Scientific Computing},\n"
11: " volume={33},\n"
12: " number={4},\n"
13: " pages={1810--1836},\n"
14: " year={2011},\n}\n";
16: typedef struct {
17: PetscReal haptol;
18: PetscReal nutol;
19: PetscBool qlp;
20: PetscReal maxxnorm;
21: PetscReal TranCond;
22: PetscBool monitor;
23: PetscViewer viewer;
24: PetscViewerFormat viewer_fmt;
25: // The following arrays are of size ksp->maxit
26: PetscScalar *e, *d;
27: PetscReal *ee, *dd; /* work space for Lanczos algorithm */
28: } KSP_MINRES;
30: static PetscErrorCode KSPSetUp_MINRES(KSP ksp)
31: {
32: PetscFunctionBegin;
33: PetscCall(KSPSetWorkVecs(ksp, 9));
34: /*
35: If user requested computations of eigenvalues then allocate
36: work space needed
37: */
38: if (ksp->calc_sings) {
39: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
40: PetscInt maxit = ksp->max_it;
41: PetscCall(PetscFree4(minres->e, minres->d, minres->ee, minres->dd));
42: PetscCall(PetscMalloc4(maxit + 1, &minres->e, maxit, &minres->d, maxit, &minres->ee, maxit, &minres->dd));
44: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_MINRES;
45: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_MINRES;
46: }
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: /* Convenience functions */
51: #define KSPMinresSwap3(V1, V2, V3) \
52: do { \
53: Vec T = V1; \
54: V1 = V2; \
55: V2 = V3; \
56: V3 = T; \
57: } while (0)
59: static inline PetscReal Norm3(PetscReal a, PetscReal b, PetscReal c)
60: {
61: return PetscSqrtReal(PetscSqr(a) + PetscSqr(b) + PetscSqr(c));
62: }
64: static inline void SymOrtho(PetscReal a, PetscReal b, PetscReal *c, PetscReal *s, PetscReal *r)
65: {
66: if (b == 0.0) {
67: if (a == 0.0) *c = 1.0;
68: else *c = PetscCopysignReal(1.0, a);
69: *s = 0.0;
70: *r = PetscAbsReal(a);
71: } else if (a == 0.0) {
72: *c = 0.0;
73: *s = PetscCopysignReal(1.0, b);
74: *r = PetscAbsReal(b);
75: } else if (PetscAbsReal(b) > PetscAbsReal(a)) {
76: PetscReal t = a / b;
78: *s = PetscCopysignReal(1.0, b) / PetscSqrtReal(1.0 + t * t);
79: *c = (*s) * t;
80: *r = b / (*s); // computationally better than d = a / c since |c| <= |s|
81: } else {
82: PetscReal t = b / a;
84: *c = PetscCopysignReal(1.0, a) / PetscSqrtReal(1.0 + t * t);
85: *s = (*c) * t;
86: *r = a / (*c); // computationally better than d = b / s since |s| <= |c|
87: }
88: }
90: /*
91: Code adapted from https://stanford.edu/group/SOL/software/minresqlp/minresqlp-matlab/CPS11.zip
92: CSP11/Algorithms/MINRESQLP/minresQLP.m
93: */
94: static PetscErrorCode KSPSolve_MINRES(KSP ksp)
95: {
96: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
97: Mat Amat;
98: Vec X, B, R1, R2, R3, V, W, WL, WL2, XL2, RN;
99: PetscReal alpha, beta, beta1, betan, betal;
100: PetscReal zero = 0.0, dbar, dltan = 0.0, dlta, cs = -1.0, sn = 0.0, epln, eplnn = 0.0, gbar, dlta_QLP;
101: PetscReal gamal3 = 0.0, gamal2 = 0.0, gamal = 0.0, gama = 0.0, gama_tmp;
102: PetscReal taul2 = 0.0, taul = 0.0, tau = 0.0, phi, phi0, phir;
103: PetscReal Axnorm, xnorm, xnorm_tmp, xl2norm = 0.0, pnorm, Anorm = 0.0, gmin = 0.0, gminl = 0.0, gminl2 = 0.0;
104: PetscReal Acond = 1.0, Acondl = 0.0, rnorml, rnorm, rootl, relAresl, relres, relresl, Arnorml, Anorml = 0.0;
105: PetscReal epsx, realmin = PETSC_REAL_MIN, eps = PETSC_MACHINE_EPSILON;
106: PetscReal veplnl2 = 0.0, veplnl = 0.0, vepln = 0.0, etal2 = 0.0, etal = 0.0, eta = 0.0;
107: PetscReal dlta_tmp, sr2 = 0.0, cr2 = -1.0, cr1 = -1.0, sr1 = 0.0;
108: PetscReal ul4 = 0.0, ul3 = 0.0, ul2 = 0.0, ul = 0.0, u = 0.0, ul_QLP = 0.0, u_QLP = 0.0;
109: PetscReal vepln_QLP = 0.0, gamal_QLP = 0.0, gama_QLP = 0.0, gamal_tmp, abs_gama;
110: PetscInt flag = -2, flag0 = -2, QLPiter = 0;
111: PetscInt stored_max_it, eigs;
112: PetscScalar *e = NULL, *d = NULL;
114: PetscFunctionBegin;
115: PetscCall(PetscCitationsRegister(QLPCitation, &QLPcite));
116: eigs = ksp->calc_sings;
117: stored_max_it = ksp->max_it;
118: if (eigs) {
119: e = minres->e;
120: d = minres->d;
121: }
123: X = ksp->vec_sol;
124: B = ksp->vec_rhs;
125: R1 = ksp->work[0];
126: R2 = ksp->work[1];
127: R3 = ksp->work[2];
128: V = ksp->work[3];
129: W = ksp->work[4];
130: WL = ksp->work[5];
131: WL2 = ksp->work[6];
132: XL2 = ksp->work[7];
133: RN = ksp->work[8];
134: PetscCall(PCGetOperators(ksp->pc, &Amat, NULL));
136: ksp->its = 0;
137: ksp->rnorm = 0.0;
138: if (!ksp->guess_zero) {
139: PetscCall(KSP_MatMult(ksp, Amat, X, R2));
140: PetscCall(VecNorm(R2, NORM_2, &Axnorm));
141: PetscCall(VecNorm(X, NORM_2, &xnorm));
142: PetscCall(VecAYPX(R2, -1.0, B));
143: } else {
144: PetscCall(VecCopy(B, R2));
145: Axnorm = 0.0;
146: xnorm = 0.0;
147: }
148: PetscCall(KSP_PCApply(ksp, R2, R3));
149: if (ksp->converged_neg_curve) PetscCall(VecCopy(R3, RN));
150: PetscCall(VecDotRealPart(R3, R2, &beta1));
151: KSPCheckDot(ksp, beta1);
152: if (beta1 < 0.0) {
153: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g", (double)beta1);
154: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
157: beta1 = PetscSqrtReal(beta1);
159: rnorm = beta1;
160: if (ksp->normtype == KSP_NORM_PRECONDITIONED) ksp->rnorm = rnorm;
161: else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) PetscCall(VecNorm(R2, NORM_2, &ksp->rnorm));
162: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
163: PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
164: PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
165: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
167: relres = rnorm / beta1;
168: betan = beta1;
169: phi0 = beta1;
170: phi = beta1;
171: betan = beta1;
172: beta = 0.0;
173: do {
174: /* Lanczos */
175: ksp->its++;
176: betal = beta;
177: beta = betan;
178: PetscCall(VecAXPBY(V, 1.0 / beta, 0.0, R3));
179: PetscCall(KSP_MatMult(ksp, Amat, V, R3));
180: if (ksp->its > 1) PetscCall(VecAXPY(R3, -beta / betal, R1));
181: PetscCall(VecDotRealPart(R3, V, &alpha));
182: PetscCall(VecAXPY(R3, -alpha / beta, R2));
183: KSPMinresSwap3(R1, R2, R3);
184: if (eigs) {
185: PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Cannot change maxit AND calculate eigenvalues");
186: d[ksp->its - 1] = alpha;
187: e[ksp->its - 1] = beta;
188: }
190: PetscCall(KSP_PCApply(ksp, R2, R3));
191: PetscCall(VecDotRealPart(R3, R2, &betan));
192: KSPCheckDot(ksp, betan);
193: if (betan < 0.0) {
194: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite preconditioner %g", (double)betan);
195: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
198: betan = PetscSqrtReal(betan);
200: pnorm = Norm3(betal, alpha, betan);
202: // Apply previous left rotation Q_{k-1}
203: dbar = dltan;
204: epln = eplnn;
205: dlta = cs * dbar + sn * alpha;
206: gbar = sn * dbar - cs * alpha;
207: eplnn = sn * betan;
208: dltan = -cs * betan;
209: dlta_QLP = dlta;
211: // Stop if negative curvature is detected and return residual
212: // This is very experimental and maybe changed in the future
213: // based on https://arxiv.org/pdf/2208.07095.pdf
214: if (ksp->converged_neg_curve) {
215: if (cs * gbar >= 0.0) {
216: PetscCall(PetscInfo(ksp, "Detected negative curvature c_nm1 %g, gbar %g\n", (double)cs, (double)gbar));
217: ksp->reason = KSP_CONVERGED_NEG_CURVE;
218: PetscCall(VecCopy(RN, X));
219: break;
220: } else {
221: PetscCall(VecAXPBY(RN, -phi * cs, PetscSqr(sn), V));
222: }
223: }
225: // Compute the current left plane rotation Q_k
226: gamal3 = gamal2;
227: gamal2 = gamal;
228: gamal = gama;
229: SymOrtho(gbar, betan, &cs, &sn, &gama);
231: // Inexactness condition from https://arxiv.org/pdf/2208.07095.pdf
232: rootl = Norm3(gbar, dltan, zero);
233: phir = PetscSqr(phi0 / phi);
234: if (ksp->its > 2 && minres->nutol > 0.0) {
235: PetscReal tmp;
237: phir = PetscSqrtReal(phir - 1.0);
238: tmp = rootl / phir;
239: PetscCall(PetscInfo(ksp, "it = %" PetscInt_FMT ": inexact check %g (%g / %g)\n", ksp->its - 2, (double)tmp, (double)rootl, (double)phir));
240: if (tmp < minres->nutol) {
241: ksp->its--;
242: ksp->reason = KSP_CONVERGED_RTOL;
243: break;
244: }
245: }
247: gama_tmp = gama;
248: taul2 = taul;
249: taul = tau;
250: tau = cs * phi;
251: Axnorm = Norm3(Axnorm, tau, zero);
252: phi = sn * phi;
254: //Apply the previous right plane rotation P{k-2,k}
255: if (ksp->its > 2) {
256: veplnl2 = veplnl;
257: etal2 = etal;
258: etal = eta;
259: dlta_tmp = sr2 * vepln - cr2 * dlta;
260: veplnl = cr2 * vepln + sr2 * dlta;
261: dlta = dlta_tmp;
262: eta = sr2 * gama;
263: gama = -cr2 * gama;
264: }
266: // Compute the current right plane rotation P{k-1,k}, P_12, P_23,...
267: if (ksp->its > 1) {
268: SymOrtho(gamal, dlta, &cr1, &sr1, &gamal);
269: vepln = sr1 * gama;
270: gama = -cr1 * gama;
271: }
273: // Update xnorm
274: ul4 = ul3;
275: ul3 = ul2;
276: if (ksp->its > 2) ul2 = (taul2 - etal2 * ul4 - veplnl2 * ul3) / gamal2;
277: if (ksp->its > 1) ul = (taul - etal * ul3 - veplnl * ul2) / gamal;
278: xnorm_tmp = Norm3(xl2norm, ul2, ul);
279: if (PetscAbsReal(gama) > realmin && xnorm_tmp < minres->maxxnorm) {
280: u = (tau - eta * ul2 - vepln * ul) / gama;
281: if (Norm3(xnorm_tmp, u, zero) > minres->maxxnorm) {
282: u = 0;
283: flag = 6;
284: }
285: } else {
286: u = 0;
287: flag = 9;
288: }
289: xl2norm = Norm3(xl2norm, ul2, zero);
290: xnorm = Norm3(xl2norm, ul, u);
292: // Update w. Update x except if it will become too big
293: //if (Acond < minres->TranCond && flag != flag0 && QLPiter == 0) { // I believe they have a typo in the MATLAB code
294: if ((Acond < minres->TranCond || !minres->qlp) && flag == flag0 && QLPiter == 0) { // MINRES
295: KSPMinresSwap3(WL2, WL, W);
296: PetscCall(VecAXPBY(W, 1.0 / gama_tmp, 0.0, V));
297: if (ksp->its > 1) {
298: Vec T[] = {WL, WL2};
299: PetscScalar alphas[] = {-dlta_QLP / gama_tmp, -epln / gama_tmp};
300: PetscInt nv = (ksp->its == 2 ? 1 : 2);
302: PetscCall(VecMAXPY(W, nv, alphas, T));
303: }
304: if (xnorm < minres->maxxnorm) {
305: PetscCall(VecAXPY(X, tau, W));
306: } else {
307: flag = 6;
308: }
309: } else if (minres->qlp) { //MINRES-QLP updates
310: QLPiter = QLPiter + 1;
311: if (QLPiter == 1) {
312: // xl2 = x - wl*ul_QLP - w*u_QLP;
313: PetscScalar maxpys[] = {1.0, -ul_QLP, -u_QLP};
314: Vec maxpyv[] = {X, WL, W};
316: PetscCall(VecSet(XL2, 0.0));
317: // construct w_{k-3}, w_{k-2}, w_{k-1}
318: if (ksp->its > 1) {
319: if (ksp->its > 3) { // w_{k-3}
320: //wl2 = gamal3*wl2 + veplnl2*wl + etal*w;
321: PetscCall(VecAXPBYPCZ(WL2, veplnl2, etal, gamal3, WL, W));
322: }
323: if (ksp->its > 2) { // w_{k-2}
324: //wl = gamal_QLP*wl + vepln_QLP*w;
325: PetscCall(VecAXPBY(WL, vepln_QLP, gamal_QLP, W));
326: }
327: // w = gama_QLP*w;
328: PetscCall(VecScale(W, gama_QLP));
329: // xl2 = x - wl*ul_QLP - w*u_QLP;
330: PetscCall(VecMAXPY(XL2, 3, maxpys, maxpyv));
331: }
332: }
333: if (ksp->its == 1) {
334: //wl2 = wl; wl = v*sr1; w = -v*cr1;
335: PetscCall(VecCopy(WL, WL2));
336: PetscCall(VecAXPBY(WL, sr1, 0, V));
337: PetscCall(VecAXPBY(W, -cr1, 0, V));
338: } else if (ksp->its == 2) {
339: //wl2 = wl;
340: //wl = w*cr1 + v*sr1;
341: //w = w*sr1 - v*cr1;
342: PetscCall(VecCopy(WL, WL2));
343: PetscCall(VecAXPBYPCZ(WL, cr1, sr1, 0.0, W, V));
344: PetscCall(VecAXPBY(W, -cr1, sr1, V));
345: } else {
346: //wl2 = wl; wl = w; w = wl2*sr2 - v*cr2;
347: //wl2 = wl2*cr2 + v*sr2; v = wl *cr1 + w*sr1;
348: //w = wl *sr1 - w*cr1; wl = v;
349: PetscCall(VecCopy(WL, WL2));
350: PetscCall(VecCopy(W, WL));
351: PetscCall(VecAXPBYPCZ(W, sr2, -cr2, 0.0, WL2, V));
352: PetscCall(VecAXPBY(WL2, sr2, cr2, V));
353: PetscCall(VecAXPBYPCZ(V, cr1, sr1, 0.0, WL, W));
354: PetscCall(VecAXPBY(W, sr1, -cr1, WL));
355: PetscCall(VecCopy(V, WL));
356: }
358: //xl2 = xl2 + wl2*ul2;
359: PetscCall(VecAXPY(XL2, ul2, WL2));
360: //x = xl2 + wl *ul + w*u;
361: PetscCall(VecCopy(XL2, X));
362: PetscCall(VecAXPBYPCZ(X, ul, u, 1.0, WL, W));
363: }
364: // Compute the next right plane rotation P{k-1,k+1}
365: gamal_tmp = gamal;
366: SymOrtho(gamal, eplnn, &cr2, &sr2, &gamal);
368: //Store quantities for transferring from MINRES to MINRES-QLP
369: gamal_QLP = gamal_tmp;
370: vepln_QLP = vepln;
371: gama_QLP = gama;
372: ul_QLP = ul;
373: u_QLP = u;
375: // Estimate various norms
376: abs_gama = PetscAbsReal(gama);
377: Anorml = Anorm;
378: Anorm = PetscMax(PetscMax(Anorm, pnorm), PetscMax(gamal, abs_gama));
379: if (ksp->its == 1) {
380: gmin = gama;
381: gminl = gmin;
382: } else {
383: gminl2 = gminl;
384: gminl = gmin;
385: gmin = PetscMin(gminl2, PetscMin(gamal, abs_gama));
386: }
387: Acondl = Acond;
388: Acond = Anorm / gmin;
389: rnorml = rnorm;
390: relresl = relres;
391: if (flag != 9) rnorm = phi;
392: relres = rnorm / (Anorm * xnorm + beta1);
393: Arnorml = rnorml * rootl;
394: relAresl = rootl / Anorm;
396: // See if any of the stopping criteria are satisfied.
397: epsx = Anorm * xnorm * eps;
398: if (flag == flag0 || flag == 9) {
399: //if (Acond >= Acondlim) flag = 7; // Huge Acond
400: if (epsx >= beta1) flag = 5; // x is an eigenvector
401: if (minres->qlp) { /* We use these indicators only if the QLP variant has been selected */
402: PetscReal t1 = 1.0 + relres;
403: PetscReal t2 = 1.0 + relAresl;
404: if (xnorm >= minres->maxxnorm) flag = 6; // xnorm exceeded its limit
405: if (t2 <= 1) flag = 4; // Accurate LS solution
406: if (t1 <= 1) flag = 3; // Accurate Ax=b solution
407: if (relAresl <= ksp->rtol) flag = 2; // Good enough LS solution
408: if (relres <= ksp->rtol) flag = 1; // Good enough Ax=b solution
409: }
410: }
412: if (flag == 2 || flag == 4 || flag == 6 || flag == 7) {
413: Acond = Acondl;
414: rnorm = rnorml;
415: relres = relresl;
416: }
418: if (minres->monitor) { /* Mimics MATLAB code with extra flag */
419: PetscCall(PetscViewerPushFormat(minres->viewer, minres->viewer_fmt));
420: if (ksp->its == 1) PetscCall(PetscViewerASCIIPrintf(minres->viewer, " flag rnorm Arnorm Compatible LS Anorm Acond xnorm\n"));
421: PetscCall(PetscViewerASCIIPrintf(minres->viewer, "%s %5" PetscInt_FMT " %2" PetscInt_FMT " %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\n", QLPiter == 1 ? "P" : " ", ksp->its - 1, flag, (double)rnorml, (double)Arnorml, (double)relresl, (double)relAresl, (double)Anorml, (double)Acondl, (double)xnorm));
422: PetscCall(PetscViewerPopFormat(minres->viewer));
423: }
425: if (ksp->normtype == KSP_NORM_PRECONDITIONED) ksp->rnorm = rnorm;
426: else if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
427: PetscCall(KSP_MatMult(ksp, Amat, X, V));
428: PetscCall(VecAYPX(V, -1.0, B));
429: PetscCall(VecNorm(V, NORM_2, &ksp->rnorm));
430: }
431: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
432: PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
433: PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
434: if (!ksp->reason) {
435: switch (flag) {
436: case 1:
437: case 2:
438: case 5: /* XXX */
439: ksp->reason = KSP_CONVERGED_RTOL;
440: break;
441: case 3:
442: case 4:
443: ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN;
444: break;
445: case 6:
446: ksp->reason = KSP_CONVERGED_STEP_LENGTH;
447: break;
448: default:
449: break;
450: }
451: }
452: if (ksp->reason) break;
453: } while (ksp->its < ksp->max_it);
455: if (minres->monitor && flag != 2 && flag != 4 && flag != 6 && flag != 7) {
456: PetscCall(VecNorm(X, NORM_2, &xnorm));
457: PetscCall(KSP_MatMult(ksp, Amat, X, R1));
458: PetscCall(VecAYPX(R1, -1.0, B));
459: PetscCall(VecNorm(R1, NORM_2, &rnorml));
460: PetscCall(KSP_MatMult(ksp, Amat, R1, R2));
461: PetscCall(VecNorm(R2, NORM_2, &Arnorml));
462: relresl = rnorml / (Anorm * xnorm + beta1);
463: relAresl = rnorml > realmin ? Arnorml / (Anorm * rnorml) : 0.0;
464: PetscCall(PetscViewerPushFormat(minres->viewer, minres->viewer_fmt));
465: PetscCall(PetscViewerASCIIPrintf(minres->viewer, "%s %5" PetscInt_FMT " %2" PetscInt_FMT " %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\n", QLPiter == 1 ? "P" : " ", ksp->its, flag, (double)rnorml, (double)Arnorml, (double)relresl, (double)relAresl, (double)Anorml, (double)Acondl, (double)xnorm));
466: PetscCall(PetscViewerPopFormat(minres->viewer));
467: }
468: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }
472: /* This was the original implementation provided by R. Scheichl */
473: static PetscErrorCode KSPSolve_MINRES_OLD(KSP ksp)
474: {
475: PetscInt i;
476: PetscScalar alpha, beta, betaold, eta, c = 1.0, ceta, cold = 1.0, coold, s = 0.0, sold = 0.0, soold;
477: PetscScalar rho0, rho1, rho2, rho3, dp = 0.0;
478: const PetscScalar none = -1.0;
479: PetscReal np;
480: Vec X, B, R, Z, U, V, W, UOLD, VOLD, WOLD, WOOLD;
481: Mat Amat;
482: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
483: PetscInt stored_max_it, eigs;
484: PetscScalar *e = NULL, *d = NULL;
486: PetscFunctionBegin;
487: X = ksp->vec_sol;
488: B = ksp->vec_rhs;
489: R = ksp->work[0];
490: Z = ksp->work[1];
491: U = ksp->work[2];
492: V = ksp->work[3];
493: W = ksp->work[4];
494: UOLD = ksp->work[5];
495: VOLD = ksp->work[6];
496: WOLD = ksp->work[7];
497: WOOLD = ksp->work[8];
499: PetscCall(PCGetOperators(ksp->pc, &Amat, NULL));
501: ksp->its = 0;
502: eigs = ksp->calc_sings;
503: stored_max_it = ksp->max_it;
504: if (eigs) {
505: e = minres->e;
506: d = minres->d;
507: }
509: if (!ksp->guess_zero) {
510: PetscCall(KSP_MatMult(ksp, Amat, X, R)); /* r <- b - A*x */
511: PetscCall(VecAYPX(R, -1.0, B));
512: } else {
513: PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
514: }
515: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */
516: PetscCall(VecNorm(Z, NORM_2, &np)); /* np <- ||z|| */
517: KSPCheckNorm(ksp, np);
518: PetscCall(VecDot(R, Z, &dp));
519: KSPCheckDot(ksp, dp);
521: if (PetscRealPart(dp) < minres->haptol && np > minres->haptol) {
522: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g tolerance %g", (double)PetscRealPart(dp), (double)minres->haptol);
523: PetscCall(PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
524: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: ksp->rnorm = 0.0;
529: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
530: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
531: PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
532: PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
533: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
535: dp = PetscAbsScalar(dp);
536: dp = PetscSqrtScalar(dp);
537: beta = dp; /* beta <- sqrt(r'*z) */
538: eta = beta;
539: PetscCall(VecAXPBY(V, 1.0 / beta, 0, R)); /* v <- r / beta */
540: PetscCall(VecAXPBY(U, 1.0 / beta, 0, Z)); /* u <- z / beta */
542: i = 0;
543: do {
544: ksp->its = i + 1;
546: /* Lanczos */
548: PetscCall(KSP_MatMult(ksp, Amat, U, R)); /* r <- A*u */
549: PetscCall(VecDot(U, R, &alpha)); /* alpha <- r'*u */
550: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- B*r */
551: if (eigs) {
552: PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Cannot change maxit AND calculate eigenvalues");
553: d[i] = alpha;
554: e[i] = beta;
555: }
557: if (ksp->its > 1) {
558: Vec T[2];
559: PetscScalar alphas[] = {-alpha, -beta};
560: /* r <- r - alpha v - beta v_old */
561: T[0] = V;
562: T[1] = VOLD;
563: PetscCall(VecMAXPY(R, 2, alphas, T));
564: /* z <- z - alpha u - beta u_old */
565: T[0] = U;
566: T[1] = UOLD;
567: PetscCall(VecMAXPY(Z, 2, alphas, T));
568: } else {
569: PetscCall(VecAXPY(R, -alpha, V)); /* r <- r - alpha v */
570: PetscCall(VecAXPY(Z, -alpha, U)); /* z <- z - alpha u */
571: }
573: betaold = beta;
575: PetscCall(VecDot(R, Z, &dp));
576: KSPCheckDot(ksp, dp);
577: dp = PetscAbsScalar(dp);
578: beta = PetscSqrtScalar(dp); /* beta <- sqrt(r'*z) */
580: /* QR factorisation */
582: coold = cold;
583: cold = c;
584: soold = sold;
585: sold = s;
587: rho0 = cold * alpha - coold * sold * betaold;
588: rho1 = PetscSqrtScalar(rho0 * rho0 + beta * beta);
589: rho2 = sold * alpha + coold * cold * betaold;
590: rho3 = soold * betaold;
592: /* Stop if negative curvature is detected */
593: if (ksp->converged_neg_curve && PetscRealPart(cold * rho0) <= 0.0) {
594: PetscCall(PetscInfo(ksp, "Detected negative curvature c_nm1=%g, gbar %g\n", (double)PetscRealPart(cold), -(double)PetscRealPart(rho0)));
595: ksp->reason = KSP_CONVERGED_NEG_CURVE;
596: break;
597: }
599: /* Givens rotation */
601: c = rho0 / rho1;
602: s = beta / rho1;
604: /* Update */
605: /* w_oold <- w_old */
606: /* w_old <- w */
607: KSPMinresSwap3(WOOLD, WOLD, W);
609: /* w <- (u - rho2 w_old - rho3 w_oold)/rho1 */
610: PetscCall(VecAXPBY(W, 1.0 / rho1, 0.0, U));
611: if (ksp->its > 1) {
612: Vec T[] = {WOLD, WOOLD};
613: PetscScalar alphas[] = {-rho2 / rho1, -rho3 / rho1};
614: PetscInt nv = (ksp->its == 2 ? 1 : 2);
616: PetscCall(VecMAXPY(W, nv, alphas, T));
617: }
619: ceta = c * eta;
620: PetscCall(VecAXPY(X, ceta, W)); /* x <- x + c eta w */
622: /*
623: when dp is really small we have either convergence or an indefinite operator so compute true
624: residual norm to check for convergence
625: */
626: if (PetscRealPart(dp) < minres->haptol) {
627: PetscCall(PetscInfo(ksp, "Possible indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
628: PetscCall(KSP_MatMult(ksp, Amat, X, VOLD));
629: PetscCall(VecAXPY(VOLD, none, B));
630: PetscCall(VecNorm(VOLD, NORM_2, &np));
631: KSPCheckNorm(ksp, np);
632: } else {
633: /* otherwise compute new residual norm via recurrence relation */
634: np *= PetscAbsScalar(s);
635: }
637: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
638: PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
639: PetscCall(KSPMonitor(ksp, i + 1, ksp->rnorm));
640: PetscCall((*ksp->converged)(ksp, i + 1, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
641: if (ksp->reason) break;
643: if (PetscRealPart(dp) < minres->haptol) {
644: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g tolerance %g", (double)PetscRealPart(dp), (double)minres->haptol);
645: PetscCall(PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
646: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
647: break;
648: }
650: eta = -s * eta;
651: KSPMinresSwap3(VOLD, V, R);
652: KSPMinresSwap3(UOLD, U, Z);
653: PetscCall(VecScale(V, 1.0 / beta)); /* v <- r / beta */
654: PetscCall(VecScale(U, 1.0 / beta)); /* u <- z / beta */
656: i++;
657: } while (i < ksp->max_it);
658: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: static PetscErrorCode KSPDestroy_MINRES(KSP ksp)
663: {
664: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
666: PetscFunctionBegin;
667: PetscCall(PetscFree4(minres->e, minres->d, minres->ee, minres->dd));
668: PetscCall(PetscViewerDestroy(&minres->viewer));
669: PetscCall(PetscFree(ksp->data));
670: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetRadius_C", NULL));
671: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetUseQLP_C", NULL));
672: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESGetUseQLP_C", NULL));
673: PetscFunctionReturn(PETSC_SUCCESS);
674: }
676: static PetscErrorCode KSPMINRESSetUseQLP_MINRES(KSP ksp, PetscBool qlp)
677: {
678: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
680: PetscFunctionBegin;
681: minres->qlp = qlp;
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: static PetscErrorCode KSPMINRESSetRadius_MINRES(KSP ksp, PetscReal radius)
686: {
687: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
689: PetscFunctionBegin;
690: minres->maxxnorm = radius;
691: PetscFunctionReturn(PETSC_SUCCESS);
692: }
694: static PetscErrorCode KSPMINRESGetUseQLP_MINRES(KSP ksp, PetscBool *qlp)
695: {
696: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
698: PetscFunctionBegin;
699: *qlp = minres->qlp;
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: static PetscErrorCode KSPSetFromOptions_MINRES(KSP ksp, PetscOptionItems PetscOptionsObject)
704: {
705: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
707: PetscFunctionBegin;
708: PetscOptionsHeadBegin(PetscOptionsObject, "KSP MINRES options");
709: { /* Allow comparing with the old code (to be removed in a few releases) */
710: PetscBool flg = PETSC_FALSE;
711: PetscCall(PetscOptionsBool("-ksp_minres_old", "Use old implementation (to be removed)", "None", flg, &flg, NULL));
712: if (flg) ksp->ops->solve = KSPSolve_MINRES_OLD;
713: else ksp->ops->solve = KSPSolve_MINRES;
714: }
715: PetscCall(PetscOptionsBool("-ksp_minres_qlp", "Solve with QLP variant", "KSPMINRESSetUseQLP", minres->qlp, &minres->qlp, NULL));
716: PetscCall(PetscOptionsReal("-ksp_minres_radius", "Maximum allowed norm of solution", "KSPMINRESSetRadius", minres->maxxnorm, &minres->maxxnorm, NULL));
717: PetscCall(PetscOptionsReal("-ksp_minres_trancond", "Threshold on condition number to dynamically switch to QLP", "None", minres->TranCond, &minres->TranCond, NULL));
718: PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)ksp), PetscOptionsObject->options, PetscOptionsObject->prefix, "-ksp_minres_monitor", &minres->viewer, &minres->viewer_fmt, &minres->monitor));
719: PetscCall(PetscOptionsReal("-ksp_minres_nutol", "Inexactness tolerance", NULL, minres->nutol, &minres->nutol, NULL));
720: PetscOptionsHeadEnd();
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }
724: /*@
725: KSPMINRESSetUseQLP - Use the QLP variant of `KSPMINRES`
727: Logically Collective
729: Input Parameters:
730: + ksp - the iterative context
731: - qlp - a Boolean indicating if the QLP variant should be used
733: Level: beginner
735: Note:
736: By default, the QLP variant is not used.
738: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESGetUseQLP()`
739: @*/
740: PetscErrorCode KSPMINRESSetUseQLP(KSP ksp, PetscBool qlp)
741: {
742: PetscFunctionBegin;
745: PetscTryMethod(ksp, "KSPMINRESSetUseQLP_C", (KSP, PetscBool), (ksp, qlp));
746: PetscFunctionReturn(PETSC_SUCCESS);
747: }
749: /*@
750: KSPMINRESSetRadius - Set the maximum solution norm allowed for use with trust region methods
752: Logically Collective
754: Input Parameters:
755: + ksp - the iterative context
756: - radius - the value
758: Level: beginner
760: Options Database Key:
761: . -ksp_minres_radius radius - maximum allowed solution norm
763: Developer Note:
764: Perhaps the KSPXXXSetRadius() should be unified
766: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESSetUseQLP()`
767: @*/
768: PetscErrorCode KSPMINRESSetRadius(KSP ksp, PetscReal radius)
769: {
770: PetscFunctionBegin;
773: PetscTryMethod(ksp, "KSPMINRESSetRadius_C", (KSP, PetscReal), (ksp, radius));
774: PetscFunctionReturn(PETSC_SUCCESS);
775: }
777: /*@
778: KSPMINRESGetUseQLP - Get the flag that indicates if the QLP variant is being used
780: Logically Collective
782: Input Parameter:
783: . ksp - the iterative context
785: Output Parameter:
786: . qlp - a Boolean indicating if the QLP variant is used
788: Level: beginner
790: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESSetUseQLP()`
791: @*/
792: PetscErrorCode KSPMINRESGetUseQLP(KSP ksp, PetscBool *qlp)
793: {
794: PetscFunctionBegin;
796: PetscAssertPointer(qlp, 2);
797: PetscUseMethod(ksp, "KSPMINRESGetUseQLP_C", (KSP, PetscBool *), (ksp, qlp));
798: PetscFunctionReturn(PETSC_SUCCESS);
799: }
801: /*MC
802: KSPMINRES - This code implements the MINRES (Minimum Residual) method and its QLP variant {cite}`paige.saunders:solution`, {cite}`choi2011minres`,
803: {cite}`liu2022newton` for solving linear systems using `KSP`.
805: Options Database Keys:
806: + -ksp_minres_qlp (true|false) - activates QLP code
807: . -ksp_minres_radius maxnorm - maximum allowed solution norm
808: . -ksp_minres_trancond condthreshold - threshold on condition number to dynamically switch to QLP iterations when QLP has been activated
809: . -ksp_minres_monitor - monitors convergence quantities
810: - -ksp_minres_nutol tol - inexactness tolerance (see https://arxiv.org/pdf/2208.07095.pdf)
812: Level: beginner
814: Notes:
815: The matrix (operator) and the preconditioner must be symmetric and the preconditioner must also be positive definite for this method.
817: `KSPMINRES` is often the best Krylov method for symmetric indefinite matrices.
819: Supports only left preconditioning.
821: Contributed by:
822: Original MINRES code - Robert Scheichl: maprs@maths.bath.ac.uk
823: QLP variant adapted from: https://stanford.edu/group/SOL/software/minresqlp/minresqlp-matlab/CPS11.zip
825: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`, `KSPCR`, `KSPMINRESGetUseQLP()`, `KSPMINRESSetUseQLP()`, `KSPMINRESSetRadius()`
826: M*/
827: PETSC_EXTERN PetscErrorCode KSPCreate_MINRES(KSP ksp)
828: {
829: KSP_MINRES *minres;
831: PetscFunctionBegin;
832: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
833: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
834: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
835: PetscCall(PetscNew(&minres));
837: /* this parameter is arbitrary and belongs to the old implementation; but e-50 didn't work for __float128 in one example */
838: #if PetscDefined(USE_REAL___FLOAT128)
839: minres->haptol = 1.e-100;
840: #elif PetscDefined(USE_REAL_SINGLE)
841: minres->haptol = 1.e-25;
842: #else
843: minres->haptol = 1.e-50;
844: #endif
845: /* those are set as 1.e7 in the MATLAB code -> use 1.0/sqrt(eps) to support single precision */
846: minres->maxxnorm = 1.0 / PETSC_SQRT_MACHINE_EPSILON;
847: minres->TranCond = 1.0 / PETSC_SQRT_MACHINE_EPSILON;
849: ksp->data = (void *)minres;
851: ksp->ops->setup = KSPSetUp_MINRES;
852: ksp->ops->solve = KSPSolve_MINRES;
853: ksp->ops->destroy = KSPDestroy_MINRES;
854: ksp->ops->setfromoptions = KSPSetFromOptions_MINRES;
855: ksp->ops->buildsolution = KSPBuildSolutionDefault;
856: ksp->ops->buildresidual = KSPBuildResidualDefault;
858: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetRadius_C", KSPMINRESSetRadius_MINRES));
859: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetUseQLP_C", KSPMINRESSetUseQLP_MINRES));
860: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESGetUseQLP_C", KSPMINRESGetUseQLP_MINRES));
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: PetscErrorCode KSPComputeEigenvalues_MINRES(KSP ksp, PetscInt nmax, PetscReal *r, PetscReal *c, PetscInt *neig)
865: {
866: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
867: PetscScalar *d, *e;
868: PetscReal *ee;
869: PetscInt n = ksp->its;
870: PetscBLASInt bn, ldz = 1;
872: PetscFunctionBegin;
873: PetscCheck(nmax >= n, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_SIZ, "Not enough room in work space r and c for eigenvalues");
874: *neig = n;
876: PetscCall(PetscArrayzero(c, nmax));
877: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
878: d = minres->d;
879: e = minres->e;
880: ee = minres->ee;
882: /* copy tridiagonal matrix to work space */
883: for (PetscInt j = 0; j < n; j++) {
884: r[j] = PetscRealPart(d[j]);
885: ee[j] = PetscRealPart(e[j + 1]);
886: }
888: PetscCall(PetscBLASIntCast(n, &bn));
889: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
890: PetscCallLAPACKInfo("LAPACKREALstev", LAPACKREALstev_("N", &bn, r, ee, NULL, &ldz, NULL, &info));
891: PetscCall(PetscFPTrapPop());
892: PetscCall(PetscSortReal(n, r));
893: PetscFunctionReturn(PETSC_SUCCESS);
894: }
896: PetscErrorCode KSPComputeExtremeSingularValues_MINRES(KSP ksp, PetscReal *emax, PetscReal *emin)
897: {
898: KSP_MINRES *minres = (KSP_MINRES *)ksp->data;
899: PetscScalar *d, *e;
900: PetscReal *dd, *ee;
901: PetscInt n = ksp->its;
902: PetscBLASInt bn, ldz = 1;
904: PetscFunctionBegin;
905: if (!n) {
906: *emax = *emin = 1.0;
907: PetscFunctionReturn(PETSC_SUCCESS);
908: }
909: d = minres->d;
910: e = minres->e;
911: dd = minres->dd;
912: ee = minres->ee;
914: /* copy tridiagonal matrix to work space */
915: for (PetscInt j = 0; j < n; j++) {
916: dd[j] = PetscRealPart(d[j]);
917: ee[j] = PetscRealPart(e[j + 1]);
918: }
920: PetscCall(PetscBLASIntCast(n, &bn));
921: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
922: PetscCallLAPACKInfo("LAPACKREALstev", LAPACKREALstev_("N", &bn, dd, ee, NULL, &ldz, NULL, &info));
923: PetscCall(PetscFPTrapPop());
924: for (PetscInt j = 0; j < n; j++) dd[j] = PetscAbsReal(dd[j]);
925: PetscCall(PetscSortReal(n, dd));
926: *emin = dd[0];
927: *emax = dd[n - 1];
928: PetscFunctionReturn(PETSC_SUCCESS);
929: }