Actual source code: gltr.c
2: #include <../src/ksp/ksp/impls/cg/gltr/gltrimpl.h>
3: #include <petscblaslapack.h>
5: #define GLTR_PRECONDITIONED_DIRECTION 0
6: #define GLTR_UNPRECONDITIONED_DIRECTION 1
7: #define GLTR_DIRECTION_TYPES 2
9: static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};
11: /*@
12: KSPGLTRGetMinEig - Get minimum eigenvalue computed by `KSPGLTR`
14: Collective
16: Input Parameter:
17: . ksp - the iterative context
19: Output Parameter:
20: . e_min - the minimum eigenvalue
22: Level: advanced
24: .seealso: [](ch_ksp), `KSPGLTR`, `KSPGLTRGetLambda()`
25: @*/
26: PetscErrorCode KSPGLTRGetMinEig(KSP ksp, PetscReal *e_min)
27: {
28: PetscFunctionBegin;
30: PetscUseMethod(ksp, "KSPGLTRGetMinEig_C", (KSP, PetscReal *), (ksp, e_min));
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
34: /*@
35: KSPGLTRGetLambda - Get the multiplier on the trust-region constraint when using `KSPGLTR`
36: t
37: Not Collective
39: Input Parameter:
40: . ksp - the iterative context
42: Output Parameter:
43: . lambda - the multiplier
45: Level: advanced
47: .seealso: [](ch_ksp), `KSPGLTR`, `KSPGLTRGetMinEig()`
48: @*/
49: PetscErrorCode KSPGLTRGetLambda(KSP ksp, PetscReal *lambda)
50: {
51: PetscFunctionBegin;
53: PetscUseMethod(ksp, "KSPGLTRGetLambda_C", (KSP, PetscReal *), (ksp, lambda));
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: static PetscErrorCode KSPCGSolve_GLTR(KSP ksp)
58: {
59: #if defined(PETSC_USE_COMPLEX)
60: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "GLTR is not available for complex systems");
61: #else
62: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
63: PetscReal *t_soln, *t_diag, *t_offd, *e_valu, *e_vect, *e_rwrk;
64: PetscBLASInt *e_iblk, *e_splt, *e_iwrk;
66: Mat Qmat, Mmat;
67: Vec r, z, p, d;
68: PC pc;
70: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
71: PetscReal alpha, beta, kappa, rz, rzm1;
72: PetscReal rr, r2, piv, step;
73: PetscReal vl, vu;
74: PetscReal coef1, coef2, coef3, root1, root2, obj1, obj2;
75: PetscReal norm_t, norm_w, pert;
77: PetscInt i, j, max_cg_its, max_lanczos_its, max_newton_its, sigma;
78: PetscBLASInt t_size = 0, l_size = 0, il, iu, info;
79: PetscBLASInt nrhs, nldb;
81: PetscBLASInt e_valus = 0, e_splts;
82: PetscBool diagonalscale;
84: PetscFunctionBegin;
85: /***************************************************************************/
86: /* Check the arguments and parameters. */
87: /***************************************************************************/
89: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
90: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
91: PetscCheck(cg->radius >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
93: /***************************************************************************/
94: /* Get the workspace vectors and initialize variables */
95: /***************************************************************************/
97: r2 = cg->radius * cg->radius;
98: r = ksp->work[0];
99: z = ksp->work[1];
100: p = ksp->work[2];
101: d = ksp->vec_sol;
102: pc = ksp->pc;
104: PetscCall(PCGetOperators(pc, &Qmat, &Mmat));
106: PetscCall(VecGetSize(d, &max_cg_its));
107: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
108: max_lanczos_its = cg->max_lanczos_its;
109: max_newton_its = cg->max_newton_its;
110: ksp->its = 0;
112: /***************************************************************************/
113: /* Initialize objective function direction, and minimum eigenvalue. */
114: /***************************************************************************/
116: cg->o_fcn = 0.0;
118: PetscCall(VecSet(d, 0.0)); /* d = 0 */
119: cg->norm_d = 0.0;
121: cg->e_min = 0.0;
122: cg->lambda = 0.0;
124: /***************************************************************************/
125: /* The first phase of GLTR performs a standard conjugate gradient method, */
126: /* but stores the values required for the Lanczos matrix. We switch to */
127: /* the Lanczos when the conjugate gradient method breaks down. Check the */
128: /* right-hand side for numerical problems. The check for not-a-number and */
129: /* infinite values need be performed only once. */
130: /***************************************************************************/
132: PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
133: PetscCall(VecDot(r, r, &rr)); /* rr = r^T r */
134: KSPCheckDot(ksp, rr);
136: /***************************************************************************/
137: /* Check the preconditioner for numerical problems and for positive */
138: /* definiteness. The check for not-a-number and infinite values need be */
139: /* performed only once. */
140: /***************************************************************************/
142: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
143: PetscCall(VecDot(r, z, &rz)); /* rz = r^T inv(M) r */
144: if (PetscIsInfOrNanScalar(rz)) {
145: /*************************************************************************/
146: /* The preconditioner contains not-a-number or an infinite value. */
147: /* Return the gradient direction intersected with the trust region. */
148: /*************************************************************************/
150: ksp->reason = KSP_DIVERGED_NANORINF;
151: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz));
153: if (cg->radius) {
154: if (r2 >= rr) {
155: alpha = 1.0;
156: cg->norm_d = PetscSqrtReal(rr);
157: } else {
158: alpha = PetscSqrtReal(r2 / rr);
159: cg->norm_d = cg->radius;
160: }
162: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
164: /***********************************************************************/
165: /* Compute objective function. */
166: /***********************************************************************/
168: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
169: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
170: PetscCall(VecDot(d, z, &cg->o_fcn));
171: cg->o_fcn = -cg->o_fcn;
172: ++ksp->its;
173: }
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: if (rz < 0.0) {
178: /*************************************************************************/
179: /* The preconditioner is indefinite. Because this is the first */
180: /* and we do not have a direction yet, we use the gradient step. Note */
181: /* that we cannot use the preconditioned norm when computing the step */
182: /* because the matrix is indefinite. */
183: /*************************************************************************/
185: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
186: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz));
188: if (cg->radius) {
189: if (r2 >= rr) {
190: alpha = 1.0;
191: cg->norm_d = PetscSqrtReal(rr);
192: } else {
193: alpha = PetscSqrtReal(r2 / rr);
194: cg->norm_d = cg->radius;
195: }
197: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
199: /***********************************************************************/
200: /* Compute objective function. */
201: /***********************************************************************/
203: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
204: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
205: PetscCall(VecDot(d, z, &cg->o_fcn));
206: cg->o_fcn = -cg->o_fcn;
207: ++ksp->its;
208: }
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: /***************************************************************************/
213: /* As far as we know, the preconditioner is positive semidefinite. */
214: /* Compute and log the residual. Check convergence because this */
215: /* initializes things, but do not terminate until at least one conjugate */
216: /* gradient iteration has been performed. */
217: /***************************************************************************/
219: cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
221: switch (ksp->normtype) {
222: case KSP_NORM_PRECONDITIONED:
223: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
224: break;
226: case KSP_NORM_UNPRECONDITIONED:
227: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
228: break;
230: case KSP_NORM_NATURAL:
231: norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
232: break;
234: default:
235: norm_r = 0.0;
236: break;
237: }
239: PetscCall(KSPLogResidualHistory(ksp, norm_r));
240: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
241: ksp->rnorm = norm_r;
243: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
245: /***************************************************************************/
246: /* Compute the first direction and update the iteration. */
247: /***************************************************************************/
249: PetscCall(VecCopy(z, p)); /* p = z */
250: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
251: ++ksp->its;
253: /***************************************************************************/
254: /* Check the matrix for numerical problems. */
255: /***************************************************************************/
257: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
258: if (PetscIsInfOrNanScalar(kappa)) {
259: /*************************************************************************/
260: /* The matrix produced not-a-number or an infinite value. In this case, */
261: /* we must stop and use the gradient direction. This condition need */
262: /* only be checked once. */
263: /*************************************************************************/
265: ksp->reason = KSP_DIVERGED_NANORINF;
266: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa));
268: if (cg->radius) {
269: if (r2 >= rr) {
270: alpha = 1.0;
271: cg->norm_d = PetscSqrtReal(rr);
272: } else {
273: alpha = PetscSqrtReal(r2 / rr);
274: cg->norm_d = cg->radius;
275: }
277: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
279: /***********************************************************************/
280: /* Compute objective function. */
281: /***********************************************************************/
283: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
284: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
285: PetscCall(VecDot(d, z, &cg->o_fcn));
286: cg->o_fcn = -cg->o_fcn;
287: ++ksp->its;
288: }
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: /***************************************************************************/
293: /* Initialize variables for calculating the norm of the direction and for */
294: /* the Lanczos tridiagonal matrix. Note that we track the diagonal value */
295: /* of the Cholesky factorization of the Lanczos matrix in order to */
296: /* determine when negative curvature is encountered. */
297: /***************************************************************************/
299: dMp = 0.0;
300: norm_d = 0.0;
301: switch (cg->dtype) {
302: case GLTR_PRECONDITIONED_DIRECTION:
303: norm_p = rz;
304: break;
306: default:
307: PetscCall(VecDot(p, p, &norm_p));
308: break;
309: }
311: cg->diag[t_size] = kappa / rz;
312: cg->offd[t_size] = 0.0;
313: ++t_size;
315: piv = 1.0;
317: /***************************************************************************/
318: /* Check for breakdown of the conjugate gradient method; this occurs when */
319: /* kappa is zero. */
320: /***************************************************************************/
322: if (PetscAbsReal(kappa) <= 0.0) {
323: /*************************************************************************/
324: /* The curvature is zero. In this case, we must stop and use follow */
325: /* the direction of negative curvature since the Lanczos matrix is zero. */
326: /*************************************************************************/
328: ksp->reason = KSP_DIVERGED_BREAKDOWN;
329: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: kappa=%g\n", (double)kappa));
331: if (cg->radius && norm_p > 0.0) {
332: /***********************************************************************/
333: /* Follow direction of negative curvature to the boundary of the */
334: /* trust region. */
335: /***********************************************************************/
337: step = PetscSqrtReal(r2 / norm_p);
338: cg->norm_d = cg->radius;
340: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
342: /***********************************************************************/
343: /* Update objective function. */
344: /***********************************************************************/
346: cg->o_fcn += step * (0.5 * step * kappa - rz);
347: } else if (cg->radius) {
348: /***********************************************************************/
349: /* The norm of the preconditioned direction is zero; use the gradient */
350: /* step. */
351: /***********************************************************************/
353: if (r2 >= rr) {
354: alpha = 1.0;
355: cg->norm_d = PetscSqrtReal(rr);
356: } else {
357: alpha = PetscSqrtReal(r2 / rr);
358: cg->norm_d = cg->radius;
359: }
361: PetscCall(VecAXPY(d, alpha, r)); /* d = d + alpha r */
363: /***********************************************************************/
364: /* Compute objective function. */
365: /***********************************************************************/
367: PetscCall(KSP_MatMult(ksp, Qmat, d, z));
368: PetscCall(VecAYPX(z, -0.5, ksp->vec_rhs));
369: PetscCall(VecDot(d, z, &cg->o_fcn));
370: cg->o_fcn = -cg->o_fcn;
371: ++ksp->its;
372: }
373: PetscFunctionReturn(PETSC_SUCCESS);
374: }
376: /***************************************************************************/
377: /* Begin the first part of the GLTR algorithm which runs the conjugate */
378: /* gradient method until either the problem is solved, we encounter the */
379: /* boundary of the trust region, or the conjugate gradient method breaks */
380: /* down. */
381: /***************************************************************************/
383: while (1) {
384: /*************************************************************************/
385: /* Know that kappa is nonzero, because we have not broken down, so we */
386: /* can compute the steplength. */
387: /*************************************************************************/
389: alpha = rz / kappa;
390: cg->alpha[l_size] = alpha;
392: /*************************************************************************/
393: /* Compute the diagonal value of the Cholesky factorization of the */
394: /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
395: /* This indicates a direction of negative curvature. */
396: /*************************************************************************/
398: piv = cg->diag[l_size] - cg->offd[l_size] * cg->offd[l_size] / piv;
399: if (piv <= 0.0) {
400: /***********************************************************************/
401: /* In this case, the matrix is indefinite and we have encountered */
402: /* a direction of negative curvature. Follow the direction to the */
403: /* boundary of the trust region. */
404: /***********************************************************************/
406: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
407: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa));
409: if (cg->radius && norm_p > 0.0) {
410: /*********************************************************************/
411: /* Follow direction of negative curvature to boundary. */
412: /*********************************************************************/
414: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
415: cg->norm_d = cg->radius;
417: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
419: /*********************************************************************/
420: /* Update objective function. */
421: /*********************************************************************/
423: cg->o_fcn += step * (0.5 * step * kappa - rz);
424: } else if (cg->radius) {
425: /*********************************************************************/
426: /* The norm of the direction is zero; there is nothing to follow. */
427: /*********************************************************************/
428: }
429: break;
430: }
432: /*************************************************************************/
433: /* Compute the steplength and check for intersection with the trust */
434: /* region. */
435: /*************************************************************************/
437: norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
438: if (cg->radius && norm_dp1 >= r2) {
439: /***********************************************************************/
440: /* In this case, the matrix is positive definite as far as we know. */
441: /* However, the full step goes beyond the trust region. */
442: /***********************************************************************/
444: ksp->reason = KSP_CONVERGED_STEP_LENGTH;
445: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius));
447: if (norm_p > 0.0) {
448: /*********************************************************************/
449: /* Follow the direction to the boundary of the trust region. */
450: /*********************************************************************/
452: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
453: cg->norm_d = cg->radius;
455: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
457: /*********************************************************************/
458: /* Update objective function. */
459: /*********************************************************************/
461: cg->o_fcn += step * (0.5 * step * kappa - rz);
462: } else {
463: /*********************************************************************/
464: /* The norm of the direction is zero; there is nothing to follow. */
465: /*********************************************************************/
466: }
467: break;
468: }
470: /*************************************************************************/
471: /* Now we can update the direction and residual. */
472: /*************************************************************************/
474: PetscCall(VecAXPY(d, alpha, p)); /* d = d + alpha p */
475: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
476: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
478: switch (cg->dtype) {
479: case GLTR_PRECONDITIONED_DIRECTION:
480: norm_d = norm_dp1;
481: break;
483: default:
484: PetscCall(VecDot(d, d, &norm_d));
485: break;
486: }
487: cg->norm_d = PetscSqrtReal(norm_d);
489: /*************************************************************************/
490: /* Update objective function. */
491: /*************************************************************************/
493: cg->o_fcn -= 0.5 * alpha * rz;
495: /*************************************************************************/
496: /* Check that the preconditioner appears positive semidefinite. */
497: /*************************************************************************/
499: rzm1 = rz;
500: PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
501: if (rz < 0.0) {
502: /***********************************************************************/
503: /* The preconditioner is indefinite. */
504: /***********************************************************************/
506: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
507: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz));
508: break;
509: }
511: /*************************************************************************/
512: /* As far as we know, the preconditioner is positive semidefinite. */
513: /* Compute the residual and check for convergence. */
514: /*************************************************************************/
516: cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
518: switch (ksp->normtype) {
519: case KSP_NORM_PRECONDITIONED:
520: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
521: break;
523: case KSP_NORM_UNPRECONDITIONED:
524: PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
525: break;
527: case KSP_NORM_NATURAL:
528: norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M */
529: break;
531: default:
532: norm_r = 0.0;
533: break;
534: }
536: PetscCall(KSPLogResidualHistory(ksp, norm_r));
537: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
538: ksp->rnorm = norm_r;
540: PetscCall((*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP));
541: if (ksp->reason) {
542: /***********************************************************************/
543: /* The method has converged. */
544: /***********************************************************************/
546: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius));
547: break;
548: }
550: /*************************************************************************/
551: /* We have not converged yet. Check for breakdown. */
552: /*************************************************************************/
554: beta = rz / rzm1;
555: if (PetscAbsReal(beta) <= 0.0) {
556: /***********************************************************************/
557: /* Conjugate gradients has broken down. */
558: /***********************************************************************/
560: ksp->reason = KSP_DIVERGED_BREAKDOWN;
561: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
562: break;
563: }
565: /*************************************************************************/
566: /* Check iteration limit. */
567: /*************************************************************************/
569: if (ksp->its >= max_cg_its) {
570: ksp->reason = KSP_DIVERGED_ITS;
571: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: iterlim: its=%" PetscInt_FMT "\n", ksp->its));
572: break;
573: }
575: /*************************************************************************/
576: /* Update p and the norms. */
577: /*************************************************************************/
579: cg->beta[l_size] = beta;
580: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
582: switch (cg->dtype) {
583: case GLTR_PRECONDITIONED_DIRECTION:
584: dMp = beta * (dMp + alpha * norm_p);
585: norm_p = beta * (rzm1 + beta * norm_p);
586: break;
588: default:
589: PetscCall(VecDot(d, p, &dMp));
590: PetscCall(VecDot(p, p, &norm_p));
591: break;
592: }
594: /*************************************************************************/
595: /* Compute the new direction and update the iteration. */
596: /*************************************************************************/
598: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
599: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
600: ++ksp->its;
602: /*************************************************************************/
603: /* Update the Lanczos tridiagonal matrix. */
604: /*************************************************************************/
606: ++l_size;
607: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
608: cg->diag[t_size] = kappa / rz + beta / alpha;
609: ++t_size;
611: /*************************************************************************/
612: /* Check for breakdown of the conjugate gradient method; this occurs */
613: /* when kappa is zero. */
614: /*************************************************************************/
616: if (PetscAbsReal(kappa) <= 0.0) {
617: /***********************************************************************/
618: /* The method breaks down; move along the direction as if the matrix */
619: /* were indefinite. */
620: /***********************************************************************/
622: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
623: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa));
625: if (cg->radius && norm_p > 0.0) {
626: /*********************************************************************/
627: /* Follow direction to boundary. */
628: /*********************************************************************/
630: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
631: cg->norm_d = cg->radius;
633: PetscCall(VecAXPY(d, step, p)); /* d = d + step p */
635: /*********************************************************************/
636: /* Update objective function. */
637: /*********************************************************************/
639: cg->o_fcn += step * (0.5 * step * kappa - rz);
640: } else if (cg->radius) {
641: /*********************************************************************/
642: /* The norm of the direction is zero; there is nothing to follow. */
643: /*********************************************************************/
644: }
645: break;
646: }
647: }
649: /***************************************************************************/
650: /* Check to see if we need to continue with the Lanczos method. */
651: /***************************************************************************/
653: if (!cg->radius) {
654: /*************************************************************************/
655: /* There is no radius. Therefore, we cannot move along the boundary. */
656: /*************************************************************************/
657: PetscFunctionReturn(PETSC_SUCCESS);
658: }
660: if (KSP_CONVERGED_NEG_CURVE != ksp->reason) {
661: /*************************************************************************/
662: /* The method either converged to an interior point, hit the boundary of */
663: /* the trust-region without encountering a direction of negative */
664: /* curvature or the iteration limit was reached. */
665: /*************************************************************************/
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: /***************************************************************************/
670: /* Switch to constructing the Lanczos basis by way of the conjugate */
671: /* directions. */
672: /***************************************************************************/
674: for (i = 0; i < max_lanczos_its; ++i) {
675: /*************************************************************************/
676: /* Check for breakdown of the conjugate gradient method; this occurs */
677: /* when kappa is zero. */
678: /*************************************************************************/
680: if (PetscAbsReal(kappa) <= 0.0) {
681: ksp->reason = KSP_DIVERGED_BREAKDOWN;
682: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa));
683: break;
684: }
686: /*************************************************************************/
687: /* Update the direction and residual. */
688: /*************************************************************************/
690: alpha = rz / kappa;
691: cg->alpha[l_size] = alpha;
693: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
694: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
696: /*************************************************************************/
697: /* Check that the preconditioner appears positive semidefinite. */
698: /*************************************************************************/
700: rzm1 = rz;
701: PetscCall(VecDot(r, z, &rz)); /* rz = r^T z */
702: if (rz < 0.0) {
703: /***********************************************************************/
704: /* The preconditioner is indefinite. */
705: /***********************************************************************/
707: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
708: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz));
709: break;
710: }
712: /*************************************************************************/
713: /* As far as we know, the preconditioner is positive definite. Compute */
714: /* the residual. Do NOT check for convergence. */
715: /*************************************************************************/
717: cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
719: switch (ksp->normtype) {
720: case KSP_NORM_PRECONDITIONED:
721: PetscCall(VecNorm(z, NORM_2, &norm_r)); /* norm_r = |z| */
722: break;
724: case KSP_NORM_UNPRECONDITIONED:
725: PetscCall(VecNorm(r, NORM_2, &norm_r)); /* norm_r = |r| */
726: break;
728: case KSP_NORM_NATURAL:
729: norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M */
730: break;
732: default:
733: norm_r = 0.0;
734: break;
735: }
737: PetscCall(KSPLogResidualHistory(ksp, norm_r));
738: PetscCall(KSPMonitor(ksp, ksp->its, norm_r));
739: ksp->rnorm = norm_r;
741: /*************************************************************************/
742: /* Check for breakdown. */
743: /*************************************************************************/
745: beta = rz / rzm1;
746: if (PetscAbsReal(beta) <= 0.0) {
747: /***********************************************************************/
748: /* Conjugate gradients has broken down. */
749: /***********************************************************************/
751: ksp->reason = KSP_DIVERGED_BREAKDOWN;
752: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta));
753: break;
754: }
756: /*************************************************************************/
757: /* Update p and the norms. */
758: /*************************************************************************/
760: cg->beta[l_size] = beta;
761: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
763: /*************************************************************************/
764: /* Compute the new direction and update the iteration. */
765: /*************************************************************************/
767: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
768: PetscCall(VecDot(p, z, &kappa)); /* kappa = p^T Q p */
769: ++ksp->its;
771: /*************************************************************************/
772: /* Update the Lanczos tridiagonal matrix. */
773: /*************************************************************************/
775: ++l_size;
776: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
777: cg->diag[t_size] = kappa / rz + beta / alpha;
778: ++t_size;
779: }
781: /***************************************************************************/
782: /* We have the Lanczos basis, solve the tridiagonal trust-region problem */
783: /* to obtain the Lanczos direction. We know that the solution lies on */
784: /* the boundary of the trust region. We start by checking that the */
785: /* workspace allocated is large enough. */
786: /***************************************************************************/
787: /* Note that the current version only computes the solution by using the */
788: /* preconditioned direction. Need to think about how to do the */
789: /* unpreconditioned direction calculation. */
790: /***************************************************************************/
792: if (t_size > cg->alloced) {
793: if (cg->alloced) {
794: PetscCall(PetscFree(cg->rwork));
795: PetscCall(PetscFree(cg->iwork));
796: cg->alloced += cg->init_alloc;
797: } else {
798: cg->alloced = cg->init_alloc;
799: }
801: while (t_size > cg->alloced) cg->alloced += cg->init_alloc;
803: cg->alloced = PetscMin(cg->alloced, t_size);
804: PetscCall(PetscMalloc2(10 * cg->alloced, &cg->rwork, 5 * cg->alloced, &cg->iwork));
805: }
807: /***************************************************************************/
808: /* Set up the required vectors. */
809: /***************************************************************************/
811: t_soln = cg->rwork + 0 * t_size; /* Solution */
812: t_diag = cg->rwork + 1 * t_size; /* Diagonal of T */
813: t_offd = cg->rwork + 2 * t_size; /* Off-diagonal of T */
814: e_valu = cg->rwork + 3 * t_size; /* Eigenvalues of T */
815: e_vect = cg->rwork + 4 * t_size; /* Eigenvector of T */
816: e_rwrk = cg->rwork + 5 * t_size; /* Eigen workspace */
818: e_iblk = cg->iwork + 0 * t_size; /* Eigen blocks */
819: e_splt = cg->iwork + 1 * t_size; /* Eigen splits */
820: e_iwrk = cg->iwork + 2 * t_size; /* Eigen workspace */
822: /***************************************************************************/
823: /* Compute the minimum eigenvalue of T. */
824: /***************************************************************************/
826: vl = 0.0;
827: vu = 0.0;
828: il = 1;
829: iu = 1;
831: 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));
833: if ((0 != info) || (1 != e_valus)) {
834: /*************************************************************************/
835: /* Calculation of the minimum eigenvalue failed. Return the */
836: /* Steihaug-Toint direction. */
837: /*************************************************************************/
839: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n"));
840: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
841: PetscFunctionReturn(PETSC_SUCCESS);
842: }
844: cg->e_min = e_valu[0];
846: /***************************************************************************/
847: /* Compute the initial value of lambda to make (T + lambda I) positive */
848: /* definite. */
849: /***************************************************************************/
851: pert = cg->init_pert;
852: if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
854: while (1) {
855: for (i = 0; i < t_size; ++i) {
856: t_diag[i] = cg->diag[i] + cg->lambda;
857: t_offd[i] = cg->offd[i];
858: }
860: PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
862: if (0 == info) break;
864: pert += pert;
865: cg->lambda = cg->lambda * (1.0 + pert) + pert;
866: }
868: /***************************************************************************/
869: /* Compute the initial step and its norm. */
870: /***************************************************************************/
872: nrhs = 1;
873: nldb = t_size;
875: t_soln[0] = -cg->norm_r[0];
876: for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
878: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
880: if (0 != info) {
881: /*************************************************************************/
882: /* Calculation of the initial step failed; return the Steihaug-Toint */
883: /* direction. */
884: /*************************************************************************/
886: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
887: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
888: PetscFunctionReturn(PETSC_SUCCESS);
889: }
891: norm_t = 0.;
892: for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
893: norm_t = PetscSqrtReal(norm_t);
895: /***************************************************************************/
896: /* Determine the case we are in. */
897: /***************************************************************************/
899: if (norm_t <= cg->radius) {
900: /*************************************************************************/
901: /* The step is within the trust region; check if we are in the hard case */
902: /* and need to move to the boundary by following a direction of negative */
903: /* curvature. */
904: /*************************************************************************/
906: if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
907: /***********************************************************************/
908: /* This is the hard case; compute the eigenvector associated with the */
909: /* minimum eigenvalue and move along this direction to the boundary. */
910: /***********************************************************************/
912: 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));
914: if (0 != info) {
915: /*********************************************************************/
916: /* Calculation of the minimum eigenvalue failed. Return the */
917: /* Steihaug-Toint direction. */
918: /*********************************************************************/
920: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n"));
921: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
922: PetscFunctionReturn(PETSC_SUCCESS);
923: }
925: coef1 = 0.0;
926: coef2 = 0.0;
927: coef3 = -cg->radius * cg->radius;
928: for (i = 0; i < t_size; ++i) {
929: coef1 += e_vect[i] * e_vect[i];
930: coef2 += e_vect[i] * t_soln[i];
931: coef3 += t_soln[i] * t_soln[i];
932: }
934: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
935: root1 = (-coef2 + coef3) / coef1;
936: root2 = (-coef2 - coef3) / coef1;
938: /***********************************************************************/
939: /* Compute objective value for (t_soln + root1 * e_vect) */
940: /***********************************************************************/
942: for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
944: obj1 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
945: 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]);
946: obj1 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
948: /***********************************************************************/
949: /* Compute objective value for (t_soln + root2 * e_vect) */
950: /***********************************************************************/
952: for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
954: obj2 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
955: 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]);
956: obj2 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
958: /***********************************************************************/
959: /* Choose the point with the best objective function value. */
960: /***********************************************************************/
962: if (obj1 <= obj2) {
963: for (i = 0; i < t_size; ++i) t_soln[i] += root1 * e_vect[i];
964: } else {
965: for (i = 0; i < t_size; ++i) t_soln[i] += root2 * e_vect[i];
966: }
967: } else {
968: /***********************************************************************/
969: /* The matrix is positive definite or there was no room to move; the */
970: /* solution is already contained in t_soln. */
971: /***********************************************************************/
972: }
973: } else {
974: /*************************************************************************/
975: /* The step is outside the trust-region. Compute the correct value for */
976: /* lambda by performing Newton's method. */
977: /*************************************************************************/
979: for (i = 0; i < max_newton_its; ++i) {
980: /***********************************************************************/
981: /* Check for convergence. */
982: /***********************************************************************/
984: if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
986: /***********************************************************************/
987: /* Compute the update. */
988: /***********************************************************************/
990: PetscCall(PetscArraycpy(e_rwrk, t_soln, t_size));
992: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
994: if (0 != info) {
995: /*********************************************************************/
996: /* Calculation of the step failed; return the Steihaug-Toint */
997: /* direction. */
998: /*********************************************************************/
1000: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
1001: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1002: PetscFunctionReturn(PETSC_SUCCESS);
1003: }
1005: /***********************************************************************/
1006: /* Modify lambda. */
1007: /***********************************************************************/
1009: norm_w = 0.;
1010: for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
1012: cg->lambda += (norm_t - cg->radius) / cg->radius * (norm_t * norm_t) / norm_w;
1014: /***********************************************************************/
1015: /* Factor T + lambda I */
1016: /***********************************************************************/
1018: for (j = 0; j < t_size; ++j) {
1019: t_diag[j] = cg->diag[j] + cg->lambda;
1020: t_offd[j] = cg->offd[j];
1021: }
1023: PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1025: if (0 != info) {
1026: /*********************************************************************/
1027: /* Calculation of factorization failed; return the Steihaug-Toint */
1028: /* direction. */
1029: /*********************************************************************/
1031: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n"));
1032: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1033: PetscFunctionReturn(PETSC_SUCCESS);
1034: }
1036: /***********************************************************************/
1037: /* Compute the new step and its norm. */
1038: /***********************************************************************/
1040: t_soln[0] = -cg->norm_r[0];
1041: for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
1043: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1045: if (0 != info) {
1046: /*********************************************************************/
1047: /* Calculation of the step failed; return the Steihaug-Toint */
1048: /* direction. */
1049: /*********************************************************************/
1051: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n"));
1052: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1053: PetscFunctionReturn(PETSC_SUCCESS);
1054: }
1056: norm_t = 0.;
1057: for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1058: norm_t = PetscSqrtReal(norm_t);
1059: }
1061: /*************************************************************************/
1062: /* Check for convergence. */
1063: /*************************************************************************/
1065: if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1066: /***********************************************************************/
1067: /* Newton method failed to converge in iteration limit. */
1068: /***********************************************************************/
1070: PetscCall(PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n"));
1071: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1072: PetscFunctionReturn(PETSC_SUCCESS);
1073: }
1074: }
1076: /***************************************************************************/
1077: /* Recover the norm of the direction and objective function value. */
1078: /***************************************************************************/
1080: cg->norm_d = norm_t;
1082: cg->o_fcn = t_soln[0] * (0.5 * (cg->diag[0] * t_soln[0] + cg->offd[1] * t_soln[1]) + cg->norm_r[0]);
1083: 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]);
1084: cg->o_fcn += 0.5 * t_soln[i] * (cg->offd[i] * t_soln[i - 1] + cg->diag[i] * t_soln[i]);
1086: /***************************************************************************/
1087: /* Recover the direction. */
1088: /***************************************************************************/
1090: sigma = -1;
1092: /***************************************************************************/
1093: /* Start conjugate gradient method from the beginning */
1094: /***************************************************************************/
1096: PetscCall(VecCopy(ksp->vec_rhs, r)); /* r = -grad */
1097: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
1099: /***************************************************************************/
1100: /* Accumulate Q * s */
1101: /***************************************************************************/
1103: PetscCall(VecCopy(z, d));
1104: PetscCall(VecScale(d, sigma * t_soln[0] / cg->norm_r[0]));
1106: /***************************************************************************/
1107: /* Compute the first direction. */
1108: /***************************************************************************/
1110: PetscCall(VecCopy(z, p)); /* p = z */
1111: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
1112: ++ksp->its;
1114: for (i = 0; i < l_size - 1; ++i) {
1115: /*************************************************************************/
1116: /* Update the residual and direction. */
1117: /*************************************************************************/
1119: alpha = cg->alpha[i];
1120: if (alpha >= 0.0) sigma = -sigma;
1122: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
1123: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
1125: /*************************************************************************/
1126: /* Accumulate Q * s */
1127: /*************************************************************************/
1129: PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));
1131: /*************************************************************************/
1132: /* Update p. */
1133: /*************************************************************************/
1135: beta = cg->beta[i];
1136: PetscCall(VecAYPX(p, beta, z)); /* p = z + beta p */
1137: PetscCall(KSP_MatMult(ksp, Qmat, p, z)); /* z = Q * p */
1138: ++ksp->its;
1139: }
1141: /***************************************************************************/
1142: /* Update the residual and direction. */
1143: /***************************************************************************/
1145: alpha = cg->alpha[i];
1146: if (alpha >= 0.0) sigma = -sigma;
1148: PetscCall(VecAXPY(r, -alpha, z)); /* r = r - alpha Q p */
1149: PetscCall(KSP_PCApply(ksp, r, z)); /* z = inv(M) r */
1151: /***************************************************************************/
1152: /* Accumulate Q * s */
1153: /***************************************************************************/
1155: PetscCall(VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z));
1157: /***************************************************************************/
1158: /* Set the termination reason. */
1159: /***************************************************************************/
1161: ksp->reason = ksp->converged_neg_curve ? KSP_CONVERGED_NEG_CURVE : KSP_DIVERGED_INDEFINITE_MAT;
1162: PetscFunctionReturn(PETSC_SUCCESS);
1163: #endif
1164: }
1166: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
1167: {
1168: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1169: PetscInt max_its;
1171: PetscFunctionBegin;
1172: /***************************************************************************/
1173: /* Determine the total maximum number of iterations. */
1174: /***************************************************************************/
1176: max_its = ksp->max_it + cg->max_lanczos_its + 1;
1178: /***************************************************************************/
1179: /* Set work vectors needed by conjugate gradient method and allocate */
1180: /* workspace for Lanczos matrix. */
1181: /***************************************************************************/
1183: PetscCall(KSPSetWorkVecs(ksp, 3));
1184: if (cg->diag) {
1185: PetscCall(PetscArrayzero(cg->diag, max_its));
1186: PetscCall(PetscArrayzero(cg->offd, max_its));
1187: PetscCall(PetscArrayzero(cg->alpha, max_its));
1188: PetscCall(PetscArrayzero(cg->beta, max_its));
1189: PetscCall(PetscArrayzero(cg->norm_r, max_its));
1190: } else {
1191: PetscCall(PetscCalloc5(max_its, &cg->diag, max_its, &cg->offd, max_its, &cg->alpha, max_its, &cg->beta, max_its, &cg->norm_r));
1192: }
1193: PetscFunctionReturn(PETSC_SUCCESS);
1194: }
1196: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
1197: {
1198: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1200: PetscFunctionBegin;
1201: /***************************************************************************/
1202: /* Free memory allocated for the data. */
1203: /***************************************************************************/
1205: PetscCall(PetscFree5(cg->diag, cg->offd, cg->alpha, cg->beta, cg->norm_r));
1206: if (cg->alloced) PetscCall(PetscFree2(cg->rwork, cg->iwork));
1208: /***************************************************************************/
1209: /* Clear composed functions */
1210: /***************************************************************************/
1212: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL));
1213: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL));
1214: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL));
1215: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", NULL));
1216: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", NULL));
1218: /***************************************************************************/
1219: /* Destroy KSP object. */
1220: /***************************************************************************/
1221: PetscCall(KSPDestroyDefault(ksp));
1222: PetscFunctionReturn(PETSC_SUCCESS);
1223: }
1225: static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
1226: {
1227: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1229: PetscFunctionBegin;
1230: cg->radius = radius;
1231: PetscFunctionReturn(PETSC_SUCCESS);
1232: }
1234: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1235: {
1236: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1238: PetscFunctionBegin;
1239: *norm_d = cg->norm_d;
1240: PetscFunctionReturn(PETSC_SUCCESS);
1241: }
1243: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1244: {
1245: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1247: PetscFunctionBegin;
1248: *o_fcn = cg->o_fcn;
1249: PetscFunctionReturn(PETSC_SUCCESS);
1250: }
1252: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1253: {
1254: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1256: PetscFunctionBegin;
1257: *e_min = cg->e_min;
1258: PetscFunctionReturn(PETSC_SUCCESS);
1259: }
1261: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1262: {
1263: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1265: PetscFunctionBegin;
1266: *lambda = cg->lambda;
1267: PetscFunctionReturn(PETSC_SUCCESS);
1268: }
1270: static PetscErrorCode KSPCGSetFromOptions_GLTR(KSP ksp, PetscOptionItems *PetscOptionsObject)
1271: {
1272: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1274: PetscFunctionBegin;
1275: PetscOptionsHeadBegin(PetscOptionsObject, "KSP GLTR options");
1277: PetscCall(PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL));
1279: PetscCall(PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL));
1281: PetscCall(PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL));
1282: PetscCall(PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL));
1283: PetscCall(PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL));
1285: PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL));
1286: PetscCall(PetscOptionsInt("-ksp_cg_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL));
1288: PetscOptionsHeadEnd();
1289: PetscFunctionReturn(PETSC_SUCCESS);
1290: }
1292: /*MC
1293: KSPGLTR - Code to run conjugate gradient method subject to a constraint
1294: on the solution norm.
1296: Options Database Keys:
1297: . -ksp_cg_radius <r> - Trust Region Radius
1299: Level: developer
1301: Notes:
1302: Uses preconditioned conjugate gradient to compute
1303: an approximate minimizer of the quadratic function
1305: q(s) = g^T * s + .5 * s^T * H * s
1307: subject to the trust region constraint
1309: || s || <= delta,
1311: where
1313: delta is the trust region radius,
1314: g is the gradient vector,
1315: H is the Hessian approximation,
1316: M is the positive definite preconditioner matrix.
1318: `KSPConvergedReason` may have the additional values
1319: .vb
1320: KSP_CONVERGED_NEG_CURVE if convergence is reached along a negative curvature direction,
1321: KSP_CONVERGED_STEP_LENGTH if convergence is reached along a constrained step.
1322: .ve
1324: The operator and the preconditioner supplied must be symmetric and positive definite.
1326: This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
1328: Reference:
1329: . * - Gould, N. and Lucidi, S. and Roma, M. and Toint, P., Solving the Trust-Region Subproblem using the Lanczos Method,
1330: SIAM Journal on Optimization, volume 9, number 2, 1999, 504-525
1332: .seealso: [](ch_ksp), `KSPQCG`, `KSPNASH`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`, `KSPGLTRGetMinEig()`, `KSPGLTRGetLambda()`, `KSPCG`
1333: M*/
1335: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1336: {
1337: KSPCG_GLTR *cg;
1339: PetscFunctionBegin;
1340: PetscCall(PetscNew(&cg));
1341: cg->radius = 0.0;
1342: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1344: cg->init_pert = 1.0e-8;
1345: cg->eigen_tol = 1.0e-10;
1346: cg->newton_tol = 1.0e-6;
1348: cg->alloced = 0;
1349: cg->init_alloc = 1024;
1351: cg->max_lanczos_its = 20;
1352: cg->max_newton_its = 10;
1354: ksp->data = (void *)cg;
1355: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3));
1356: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
1357: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));
1358: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
1359: PetscCall(KSPSetConvergedNegativeCurvature(ksp, PETSC_TRUE));
1361: /***************************************************************************/
1362: /* Sets the functions that are associated with this data structure */
1363: /* (in C++ this is the same as defining virtual functions). */
1364: /***************************************************************************/
1366: ksp->ops->setup = KSPCGSetUp_GLTR;
1367: ksp->ops->solve = KSPCGSolve_GLTR;
1368: ksp->ops->destroy = KSPCGDestroy_GLTR;
1369: ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1370: ksp->ops->buildsolution = KSPBuildSolutionDefault;
1371: ksp->ops->buildresidual = KSPBuildResidualDefault;
1372: ksp->ops->view = NULL;
1374: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_GLTR));
1375: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_GLTR));
1376: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_GLTR));
1377: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", KSPGLTRGetMinEig_GLTR));
1378: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", KSPGLTRGetLambda_GLTR));
1379: PetscFunctionReturn(PETSC_SUCCESS);
1380: }