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(PetscFree2(cg->rwork, cg->iwork));
626: cg->alloced += cg->init_alloc;
627: } else {
628: cg->alloced = cg->init_alloc;
629: }
631: while (t_size > cg->alloced) cg->alloced += cg->init_alloc;
633: cg->alloced = PetscMin(cg->alloced, t_size);
634: PetscCall(PetscMalloc2(10 * cg->alloced, &cg->rwork, 5 * cg->alloced, &cg->iwork));
635: }
637: /* Set up the required vectors. */
638: t_soln = cg->rwork + 0 * t_size; /* Solution */
639: t_diag = cg->rwork + 1 * t_size; /* Diagonal of T */
640: t_offd = cg->rwork + 2 * t_size; /* Off-diagonal of T */
641: e_valu = cg->rwork + 3 * t_size; /* Eigenvalues of T */
642: e_vect = cg->rwork + 4 * t_size; /* Eigenvector of T */
643: e_rwrk = cg->rwork + 5 * t_size; /* Eigen workspace */
645: e_iblk = cg->iwork + 0 * t_size; /* Eigen blocks */
646: e_splt = cg->iwork + 1 * t_size; /* Eigen splits */
647: e_iwrk = cg->iwork + 2 * t_size; /* Eigen workspace */
649: /* Compute the minimum eigenvalue of T. */
650: vl = 0.0;
651: vu = 0.0;
652: il = 1;
653: iu = 1;
655: 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));
657: if ((0 != info) || (1 != e_valus)) {
658: /* Calculation of the minimum eigenvalue failed. Return the */
659: /* Steihaug-Toint direction. */
660: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n"));
661: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }
665: cg->e_min = e_valu[0];
667: /* Compute the initial value of lambda to make (T + lambda I) positive */
668: /* definite. */
669: pert = cg->init_pert;
670: if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
672: while (1) {
673: for (i = 0; i < t_size; ++i) {
674: t_diag[i] = cg->diag[i] + cg->lambda;
675: t_offd[i] = cg->offd[i];
676: }
678: PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
680: if (0 == info) break;
682: pert += pert;
683: cg->lambda = cg->lambda * (1.0 + pert) + pert;
684: }
686: /* Compute the initial step and its norm. */
687: nrhs = 1;
688: nldb = t_size;
690: t_soln[0] = -cg->norm_r[0];
691: for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
693: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
695: if (0 != info) {
696: /* Calculation of the initial step failed; return the Steihaug-Toint */
697: /* direction. */
698: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
699: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: norm_t = 0.;
704: for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
705: norm_t = PetscSqrtReal(norm_t);
707: /* Determine the case we are in. */
708: if (norm_t <= cg->radius) {
709: /* The step is within the trust region; check if we are in the hard case */
710: /* and need to move to the boundary by following a direction of negative */
711: /* curvature. */
712: if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
713: /* This is the hard case; compute the eigenvector associated with the */
714: /* minimum eigenvalue and move along this direction to the boundary. */
715: 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));
717: if (0 != info) {
718: /* Calculation of the minimum eigenvalue failed. Return the */
719: /* Steihaug-Toint direction. */
720: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n"));
721: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
725: coef1 = 0.0;
726: coef2 = 0.0;
727: coef3 = -cg->radius * cg->radius;
728: for (i = 0; i < t_size; ++i) {
729: coef1 += e_vect[i] * e_vect[i];
730: coef2 += e_vect[i] * t_soln[i];
731: coef3 += t_soln[i] * t_soln[i];
732: }
734: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
735: root1 = (-coef2 + coef3) / coef1;
736: root2 = (-coef2 - coef3) / coef1;
738: /* Compute objective value for (t_soln + root1 * e_vect) */
739: for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
741: obj1 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
742: 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]);
743: obj1 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
745: /* Compute objective value for (t_soln + root2 * e_vect) */
746: for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
748: obj2 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
749: 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]);
750: obj2 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
752: /* Choose the point with the best objective function value. */
753: if (obj1 <= obj2) {
754: for (i = 0; i < t_size; ++i) t_soln[i] += root1 * e_vect[i];
755: } else {
756: for (i = 0; i < t_size; ++i) t_soln[i] += root2 * e_vect[i];
757: }
758: } else {
759: /* The matrix is positive definite or there was no room to move; the */
760: /* solution is already contained in t_soln. */
761: }
762: } else {
763: /* The step is outside the trust-region. Compute the correct value for */
764: /* lambda by performing Newton's method. */
766: for (i = 0; i < max_newton_its; ++i) {
767: /* Check for convergence. */
768: if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
770: /* Compute the update. */
771: PetscCall(PetscArraycpy(e_rwrk, t_soln, t_size));
773: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
775: if (0 != info) {
776: /* Calculation of the step failed; return the Steihaug-Toint */
777: /* direction. */
778: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
779: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
780: PetscFunctionReturn(PETSC_SUCCESS);
781: }
783: /* Modify lambda. */
784: norm_w = 0.;
785: for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
787: cg->lambda += (norm_t - cg->radius) / cg->radius * (norm_t * norm_t) / norm_w;
789: /* Factor T + lambda I */
790: for (j = 0; j < t_size; ++j) {
791: t_diag[j] = cg->diag[j] + cg->lambda;
792: t_offd[j] = cg->offd[j];
793: }
795: PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
797: if (0 != info) {
798: /* Calculation of factorization failed; return the Steihaug-Toint */
799: /* direction. */
800: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n"));
801: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
802: PetscFunctionReturn(PETSC_SUCCESS);
803: }
805: /* Compute the new step and its norm. */
806: t_soln[0] = -cg->norm_r[0];
807: for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
809: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
811: if (0 != info) {
812: /* Calculation of the step failed; return the Steihaug-Toint */
813: /* direction. */
814: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
815: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
816: PetscFunctionReturn(PETSC_SUCCESS);
817: }
819: norm_t = 0.;
820: for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
821: norm_t = PetscSqrtReal(norm_t);
822: }
824: /* Check for convergence. */
825: if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
826: /* Newton method failed to converge in iteration limit. */
827: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n"));
828: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
831: }
833: /* Recover the norm of the direction and objective function value. */
834: cg->norm_d = norm_t;
836: cg->o_fcn = t_soln[0] * (0.5 * (cg->diag[0] * t_soln[0] + cg->offd[1] * t_soln[1]) + cg->norm_r[0]);
837: 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]);
838: cg->o_fcn += 0.5 * t_soln[i] * (cg->offd[i] * t_soln[i - 1] + cg->diag[i] * t_soln[i]);
840: /* Recover the direction. */
841: sigma = -1;
843: /* Start conjugate gradient method from the beginning */
844: PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
845: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
847: /* Accumulate Q * s */
848: PetscCall(VecCopy(z, d));
849: PetscCall(VecScale(d, sigma * t_soln[0] / cg->norm_r[0]));
851: /* Compute the first direction. */
852: PetscCall(VecCopy(z, p)); /* p = z */
853: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
854: ++ksp->its;
856: for (i = 0; i < l_size - 1; ++i) {
857: /* Update the residual and direction. */
858: alpha = cg->alpha[i];
859: if (alpha >= 0.0) sigma = -sigma;
861: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
862: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
864: /* Accumulate Q * s */
865: PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));
867: /* Update p. */
868: beta = cg->beta[i];
869: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
870: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
871: ++ksp->its;
872: }
874: /* Update the residual and direction. */
875: alpha = cg->alpha[i];
876: if (alpha >= 0.0) sigma = -sigma;
878: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
879: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
881: /* Accumulate Q * s */
882: PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));
884: /* Set the termination reason. */
885: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
886: PetscFunctionReturn(PETSC_SUCCESS);
887: #endif
888: }
890: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
891: {
892: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
893: PetscInt max_its;
895: PetscFunctionBegin;
896: /* Determine the total maximum number of iterations. */
897: max_its = ksp->max_it + cg->max_lanczos_its + 1;
899: /* Set work vectors needed by conjugate gradient method and allocate */
900: /* workspace for Lanczos matrix. */
901: PetscCall(KSPSetWorkVecs(ksp, 3));
902: if (cg->diag) {
903: PetscCall(PetscArrayzero(cg->diag, max_its));
904: PetscCall(PetscArrayzero(cg->offd, max_its));
905: PetscCall(PetscArrayzero(cg->alpha, max_its));
906: PetscCall(PetscArrayzero(cg->beta, max_its));
907: PetscCall(PetscArrayzero(cg->norm_r, max_its));
908: } else {
909: PetscCall(PetscCalloc5(max_its, &cg->diag, max_its, &cg->offd, max_its, &cg->alpha, max_its, &cg->beta, max_its, &cg->norm_r));
910: }
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }
914: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
915: {
916: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
918: PetscFunctionBegin;
919: PetscCall(PetscFree5(cg->diag, cg->offd, cg->alpha, cg->beta, cg->norm_r));
920: if (cg->alloced) PetscCall(PetscFree2(cg->rwork, cg->iwork));
921: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
922: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
923: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
924: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", NULL));
925: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", NULL));
926: PetscCall(KSPDestroyDefault(ksp));
927: PetscFunctionReturn(PETSC_SUCCESS);
928: }
930: static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
931: {
932: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
934: PetscFunctionBegin;
935: cg->radius = radius;
936: PetscFunctionReturn(PETSC_SUCCESS);
937: }
939: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
940: {
941: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
943: PetscFunctionBegin;
944: *norm_d = cg->norm_d;
945: PetscFunctionReturn(PETSC_SUCCESS);
946: }
948: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
949: {
950: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
952: PetscFunctionBegin;
953: *o_fcn = cg->o_fcn;
954: PetscFunctionReturn(PETSC_SUCCESS);
955: }
957: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
958: {
959: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
961: PetscFunctionBegin;
962: *e_min = cg->e_min;
963: PetscFunctionReturn(PETSC_SUCCESS);
964: }
966: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
967: {
968: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
970: PetscFunctionBegin;
971: *lambda = cg->lambda;
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: static PetscErrorCode KSPCGSetFromOptions_GLTR(KSP ksp, PetscOptionItems PetscOptionsObject)
976: {
977: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
979: PetscFunctionBegin;
980: PetscOptionsHeadBegin(PetscOptionsObject, "KSP GLTR options");
982: PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));
984: PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));
986: PetscCall(PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL));
987: PetscCall(PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL));
988: PetscCall(PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL));
990: PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL));
991: PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL));
993: PetscOptionsHeadEnd();
994: PetscFunctionReturn(PETSC_SUCCESS);
995: }
997: /*MC
998: KSPGLTR - Code to run conjugate gradient method subject to a constraint on the solution norm, used within trust region methods {cite}`gould1999solving`
1000: Options Database Key:
1001: . -ksp_cg_radius radius - Trust Region Radius
1003: Level: developer
1005: Notes:
1006: Uses preconditioned conjugate gradient to compute an approximate minimizer of the quadratic function
1008: $$
1009: q(s) = g^T * s + .5 * s^T * H * s
1010: $$
1012: subject to the trust region constraint
1014: $$
1015: || s || \le delta,
1016: $$
1018: where
1019: .vb
1020: delta is the trust region radius,
1021: g is the gradient vector,
1022: H is the Hessian approximation,
1023: .ve
1025: `KSPConvergedReason` may have the additional values
1026: + `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
1027: - `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step.
1029: The operator and the preconditioner supplied must be symmetric and positive definite.
1031: This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
1033: .seealso: [](ch_ksp), `KSPQCG`, `KSPNASH`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`, `KSPGLTRGetMinEig()`, `KSPGLTRGetLambda()`, `KSPCG`
1034: M*/
1036: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1037: {
1038: KSPCG_GLTR *cg;
1040: PetscFunctionBegin;
1041: PetscCall(PetscNew(&cg));
1042: cg->radius = 0.0;
1043: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1045: cg->init_pert = 1.0e-8;
1046: cg->eigen_tol = 1.0e-10;
1047: cg->newton_tol = 1.0e-6;
1049: cg->alloced = 0;
1050: cg->init_alloc = 1024;
1052: cg->max_lanczos_its = 20;
1053: cg->max_newton_its = 10;
1055: ksp->data = (void *)cg;
1056: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
1057: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
1058: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
1059: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
1060: PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
1062: /* Sets the functions that are associated with this data structure */
1063: /* (in C++ this is the same as defining virtual functions). */
1065: ksp->ops->setup = KSPCGSetUp_GLTR;
1066: ksp->ops->solve = KSPCGSolve_GLTR;
1067: ksp->ops->destroy = KSPCGDestroy_GLTR;
1068: ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1069: ksp->ops->buildsolution = KSPBuildSolutionDefault;
1070: ksp->ops->buildresidual = KSPBuildResidualDefault;
1071: ksp->ops->view = NULL;
1073: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_GLTR));
1074: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_GLTR));
1075: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_GLTR));
1076: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", KSPGLTRGetMinEig_GLTR));
1077: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", KSPGLTRGetLambda_GLTR));
1078: PetscFunctionReturn(PETSC_SUCCESS);
1079: }