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: /***************************************************************************/
85: /* Check the arguments and parameters. */
86: /***************************************************************************/
88: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
89: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
90: PetscCheck(cg->radius >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
92: /***************************************************************************/
93: /* Get the workspace vectors and initialize variables */
94: /***************************************************************************/
96: r2 = cg->radius * cg->radius;
97: r = ksp->work[0];
98: z = ksp->work[1];
99: p = ksp->work[2];
100: d = ksp->vec_sol;
101: pc = ksp->pc;
103: PetscCall(PCGetOperators(pc, &Qmat, &Mmat));
105: PetscCall(VecGetSize(d, &max_cg_its));
106: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
107: max_lanczos_its = cg->max_lanczos_its;
108: max_newton_its = cg->max_newton_its;
109: ksp->its = 0;
111: /***************************************************************************/
112: /* Initialize objective function direction, and minimum eigenvalue. */
113: /***************************************************************************/
115: cg->o_fcn = 0.0;
117: PetscCall(VecSet(d, 0.0)); /* d = 0 */
118: cg->norm_d = 0.0;
120: cg->e_min = 0.0;
121: cg->lambda = 0.0;
123: /***************************************************************************/
124: /* The first phase of GLTR performs a standard conjugate gradient method, */
125: /* but stores the values required for the Lanczos matrix. We switch to */
126: /* the Lanczos when the conjugate gradient method breaks down. Check the */
127: /* right-hand side for numerical problems. The check for not-a-number and */
128: /* infinite values need be performed only once. */
129: /***************************************************************************/
131: PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
132: PetscCall(VecDot(r, r, &rr)); /* rr = r^T r */
133: KSPCheckDot(ksp, rr);
135: /***************************************************************************/
136: /* Check the preconditioner for numerical problems and for positive */
137: /* definiteness. The check for not-a-number and infinite values need be */
138: /* performed only once. */
139: /***************************************************************************/
141: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
142: PetscCall(VecDot(r, z, &rz)); /* rz = r^T inv(M) r */
143: if (PetscIsInfOrNanScalar(rz)) {
144: /*************************************************************************/
145: /* The preconditioner contains not-a-number or an infinite value. */
146: /* Return the gradient direction intersected with the trust region. */
147: /*************************************************************************/
149: ksp->reason = KSP_DIVERGED_NANORINF;
150: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz));
152: if (cg->radius) {
153: if (r2 >= rr) {
154: alpha = 1.0;
155: cg->norm_d = PetscSqrtReal(rr);
156: } else {
157: alpha = PetscSqrtReal(r2 / rr);
158: cg->norm_d = cg->radius;
159: }
161: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
163: /***********************************************************************/
164: /* Compute objective function. */
165: /***********************************************************************/
167: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
168: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
169: PetscCall(VecDot(d, z, &cg->o_fcn));
170: cg->o_fcn = -cg->o_fcn;
171: ++ksp->its;
172: }
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: if (rz < 0.0) {
177: /*************************************************************************/
178: /* The preconditioner is indefinite. Because this is the first */
179: /* and we do not have a direction yet, we use the gradient step. Note */
180: /* that we cannot use the preconditioned norm when computing the step */
181: /* because the matrix is indefinite. */
182: /*************************************************************************/
184: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
185: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz));
187: if (cg->radius) {
188: if (r2 >= rr) {
189: alpha = 1.0;
190: cg->norm_d = PetscSqrtReal(rr);
191: } else {
192: alpha = PetscSqrtReal(r2 / rr);
193: cg->norm_d = cg->radius;
194: }
196: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
198: /***********************************************************************/
199: /* Compute objective function. */
200: /***********************************************************************/
202: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
203: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
204: PetscCall(VecDot(d, z, &cg->o_fcn));
205: cg->o_fcn = -cg->o_fcn;
206: ++ksp->its;
207: }
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: /***************************************************************************/
212: /* As far as we know, the preconditioner is positive semidefinite. */
213: /* Compute and log the residual. Check convergence because this */
214: /* initializes things, but do not terminate until at least one conjugate */
215: /* gradient iteration has been performed. */
216: /***************************************************************************/
218: cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
220: switch (ksp->normtype) {
221: case KSP_NORM_PRECONDITIONED:
222: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
223: break;
225: case KSP_NORM_UNPRECONDITIONED:
226: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
227: break;
229: case KSP_NORM_NATURAL:
230: norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
231: break;
233: default:
234: norm_r = 0.0;
235: break;
236: }
238: PetscCall(KSPLogResidualHistory(ksp, norm_r));
239: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
240: ksp->rnorm = norm_r;
242: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
244: /***************************************************************************/
245: /* Compute the first direction and update the iteration. */
246: /***************************************************************************/
248: PetscCall(VecCopy(z, p)); /* p = z */
249: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
250: ++ksp->its;
252: /***************************************************************************/
253: /* Check the matrix for numerical problems. */
254: /***************************************************************************/
256: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
257: if (PetscIsInfOrNanScalar(kappa)) {
258: /*************************************************************************/
259: /* The matrix produced not-a-number or an infinite value. In this case, */
260: /* we must stop and use the gradient direction. This condition need */
261: /* only be checked once. */
262: /*************************************************************************/
264: ksp->reason = KSP_DIVERGED_NANORINF;
265: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa));
267: if (cg->radius) {
268: if (r2 >= rr) {
269: alpha = 1.0;
270: cg->norm_d = PetscSqrtReal(rr);
271: } else {
272: alpha = PetscSqrtReal(r2 / rr);
273: cg->norm_d = cg->radius;
274: }
276: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
278: /***********************************************************************/
279: /* Compute objective function. */
280: /***********************************************************************/
282: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
283: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
284: PetscCall(VecDot(d, z, &cg->o_fcn));
285: cg->o_fcn = -cg->o_fcn;
286: ++ksp->its;
287: }
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: /***************************************************************************/
292: /* Initialize variables for calculating the norm of the direction and for */
293: /* the Lanczos tridiagonal matrix. Note that we track the diagonal value */
294: /* of the Cholesky factorization of the Lanczos matrix in order to */
295: /* determine when negative curvature is encountered. */
296: /***************************************************************************/
298: dMp = 0.0;
299: norm_d = 0.0;
300: switch (cg->dtype) {
301: case GLTR_PRECONDITIONED_DIRECTION:
302: norm_p = rz;
303: break;
305: default:
306: PetscCall(VecDot(p, p, &norm_p));
307: break;
308: }
310: cg->diag[t_size] = kappa / rz;
311: cg->offd[t_size] = 0.0;
312: ++t_size;
314: piv = 1.0;
316: /***************************************************************************/
317: /* Check for breakdown of the conjugate gradient method; this occurs when */
318: /* kappa is zero. */
319: /***************************************************************************/
321: if (PetscAbsReal(kappa) <= 0.0) {
322: /*************************************************************************/
323: /* The curvature is zero. In this case, we must stop and use follow */
324: /* the direction of negative curvature since the Lanczos matrix is zero. */
325: /*************************************************************************/
327: ksp->reason = KSP_DIVERGED_BREAKDOWN;
328: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: kappa=%g\n", (double)kappa));
330: if (cg->radius && norm_p > 0.0) {
331: /***********************************************************************/
332: /* Follow direction of negative curvature to the boundary of the */
333: /* trust region. */
334: /***********************************************************************/
336: step = PetscSqrtReal(r2 / norm_p);
337: cg->norm_d = cg->radius;
339: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
341: /***********************************************************************/
342: /* Update objective function. */
343: /***********************************************************************/
345: cg->o_fcn += step * (0.5 * step * kappa - rz);
346: } else if (cg->radius) {
347: /***********************************************************************/
348: /* The norm of the preconditioned direction is zero; use the gradient */
349: /* step. */
350: /***********************************************************************/
352: if (r2 >= rr) {
353: alpha = 1.0;
354: cg->norm_d = PetscSqrtReal(rr);
355: } else {
356: alpha = PetscSqrtReal(r2 / rr);
357: cg->norm_d = cg->radius;
358: }
360: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
362: /***********************************************************************/
363: /* Compute objective function. */
364: /***********************************************************************/
366: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
367: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
368: PetscCall(VecDot(d, z, &cg->o_fcn));
369: cg->o_fcn = -cg->o_fcn;
370: ++ksp->its;
371: }
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: /***************************************************************************/
376: /* Begin the first part of the GLTR algorithm which runs the conjugate */
377: /* gradient method until either the problem is solved, we encounter the */
378: /* boundary of the trust region, or the conjugate gradient method breaks */
379: /* down. */
380: /***************************************************************************/
382: while (1) {
383: /*************************************************************************/
384: /* Know that kappa is nonzero, because we have not broken down, so we */
385: /* can compute the steplength. */
386: /*************************************************************************/
388: alpha = rz / kappa;
389: cg->alpha[l_size] = alpha;
391: /*************************************************************************/
392: /* Compute the diagonal value of the Cholesky factorization of the */
393: /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
394: /* This indicates a direction of negative curvature. */
395: /*************************************************************************/
397: piv = cg->diag[l_size] - cg->offd[l_size] * cg->offd[l_size] / piv;
398: if (piv <= 0.0) {
399: /***********************************************************************/
400: /* In this case, the matrix is indefinite and we have encountered */
401: /* a direction of negative curvature. Follow the direction to the */
402: /* boundary of the trust region. */
403: /***********************************************************************/
405: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
406: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa));
408: if (cg->radius && norm_p > 0.0) {
409: /*********************************************************************/
410: /* Follow direction of negative curvature to boundary. */
411: /*********************************************************************/
413: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
414: cg->norm_d = cg->radius;
416: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
418: /*********************************************************************/
419: /* Update objective function. */
420: /*********************************************************************/
422: cg->o_fcn += step * (0.5 * step * kappa - rz);
423: } else if (cg->radius) {
424: /*********************************************************************/
425: /* The norm of the direction is zero; there is nothing to follow. */
426: /*********************************************************************/
427: }
428: break;
429: }
431: /*************************************************************************/
432: /* Compute the steplength and check for intersection with the trust */
433: /* region. */
434: /*************************************************************************/
436: norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
437: if (cg->radius && norm_dp1 >= r2) {
438: /***********************************************************************/
439: /* In this case, the matrix is positive definite as far as we know. */
440: /* However, the full step goes beyond the trust region. */
441: /***********************************************************************/
443: ksp->reason = KSP_CONVERGED_STEP_LENGTH;
444: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius));
446: if (norm_p > 0.0) {
447: /*********************************************************************/
448: /* Follow the direction to the boundary of the trust region. */
449: /*********************************************************************/
451: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
452: cg->norm_d = cg->radius;
454: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
456: /*********************************************************************/
457: /* Update objective function. */
458: /*********************************************************************/
460: cg->o_fcn += step * (0.5 * step * kappa - rz);
461: } else {
462: /*********************************************************************/
463: /* The norm of the direction is zero; there is nothing to follow. */
464: /*********************************************************************/
465: }
466: break;
467: }
469: /*************************************************************************/
470: /* Now we can update the direction and residual. */
471: /*************************************************************************/
473: PetscCall(VecAXPY(d, alpha, p)); /* d = d + alpha p */
474: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
475: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
477: switch (cg->dtype) {
478: case GLTR_PRECONDITIONED_DIRECTION:
479: norm_d = norm_dp1;
480: break;
482: default:
483: PetscCall(VecDot(d, d, &norm_d));
484: break;
485: }
486: cg->norm_d = PetscSqrtReal(norm_d);
488: /*************************************************************************/
489: /* Update objective function. */
490: /*************************************************************************/
492: cg->o_fcn -= 0.5 * alpha * rz;
494: /*************************************************************************/
495: /* Check that the preconditioner appears positive semidefinite. */
496: /*************************************************************************/
498: rzm1 = rz;
499: PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
500: if (rz < 0.0) {
501: /***********************************************************************/
502: /* The preconditioner is indefinite. */
503: /***********************************************************************/
505: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
506: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz));
507: break;
508: }
510: /*************************************************************************/
511: /* As far as we know, the preconditioner is positive semidefinite. */
512: /* Compute the residual and check for convergence. */
513: /*************************************************************************/
515: cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
517: switch (ksp->normtype) {
518: case KSP_NORM_PRECONDITIONED:
519: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
520: break;
522: case KSP_NORM_UNPRECONDITIONED:
523: PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
524: break;
526: case KSP_NORM_NATURAL:
527: norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M */
528: break;
530: default:
531: norm_r = 0.0;
532: break;
533: }
535: PetscCall(KSPLogResidualHistory(ksp, norm_r));
536: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
537: ksp->rnorm = norm_r;
539: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
540: if (ksp->reason) {
541: /***********************************************************************/
542: /* The method has converged. */
543: /***********************************************************************/
545: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius));
546: break;
547: }
549: /*************************************************************************/
550: /* We have not converged yet. Check for breakdown. */
551: /*************************************************************************/
553: beta = rz / rzm1;
554: if (PetscAbsReal(beta) <= 0.0) {
555: /***********************************************************************/
556: /* Conjugate gradients has broken down. */
557: /***********************************************************************/
559: ksp->reason = KSP_DIVERGED_BREAKDOWN;
560: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
561: break;
562: }
564: /*************************************************************************/
565: /* Check iteration limit. */
566: /*************************************************************************/
568: if (ksp->its >= max_cg_its) {
569: ksp->reason = KSP_DIVERGED_ITS;
570: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: iterlim: its=%" PetscInt_FMT "\n", ksp->its));
571: break;
572: }
574: /*************************************************************************/
575: /* Update p and the norms. */
576: /*************************************************************************/
578: cg->beta[l_size] = beta;
579: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
581: switch (cg->dtype) {
582: case GLTR_PRECONDITIONED_DIRECTION:
583: dMp = beta * (dMp + alpha * norm_p);
584: norm_p = beta * (rzm1 + beta * norm_p);
585: break;
587: default:
588: PetscCall(VecDot(d, p, &dMp));
589: PetscCall(VecDot(p, p, &norm_p));
590: break;
591: }
593: /*************************************************************************/
594: /* Compute the new direction and update the iteration. */
595: /*************************************************************************/
597: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
598: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
599: ++ksp->its;
601: /*************************************************************************/
602: /* Update the Lanczos tridiagonal matrix. */
603: /*************************************************************************/
605: ++l_size;
606: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
607: cg->diag[t_size] = kappa / rz + beta / alpha;
608: ++t_size;
610: /*************************************************************************/
611: /* Check for breakdown of the conjugate gradient method; this occurs */
612: /* when kappa is zero. */
613: /*************************************************************************/
615: if (PetscAbsReal(kappa) <= 0.0) {
616: /***********************************************************************/
617: /* The method breaks down; move along the direction as if the matrix */
618: /* were indefinite. */
619: /***********************************************************************/
621: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
622: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa));
624: if (cg->radius && norm_p > 0.0) {
625: /*********************************************************************/
626: /* Follow direction to boundary. */
627: /*********************************************************************/
629: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
630: cg->norm_d = cg->radius;
632: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
634: /*********************************************************************/
635: /* Update objective function. */
636: /*********************************************************************/
638: cg->o_fcn += step * (0.5 * step * kappa - rz);
639: } else if (cg->radius) {
640: /*********************************************************************/
641: /* The norm of the direction is zero; there is nothing to follow. */
642: /*********************************************************************/
643: }
644: break;
645: }
646: }
648: /***************************************************************************/
649: /* Check to see if we need to continue with the Lanczos method. */
650: /***************************************************************************/
652: if (!cg->radius) {
653: /*************************************************************************/
654: /* There is no radius. Therefore, we cannot move along the boundary. */
655: /*************************************************************************/
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: if (KSP_CONVERGED_NEG_CURVE != ksp->reason) {
660: /*************************************************************************/
661: /* The method either converged to an interior point, hit the boundary of */
662: /* the trust-region without encountering a direction of negative */
663: /* curvature or the iteration limit was reached. */
664: /*************************************************************************/
665: PetscFunctionReturn(PETSC_SUCCESS);
666: }
668: /***************************************************************************/
669: /* Switch to constructing the Lanczos basis by way of the conjugate */
670: /* directions. */
671: /***************************************************************************/
673: for (i = 0; i < max_lanczos_its; ++i) {
674: /*************************************************************************/
675: /* Check for breakdown of the conjugate gradient method; this occurs */
676: /* when kappa is zero. */
677: /*************************************************************************/
679: if (PetscAbsReal(kappa) <= 0.0) {
680: ksp->reason = KSP_DIVERGED_BREAKDOWN;
681: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa));
682: break;
683: }
685: /*************************************************************************/
686: /* Update the direction and residual. */
687: /*************************************************************************/
689: alpha = rz / kappa;
690: cg->alpha[l_size] = alpha;
692: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
693: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
695: /*************************************************************************/
696: /* Check that the preconditioner appears positive semidefinite. */
697: /*************************************************************************/
699: rzm1 = rz;
700: PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
701: if (rz < 0.0) {
702: /***********************************************************************/
703: /* The preconditioner is indefinite. */
704: /***********************************************************************/
706: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
707: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz));
708: break;
709: }
711: /*************************************************************************/
712: /* As far as we know, the preconditioner is positive definite. Compute */
713: /* the residual. Do NOT check for convergence. */
714: /*************************************************************************/
716: cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
718: switch (ksp->normtype) {
719: case KSP_NORM_PRECONDITIONED:
720: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
721: break;
723: case KSP_NORM_UNPRECONDITIONED:
724: PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
725: break;
727: case KSP_NORM_NATURAL:
728: norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M */
729: break;
731: default:
732: norm_r = 0.0;
733: break;
734: }
736: PetscCall(KSPLogResidualHistory(ksp, norm_r));
737: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
738: ksp->rnorm = norm_r;
740: /*************************************************************************/
741: /* Check for breakdown. */
742: /*************************************************************************/
744: beta = rz / rzm1;
745: if (PetscAbsReal(beta) <= 0.0) {
746: /***********************************************************************/
747: /* Conjugate gradients has broken down. */
748: /***********************************************************************/
750: ksp->reason = KSP_DIVERGED_BREAKDOWN;
751: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
752: break;
753: }
755: /*************************************************************************/
756: /* Update p and the norms. */
757: /*************************************************************************/
759: cg->beta[l_size] = beta;
760: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
762: /*************************************************************************/
763: /* Compute the new direction and update the iteration. */
764: /*************************************************************************/
766: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
767: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
768: ++ksp->its;
770: /*************************************************************************/
771: /* Update the Lanczos tridiagonal matrix. */
772: /*************************************************************************/
774: ++l_size;
775: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
776: cg->diag[t_size] = kappa / rz + beta / alpha;
777: ++t_size;
778: }
780: /***************************************************************************/
781: /* We have the Lanczos basis, solve the tridiagonal trust-region problem */
782: /* to obtain the Lanczos direction. We know that the solution lies on */
783: /* the boundary of the trust region. We start by checking that the */
784: /* workspace allocated is large enough. */
785: /***************************************************************************/
786: /* Note that the current version only computes the solution by using the */
787: /* preconditioned direction. Need to think about how to do the */
788: /* unpreconditioned direction calculation. */
789: /***************************************************************************/
791: if (t_size > cg->alloced) {
792: if (cg->alloced) {
793: PetscCall(PetscFree(cg->rwork));
794: PetscCall(PetscFree(cg->iwork));
795: cg->alloced += cg->init_alloc;
796: } else {
797: cg->alloced = cg->init_alloc;
798: }
800: while (t_size > cg->alloced) cg->alloced += cg->init_alloc;
802: cg->alloced = PetscMin(cg->alloced, t_size);
803: PetscCall(PetscMalloc2(10 * cg->alloced, &cg->rwork, 5 * cg->alloced, &cg->iwork));
804: }
806: /***************************************************************************/
807: /* Set up the required vectors. */
808: /***************************************************************************/
810: t_soln = cg->rwork + 0 * t_size; /* Solution */
811: t_diag = cg->rwork + 1 * t_size; /* Diagonal of T */
812: t_offd = cg->rwork + 2 * t_size; /* Off-diagonal of T */
813: e_valu = cg->rwork + 3 * t_size; /* Eigenvalues of T */
814: e_vect = cg->rwork + 4 * t_size; /* Eigenvector of T */
815: e_rwrk = cg->rwork + 5 * t_size; /* Eigen workspace */
817: e_iblk = cg->iwork + 0 * t_size; /* Eigen blocks */
818: e_splt = cg->iwork + 1 * t_size; /* Eigen splits */
819: e_iwrk = cg->iwork + 2 * t_size; /* Eigen workspace */
821: /***************************************************************************/
822: /* Compute the minimum eigenvalue of T. */
823: /***************************************************************************/
825: vl = 0.0;
826: vu = 0.0;
827: il = 1;
828: iu = 1;
830: 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));
832: if ((0 != info) || (1 != e_valus)) {
833: /*************************************************************************/
834: /* Calculation of the minimum eigenvalue failed. Return the */
835: /* Steihaug-Toint direction. */
836: /*************************************************************************/
838: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n"));
839: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
840: PetscFunctionReturn(PETSC_SUCCESS);
841: }
843: cg->e_min = e_valu[0];
845: /***************************************************************************/
846: /* Compute the initial value of lambda to make (T + lambda I) positive */
847: /* definite. */
848: /***************************************************************************/
850: pert = cg->init_pert;
851: if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
853: while (1) {
854: for (i = 0; i < t_size; ++i) {
855: t_diag[i] = cg->diag[i] + cg->lambda;
856: t_offd[i] = cg->offd[i];
857: }
859: PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
861: if (0 == info) break;
863: pert += pert;
864: cg->lambda = cg->lambda * (1.0 + pert) + pert;
865: }
867: /***************************************************************************/
868: /* Compute the initial step and its norm. */
869: /***************************************************************************/
871: nrhs = 1;
872: nldb = t_size;
874: t_soln[0] = -cg->norm_r[0];
875: for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
877: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
879: if (0 != info) {
880: /*************************************************************************/
881: /* Calculation of the initial step failed; return the Steihaug-Toint */
882: /* direction. */
883: /*************************************************************************/
885: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
886: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
887: PetscFunctionReturn(PETSC_SUCCESS);
888: }
890: norm_t = 0.;
891: for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
892: norm_t = PetscSqrtReal(norm_t);
894: /***************************************************************************/
895: /* Determine the case we are in. */
896: /***************************************************************************/
898: if (norm_t <= cg->radius) {
899: /*************************************************************************/
900: /* The step is within the trust region; check if we are in the hard case */
901: /* and need to move to the boundary by following a direction of negative */
902: /* curvature. */
903: /*************************************************************************/
905: if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
906: /***********************************************************************/
907: /* This is the hard case; compute the eigenvector associated with the */
908: /* minimum eigenvalue and move along this direction to the boundary. */
909: /***********************************************************************/
911: 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));
913: if (0 != info) {
914: /*********************************************************************/
915: /* Calculation of the minimum eigenvalue failed. Return the */
916: /* Steihaug-Toint direction. */
917: /*********************************************************************/
919: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n"));
920: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
921: PetscFunctionReturn(PETSC_SUCCESS);
922: }
924: coef1 = 0.0;
925: coef2 = 0.0;
926: coef3 = -cg->radius * cg->radius;
927: for (i = 0; i < t_size; ++i) {
928: coef1 += e_vect[i] * e_vect[i];
929: coef2 += e_vect[i] * t_soln[i];
930: coef3 += t_soln[i] * t_soln[i];
931: }
933: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
934: root1 = (-coef2 + coef3) / coef1;
935: root2 = (-coef2 - coef3) / coef1;
937: /***********************************************************************/
938: /* Compute objective value for (t_soln + root1 * e_vect) */
939: /***********************************************************************/
941: for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
943: obj1 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
944: 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]);
945: obj1 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
947: /***********************************************************************/
948: /* Compute objective value for (t_soln + root2 * e_vect) */
949: /***********************************************************************/
951: for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
953: obj2 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
954: 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]);
955: obj2 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
957: /***********************************************************************/
958: /* Choose the point with the best objective function value. */
959: /***********************************************************************/
961: if (obj1 <= obj2) {
962: for (i = 0; i < t_size; ++i) t_soln[i] += root1 * e_vect[i];
963: } else {
964: for (i = 0; i < t_size; ++i) t_soln[i] += root2 * e_vect[i];
965: }
966: } else {
967: /***********************************************************************/
968: /* The matrix is positive definite or there was no room to move; the */
969: /* solution is already contained in t_soln. */
970: /***********************************************************************/
971: }
972: } else {
973: /*************************************************************************/
974: /* The step is outside the trust-region. Compute the correct value for */
975: /* lambda by performing Newton's method. */
976: /*************************************************************************/
978: for (i = 0; i < max_newton_its; ++i) {
979: /***********************************************************************/
980: /* Check for convergence. */
981: /***********************************************************************/
983: if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
985: /***********************************************************************/
986: /* Compute the update. */
987: /***********************************************************************/
989: PetscCall(PetscArraycpy(e_rwrk, t_soln, t_size));
991: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
993: if (0 != info) {
994: /*********************************************************************/
995: /* Calculation of the step failed; return the Steihaug-Toint */
996: /* direction. */
997: /*********************************************************************/
999: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
1000: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1001: PetscFunctionReturn(PETSC_SUCCESS);
1002: }
1004: /***********************************************************************/
1005: /* Modify lambda. */
1006: /***********************************************************************/
1008: norm_w = 0.;
1009: for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
1011: cg->lambda += (norm_t - cg->radius) / cg->radius * (norm_t * norm_t) / norm_w;
1013: /***********************************************************************/
1014: /* Factor T + lambda I */
1015: /***********************************************************************/
1017: for (j = 0; j < t_size; ++j) {
1018: t_diag[j] = cg->diag[j] + cg->lambda;
1019: t_offd[j] = cg->offd[j];
1020: }
1022: PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1024: if (0 != info) {
1025: /*********************************************************************/
1026: /* Calculation of factorization failed; return the Steihaug-Toint */
1027: /* direction. */
1028: /*********************************************************************/
1030: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n"));
1031: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1032: PetscFunctionReturn(PETSC_SUCCESS);
1033: }
1035: /***********************************************************************/
1036: /* Compute the new step and its norm. */
1037: /***********************************************************************/
1039: t_soln[0] = -cg->norm_r[0];
1040: for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
1042: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1044: if (0 != info) {
1045: /*********************************************************************/
1046: /* Calculation of the step failed; return the Steihaug-Toint */
1047: /* direction. */
1048: /*********************************************************************/
1050: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
1051: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1052: PetscFunctionReturn(PETSC_SUCCESS);
1053: }
1055: norm_t = 0.;
1056: for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1057: norm_t = PetscSqrtReal(norm_t);
1058: }
1060: /*************************************************************************/
1061: /* Check for convergence. */
1062: /*************************************************************************/
1064: if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1065: /***********************************************************************/
1066: /* Newton method failed to converge in iteration limit. */
1067: /***********************************************************************/
1069: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n"));
1070: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }
1073: }
1075: /***************************************************************************/
1076: /* Recover the norm of the direction and objective function value. */
1077: /***************************************************************************/
1079: cg->norm_d = norm_t;
1081: cg->o_fcn = t_soln[0] * (0.5 * (cg->diag[0] * t_soln[0] + cg->offd[1] * t_soln[1]) + cg->norm_r[0]);
1082: 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]);
1083: cg->o_fcn += 0.5 * t_soln[i] * (cg->offd[i] * t_soln[i - 1] + cg->diag[i] * t_soln[i]);
1085: /***************************************************************************/
1086: /* Recover the direction. */
1087: /***************************************************************************/
1089: sigma = -1;
1091: /***************************************************************************/
1092: /* Start conjugate gradient method from the beginning */
1093: /***************************************************************************/
1095: PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
1096: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
1098: /***************************************************************************/
1099: /* Accumulate Q * s */
1100: /***************************************************************************/
1102: PetscCall(VecCopy(z, d));
1103: PetscCall(VecScale(d, sigma * t_soln[0] / cg->norm_r[0]));
1105: /***************************************************************************/
1106: /* Compute the first direction. */
1107: /***************************************************************************/
1109: PetscCall(VecCopy(z, p)); /* p = z */
1110: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
1111: ++ksp->its;
1113: for (i = 0; i < l_size - 1; ++i) {
1114: /*************************************************************************/
1115: /* Update the residual and direction. */
1116: /*************************************************************************/
1118: alpha = cg->alpha[i];
1119: if (alpha >= 0.0) sigma = -sigma;
1121: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
1122: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
1124: /*************************************************************************/
1125: /* Accumulate Q * s */
1126: /*************************************************************************/
1128: PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));
1130: /*************************************************************************/
1131: /* Update p. */
1132: /*************************************************************************/
1134: beta = cg->beta[i];
1135: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
1136: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
1137: ++ksp->its;
1138: }
1140: /***************************************************************************/
1141: /* Update the residual and direction. */
1142: /***************************************************************************/
1144: alpha = cg->alpha[i];
1145: if (alpha >= 0.0) sigma = -sigma;
1147: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
1148: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
1150: /***************************************************************************/
1151: /* Accumulate Q * s */
1152: /***************************************************************************/
1154: PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));
1156: /***************************************************************************/
1157: /* Set the termination reason. */
1158: /***************************************************************************/
1160: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1161: PetscFunctionReturn(PETSC_SUCCESS);
1162: #endif
1163: }
1165: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
1166: {
1167: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1168: PetscInt max_its;
1170: PetscFunctionBegin;
1171: /***************************************************************************/
1172: /* Determine the total maximum number of iterations. */
1173: /***************************************************************************/
1175: max_its = ksp->max_it + cg->max_lanczos_its + 1;
1177: /***************************************************************************/
1178: /* Set work vectors needed by conjugate gradient method and allocate */
1179: /* workspace for Lanczos matrix. */
1180: /***************************************************************************/
1182: PetscCall(KSPSetWorkVecs(ksp, 3));
1183: if (cg->diag) {
1184: PetscCall(PetscArrayzero(cg->diag, max_its));
1185: PetscCall(PetscArrayzero(cg->offd, max_its));
1186: PetscCall(PetscArrayzero(cg->alpha, max_its));
1187: PetscCall(PetscArrayzero(cg->beta, max_its));
1188: PetscCall(PetscArrayzero(cg->norm_r, max_its));
1189: } else {
1190: PetscCall(PetscCalloc5(max_its, &cg->diag, max_its, &cg->offd, max_its, &cg->alpha, max_its, &cg->beta, max_its, &cg->norm_r));
1191: }
1192: PetscFunctionReturn(PETSC_SUCCESS);
1193: }
1195: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
1196: {
1197: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1199: PetscFunctionBegin;
1200: PetscCall(PetscFree5(cg->diag, cg->offd, cg->alpha, cg->beta, cg->norm_r));
1201: if (cg->alloced) PetscCall(PetscFree2(cg->rwork, cg->iwork));
1202: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
1203: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
1204: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
1205: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", NULL));
1206: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", NULL));
1207: PetscCall(KSPDestroyDefault(ksp));
1208: PetscFunctionReturn(PETSC_SUCCESS);
1209: }
1211: static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
1212: {
1213: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1215: PetscFunctionBegin;
1216: cg->radius = radius;
1217: PetscFunctionReturn(PETSC_SUCCESS);
1218: }
1220: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1221: {
1222: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1224: PetscFunctionBegin;
1225: *norm_d = cg->norm_d;
1226: PetscFunctionReturn(PETSC_SUCCESS);
1227: }
1229: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1230: {
1231: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1233: PetscFunctionBegin;
1234: *o_fcn = cg->o_fcn;
1235: PetscFunctionReturn(PETSC_SUCCESS);
1236: }
1238: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1239: {
1240: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1242: PetscFunctionBegin;
1243: *e_min = cg->e_min;
1244: PetscFunctionReturn(PETSC_SUCCESS);
1245: }
1247: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1248: {
1249: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1251: PetscFunctionBegin;
1252: *lambda = cg->lambda;
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: static PetscErrorCode KSPCGSetFromOptions_GLTR(KSP ksp, PetscOptionItems *PetscOptionsObject)
1257: {
1258: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1260: PetscFunctionBegin;
1261: PetscOptionsHeadBegin(PetscOptionsObject, "KSP GLTR options");
1263: PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));
1265: PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));
1267: PetscCall(PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL));
1268: PetscCall(PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL));
1269: PetscCall(PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL));
1271: PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL));
1272: PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL));
1274: PetscOptionsHeadEnd();
1275: PetscFunctionReturn(PETSC_SUCCESS);
1276: }
1278: /*MC
1279: KSPGLTR - Code to run conjugate gradient method subject to a constraint on the solution norm, used within trust region methods {cite}`gould1999solving`
1281: Options Database Key:
1282: . -ksp_cg_radius <r> - Trust Region Radius
1284: Level: developer
1286: Notes:
1287: Uses preconditioned conjugate gradient to compute an approximate minimizer of the quadratic function
1289: $$
1290: q(s) = g^T * s + .5 * s^T * H * s
1291: $$
1293: subject to the trust region constraint
1295: $$
1296: || s || \le delta,
1297: $$
1299: where
1300: .vb
1301: delta is the trust region radius,
1302: g is the gradient vector,
1303: H is the Hessian approximation,
1304: M is the positive definite preconditioner matrix.
1305: .ve
1307: `KSPConvergedReason` may have the additional values
1308: + `KSP_CONVERGED_NEG_CURVE` - if convergence is reached along a negative curvature direction,
1309: - `KSP_CONVERGED_STEP_LENGTH` - if convergence is reached along a constrained step.
1311: The operator and the preconditioner supplied must be symmetric and positive definite.
1313: This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
1315: .seealso: [](ch_ksp), `KSPQCG`, `KSPNASH`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`, `KSPGLTRGetMinEig()`, `KSPGLTRGetLambda()`, `KSPCG`
1316: M*/
1318: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1319: {
1320: KSPCG_GLTR *cg;
1322: PetscFunctionBegin;
1323: PetscCall(PetscNew(&cg));
1324: cg->radius = 0.0;
1325: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1327: cg->init_pert = 1.0e-8;
1328: cg->eigen_tol = 1.0e-10;
1329: cg->newton_tol = 1.0e-6;
1331: cg->alloced = 0;
1332: cg->init_alloc = 1024;
1334: cg->max_lanczos_its = 20;
1335: cg->max_newton_its = 10;
1337: ksp->data = (void *)cg;
1338: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
1339: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
1340: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
1341: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
1342: PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
1344: /***************************************************************************/
1345: /* Sets the functions that are associated with this data structure */
1346: /* (in C++ this is the same as defining virtual functions). */
1347: /***************************************************************************/
1349: ksp->ops->setup = KSPCGSetUp_GLTR;
1350: ksp->ops->solve = KSPCGSolve_GLTR;
1351: ksp->ops->destroy = KSPCGDestroy_GLTR;
1352: ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1353: ksp->ops->buildsolution = KSPBuildSolutionDefault;
1354: ksp->ops->buildresidual = KSPBuildResidualDefault;
1355: ksp->ops->view = NULL;
1357: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_GLTR));
1358: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_GLTR));
1359: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_GLTR));
1360: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", KSPGLTRGetMinEig_GLTR));
1361: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", KSPGLTRGetLambda_GLTR));
1362: PetscFunctionReturn(PETSC_SUCCESS);
1363: }