Actual source code: gltr.c
1: #include <../src/ksp/ksp/impls/cg/gltr/gltrimpl.h>
2: #include <petscblaslapack.h>
4: #define GLTR_PRECONDITIONED_DIRECTION 0
5: #define GLTR_UNPRECONDITIONED_DIRECTION 1
6: #define GLTR_DIRECTION_TYPES 2
8: static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};
10: /*@
11: KSPGLTRGetMinEig - Get minimum eigenvalue computed by `KSPGLTR`
13: Collective
15: Input Parameter:
16: . ksp - the iterative context
18: Output Parameter:
19: . e_min - the minimum eigenvalue
21: Level: advanced
23: .seealso: [](ch_ksp), `KSP`, `KSPGLTR`, `KSPGLTRGetLambda()`
24: @*/
25: PetscErrorCode KSPGLTRGetMinEig(KSP ksp, PetscReal *e_min)
26: {
27: PetscFunctionBegin;
29: PetscUseMethod(ksp, "KSPGLTRGetMinEig_C", (KSP, PetscReal *), (ksp, e_min));
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: /*@
34: KSPGLTRGetLambda - Get the multiplier on the trust-region constraint when using `KSPGLTR`
36: Not Collective
38: Input Parameter:
39: . ksp - the iterative context
41: Output Parameter:
42: . lambda - the multiplier
44: Level: advanced
46: .seealso: [](ch_ksp), `KSP`, `KSPGLTR`, `KSPGLTRGetMinEig()`
47: @*/
48: PetscErrorCode KSPGLTRGetLambda(KSP ksp, PetscReal *lambda)
49: {
50: PetscFunctionBegin;
52: PetscUseMethod(ksp, "KSPGLTRGetLambda_C", (KSP, PetscReal *), (ksp, lambda));
53: PetscFunctionReturn(PETSC_SUCCESS);
54: }
56: static PetscErrorCode KSPCGSolve_GLTR(KSP ksp)
57: {
58: #if PetscDefined(USE_COMPLEX)
59: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "GLTR is not available for complex systems");
60: #else
61: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
62: PetscReal *t_soln, *t_diag, *t_offd, *e_valu, *e_vect, *e_rwrk;
63: PetscBLASInt *e_iblk, *e_splt, *e_iwrk;
65: Mat Qmat, Mmat;
66: Vec r, z, p, d;
67: PC pc;
69: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
70: PetscReal alpha, beta, kappa, rz, rzm1;
71: PetscReal rr, r2, piv, step;
72: PetscReal vl, vu;
73: PetscReal coef1, coef2, coef3, root1, root2, obj1, obj2;
74: PetscReal norm_t, norm_w, pert;
76: PetscInt i, j, max_cg_its, max_lanczos_its, max_newton_its, sigma;
77: PetscBLASInt t_size = 0, l_size = 0, il, iu, info;
78: PetscBLASInt nrhs, nldb;
80: PetscBLASInt e_valus = 0, e_splts;
82: PetscFunctionBegin;
83: /* Check the arguments and parameters. */
84: PetscCheck(cg->radius >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
86: /* Get the workspace vectors and initialize variables */
87: r2 = cg->radius * cg->radius;
88: r = ksp->work[0];
89: z = ksp->work[1];
90: p = ksp->work[2];
91: d = ksp->vec_sol;
92: pc = ksp->pc;
94: PetscCall(PCGetOperators(pc, &Qmat, &Mmat));
96: PetscCall(VecGetSize(d, &max_cg_its));
97: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
98: max_lanczos_its = cg->max_lanczos_its;
99: max_newton_its = cg->max_newton_its;
100: ksp->its = 0;
102: /* Initialize objective function direction, and minimum eigenvalue. */
103: cg->o_fcn = 0.0;
105: PetscCall(VecSet(d, 0.0)); /* d = 0 */
106: cg->norm_d = 0.0;
108: cg->e_min = 0.0;
109: cg->lambda = 0.0;
111: /*
112: The first phase of GLTR performs a standard conjugate gradient method,
113: but stores the values required for the Lanczos matrix. We switch to
114: the Lanczos when the conjugate gradient method breaks down. Check the
115: right-hand side for numerical problems. The check for not-a-number and
116: infinite values need be performed only once.
117: */
118: PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
119: PetscCall(VecDot(r, r, &rr)); /* rr = r^T r */
120: KSPCheckDot(ksp, rr);
122: /*
123: Check the preconditioner for numerical problems and for positive
124: definiteness. The check for not-a-number and infinite values need be
125: performed only once.
126: */
127: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
128: PetscCall(VecDot(r, z, &rz)); /* rz = r^T inv(M) r */
129: if (PetscIsInfOrNanScalar(rz)) {
130: /*
131: The preconditioner contains not-a-number or an infinite value.
132: Return the gradient direction intersected with the trust region.
133: */
134: ksp->reason = KSP_DIVERGED_NANORINF;
135: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz));
137: if (cg->radius) {
138: if (r2 >= rr) {
139: alpha = 1.0;
140: cg->norm_d = PetscSqrtReal(rr);
141: } else {
142: alpha = PetscSqrtReal(r2 / rr);
143: cg->norm_d = cg->radius;
144: }
146: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
148: /* Compute objective function. */
149: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
150: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
151: PetscCall(VecDot(d, z, &cg->o_fcn));
152: cg->o_fcn = -cg->o_fcn;
153: ++ksp->its;
154: }
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: if (rz < 0.0) {
159: /*
160: The preconditioner is indefinite. Because this is the first
161: and we do not have a direction yet, we use the gradient step. Note
162: that we cannot use the preconditioned norm when computing the step
163: because the matrix is indefinite.
164: */
165: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
166: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz));
168: if (cg->radius) {
169: if (r2 >= rr) {
170: alpha = 1.0;
171: cg->norm_d = PetscSqrtReal(rr);
172: } else {
173: alpha = PetscSqrtReal(r2 / rr);
174: cg->norm_d = cg->radius;
175: }
177: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
179: /* Compute objective function. */
180: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
181: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
182: PetscCall(VecDot(d, z, &cg->o_fcn));
183: cg->o_fcn = -cg->o_fcn;
184: ++ksp->its;
185: }
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: /*
190: As far as we know, the preconditioner is positive semidefinite.
191: Compute and log the residual. Check convergence because this
192: initializes things, but do not terminate until at least one conjugate
193: gradient iteration has been performed.
194: */
195: cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
197: switch (ksp->normtype) {
198: case KSP_NORM_PRECONDITIONED:
199: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
200: break;
202: case KSP_NORM_UNPRECONDITIONED:
203: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
204: break;
206: case KSP_NORM_NATURAL:
207: norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
208: break;
210: default:
211: norm_r = 0.0;
212: break;
213: }
215: PetscCall(KSPLogResidualHistory(ksp, norm_r));
216: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
217: ksp->rnorm = norm_r;
219: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
221: /* Compute the first direction and update the iteration. */
222: PetscCall(VecCopy(z, p)); /* p = z */
223: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
224: ++ksp->its;
226: /* Check the matrix for numerical problems. */
227: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
228: if (PetscIsInfOrNanScalar(kappa)) {
229: /*
230: The matrix produced not-a-number or an infinite value. In this case
231: we must stop and use the gradient direction. This condition need
232: only be checked once.
233: */
234: ksp->reason = KSP_DIVERGED_NANORINF;
235: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa));
237: if (cg->radius) {
238: if (r2 >= rr) {
239: alpha = 1.0;
240: cg->norm_d = PetscSqrtReal(rr);
241: } else {
242: alpha = PetscSqrtReal(r2 / rr);
243: cg->norm_d = cg->radius;
244: }
246: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
248: /* Compute objective function. */
249: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
250: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
251: PetscCall(VecDot(d, z, &cg->o_fcn));
252: cg->o_fcn = -cg->o_fcn;
253: ++ksp->its;
254: }
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: /*
259: Initialize variables for calculating the norm of the direction and for
260: the Lanczos tridiagonal matrix. Note that we track the diagonal value
261: of the Cholesky factorization of the Lanczos matrix in order to
262: determine when negative curvature is encountered.
263: */
264: dMp = 0.0;
265: norm_d = 0.0;
266: switch (cg->dtype) {
267: case GLTR_PRECONDITIONED_DIRECTION:
268: norm_p = rz;
269: break;
271: default:
272: PetscCall(VecDot(p, p, &norm_p));
273: break;
274: }
276: cg->diag[t_size] = kappa / rz;
277: cg->offd[t_size] = 0.0;
278: ++t_size;
280: piv = 1.0;
282: /*
283: Check for breakdown of the conjugate gradient method; this occurs when
284: kappa is zero.
285: */
286: if (PetscAbsReal(kappa) <= 0.0) {
287: /* The curvature is zero. In this case, we must stop and use follow
288: the direction of negative curvature since the Lanczos matrix is zero. */
289: ksp->reason = KSP_DIVERGED_BREAKDOWN;
290: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: kappa=%g\n", (double)kappa));
292: if (cg->radius && norm_p > 0.0) {
293: /* Follow direction of negative curvature to the boundary of the
294: trust region. */
295: step = PetscSqrtReal(r2 / norm_p);
296: cg->norm_d = cg->radius;
298: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
300: /* Update objective function. */
301: cg->o_fcn += step * (0.5 * step * kappa - rz);
302: } else if (cg->radius) {
303: /* The norm of the preconditioned direction is zero; use the gradient
304: step. */
305: if (r2 >= rr) {
306: alpha = 1.0;
307: cg->norm_d = PetscSqrtReal(rr);
308: } else {
309: alpha = PetscSqrtReal(r2 / rr);
310: cg->norm_d = cg->radius;
311: }
313: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
315: /* Compute objective function. */
316: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
317: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
318: PetscCall(VecDot(d, z, &cg->o_fcn));
319: cg->o_fcn = -cg->o_fcn;
320: ++ksp->its;
321: }
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: /*
326: Begin the first part of the GLTR algorithm which runs the conjugate
327: gradient method until either the problem is solved, we encounter the
328: boundary of the trust region, or the conjugate gradient method breaks
329: down.
330: */
331: while (1) {
332: /* Know that kappa is nonzero, because we have not broken down, so we */
333: /* can compute the steplength. */
334: alpha = rz / kappa;
335: cg->alpha[l_size] = alpha;
337: /* Compute the diagonal value of the Cholesky factorization of the */
338: /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
339: /* This indicates a direction of negative curvature. */
340: piv = cg->diag[l_size] - cg->offd[l_size] * cg->offd[l_size] / piv;
341: if (piv <= 0.0) {
342: /* In this case, the matrix is indefinite and we have encountered */
343: /* a direction of negative curvature. Follow the direction to the */
344: /* boundary of the trust region. */
345: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
346: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa));
348: if (cg->radius && norm_p > 0.0) {
349: /* Follow direction of negative curvature to boundary. */
350: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
351: cg->norm_d = cg->radius;
353: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
355: /* Update objective function. */
356: cg->o_fcn += step * (0.5 * step * kappa - rz);
357: } else if (cg->radius) {
358: /* The norm of the direction is zero; there is nothing to follow. */
359: }
360: break;
361: }
363: /* Compute the steplength and check for intersection with the trust */
364: /* region. */
365: norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
366: if (cg->radius && norm_dp1 >= r2) {
367: /* In this case, the matrix is positive definite as far as we know. */
368: /* However, the full step goes beyond the trust region. */
369: ksp->reason = KSP_CONVERGED_STEP_LENGTH;
370: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius));
372: if (norm_p > 0.0) {
373: /* Follow the direction to the boundary of the trust region. */
375: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
376: cg->norm_d = cg->radius;
378: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
380: /* Update objective function. */
381: cg->o_fcn += step * (0.5 * step * kappa - rz);
382: } else {
383: /* The norm of the direction is zero; there is nothing to follow. */
384: }
385: break;
386: }
388: /* Now we can update the direction and residual. */
389: PetscCall(VecAXPY(d, alpha, p)); /* d = d + alpha p */
390: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
391: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
393: switch (cg->dtype) {
394: case GLTR_PRECONDITIONED_DIRECTION:
395: norm_d = norm_dp1;
396: break;
398: default:
399: PetscCall(VecDot(d, d, &norm_d));
400: break;
401: }
402: cg->norm_d = PetscSqrtReal(norm_d);
404: /* Update objective function. */
405: cg->o_fcn -= 0.5 * alpha * rz;
407: /* Check that the preconditioner appears positive semidefinite. */
408: rzm1 = rz;
409: PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
410: if (rz < 0.0) {
411: /* The preconditioner is indefinite. */
412: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
413: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz));
414: break;
415: }
417: /* As far as we know, the preconditioner is positive semidefinite. */
418: /* Compute the residual and check for convergence. */
419: cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
421: switch (ksp->normtype) {
422: case KSP_NORM_PRECONDITIONED:
423: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
424: break;
426: case KSP_NORM_UNPRECONDITIONED:
427: PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
428: break;
430: case KSP_NORM_NATURAL:
431: norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M */
432: break;
434: default:
435: norm_r = 0.0;
436: break;
437: }
439: PetscCall(KSPLogResidualHistory(ksp, norm_r));
440: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
441: ksp->rnorm = norm_r;
443: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
444: if (ksp->reason) {
445: /* The method has converged. */
446: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius));
447: break;
448: }
450: /* We have not converged yet. Check for breakdown. */
451: beta = rz / rzm1;
452: if (PetscAbsReal(beta) <= 0.0) {
453: /* Conjugate gradients has broken down. */
454: ksp->reason = KSP_DIVERGED_BREAKDOWN;
455: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
456: break;
457: }
459: /* Check iteration limit. */
460: if (ksp->its >= max_cg_its) {
461: ksp->reason = KSP_DIVERGED_ITS;
462: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: iterlim: its=%" PetscInt_FMT "\n", ksp->its));
463: break;
464: }
466: /* Update p and the norms. */
467: cg->beta[l_size] = beta;
468: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
470: switch (cg->dtype) {
471: case GLTR_PRECONDITIONED_DIRECTION:
472: dMp = beta * (dMp + alpha * norm_p);
473: norm_p = beta * (rzm1 + beta * norm_p);
474: break;
476: default:
477: PetscCall(VecDot(d, p, &dMp));
478: PetscCall(VecDot(p, p, &norm_p));
479: break;
480: }
482: /* Compute the new direction and update the iteration. */
483: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
484: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
485: ++ksp->its;
487: /* Update the Lanczos tridiagonal matrix. */
488: ++l_size;
489: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
490: cg->diag[t_size] = kappa / rz + beta / alpha;
491: ++t_size;
493: /* Check for breakdown of the conjugate gradient method; this occurs */
494: /* when kappa is zero. */
495: if (PetscAbsReal(kappa) <= 0.0) {
496: /* The method breaks down; move along the direction as if the matrix */
497: /* were indefinite. */
498: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
499: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa));
501: if (cg->radius && norm_p > 0.0) {
502: /* Follow direction to boundary. */
503: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
504: cg->norm_d = cg->radius;
506: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
508: /* Update objective function. */
509: cg->o_fcn += step * (0.5 * step * kappa - rz);
510: } else if (cg->radius) {
511: /* The norm of the direction is zero; there is nothing to follow. */
512: }
513: break;
514: }
515: }
517: /* Check to see if we need to continue with the Lanczos method. */
518: if (!cg->radius) {
519: /* There is no radius. Therefore, we cannot move along the boundary. */
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: if (KSP_CONVERGED_NEG_CURVE != ksp->reason) {
524: /* The method either converged to an interior point, hit the boundary of */
525: /* the trust-region without encountering a direction of negative */
526: /* curvature or the iteration limit was reached. */
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: /* Switch to constructing the Lanczos basis by way of the conjugate */
531: /* directions. */
532: for (i = 0; i < max_lanczos_its; ++i) {
533: /* Check for breakdown of the conjugate gradient method; this occurs */
534: /* when kappa is zero. */
535: if (PetscAbsReal(kappa) <= 0.0) {
536: ksp->reason = KSP_DIVERGED_BREAKDOWN;
537: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa));
538: break;
539: }
541: /* Update the direction and residual. */
542: alpha = rz / kappa;
543: cg->alpha[l_size] = alpha;
545: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
546: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
548: /* Check that the preconditioner appears positive semidefinite. */
549: rzm1 = rz;
550: PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
551: if (rz < 0.0) {
552: /* The preconditioner is indefinite. */
553: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
554: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz));
555: break;
556: }
558: /* As far as we know, the preconditioner is positive definite. Compute */
559: /* the residual. Do NOT check for convergence. */
560: cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
562: switch (ksp->normtype) {
563: case KSP_NORM_PRECONDITIONED:
564: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
565: break;
567: case KSP_NORM_UNPRECONDITIONED:
568: PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
569: break;
571: case KSP_NORM_NATURAL:
572: norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M */
573: break;
575: default:
576: norm_r = 0.0;
577: break;
578: }
580: PetscCall(KSPLogResidualHistory(ksp, norm_r));
581: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
582: ksp->rnorm = norm_r;
584: /* Check for breakdown. */
585: beta = rz / rzm1;
586: if (PetscAbsReal(beta) <= 0.0) {
587: /* Conjugate gradients has broken down. */
588: ksp->reason = KSP_DIVERGED_BREAKDOWN;
589: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
590: break;
591: }
593: /* Update p and the norms. */
594: cg->beta[l_size] = beta;
595: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
597: /* Compute the new direction and update the iteration. */
598: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
599: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
600: ++ksp->its;
602: /* Update the Lanczos tridiagonal matrix. */
603: ++l_size;
604: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
605: cg->diag[t_size] = kappa / rz + beta / alpha;
606: ++t_size;
607: }
609: /*
610: We have the Lanczos basis, solve the tridiagonal trust-region problem
611: to obtain the Lanczos direction. We know that the solution lies on
612: the boundary of the trust region. We start by checking that the
613: workspace allocated is large enough.
615: Note that the current version only computes the solution by using the
616: preconditioned direction. Need to think about how to do the
617: unpreconditioned direction calculation.
618: */
620: if (t_size > cg->alloced) {
621: if (cg->alloced) {
622: PetscCall(PetscFree2(cg->rwork, cg->iwork));
623: cg->alloced += cg->init_alloc;
624: } else {
625: cg->alloced = cg->init_alloc;
626: }
628: while (t_size > cg->alloced) cg->alloced += cg->init_alloc;
630: cg->alloced = PetscMin(cg->alloced, t_size);
631: PetscCall(PetscMalloc2(10 * cg->alloced, &cg->rwork, 5 * cg->alloced, &cg->iwork));
632: }
634: /* Set up the required vectors. */
635: t_soln = cg->rwork + 0 * t_size; /* Solution */
636: t_diag = cg->rwork + 1 * t_size; /* Diagonal of T */
637: t_offd = cg->rwork + 2 * t_size; /* Off-diagonal of T */
638: e_valu = cg->rwork + 3 * t_size; /* Eigenvalues of T */
639: e_vect = cg->rwork + 4 * t_size; /* Eigenvector of T */
640: e_rwrk = cg->rwork + 5 * t_size; /* Eigen workspace */
642: e_iblk = cg->iwork + 0 * t_size; /* Eigen blocks */
643: e_splt = cg->iwork + 1 * t_size; /* Eigen splits */
644: e_iwrk = cg->iwork + 2 * t_size; /* Eigen workspace */
646: /* Compute the minimum eigenvalue of T. */
647: vl = 0.0;
648: vu = 0.0;
649: il = 1;
650: iu = 1;
652: PetscCallBLAS("LAPACKstebz", LAPACKstebz_("I", "E", &t_size, &vl, &vu, &il, &iu, &cg->eigen_tol, cg->diag, cg->offd + 1, &e_valus, &e_splts, e_valu, e_iblk, e_splt, e_rwrk, e_iwrk, &info));
654: if (0 != info || 1 != e_valus) {
655: /* Calculation of the minimum eigenvalue failed. Return the */
656: /* Steihaug-Toint direction. */
657: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n"));
658: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: cg->e_min = e_valu[0];
664: /* Compute the initial value of lambda to make (T + lambda I) positive */
665: /* definite. */
666: pert = cg->init_pert;
667: if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
669: while (1) {
670: for (i = 0; i < t_size; ++i) {
671: t_diag[i] = cg->diag[i] + cg->lambda;
672: t_offd[i] = cg->offd[i];
673: }
675: PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
676: if (0 == info) break;
678: pert += pert;
679: cg->lambda = cg->lambda * (1.0 + pert) + pert;
680: }
682: /* Compute the initial step and its norm. */
683: nrhs = 1;
684: nldb = t_size;
686: t_soln[0] = -cg->norm_r[0];
687: for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
689: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
690: if (0 != info) {
691: /* Calculation of the initial step failed; return the Steihaug-Toint */
692: /* direction. */
693: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
694: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
695: PetscFunctionReturn(PETSC_SUCCESS);
696: }
698: norm_t = 0.;
699: for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
700: norm_t = PetscSqrtReal(norm_t);
702: /* Determine the case we are in. */
703: if (norm_t <= cg->radius) {
704: /* The step is within the trust region; check if we are in the hard case */
705: /* and need to move to the boundary by following a direction of negative */
706: /* curvature. */
707: if (e_valu[0] <= 0.0 && norm_t < cg->radius) {
708: /* This is the hard case; compute the eigenvector associated with the */
709: /* minimum eigenvalue and move along this direction to the boundary. */
710: PetscCallBLAS("LAPACKstein", LAPACKstein_(&t_size, cg->diag, cg->offd + 1, &e_valus, e_valu, e_iblk, e_splt, e_vect, &nldb, e_rwrk, e_iwrk, e_iwrk + t_size, &info));
711: if (0 != info) {
712: /* Calculation of the minimum eigenvalue failed. Return the */
713: /* Steihaug-Toint direction. */
714: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n"));
715: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: coef1 = 0.0;
720: coef2 = 0.0;
721: coef3 = -cg->radius * cg->radius;
722: for (i = 0; i < t_size; ++i) {
723: coef1 += e_vect[i] * e_vect[i];
724: coef2 += e_vect[i] * t_soln[i];
725: coef3 += t_soln[i] * t_soln[i];
726: }
728: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
729: root1 = (-coef2 + coef3) / coef1;
730: root2 = (-coef2 - coef3) / coef1;
732: /* Compute objective value for (t_soln + root1 * e_vect) */
733: for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
735: obj1 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
736: for (i = 1; i < t_size - 1; ++i) obj1 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i] + cg->offd[i + 1] * e_rwrk[i + 1]);
737: obj1 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
739: /* Compute objective value for (t_soln + root2 * e_vect) */
740: for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
742: obj2 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
743: for (i = 1; i < t_size - 1; ++i) obj2 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i] + cg->offd[i + 1] * e_rwrk[i + 1]);
744: obj2 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
746: /* Choose the point with the best objective function value. */
747: if (obj1 <= obj2) {
748: for (i = 0; i < t_size; ++i) t_soln[i] += root1 * e_vect[i];
749: } else {
750: for (i = 0; i < t_size; ++i) t_soln[i] += root2 * e_vect[i];
751: }
752: } else {
753: /* The matrix is positive definite or there was no room to move; the */
754: /* solution is already contained in t_soln. */
755: }
756: } else {
757: /* The step is outside the trust-region. Compute the correct value for */
758: /* lambda by performing Newton's method. */
760: for (i = 0; i < max_newton_its; ++i) {
761: /* Check for convergence. */
762: if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
764: /* Compute the update. */
765: PetscCall(PetscArraycpy(e_rwrk, t_soln, t_size));
767: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
768: if (0 != info) {
769: /* Calculation of the step failed; return the Steihaug-Toint */
770: /* direction. */
771: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
772: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
773: PetscFunctionReturn(PETSC_SUCCESS);
774: }
776: /* Modify lambda. */
777: norm_w = 0.;
778: for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
780: cg->lambda += (norm_t - cg->radius) / cg->radius * (norm_t * norm_t) / norm_w;
782: /* Factor T + lambda I */
783: for (j = 0; j < t_size; ++j) {
784: t_diag[j] = cg->diag[j] + cg->lambda;
785: t_offd[j] = cg->offd[j];
786: }
788: PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
789: if (0 != info) {
790: /* Calculation of factorization failed; return the Steihaug-Toint */
791: /* direction. */
792: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n"));
793: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: /* Compute the new step and its norm. */
798: t_soln[0] = -cg->norm_r[0];
799: for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
801: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
802: if (0 != info) {
803: /* Calculation of the step failed; return the Steihaug-Toint */
804: /* direction. */
805: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
806: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: norm_t = 0.;
811: for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
812: norm_t = PetscSqrtReal(norm_t);
813: }
815: /* Check for convergence. */
816: if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
817: /* Newton method failed to converge in iteration limit. */
818: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n"));
819: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
820: PetscFunctionReturn(PETSC_SUCCESS);
821: }
822: }
824: /* Recover the norm of the direction and objective function value. */
825: cg->norm_d = norm_t;
827: cg->o_fcn = t_soln[0] * (0.5 * (cg->diag[0] * t_soln[0] + cg->offd[1] * t_soln[1]) + cg->norm_r[0]);
828: for (i = 1; i < t_size - 1; ++i) cg->o_fcn += 0.5 * t_soln[i] * (cg->offd[i] * t_soln[i - 1] + cg->diag[i] * t_soln[i] + cg->offd[i + 1] * t_soln[i + 1]);
829: cg->o_fcn += 0.5 * t_soln[i] * (cg->offd[i] * t_soln[i - 1] + cg->diag[i] * t_soln[i]);
831: /* Recover the direction. */
832: sigma = -1;
834: /* Start conjugate gradient method from the beginning */
835: PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
836: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
838: /* Accumulate Q * s */
839: PetscCall(VecCopy(z, d));
840: PetscCall(VecScale(d, sigma * t_soln[0] / cg->norm_r[0]));
842: /* Compute the first direction. */
843: PetscCall(VecCopy(z, p)); /* p = z */
844: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
845: ++ksp->its;
847: for (i = 0; i < l_size - 1; ++i) {
848: /* Update the residual and direction. */
849: alpha = cg->alpha[i];
850: if (alpha >= 0.0) sigma = -sigma;
852: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
853: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
855: /* Accumulate Q * s */
856: PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));
858: /* Update p. */
859: beta = cg->beta[i];
860: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
861: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
862: ++ksp->its;
863: }
865: /* Update the residual and direction. */
866: alpha = cg->alpha[i];
867: if (alpha >= 0.0) sigma = -sigma;
869: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
870: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
872: /* Accumulate Q * s */
873: PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));
875: /* Set the termination reason. */
876: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
877: PetscFunctionReturn(PETSC_SUCCESS);
878: #endif
879: }
881: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
882: {
883: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
884: PetscInt max_its;
886: PetscFunctionBegin;
887: /* Determine the total maximum number of iterations. */
888: max_its = ksp->max_it + cg->max_lanczos_its + 1;
890: /* Set work vectors needed by conjugate gradient method and allocate */
891: /* workspace for Lanczos matrix. */
892: PetscCall(KSPSetWorkVecs(ksp, 3));
893: if (cg->diag) {
894: PetscCall(PetscArrayzero(cg->diag, max_its));
895: PetscCall(PetscArrayzero(cg->offd, max_its));
896: PetscCall(PetscArrayzero(cg->alpha, max_its));
897: PetscCall(PetscArrayzero(cg->beta, max_its));
898: PetscCall(PetscArrayzero(cg->norm_r, max_its));
899: } else {
900: PetscCall(PetscCalloc5(max_its, &cg->diag, max_its, &cg->offd, max_its, &cg->alpha, max_its, &cg->beta, max_its, &cg->norm_r));
901: }
902: PetscFunctionReturn(PETSC_SUCCESS);
903: }
905: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
906: {
907: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
909: PetscFunctionBegin;
910: PetscCall(PetscFree5(cg->diag, cg->offd, cg->alpha, cg->beta, cg->norm_r));
911: if (cg->alloced) PetscCall(PetscFree2(cg->rwork, cg->iwork));
912: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
913: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
914: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
915: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", NULL));
916: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", NULL));
917: PetscCall(KSPDestroyDefault(ksp));
918: PetscFunctionReturn(PETSC_SUCCESS);
919: }
921: static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
922: {
923: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
925: PetscFunctionBegin;
926: cg->radius = radius;
927: PetscFunctionReturn(PETSC_SUCCESS);
928: }
930: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
931: {
932: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
934: PetscFunctionBegin;
935: *norm_d = cg->norm_d;
936: PetscFunctionReturn(PETSC_SUCCESS);
937: }
939: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
940: {
941: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
943: PetscFunctionBegin;
944: *o_fcn = cg->o_fcn;
945: PetscFunctionReturn(PETSC_SUCCESS);
946: }
948: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
949: {
950: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
952: PetscFunctionBegin;
953: *e_min = cg->e_min;
954: PetscFunctionReturn(PETSC_SUCCESS);
955: }
957: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
958: {
959: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
961: PetscFunctionBegin;
962: *lambda = cg->lambda;
963: PetscFunctionReturn(PETSC_SUCCESS);
964: }
966: static PetscErrorCode KSPCGSetFromOptions_GLTR(KSP ksp, PetscOptionItems PetscOptionsObject)
967: {
968: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
970: PetscFunctionBegin;
971: PetscOptionsHeadBegin(PetscOptionsObject, "KSP GLTR options");
973: PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));
975: PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));
977: PetscCall(PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL));
978: PetscCall(PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL));
979: PetscCall(PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL));
981: PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL));
982: PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL));
984: PetscOptionsHeadEnd();
985: PetscFunctionReturn(PETSC_SUCCESS);
986: }
988: /*MC
989: KSPGLTR - Code to run conjugate gradient method subject to a constraint on the solution norm, used within trust region methods {cite}`gould1999solving`
991: Options Database Key:
992: . -ksp_cg_radius radius - Trust Region Radius
994: Level: developer
996: Notes:
997: Uses preconditioned conjugate gradient to compute an approximate minimizer of the quadratic function
999: $$
1000: q(s) = g^T * s + .5 * s^T * H * s
1001: $$
1003: subject to the trust region constraint
1005: $$
1006: || s || \le delta,
1007: $$
1009: where
1010: .vb
1011: delta is the trust region radius,
1012: g is the gradient vector,
1013: H is the Hessian approximation,
1014: .ve
1016: `KSPConvergedReason` may have the additional values
1017: + `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
1018: - `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step.
1020: The operator and the preconditioner supplied must be symmetric and positive definite.
1022: This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
1024: .seealso: [](ch_ksp), `KSPQCG`, `KSPNASH`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`, `KSPGLTRGetMinEig()`, `KSPGLTRGetLambda()`, `KSPCG`
1025: M*/
1027: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1028: {
1029: KSPCG_GLTR *cg;
1031: PetscFunctionBegin;
1032: PetscCall(PetscNew(&cg));
1033: cg->radius = 0.0;
1034: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1036: cg->init_pert = 1.0e-8;
1037: cg->eigen_tol = 1.0e-10;
1038: cg->newton_tol = 1.0e-6;
1040: cg->alloced = 0;
1041: cg->init_alloc = 1024;
1043: cg->max_lanczos_its = 20;
1044: cg->max_newton_its = 10;
1046: ksp->data = (void *)cg;
1047: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
1048: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
1049: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
1050: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
1051: PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
1053: /* Sets the functions that are associated with this data structure */
1054: /* (in C++ this is the same as defining virtual functions). */
1056: ksp->ops->setup = KSPCGSetUp_GLTR;
1057: ksp->ops->solve = KSPCGSolve_GLTR;
1058: ksp->ops->destroy = KSPCGDestroy_GLTR;
1059: ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1060: ksp->ops->buildsolution = KSPBuildSolutionDefault;
1061: ksp->ops->buildresidual = KSPBuildResidualDefault;
1062: ksp->ops->view = NULL;
1064: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_GLTR));
1065: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_GLTR));
1066: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_GLTR));
1067: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", KSPGLTRGetMinEig_GLTR));
1068: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", KSPGLTRGetLambda_GLTR));
1069: PetscFunctionReturn(PETSC_SUCCESS);
1070: }