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 defined(PETSC_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;
81: PetscBool diagonalscale;
83: PetscFunctionBegin;
84: /* Check the arguments and parameters. */
85: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
86: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
87: PetscCheck(cg->radius >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
89: /* Get the workspace vectors and initialize variables */
90: r2 = cg->radius * cg->radius;
91: r = ksp->work[0];
92: z = ksp->work[1];
93: p = ksp->work[2];
94: d = ksp->vec_sol;
95: pc = ksp->pc;
97: PetscCall(PCGetOperators(pc, &Qmat, &Mmat));
99: PetscCall(VecGetSize(d, &max_cg_its));
100: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
101: max_lanczos_its = cg->max_lanczos_its;
102: max_newton_its = cg->max_newton_its;
103: ksp->its = 0;
105: /* Initialize objective function direction, and minimum eigenvalue. */
106: cg->o_fcn = 0.0;
108: PetscCall(VecSet(d, 0.0)); /* d = 0 */
109: cg->norm_d = 0.0;
111: cg->e_min = 0.0;
112: cg->lambda = 0.0;
114: /*
115: The first phase of GLTR performs a standard conjugate gradient method,
116: but stores the values required for the Lanczos matrix. We switch to
117: the Lanczos when the conjugate gradient method breaks down. Check the
118: right-hand side for numerical problems. The check for not-a-number and
119: infinite values need be performed only once.
120: */
121: PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
122: PetscCall(VecDot(r, r, &rr)); /* rr = r^T r */
123: KSPCheckDot(ksp, rr);
125: /*
126: Check the preconditioner for numerical problems and for positive
127: definiteness. The check for not-a-number and infinite values need be
128: performed only once.
129: */
130: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
131: PetscCall(VecDot(r, z, &rz)); /* rz = r^T inv(M) r */
132: if (PetscIsInfOrNanScalar(rz)) {
133: /*
134: The preconditioner contains not-a-number or an infinite value.
135: Return the gradient direction intersected with the trust region.
136: */
137: ksp->reason = KSP_DIVERGED_NANORINF;
138: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz));
140: if (cg->radius) {
141: if (r2 >= rr) {
142: alpha = 1.0;
143: cg->norm_d = PetscSqrtReal(rr);
144: } else {
145: alpha = PetscSqrtReal(r2 / rr);
146: cg->norm_d = cg->radius;
147: }
149: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
151: /* Compute objective function. */
152: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
153: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
154: PetscCall(VecDot(d, z, &cg->o_fcn));
155: cg->o_fcn = -cg->o_fcn;
156: ++ksp->its;
157: }
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: if (rz < 0.0) {
162: /*
163: The preconditioner is indefinite. Because this is the first
164: and we do not have a direction yet, we use the gradient step. Note
165: that we cannot use the preconditioned norm when computing the step
166: because the matrix is indefinite.
167: */
168: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
169: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz));
171: if (cg->radius) {
172: if (r2 >= rr) {
173: alpha = 1.0;
174: cg->norm_d = PetscSqrtReal(rr);
175: } else {
176: alpha = PetscSqrtReal(r2 / rr);
177: cg->norm_d = cg->radius;
178: }
180: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
182: /* Compute objective function. */
183: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
184: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
185: PetscCall(VecDot(d, z, &cg->o_fcn));
186: cg->o_fcn = -cg->o_fcn;
187: ++ksp->its;
188: }
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: /*
193: As far as we know, the preconditioner is positive semidefinite.
194: Compute and log the residual. Check convergence because this
195: initializes things, but do not terminate until at least one conjugate
196: gradient iteration has been performed.
197: */
198: cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
200: switch (ksp->normtype) {
201: case KSP_NORM_PRECONDITIONED:
202: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
203: break;
205: case KSP_NORM_UNPRECONDITIONED:
206: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
207: break;
209: case KSP_NORM_NATURAL:
210: norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
211: break;
213: default:
214: norm_r = 0.0;
215: break;
216: }
218: PetscCall(KSPLogResidualHistory(ksp, norm_r));
219: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
220: ksp->rnorm = norm_r;
222: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
224: /* Compute the first direction and update the iteration. */
225: PetscCall(VecCopy(z, p)); /* p = z */
226: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
227: ++ksp->its;
229: /* Check the matrix for numerical problems. */
230: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
231: if (PetscIsInfOrNanScalar(kappa)) {
232: /*
233: The matrix produced not-a-number or an infinite value. In this case
234: we must stop and use the gradient direction. This condition need
235: only be checked once.
236: */
237: ksp->reason = KSP_DIVERGED_NANORINF;
238: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa));
240: if (cg->radius) {
241: if (r2 >= rr) {
242: alpha = 1.0;
243: cg->norm_d = PetscSqrtReal(rr);
244: } else {
245: alpha = PetscSqrtReal(r2 / rr);
246: cg->norm_d = cg->radius;
247: }
249: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
251: /* Compute objective function. */
252: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
253: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
254: PetscCall(VecDot(d, z, &cg->o_fcn));
255: cg->o_fcn = -cg->o_fcn;
256: ++ksp->its;
257: }
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: /*
262: Initialize variables for calculating the norm of the direction and for
263: the Lanczos tridiagonal matrix. Note that we track the diagonal value
264: of the Cholesky factorization of the Lanczos matrix in order to
265: determine when negative curvature is encountered.
266: */
267: dMp = 0.0;
268: norm_d = 0.0;
269: switch (cg->dtype) {
270: case GLTR_PRECONDITIONED_DIRECTION:
271: norm_p = rz;
272: break;
274: default:
275: PetscCall(VecDot(p, p, &norm_p));
276: break;
277: }
279: cg->diag[t_size] = kappa / rz;
280: cg->offd[t_size] = 0.0;
281: ++t_size;
283: piv = 1.0;
285: /*
286: Check for breakdown of the conjugate gradient method; this occurs when
287: kappa is zero.
288: */
289: if (PetscAbsReal(kappa) <= 0.0) {
290: /* The curvature is zero. In this case, we must stop and use follow
291: the direction of negative curvature since the Lanczos matrix is zero. */
292: ksp->reason = KSP_DIVERGED_BREAKDOWN;
293: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: kappa=%g\n", (double)kappa));
295: if (cg->radius && norm_p > 0.0) {
296: /* Follow direction of negative curvature to the boundary of the
297: trust region. */
298: step = PetscSqrtReal(r2 / norm_p);
299: cg->norm_d = cg->radius;
301: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
303: /* Update objective function. */
304: cg->o_fcn += step * (0.5 * step * kappa - rz);
305: } else if (cg->radius) {
306: /* The norm of the preconditioned direction is zero; use the gradient
307: step. */
308: if (r2 >= rr) {
309: alpha = 1.0;
310: cg->norm_d = PetscSqrtReal(rr);
311: } else {
312: alpha = PetscSqrtReal(r2 / rr);
313: cg->norm_d = cg->radius;
314: }
316: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
318: /* Compute objective function. */
319: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
320: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
321: PetscCall(VecDot(d, z, &cg->o_fcn));
322: cg->o_fcn = -cg->o_fcn;
323: ++ksp->its;
324: }
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: /*
329: Begin the first part of the GLTR algorithm which runs the conjugate
330: gradient method until either the problem is solved, we encounter the
331: boundary of the trust region, or the conjugate gradient method breaks
332: down.
333: */
334: while (1) {
335: /* Know that kappa is nonzero, because we have not broken down, so we */
336: /* can compute the steplength. */
337: alpha = rz / kappa;
338: cg->alpha[l_size] = alpha;
340: /* Compute the diagonal value of the Cholesky factorization of the */
341: /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
342: /* This indicates a direction of negative curvature. */
343: piv = cg->diag[l_size] - cg->offd[l_size] * cg->offd[l_size] / piv;
344: if (piv <= 0.0) {
345: /* In this case, the matrix is indefinite and we have encountered */
346: /* a direction of negative curvature. Follow the direction to the */
347: /* boundary of the trust region. */
348: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
349: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa));
351: if (cg->radius && norm_p > 0.0) {
352: /* Follow direction of negative curvature to boundary. */
353: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
354: cg->norm_d = cg->radius;
356: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
358: /* Update objective function. */
359: cg->o_fcn += step * (0.5 * step * kappa - rz);
360: } else if (cg->radius) {
361: /* The norm of the direction is zero; there is nothing to follow. */
362: }
363: break;
364: }
366: /* Compute the steplength and check for intersection with the trust */
367: /* region. */
368: norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
369: if (cg->radius && norm_dp1 >= r2) {
370: /* In this case, the matrix is positive definite as far as we know. */
371: /* However, the full step goes beyond the trust region. */
372: ksp->reason = KSP_CONVERGED_STEP_LENGTH;
373: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius));
375: if (norm_p > 0.0) {
376: /* Follow the direction to the boundary of the trust region. */
378: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
379: cg->norm_d = cg->radius;
381: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
383: /* Update objective function. */
384: cg->o_fcn += step * (0.5 * step * kappa - rz);
385: } else {
386: /* The norm of the direction is zero; there is nothing to follow. */
387: }
388: break;
389: }
391: /* Now we can update the direction and residual. */
392: PetscCall(VecAXPY(d, alpha, p)); /* d = d + alpha p */
393: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
394: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
396: switch (cg->dtype) {
397: case GLTR_PRECONDITIONED_DIRECTION:
398: norm_d = norm_dp1;
399: break;
401: default:
402: PetscCall(VecDot(d, d, &norm_d));
403: break;
404: }
405: cg->norm_d = PetscSqrtReal(norm_d);
407: /* Update objective function. */
408: cg->o_fcn -= 0.5 * alpha * rz;
410: /* Check that the preconditioner appears positive semidefinite. */
411: rzm1 = rz;
412: PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
413: if (rz < 0.0) {
414: /* The preconditioner is indefinite. */
415: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
416: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz));
417: break;
418: }
420: /* As far as we know, the preconditioner is positive semidefinite. */
421: /* Compute the residual and check for convergence. */
422: cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
424: switch (ksp->normtype) {
425: case KSP_NORM_PRECONDITIONED:
426: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
427: break;
429: case KSP_NORM_UNPRECONDITIONED:
430: PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
431: break;
433: case KSP_NORM_NATURAL:
434: norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M */
435: break;
437: default:
438: norm_r = 0.0;
439: break;
440: }
442: PetscCall(KSPLogResidualHistory(ksp, norm_r));
443: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
444: ksp->rnorm = norm_r;
446: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
447: if (ksp->reason) {
448: /* The method has converged. */
449: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius));
450: break;
451: }
453: /* We have not converged yet. Check for breakdown. */
454: beta = rz / rzm1;
455: if (PetscAbsReal(beta) <= 0.0) {
456: /* Conjugate gradients has broken down. */
457: ksp->reason = KSP_DIVERGED_BREAKDOWN;
458: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
459: break;
460: }
462: /* Check iteration limit. */
463: if (ksp->its >= max_cg_its) {
464: ksp->reason = KSP_DIVERGED_ITS;
465: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: iterlim: its=%" PetscInt_FMT "\n", ksp->its));
466: break;
467: }
469: /* Update p and the norms. */
470: cg->beta[l_size] = beta;
471: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
473: switch (cg->dtype) {
474: case GLTR_PRECONDITIONED_DIRECTION:
475: dMp = beta * (dMp + alpha * norm_p);
476: norm_p = beta * (rzm1 + beta * norm_p);
477: break;
479: default:
480: PetscCall(VecDot(d, p, &dMp));
481: PetscCall(VecDot(p, p, &norm_p));
482: break;
483: }
485: /* Compute the new direction and update the iteration. */
486: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
487: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
488: ++ksp->its;
490: /* Update the Lanczos tridiagonal matrix. */
491: ++l_size;
492: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
493: cg->diag[t_size] = kappa / rz + beta / alpha;
494: ++t_size;
496: /* Check for breakdown of the conjugate gradient method; this occurs */
497: /* when kappa is zero. */
498: if (PetscAbsReal(kappa) <= 0.0) {
499: /* The method breaks down; move along the direction as if the matrix */
500: /* were indefinite. */
501: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
502: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa));
504: if (cg->radius && norm_p > 0.0) {
505: /* Follow direction to boundary. */
506: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
507: cg->norm_d = cg->radius;
509: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
511: /* Update objective function. */
512: cg->o_fcn += step * (0.5 * step * kappa - rz);
513: } else if (cg->radius) {
514: /* The norm of the direction is zero; there is nothing to follow. */
515: }
516: break;
517: }
518: }
520: /* Check to see if we need to continue with the Lanczos method. */
521: if (!cg->radius) {
522: /* There is no radius. Therefore, we cannot move along the boundary. */
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: if (KSP_CONVERGED_NEG_CURVE != ksp->reason) {
527: /* The method either converged to an interior point, hit the boundary of */
528: /* the trust-region without encountering a direction of negative */
529: /* curvature or the iteration limit was reached. */
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: /* Switch to constructing the Lanczos basis by way of the conjugate */
534: /* directions. */
535: for (i = 0; i < max_lanczos_its; ++i) {
536: /* Check for breakdown of the conjugate gradient method; this occurs */
537: /* when kappa is zero. */
538: if (PetscAbsReal(kappa) <= 0.0) {
539: ksp->reason = KSP_DIVERGED_BREAKDOWN;
540: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa));
541: break;
542: }
544: /* Update the direction and residual. */
545: alpha = rz / kappa;
546: cg->alpha[l_size] = alpha;
548: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
549: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
551: /* Check that the preconditioner appears positive semidefinite. */
552: rzm1 = rz;
553: PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
554: if (rz < 0.0) {
555: /* The preconditioner is indefinite. */
556: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
557: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz));
558: break;
559: }
561: /* As far as we know, the preconditioner is positive definite. Compute */
562: /* the residual. Do NOT check for convergence. */
563: cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
565: switch (ksp->normtype) {
566: case KSP_NORM_PRECONDITIONED:
567: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
568: break;
570: case KSP_NORM_UNPRECONDITIONED:
571: PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
572: break;
574: case KSP_NORM_NATURAL:
575: norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M */
576: break;
578: default:
579: norm_r = 0.0;
580: break;
581: }
583: PetscCall(KSPLogResidualHistory(ksp, norm_r));
584: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
585: ksp->rnorm = norm_r;
587: /* Check for breakdown. */
588: beta = rz / rzm1;
589: if (PetscAbsReal(beta) <= 0.0) {
590: /* Conjugate gradients has broken down. */
591: ksp->reason = KSP_DIVERGED_BREAKDOWN;
592: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
593: break;
594: }
596: /* Update p and the norms. */
597: cg->beta[l_size] = beta;
598: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
600: /* Compute the new direction and update the iteration. */
601: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
602: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
603: ++ksp->its;
605: /* Update the Lanczos tridiagonal matrix. */
606: ++l_size;
607: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
608: cg->diag[t_size] = kappa / rz + beta / alpha;
609: ++t_size;
610: }
612: /*
613: We have the Lanczos basis, solve the tridiagonal trust-region problem
614: to obtain the Lanczos direction. We know that the solution lies on
615: the boundary of the trust region. We start by checking that the
616: workspace allocated is large enough.
618: Note that the current version only computes the solution by using the
619: preconditioned direction. Need to think about how to do the
620: unpreconditioned direction calculation.
621: */
623: if (t_size > cg->alloced) {
624: if (cg->alloced) {
625: PetscCall(PetscFree(cg->rwork));
626: PetscCall(PetscFree(cg->iwork));
627: cg->alloced += cg->init_alloc;
628: } else {
629: cg->alloced = cg->init_alloc;
630: }
632: while (t_size > cg->alloced) cg->alloced += cg->init_alloc;
634: cg->alloced = PetscMin(cg->alloced, t_size);
635: PetscCall(PetscMalloc2(10 * cg->alloced, &cg->rwork, 5 * cg->alloced, &cg->iwork));
636: }
638: /* Set up the required vectors. */
639: t_soln = cg->rwork + 0 * t_size; /* Solution */
640: t_diag = cg->rwork + 1 * t_size; /* Diagonal of T */
641: t_offd = cg->rwork + 2 * t_size; /* Off-diagonal of T */
642: e_valu = cg->rwork + 3 * t_size; /* Eigenvalues of T */
643: e_vect = cg->rwork + 4 * t_size; /* Eigenvector of T */
644: e_rwrk = cg->rwork + 5 * t_size; /* Eigen workspace */
646: e_iblk = cg->iwork + 0 * t_size; /* Eigen blocks */
647: e_splt = cg->iwork + 1 * t_size; /* Eigen splits */
648: e_iwrk = cg->iwork + 2 * t_size; /* Eigen workspace */
650: /* Compute the minimum eigenvalue of T. */
651: vl = 0.0;
652: vu = 0.0;
653: il = 1;
654: iu = 1;
656: 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));
658: if ((0 != info) || (1 != e_valus)) {
659: /* Calculation of the minimum eigenvalue failed. Return the */
660: /* Steihaug-Toint direction. */
661: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n"));
662: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: cg->e_min = e_valu[0];
668: /* Compute the initial value of lambda to make (T + lambda I) positive */
669: /* definite. */
670: pert = cg->init_pert;
671: if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
673: while (1) {
674: for (i = 0; i < t_size; ++i) {
675: t_diag[i] = cg->diag[i] + cg->lambda;
676: t_offd[i] = cg->offd[i];
677: }
679: PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
681: if (0 == info) break;
683: pert += pert;
684: cg->lambda = cg->lambda * (1.0 + pert) + pert;
685: }
687: /* Compute the initial step and its norm. */
688: nrhs = 1;
689: nldb = t_size;
691: t_soln[0] = -cg->norm_r[0];
692: for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
694: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
696: if (0 != info) {
697: /* Calculation of the initial step failed; return the Steihaug-Toint */
698: /* direction. */
699: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
700: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
701: PetscFunctionReturn(PETSC_SUCCESS);
702: }
704: norm_t = 0.;
705: for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
706: norm_t = PetscSqrtReal(norm_t);
708: /* Determine the case we are in. */
709: if (norm_t <= cg->radius) {
710: /* The step is within the trust region; check if we are in the hard case */
711: /* and need to move to the boundary by following a direction of negative */
712: /* curvature. */
713: if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
714: /* This is the hard case; compute the eigenvector associated with the */
715: /* minimum eigenvalue and move along this direction to the boundary. */
716: 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));
718: if (0 != info) {
719: /* Calculation of the minimum eigenvalue failed. Return the */
720: /* Steihaug-Toint direction. */
721: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n"));
722: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
723: PetscFunctionReturn(PETSC_SUCCESS);
724: }
726: coef1 = 0.0;
727: coef2 = 0.0;
728: coef3 = -cg->radius * cg->radius;
729: for (i = 0; i < t_size; ++i) {
730: coef1 += e_vect[i] * e_vect[i];
731: coef2 += e_vect[i] * t_soln[i];
732: coef3 += t_soln[i] * t_soln[i];
733: }
735: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
736: root1 = (-coef2 + coef3) / coef1;
737: root2 = (-coef2 - coef3) / coef1;
739: /* Compute objective value for (t_soln + root1 * e_vect) */
740: for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
742: obj1 = 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) 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]);
744: obj1 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
746: /* Compute objective value for (t_soln + root2 * e_vect) */
747: for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
749: obj2 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
750: 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]);
751: obj2 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
753: /* Choose the point with the best objective function value. */
754: if (obj1 <= obj2) {
755: for (i = 0; i < t_size; ++i) t_soln[i] += root1 * e_vect[i];
756: } else {
757: for (i = 0; i < t_size; ++i) t_soln[i] += root2 * e_vect[i];
758: }
759: } else {
760: /* The matrix is positive definite or there was no room to move; the */
761: /* solution is already contained in t_soln. */
762: }
763: } else {
764: /* The step is outside the trust-region. Compute the correct value for */
765: /* lambda by performing Newton's method. */
767: for (i = 0; i < max_newton_its; ++i) {
768: /* Check for convergence. */
769: if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
771: /* Compute the update. */
772: PetscCall(PetscArraycpy(e_rwrk, t_soln, t_size));
774: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
776: if (0 != info) {
777: /* Calculation of the step failed; return the Steihaug-Toint */
778: /* direction. */
779: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
780: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
781: PetscFunctionReturn(PETSC_SUCCESS);
782: }
784: /* Modify lambda. */
785: norm_w = 0.;
786: for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
788: cg->lambda += (norm_t - cg->radius) / cg->radius * (norm_t * norm_t) / norm_w;
790: /* Factor T + lambda I */
791: for (j = 0; j < t_size; ++j) {
792: t_diag[j] = cg->diag[j] + cg->lambda;
793: t_offd[j] = cg->offd[j];
794: }
796: PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
798: if (0 != info) {
799: /* Calculation of factorization failed; return the Steihaug-Toint */
800: /* direction. */
801: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n"));
802: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
803: PetscFunctionReturn(PETSC_SUCCESS);
804: }
806: /* Compute the new step and its norm. */
807: t_soln[0] = -cg->norm_r[0];
808: for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
810: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
812: if (0 != info) {
813: /* Calculation of the step failed; return the Steihaug-Toint */
814: /* direction. */
815: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
816: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
817: PetscFunctionReturn(PETSC_SUCCESS);
818: }
820: norm_t = 0.;
821: for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
822: norm_t = PetscSqrtReal(norm_t);
823: }
825: /* Check for convergence. */
826: if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
827: /* Newton method failed to converge in iteration limit. */
828: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n"));
829: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
830: PetscFunctionReturn(PETSC_SUCCESS);
831: }
832: }
834: /* Recover the norm of the direction and objective function value. */
835: cg->norm_d = norm_t;
837: cg->o_fcn = t_soln[0] * (0.5 * (cg->diag[0] * t_soln[0] + cg->offd[1] * t_soln[1]) + cg->norm_r[0]);
838: 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]);
839: cg->o_fcn += 0.5 * t_soln[i] * (cg->offd[i] * t_soln[i - 1] + cg->diag[i] * t_soln[i]);
841: /* Recover the direction. */
842: sigma = -1;
844: /* Start conjugate gradient method from the beginning */
845: PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
846: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
848: /* Accumulate Q * s */
849: PetscCall(VecCopy(z, d));
850: PetscCall(VecScale(d, sigma * t_soln[0] / cg->norm_r[0]));
852: /* Compute the first direction. */
853: PetscCall(VecCopy(z, p)); /* p = z */
854: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
855: ++ksp->its;
857: for (i = 0; i < l_size - 1; ++i) {
858: /* Update the residual and direction. */
859: alpha = cg->alpha[i];
860: if (alpha >= 0.0) sigma = -sigma;
862: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
863: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
865: /* Accumulate Q * s */
866: PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));
868: /* Update p. */
869: beta = cg->beta[i];
870: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
871: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
872: ++ksp->its;
873: }
875: /* Update the residual and direction. */
876: alpha = cg->alpha[i];
877: if (alpha >= 0.0) sigma = -sigma;
879: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
880: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
882: /* Accumulate Q * s */
883: PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));
885: /* Set the termination reason. */
886: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
887: PetscFunctionReturn(PETSC_SUCCESS);
888: #endif
889: }
891: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
892: {
893: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
894: PetscInt max_its;
896: PetscFunctionBegin;
897: /* Determine the total maximum number of iterations. */
898: max_its = ksp->max_it + cg->max_lanczos_its + 1;
900: /* Set work vectors needed by conjugate gradient method and allocate */
901: /* workspace for Lanczos matrix. */
902: PetscCall(KSPSetWorkVecs(ksp, 3));
903: if (cg->diag) {
904: PetscCall(PetscArrayzero(cg->diag, max_its));
905: PetscCall(PetscArrayzero(cg->offd, max_its));
906: PetscCall(PetscArrayzero(cg->alpha, max_its));
907: PetscCall(PetscArrayzero(cg->beta, max_its));
908: PetscCall(PetscArrayzero(cg->norm_r, max_its));
909: } else {
910: PetscCall(PetscCalloc5(max_its, &cg->diag, max_its, &cg->offd, max_its, &cg->alpha, max_its, &cg->beta, max_its, &cg->norm_r));
911: }
912: PetscFunctionReturn(PETSC_SUCCESS);
913: }
915: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
916: {
917: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
919: PetscFunctionBegin;
920: PetscCall(PetscFree5(cg->diag, cg->offd, cg->alpha, cg->beta, cg->norm_r));
921: if (cg->alloced) PetscCall(PetscFree2(cg->rwork, cg->iwork));
922: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
923: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
924: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
925: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", NULL));
926: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", NULL));
927: PetscCall(KSPDestroyDefault(ksp));
928: PetscFunctionReturn(PETSC_SUCCESS);
929: }
931: static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
932: {
933: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
935: PetscFunctionBegin;
936: cg->radius = radius;
937: PetscFunctionReturn(PETSC_SUCCESS);
938: }
940: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
941: {
942: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
944: PetscFunctionBegin;
945: *norm_d = cg->norm_d;
946: PetscFunctionReturn(PETSC_SUCCESS);
947: }
949: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
950: {
951: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
953: PetscFunctionBegin;
954: *o_fcn = cg->o_fcn;
955: PetscFunctionReturn(PETSC_SUCCESS);
956: }
958: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
959: {
960: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
962: PetscFunctionBegin;
963: *e_min = cg->e_min;
964: PetscFunctionReturn(PETSC_SUCCESS);
965: }
967: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
968: {
969: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
971: PetscFunctionBegin;
972: *lambda = cg->lambda;
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: static PetscErrorCode KSPCGSetFromOptions_GLTR(KSP ksp, PetscOptionItems *PetscOptionsObject)
977: {
978: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
980: PetscFunctionBegin;
981: PetscOptionsHeadBegin(PetscOptionsObject, "KSP GLTR options");
983: PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));
985: PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));
987: PetscCall(PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL));
988: PetscCall(PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL));
989: PetscCall(PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL));
991: PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL));
992: PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL));
994: PetscOptionsHeadEnd();
995: PetscFunctionReturn(PETSC_SUCCESS);
996: }
998: /*MC
999: KSPGLTR - Code to run conjugate gradient method subject to a constraint on the solution norm, used within trust region methods {cite}`gould1999solving`
1001: Options Database Key:
1002: . -ksp_cg_radius <r> - Trust Region Radius
1004: Level: developer
1006: Notes:
1007: Uses preconditioned conjugate gradient to compute an approximate minimizer of the quadratic function
1009: $$
1010: q(s) = g^T * s + .5 * s^T * H * s
1011: $$
1013: subject to the trust region constraint
1015: $$
1016: || s || \le delta,
1017: $$
1019: where
1020: .vb
1021: delta is the trust region radius,
1022: g is the gradient vector,
1023: H is the Hessian approximation,
1024: M is the positive definite preconditioner matrix.
1025: .ve
1027: `KSPConvergedReason` may have the additional values
1028: + `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
1029: - `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step.
1031: The operator and the preconditioner supplied must be symmetric and positive definite.
1033: This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
1035: .seealso: [](ch_ksp), `KSPQCG`, `KSPNASH`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`, `KSPGLTRGetMinEig()`, `KSPGLTRGetLambda()`, `KSPCG`
1036: M*/
1038: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1039: {
1040: KSPCG_GLTR *cg;
1042: PetscFunctionBegin;
1043: PetscCall(PetscNew(&cg));
1044: cg->radius = 0.0;
1045: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1047: cg->init_pert = 1.0e-8;
1048: cg->eigen_tol = 1.0e-10;
1049: cg->newton_tol = 1.0e-6;
1051: cg->alloced = 0;
1052: cg->init_alloc = 1024;
1054: cg->max_lanczos_its = 20;
1055: cg->max_newton_its = 10;
1057: ksp->data = (void *)cg;
1058: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
1059: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
1060: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
1061: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
1062: PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
1064: /* Sets the functions that are associated with this data structure */
1065: /* (in C++ this is the same as defining virtual functions). */
1067: ksp->ops->setup = KSPCGSetUp_GLTR;
1068: ksp->ops->solve = KSPCGSolve_GLTR;
1069: ksp->ops->destroy = KSPCGDestroy_GLTR;
1070: ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1071: ksp->ops->buildsolution = KSPBuildSolutionDefault;
1072: ksp->ops->buildresidual = KSPBuildResidualDefault;
1073: ksp->ops->view = NULL;
1075: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_GLTR));
1076: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_GLTR));
1077: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_GLTR));
1078: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", KSPGLTRGetMinEig_GLTR));
1079: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", KSPGLTRGetLambda_GLTR));
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }