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: }