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, &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:   PetscBool    diagonalscale;
101:   PetscReal    zero = 0.0, dbar, dltan = 0.0, dlta, cs = -1.0, sn = 0.0, epln, eplnn = 0.0, gbar, dlta_QLP;
102:   PetscReal    gamal3 = 0.0, gamal2 = 0.0, gamal = 0.0, gama = 0.0, gama_tmp;
103:   PetscReal    taul2 = 0.0, taul = 0.0, tau = 0.0, phi, phi0, phir;
104:   PetscReal    Axnorm, xnorm, xnorm_tmp, xl2norm = 0.0, pnorm, Anorm = 0.0, gmin = 0.0, gminl = 0.0, gminl2 = 0.0;
105:   PetscReal    Acond = 1.0, Acondl = 0.0, rnorml, rnorm, rootl, relAresl, relres, relresl, Arnorml, Anorml = 0.0, xnorml = 0.0;
106:   PetscReal    epsx, realmin = PETSC_REAL_MIN, eps = PETSC_MACHINE_EPSILON;
107:   PetscReal    veplnl2 = 0.0, veplnl = 0.0, vepln = 0.0, etal2 = 0.0, etal = 0.0, eta = 0.0;
108:   PetscReal    dlta_tmp, sr2 = 0.0, cr2 = -1.0, cr1 = -1.0, sr1 = 0.0;
109:   PetscReal    ul4 = 0.0, ul3 = 0.0, ul2 = 0.0, ul = 0.0, u = 0.0, ul_QLP = 0.0, u_QLP = 0.0;
110:   PetscReal    vepln_QLP = 0.0, gamal_QLP = 0.0, gama_QLP = 0.0, gamal_tmp, abs_gama;
111:   PetscInt     flag = -2, flag0 = -2, QLPiter = 0;
112:   PetscInt     stored_max_it, eigs;
113:   PetscScalar *e = NULL, *d = NULL;

115:   PetscFunctionBegin;
116:   PetscCall(PetscCitationsRegister(QLPCitation, &QLPcite));
117:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
118:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

120:   eigs          = ksp->calc_sings;
121:   stored_max_it = ksp->max_it;
122:   if (eigs) {
123:     e = minres->e;
124:     d = minres->d;
125:   }

127:   X   = ksp->vec_sol;
128:   B   = ksp->vec_rhs;
129:   R1  = ksp->work[0];
130:   R2  = ksp->work[1];
131:   R3  = ksp->work[2];
132:   V   = ksp->work[3];
133:   W   = ksp->work[4];
134:   WL  = ksp->work[5];
135:   WL2 = ksp->work[6];
136:   XL2 = ksp->work[7];
137:   RN  = ksp->work[8];
138:   PetscCall(PCGetOperators(ksp->pc, &Amat, NULL));

140:   ksp->its   = 0;
141:   ksp->rnorm = 0.0;
142:   if (!ksp->guess_zero) {
143:     PetscCall(KSP_MatMult(ksp, Amat, X, R2));
144:     PetscCall(VecNorm(R2, NORM_2, &Axnorm));
145:     PetscCall(VecNorm(X, NORM_2, &xnorm));
146:     PetscCall(VecAYPX(R2, -1.0, B));
147:   } else {
148:     PetscCall(VecCopy(B, R2));
149:     Axnorm = 0.0;
150:     xnorm  = 0.0;
151:   }
152:   PetscCall(KSP_PCApply(ksp, R2, R3));
153:   if (ksp->converged_neg_curve) PetscCall(VecCopy(R3, RN));
154:   PetscCall(VecDotRealPart(R3, R2, &beta1));
155:   KSPCheckDot(ksp, beta1);
156:   if (beta1 < 0.0) {
157:     PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g", (double)beta1);
158:     ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
159:     PetscFunctionReturn(PETSC_SUCCESS);
160:   }
161:   beta1 = PetscSqrtReal(beta1);

163:   rnorm = beta1;
164:   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
165:   PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
166:   PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
167:   PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
168:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

170:   relres = rnorm / beta1;
171:   betan  = beta1;
172:   phi0   = beta1;
173:   phi    = beta1;
174:   betan  = beta1;
175:   beta   = 0.0;
176:   do {
177:     /* Lanczos */
178:     ksp->its++;
179:     betal = beta;
180:     beta  = betan;
181:     PetscCall(VecAXPBY(V, 1.0 / beta, 0.0, R3));
182:     PetscCall(KSP_MatMult(ksp, Amat, V, R3));
183:     if (ksp->its > 1) PetscCall(VecAXPY(R3, -beta / betal, R1));
184:     PetscCall(VecDotRealPart(R3, V, &alpha));
185:     PetscCall(VecAXPY(R3, -alpha / beta, R2));
186:     KSPMinresSwap3(R1, R2, R3);
187:     if (eigs) {
188:       PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Cannot change maxit AND calculate eigenvalues");
189:       d[ksp->its - 1] = alpha;
190:       e[ksp->its - 1] = beta;
191:     }

193:     PetscCall(KSP_PCApply(ksp, R2, R3));
194:     PetscCall(VecDotRealPart(R3, R2, &betan));
195:     KSPCheckDot(ksp, betan);
196:     if (betan < 0.0) {
197:       PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite preconditioner %g", (double)betan);
198:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
199:       PetscFunctionReturn(PETSC_SUCCESS);
200:     }
201:     betan = PetscSqrtReal(betan);

203:     pnorm = Norm3(betal, alpha, betan);

205:     // Apply previous left rotation Q_{k-1}
206:     dbar     = dltan;
207:     epln     = eplnn;
208:     dlta     = cs * dbar + sn * alpha;
209:     gbar     = sn * dbar - cs * alpha;
210:     eplnn    = sn * betan;
211:     dltan    = -cs * betan;
212:     dlta_QLP = dlta;

214:     // Stop if negative curvature is detected and return residual
215:     // This is very experimental and maybe changed in the future
216:     // based on https://arxiv.org/pdf/2208.07095.pdf
217:     if (ksp->converged_neg_curve) {
218:       if (cs * gbar >= 0.0) {
219:         PetscCall(PetscInfo(ksp, "Detected negative curvature c_nm1 %g, gbar %g\n", (double)cs, (double)gbar));
220:         ksp->reason = KSP_CONVERGED_NEG_CURVE;
221:         PetscCall(VecCopy(RN, X));
222:         break;
223:       } else {
224:         PetscCall(VecAXPBY(RN, -phi * cs, PetscSqr(sn), V));
225:       }
226:     }

228:     // Compute the current left plane rotation Q_k
229:     gamal3 = gamal2;
230:     gamal2 = gamal;
231:     gamal  = gama;
232:     SymOrtho(gbar, betan, &cs, &sn, &gama);

234:     // Inexactness condition from https://arxiv.org/pdf/2208.07095.pdf
235:     rootl = Norm3(gbar, dltan, zero);
236:     phir  = PetscSqr(phi0 / phi);
237:     if (ksp->its > 2 && minres->nutol > 0.0) {
238:       PetscReal tmp;

240:       phir = PetscSqrtReal(phir - 1.0);
241:       tmp  = rootl / phir;
242:       PetscCall(PetscInfo(ksp, "it = %" PetscInt_FMT ": inexact check %g (%g / %g)\n", ksp->its - 2, (double)tmp, (double)rootl, (double)phir));
243:       if (tmp < minres->nutol) {
244:         ksp->its--;
245:         ksp->reason = KSP_CONVERGED_RTOL;
246:         break;
247:       }
248:     }

250:     gama_tmp = gama;
251:     taul2    = taul;
252:     taul     = tau;
253:     tau      = cs * phi;
254:     Axnorm   = Norm3(Axnorm, tau, zero);
255:     phi      = sn * phi;

257:     //Apply the previous right plane rotation P{k-2,k}
258:     if (ksp->its > 2) {
259:       veplnl2  = veplnl;
260:       etal2    = etal;
261:       etal     = eta;
262:       dlta_tmp = sr2 * vepln - cr2 * dlta;
263:       veplnl   = cr2 * vepln + sr2 * dlta;
264:       dlta     = dlta_tmp;
265:       eta      = sr2 * gama;
266:       gama     = -cr2 * gama;
267:     }

269:     // Compute the current right plane rotation P{k-1,k}, P_12, P_23,...
270:     if (ksp->its > 1) {
271:       SymOrtho(gamal, dlta, &cr1, &sr1, &gamal);
272:       vepln = sr1 * gama;
273:       gama  = -cr1 * gama;
274:     }

276:     // Update xnorm
277:     xnorml = xnorm;
278:     ul4    = ul3;
279:     ul3    = ul2;
280:     if (ksp->its > 2) ul2 = (taul2 - etal2 * ul4 - veplnl2 * ul3) / gamal2;
281:     if (ksp->its > 1) ul = (taul - etal * ul3 - veplnl * ul2) / gamal;
282:     xnorm_tmp = Norm3(xl2norm, ul2, ul);
283:     if (PetscAbsReal(gama) > realmin && xnorm_tmp < minres->maxxnorm) {
284:       u = (tau - eta * ul2 - vepln * ul) / gama;
285:       if (Norm3(xnorm_tmp, u, zero) > minres->maxxnorm) {
286:         u    = 0;
287:         flag = 6;
288:       }
289:     } else {
290:       u    = 0;
291:       flag = 9;
292:     }
293:     xl2norm = Norm3(xl2norm, ul2, zero);
294:     xnorm   = Norm3(xl2norm, ul, u);

296:     // Update w. Update x except if it will become too big
297:     //if (Acond < minres->TranCond && flag != flag0 && QLPiter == 0) { // I believe they have a typo in the MATLAB code
298:     if ((Acond < minres->TranCond || !minres->qlp) && flag == flag0 && QLPiter == 0) { // MINRES
299:       KSPMinresSwap3(WL2, WL, W);
300:       PetscCall(VecAXPBY(W, 1.0 / gama_tmp, 0.0, V));
301:       if (ksp->its > 1) {
302:         Vec         T[]      = {WL, WL2};
303:         PetscScalar alphas[] = {-dlta_QLP / gama_tmp, -epln / gama_tmp};
304:         PetscInt    nv       = (ksp->its == 2 ? 1 : 2);

306:         PetscCall(VecMAXPY(W, nv, alphas, T));
307:       }
308:       if (xnorm < minres->maxxnorm) {
309:         PetscCall(VecAXPY(X, tau, W));
310:       } else {
311:         flag = 6;
312:       }
313:     } else if (minres->qlp) { //MINRES-QLP updates
314:       QLPiter = QLPiter + 1;
315:       if (QLPiter == 1) {
316:         // xl2 = x - wl*ul_QLP - w*u_QLP;
317:         PetscScalar maxpys[] = {1.0, -ul_QLP, -u_QLP};
318:         Vec         maxpyv[] = {X, WL, W};

320:         PetscCall(VecSet(XL2, 0.0));
321:         // construct w_{k-3}, w_{k-2}, w_{k-1}
322:         if (ksp->its > 1) {
323:           if (ksp->its > 3) { // w_{k-3}
324:             //wl2 = gamal3*wl2 + veplnl2*wl + etal*w;
325:             PetscCall(VecAXPBYPCZ(WL2, veplnl2, etal, gamal3, WL, W));
326:           }
327:           if (ksp->its > 2) { // w_{k-2}
328:             //wl = gamal_QLP*wl + vepln_QLP*w;
329:             PetscCall(VecAXPBY(WL, vepln_QLP, gamal_QLP, W));
330:           }
331:           // w = gama_QLP*w;
332:           PetscCall(VecScale(W, gama_QLP));
333:           // xl2 = x - wl*ul_QLP - w*u_QLP;
334:           PetscCall(VecMAXPY(XL2, 3, maxpys, maxpyv));
335:         }
336:       }
337:       if (ksp->its == 1) {
338:         //wl2 = wl;      wl = v*sr1;     w  = -v*cr1;
339:         PetscCall(VecCopy(WL, WL2));
340:         PetscCall(VecAXPBY(WL, sr1, 0, V));
341:         PetscCall(VecAXPBY(W, -cr1, 0, V));
342:       } else if (ksp->its == 2) {
343:         //wl2 = wl;
344:         //wl  = w*cr1 + v*sr1;
345:         //w   = w*sr1 - v*cr1;
346:         PetscCall(VecCopy(WL, WL2));
347:         PetscCall(VecAXPBYPCZ(WL, cr1, sr1, 0.0, W, V));
348:         PetscCall(VecAXPBY(W, -cr1, sr1, V));
349:       } else {
350:         //wl2 = wl;      wl = w;         w  = wl2*sr2 - v*cr2;
351:         //wl2 = wl2*cr2 + v*sr2;         v  = wl *cr1 + w*sr1;
352:         //w   = wl *sr1 - w*cr1;         wl = v;
353:         PetscCall(VecCopy(WL, WL2));
354:         PetscCall(VecCopy(W, WL));
355:         PetscCall(VecAXPBYPCZ(W, sr2, -cr2, 0.0, WL2, V));
356:         PetscCall(VecAXPBY(WL2, sr2, cr2, V));
357:         PetscCall(VecAXPBYPCZ(V, cr1, sr1, 0.0, WL, W));
358:         PetscCall(VecAXPBY(W, sr1, -cr1, WL));
359:         PetscCall(VecCopy(V, WL));
360:       }

362:       //xl2 = xl2 + wl2*ul2;
363:       PetscCall(VecAXPY(XL2, ul2, WL2));
364:       //x   = xl2 + wl *ul + w*u;
365:       PetscCall(VecCopy(XL2, X));
366:       PetscCall(VecAXPBYPCZ(X, ul, u, 1.0, WL, W));
367:     }
368:     // Compute the next right plane rotation P{k-1,k+1}
369:     gamal_tmp = gamal;
370:     SymOrtho(gamal, eplnn, &cr2, &sr2, &gamal);

372:     //Store quantities for transferring from MINRES to MINRES-QLP
373:     gamal_QLP = gamal_tmp;
374:     vepln_QLP = vepln;
375:     gama_QLP  = gama;
376:     ul_QLP    = ul;
377:     u_QLP     = u;

379:     // Estimate various norms
380:     abs_gama = PetscAbsReal(gama);
381:     Anorml   = Anorm;
382:     Anorm    = PetscMax(PetscMax(Anorm, pnorm), PetscMax(gamal, abs_gama));
383:     if (ksp->its == 1) {
384:       gmin  = gama;
385:       gminl = gmin;
386:     } else {
387:       gminl2 = gminl;
388:       gminl  = gmin;
389:       gmin   = PetscMin(gminl2, PetscMin(gamal, abs_gama));
390:     }
391:     Acondl  = Acond;
392:     Acond   = Anorm / gmin;
393:     rnorml  = rnorm;
394:     relresl = relres;
395:     if (flag != 9) rnorm = phi;
396:     relres   = rnorm / (Anorm * xnorm + beta1);
397:     Arnorml  = rnorml * rootl;
398:     relAresl = rootl / Anorm;

400:     // See if any of the stopping criteria are satisfied.
401:     epsx = Anorm * xnorm * eps;
402:     if (flag == flag0 || flag == 9) {
403:       //if (Acond >= Acondlim) flag = 7; // Huge Acond
404:       if (epsx >= beta1) flag = 5; // x is an eigenvector
405:       if (minres->qlp) {           /* We use these indicators only if the QLP variant has been selected */
406:         PetscReal t1 = 1.0 + relres;
407:         PetscReal t2 = 1.0 + relAresl;
408:         if (xnorm >= minres->maxxnorm) flag = 6; // xnorm exceeded its limit
409:         if (t2 <= 1) flag = 4;                   // Accurate LS solution
410:         if (t1 <= 1) flag = 3;                   // Accurate Ax=b solution
411:         if (relAresl <= ksp->rtol) flag = 2;     // Good enough LS solution
412:         if (relres <= ksp->rtol) flag = 1;       // Good enough Ax=b solution
413:       }
414:     }

416:     if (flag == 2 || flag == 4 || flag == 6 || flag == 7) {
417:       Acond  = Acondl;
418:       rnorm  = rnorml;
419:       relres = relresl;
420:     }

422:     if (minres->monitor) { /* Mimics MATLAB code with extra flag */
423:       PetscCall(PetscViewerPushFormat(minres->viewer, minres->viewer_fmt));
424:       if (ksp->its == 1) PetscCall(PetscViewerASCIIPrintf(minres->viewer, "        flag      rnorm     Arnorm   Compatible         LS      Anorm      Acond      xnorm\n"));
425:       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));
426:       PetscCall(PetscViewerPopFormat(minres->viewer));
427:     }

429:     if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
430:     PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
431:     PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
432:     PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
433:     if (!ksp->reason) {
434:       switch (flag) {
435:       case 1:
436:       case 2:
437:       case 5: /* XXX */
438:         ksp->reason = KSP_CONVERGED_RTOL;
439:         break;
440:       case 3:
441:       case 4:
442:         ksp->reason = KSP_CONVERGED_HAPPY_BREAKDOWN;
443:         break;
444:       case 6:
445:         ksp->reason = KSP_CONVERGED_STEP_LENGTH;
446:         break;
447:       default:
448:         break;
449:       }
450:     }
451:     if (ksp->reason) break;
452:   } while (ksp->its < ksp->max_it);

454:   if (minres->monitor && flag != 2 && flag != 4 && flag != 6 && flag != 7) {
455:     PetscCall(VecNorm(X, NORM_2, &xnorm));
456:     PetscCall(KSP_MatMult(ksp, Amat, X, R1));
457:     PetscCall(VecAYPX(R1, -1.0, B));
458:     PetscCall(VecNorm(R1, NORM_2, &rnorml));
459:     PetscCall(KSP_MatMult(ksp, Amat, R1, R2));
460:     PetscCall(VecNorm(R2, NORM_2, &Arnorml));
461:     relresl  = rnorml / (Anorm * xnorm + beta1);
462:     relAresl = rnorml > realmin ? Arnorml / (Anorm * rnorml) : 0.0;
463:     PetscCall(PetscViewerPushFormat(minres->viewer, minres->viewer_fmt));
464:     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));
465:     PetscCall(PetscViewerPopFormat(minres->viewer));
466:   }
467:   if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
468:   PetscFunctionReturn(PETSC_SUCCESS);
469: }

471: /* This was the original implementation provided by R. Scheichl */
472: static PetscErrorCode KSPSolve_MINRES_OLD(KSP ksp)
473: {
474:   PetscInt          i;
475:   PetscScalar       alpha, beta, betaold, eta, c = 1.0, ceta, cold = 1.0, coold, s = 0.0, sold = 0.0, soold;
476:   PetscScalar       rho0, rho1, rho2, rho3, dp = 0.0;
477:   const PetscScalar none = -1.0;
478:   PetscReal         np;
479:   Vec               X, B, R, Z, U, V, W, UOLD, VOLD, WOLD, WOOLD;
480:   Mat               Amat;
481:   KSP_MINRES       *minres = (KSP_MINRES *)ksp->data;
482:   PetscBool         diagonalscale;
483:   PetscInt          stored_max_it, eigs;
484:   PetscScalar      *e = NULL, *d = NULL;

486:   PetscFunctionBegin;
487:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
488:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

490:   X     = ksp->vec_sol;
491:   B     = ksp->vec_rhs;
492:   R     = ksp->work[0];
493:   Z     = ksp->work[1];
494:   U     = ksp->work[2];
495:   V     = ksp->work[3];
496:   W     = ksp->work[4];
497:   UOLD  = ksp->work[5];
498:   VOLD  = ksp->work[6];
499:   WOLD  = ksp->work[7];
500:   WOOLD = ksp->work[8];

502:   PetscCall(PCGetOperators(ksp->pc, &Amat, NULL));

504:   ksp->its      = 0;
505:   eigs          = ksp->calc_sings;
506:   stored_max_it = ksp->max_it;
507:   if (eigs) {
508:     e = minres->e;
509:     d = minres->d;
510:   }

512:   if (!ksp->guess_zero) {
513:     PetscCall(KSP_MatMult(ksp, Amat, X, R)); /*     r <- b - A*x    */
514:     PetscCall(VecAYPX(R, -1.0, B));
515:   } else {
516:     PetscCall(VecCopy(B, R)); /*     r <- b (x is 0) */
517:   }
518:   PetscCall(KSP_PCApply(ksp, R, Z));  /*     z  <- B*r       */
519:   PetscCall(VecNorm(Z, NORM_2, &np)); /*   np <- ||z||        */
520:   KSPCheckNorm(ksp, np);
521:   PetscCall(VecDot(R, Z, &dp));
522:   KSPCheckDot(ksp, dp);

524:   if (PetscRealPart(dp) < minres->haptol && np > minres->haptol) {
525:     PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g tolerance %g", (double)PetscRealPart(dp), (double)minres->haptol);
526:     PetscCall(PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
527:     ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
528:     PetscFunctionReturn(PETSC_SUCCESS);
529:   }

531:   ksp->rnorm = 0.0;
532:   if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
533:   PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
534:   PetscCall(KSPMonitor(ksp, 0, ksp->rnorm));
535:   PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
536:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

538:   dp   = PetscAbsScalar(dp);
539:   dp   = PetscSqrtScalar(dp);
540:   beta = dp; /*  beta <- sqrt(r'*z)  */
541:   eta  = beta;
542:   PetscCall(VecAXPBY(V, 1.0 / beta, 0, R)); /* v <- r / beta */
543:   PetscCall(VecAXPBY(U, 1.0 / beta, 0, Z)); /* u <- z / beta */

545:   i = 0;
546:   do {
547:     ksp->its = i + 1;

549:     /*   Lanczos  */

551:     PetscCall(KSP_MatMult(ksp, Amat, U, R)); /*      r <- A*u   */
552:     PetscCall(VecDot(U, R, &alpha));         /*  alpha <- r'*u  */
553:     PetscCall(KSP_PCApply(ksp, R, Z));       /*      z <- B*r   */
554:     if (eigs) {
555:       PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Cannot change maxit AND calculate eigenvalues");
556:       d[i] = alpha;
557:       e[i] = beta;
558:     }

560:     if (ksp->its > 1) {
561:       Vec         T[2];
562:       PetscScalar alphas[] = {-alpha, -beta};
563:       /*  r <- r - alpha v - beta v_old    */
564:       T[0] = V;
565:       T[1] = VOLD;
566:       PetscCall(VecMAXPY(R, 2, alphas, T));
567:       /*  z <- z - alpha u - beta u_old    */
568:       T[0] = U;
569:       T[1] = UOLD;
570:       PetscCall(VecMAXPY(Z, 2, alphas, T));
571:     } else {
572:       PetscCall(VecAXPY(R, -alpha, V)); /*  r <- r - alpha v     */
573:       PetscCall(VecAXPY(Z, -alpha, U)); /*  z <- z - alpha u     */
574:     }

576:     betaold = beta;

578:     PetscCall(VecDot(R, Z, &dp));
579:     KSPCheckDot(ksp, dp);
580:     dp   = PetscAbsScalar(dp);
581:     beta = PetscSqrtScalar(dp); /*  beta <- sqrt(r'*z)   */

583:     /*    QR factorisation    */

585:     coold = cold;
586:     cold  = c;
587:     soold = sold;
588:     sold  = s;

590:     rho0 = cold * alpha - coold * sold * betaold;
591:     rho1 = PetscSqrtScalar(rho0 * rho0 + beta * beta);
592:     rho2 = sold * alpha + coold * cold * betaold;
593:     rho3 = soold * betaold;

595:     /* Stop if negative curvature is detected */
596:     if (ksp->converged_neg_curve && PetscRealPart(cold * rho0) <= 0.0) {
597:       PetscCall(PetscInfo(ksp, "Detected negative curvature c_nm1=%g, gbar %g\n", (double)PetscRealPart(cold), -(double)PetscRealPart(rho0)));
598:       ksp->reason = KSP_CONVERGED_NEG_CURVE;
599:       break;
600:     }

602:     /*     Givens rotation    */

604:     c = rho0 / rho1;
605:     s = beta / rho1;

607:     /* Update */
608:     /*  w_oold <- w_old */
609:     /*  w_old  <- w     */
610:     KSPMinresSwap3(WOOLD, WOLD, W);

612:     /* w <- (u - rho2 w_old - rho3 w_oold)/rho1 */
613:     PetscCall(VecAXPBY(W, 1.0 / rho1, 0.0, U));
614:     if (ksp->its > 1) {
615:       Vec         T[]      = {WOLD, WOOLD};
616:       PetscScalar alphas[] = {-rho2 / rho1, -rho3 / rho1};
617:       PetscInt    nv       = (ksp->its == 2 ? 1 : 2);

619:       PetscCall(VecMAXPY(W, nv, alphas, T));
620:     }

622:     ceta = c * eta;
623:     PetscCall(VecAXPY(X, ceta, W)); /*  x <- x + c eta w     */

625:     /*
626:         when dp is really small we have either convergence or an indefinite operator so compute true
627:         residual norm to check for convergence
628:     */
629:     if (PetscRealPart(dp) < minres->haptol) {
630:       PetscCall(PetscInfo(ksp, "Possible indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
631:       PetscCall(KSP_MatMult(ksp, Amat, X, VOLD));
632:       PetscCall(VecAXPY(VOLD, none, B));
633:       PetscCall(VecNorm(VOLD, NORM_2, &np));
634:       KSPCheckNorm(ksp, np);
635:     } else {
636:       /* otherwise compute new residual norm via recurrence relation */
637:       np *= PetscAbsScalar(s);
638:     }

640:     if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = np;
641:     PetscCall(KSPLogResidualHistory(ksp, ksp->rnorm));
642:     PetscCall(KSPMonitor(ksp, i + 1, ksp->rnorm));
643:     PetscCall((*ksp->converged)(ksp, i + 1, ksp->rnorm, &ksp->reason, ksp->cnvP)); /* test for convergence */
644:     if (ksp->reason) break;

646:     if (PetscRealPart(dp) < minres->haptol) {
647:       PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_CONV_FAILED, "Detected indefinite operator %g tolerance %g", (double)PetscRealPart(dp), (double)minres->haptol);
648:       PetscCall(PetscInfo(ksp, "Detected indefinite operator %g tolerance %g\n", (double)PetscRealPart(dp), (double)minres->haptol));
649:       ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
650:       break;
651:     }

653:     eta = -s * eta;
654:     KSPMinresSwap3(VOLD, V, R);
655:     KSPMinresSwap3(UOLD, U, Z);
656:     PetscCall(VecScale(V, 1.0 / beta)); /* v <- r / beta */
657:     PetscCall(VecScale(U, 1.0 / beta)); /* u <- z / beta */

659:     i++;
660:   } while (i < ksp->max_it);
661:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
662:   PetscFunctionReturn(PETSC_SUCCESS);
663: }

665: static PetscErrorCode KSPDestroy_MINRES(KSP ksp)
666: {
667:   KSP_MINRES *minres = (KSP_MINRES *)ksp->data;

669:   PetscFunctionBegin;
670:   PetscCall(PetscFree4(minres->e, minres->d, minres->ee, minres->dd));
671:   PetscCall(PetscViewerDestroy(&minres->viewer));
672:   PetscCall(PetscFree(ksp->data));
673:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetRadius_C", NULL));
674:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetUseQLP_C", NULL));
675:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESGetUseQLP_C", NULL));
676:   PetscFunctionReturn(PETSC_SUCCESS);
677: }

679: static PetscErrorCode KSPMINRESSetUseQLP_MINRES(KSP ksp, PetscBool qlp)
680: {
681:   KSP_MINRES *minres = (KSP_MINRES *)ksp->data;

683:   PetscFunctionBegin;
684:   minres->qlp = qlp;
685:   PetscFunctionReturn(PETSC_SUCCESS);
686: }

688: static PetscErrorCode KSPMINRESSetRadius_MINRES(KSP ksp, PetscReal radius)
689: {
690:   KSP_MINRES *minres = (KSP_MINRES *)ksp->data;

692:   PetscFunctionBegin;
693:   minres->maxxnorm = radius;
694:   PetscFunctionReturn(PETSC_SUCCESS);
695: }

697: static PetscErrorCode KSPMINRESGetUseQLP_MINRES(KSP ksp, PetscBool *qlp)
698: {
699:   KSP_MINRES *minres = (KSP_MINRES *)ksp->data;

701:   PetscFunctionBegin;
702:   *qlp = minres->qlp;
703:   PetscFunctionReturn(PETSC_SUCCESS);
704: }

706: static PetscErrorCode KSPSetFromOptions_MINRES(KSP ksp, PetscOptionItems *PetscOptionsObject)
707: {
708:   KSP_MINRES *minres = (KSP_MINRES *)ksp->data;

710:   PetscFunctionBegin;
711:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP MINRES options");
712:   { /* Allow comparing with the old code (to be removed in a few releases) */
713:     PetscBool flg = PETSC_FALSE;
714:     PetscCall(PetscOptionsBool("-ksp_minres_old", "Use old implementation (to be removed)", "None", flg, &flg, NULL));
715:     if (flg) ksp->ops->solve = KSPSolve_MINRES_OLD;
716:     else ksp->ops->solve = KSPSolve_MINRES;
717:   }
718:   PetscCall(PetscOptionsBool("-ksp_minres_qlp", "Solve with QLP variant", "KSPMINRESSetUseQLP", minres->qlp, &minres->qlp, NULL));
719:   PetscCall(PetscOptionsReal("-ksp_minres_radius", "Maximum allowed norm of solution", "KSPMINRESSetRadius", minres->maxxnorm, &minres->maxxnorm, NULL));
720:   PetscCall(PetscOptionsReal("-ksp_minres_trancond", "Threshold on condition number to dynamically switch to QLP", "None", minres->TranCond, &minres->TranCond, NULL));
721:   PetscCall(PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp), PetscOptionsObject->options, PetscOptionsObject->prefix, "-ksp_minres_monitor", &minres->viewer, &minres->viewer_fmt, &minres->monitor));
722:   PetscCall(PetscOptionsReal("-ksp_minres_nutol", "Inexactness tolerance", NULL, minres->nutol, &minres->nutol, NULL));
723:   PetscOptionsHeadEnd();
724:   PetscFunctionReturn(PETSC_SUCCESS);
725: }

727: /*@
728:   KSPMINRESSetUseQLP - Use the QLP variant of `KSPMINRES`

730:   Logically Collective

732:   Input Parameters:
733: + ksp - the iterative context
734: - qlp - a Boolean indicating if the QLP variant should be used

736:   Level: beginner

738:   Note:
739:   By default, the QLP variant is not used.

741: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESGetUseQLP()`
742: @*/
743: PetscErrorCode KSPMINRESSetUseQLP(KSP ksp, PetscBool qlp)
744: {
745:   PetscFunctionBegin;
748:   PetscTryMethod(ksp, "KSPMINRESSetUseQLP_C", (KSP, PetscBool), (ksp, qlp));
749:   PetscFunctionReturn(PETSC_SUCCESS);
750: }

752: /*@
753:   KSPMINRESSetRadius - Set the maximum solution norm allowed for use with trust region methods

755:   Logically Collective

757:   Input Parameters:
758: + ksp    - the iterative context
759: - radius - the value

761:   Level: beginner

763:   Options Database Key:
764: . -ksp_minres_radius <real> - maximum allowed solution norm

766:   Developer Note:
767:   Perhaps the KSPXXXSetRadius() should be unified

769: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESSetUseQLP()`
770: @*/
771: PetscErrorCode KSPMINRESSetRadius(KSP ksp, PetscReal radius)
772: {
773:   PetscFunctionBegin;
776:   PetscTryMethod(ksp, "KSPMINRESSetRadius_C", (KSP, PetscReal), (ksp, radius));
777:   PetscFunctionReturn(PETSC_SUCCESS);
778: }

780: /*@
781:   KSPMINRESGetUseQLP - Get the flag that indicates if the QLP variant is being used

783:   Logically Collective

785:   Input Parameter:
786: . ksp - the iterative context

788:   Output Parameter:
789: . qlp - a Boolean indicating if the QLP variant is used

791:   Level: beginner

793: .seealso: [](ch_ksp), `KSP`, `KSPMINRES`, `KSPMINRESSetUseQLP()`
794: @*/
795: PetscErrorCode KSPMINRESGetUseQLP(KSP ksp, PetscBool *qlp)
796: {
797:   PetscFunctionBegin;
799:   PetscAssertPointer(qlp, 2);
800:   PetscUseMethod(ksp, "KSPMINRESGetUseQLP_C", (KSP, PetscBool *), (ksp, qlp));
801:   PetscFunctionReturn(PETSC_SUCCESS);
802: }

804: /*MC
805:    KSPMINRES - This code implements the MINRES (Minimum Residual) method and its QLP variant {cite}`paige.saunders:solution`, {cite}`choi2011minres`,
806:    {cite}`liu2022newton`.

808:    Options Database Keys:
809: +   -ksp_minres_qlp <bool> - activates QLP code
810: .   -ksp_minres_radius <real> - maximum allowed solution norm
811: .   -ksp_minres_trancond <real> - threshold on condition number to dynamically switch to QLP iterations when QLP has been activated
812: .   -ksp_minres_monitor - monitors convergence quantities
813: -   -ksp_minres_nutol <real> - inexactness tolerance (see https://arxiv.org/pdf/2208.07095.pdf)

815:    Level: beginner

817:    Notes:
818:    The operator and the preconditioner must be symmetric and the preconditioner must be positive definite for this method.

820:    Supports only left preconditioning.

822:    Contributed by:
823:    Original MINRES code - Robert Scheichl: maprs@maths.bath.ac.uk
824:    QLP variant adapted from: https://stanford.edu/group/SOL/software/minresqlp/minresqlp-matlab/CPS11.zip

826: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCG`, `KSPCR`, `KSPMINRESGetUseQLP()`, `KSPMINRESSetUseQLP()`, `KSPMINRESSetRadius()`
827:           `KSPMINRESGetRadius()`
828: M*/
829: PETSC_EXTERN PetscErrorCode KSPCreate_MINRES(KSP ksp)
830: {
831:   KSP_MINRES *minres;

833:   PetscFunctionBegin;
834:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
835:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
836:   PetscCall(PetscNew(&minres));

838:   /* this parameter is arbitrary and belongs to the old implementation; but e-50 didn't work for __float128 in one example */
839: #if defined(PETSC_USE_REAL___FLOAT128)
840:   minres->haptol = 1.e-100;
841: #elif defined(PETSC_USE_REAL_SINGLE)
842:   minres->haptol = 1.e-25;
843: #else
844:   minres->haptol = 1.e-50;
845: #endif
846:   /* those are set as 1.e7 in the MATLAB code -> use 1.0/sqrt(eps) to support single precision */
847:   minres->maxxnorm = 1.0 / PETSC_SQRT_MACHINE_EPSILON;
848:   minres->TranCond = 1.0 / PETSC_SQRT_MACHINE_EPSILON;

850:   ksp->data = (void *)minres;

852:   ksp->ops->setup          = KSPSetUp_MINRES;
853:   ksp->ops->solve          = KSPSolve_MINRES;
854:   ksp->ops->destroy        = KSPDestroy_MINRES;
855:   ksp->ops->setfromoptions = KSPSetFromOptions_MINRES;
856:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
857:   ksp->ops->buildresidual  = KSPBuildResidualDefault;

859:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetRadius_C", KSPMINRESSetRadius_MINRES));
860:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESSetUseQLP_C", KSPMINRESSetUseQLP_MINRES));
861:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPMINRESGetUseQLP_C", KSPMINRESGetUseQLP_MINRES));
862:   PetscFunctionReturn(PETSC_SUCCESS);
863: }

865: PetscErrorCode KSPComputeEigenvalues_MINRES(KSP ksp, PetscInt nmax, PetscReal *r, PetscReal *c, PetscInt *neig)
866: {
867:   KSP_MINRES  *minres = (KSP_MINRES *)ksp->data;
868:   PetscScalar *d, *e;
869:   PetscReal   *ee;
870:   PetscInt     n = ksp->its;
871:   PetscBLASInt bn, lierr = 0, ldz = 1;

873:   PetscFunctionBegin;
874:   PetscCheck(nmax >= n, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_SIZ, "Not enough room in work space r and c for eigenvalues");
875:   *neig = n;

877:   PetscCall(PetscArrayzero(c, nmax));
878:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);
879:   d  = minres->d;
880:   e  = minres->e;
881:   ee = minres->ee;

883:   /* copy tridiagonal matrix to work space */
884:   for (PetscInt j = 0; j < n; j++) {
885:     r[j]  = PetscRealPart(d[j]);
886:     ee[j] = PetscRealPart(e[j]);
887:   }

889:   PetscCall(PetscBLASIntCast(n, &bn));
890:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
891:   PetscCallBLAS("LAPACKREALstev", LAPACKREALstev_("N", &bn, r, &ee[1], NULL, &ldz, NULL, &lierr));
892:   PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEV error");
893:   PetscCall(PetscFPTrapPop());
894:   PetscCall(PetscSortReal(n, r));
895:   PetscFunctionReturn(PETSC_SUCCESS);
896: }

898: PetscErrorCode KSPComputeExtremeSingularValues_MINRES(KSP ksp, PetscReal *emax, PetscReal *emin)
899: {
900:   KSP_MINRES  *minres = (KSP_MINRES *)ksp->data;
901:   PetscScalar *d, *e;
902:   PetscReal   *dd, *ee;
903:   PetscInt     n = ksp->its;
904:   PetscBLASInt bn, lierr = 0, ldz = 1;

906:   PetscFunctionBegin;
907:   if (!n) {
908:     *emax = *emin = 1.0;
909:     PetscFunctionReturn(PETSC_SUCCESS);
910:   }
911:   d  = minres->d;
912:   e  = minres->e;
913:   dd = minres->dd;
914:   ee = minres->ee;

916:   /* copy tridiagonal matrix to work space */
917:   for (PetscInt j = 0; j < n; j++) {
918:     dd[j] = PetscRealPart(d[j]);
919:     ee[j] = PetscRealPart(e[j]);
920:   }

922:   PetscCall(PetscBLASIntCast(n, &bn));
923:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
924:   PetscCallBLAS("LAPACKREALstev", LAPACKREALstev_("N", &bn, dd, &ee[1], NULL, &ldz, NULL, &lierr));
925:   PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_PLIB, "xSTEV error");
926:   PetscCall(PetscFPTrapPop());
927:   for (PetscInt j = 0; j < n; j++) dd[j] = PetscAbsReal(dd[j]);
928:   PetscCall(PetscSortReal(n, dd));
929:   *emin = dd[0];
930:   *emax = dd[n - 1];
931:   PetscFunctionReturn(PETSC_SUCCESS);
932: }