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: [](chapter_ksp), `KSPGLTR`, `KSPGLTRGetLambda()`
25: @*/
26: PetscErrorCode KSPGLTRGetMinEig(KSP ksp, PetscReal *e_min)
27: {
29: PetscUseMethod(ksp, "KSPGLTRGetMinEig_C", (KSP, PetscReal *), (ksp, e_min));
30: return 0;
31: }
33: /*@
34: KSPGLTRGetLambda - Get the multiplier on the trust-region constraint when using `KSPGLTR`
35: t
36: Not Collective
38: Input Parameter:
39: . ksp - the iterative context
41: Output Parameter:
42: . lambda - the multiplier
44: Level: advanced
46: .seealso: [](chapter_ksp), `KSPGLTR`, `KSPGLTRGetMinEig()`
47: @*/
48: PetscErrorCode KSPGLTRGetLambda(KSP ksp, PetscReal *lambda)
49: {
51: PetscUseMethod(ksp, "KSPGLTRGetLambda_C", (KSP, PetscReal *), (ksp, lambda));
52: return 0;
53: }
55: static PetscErrorCode KSPCGSolve_GLTR(KSP ksp)
56: {
57: #if defined(PETSC_USE_COMPLEX)
58: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "GLTR is not available for complex systems");
59: #else
60: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
61: PetscReal *t_soln, *t_diag, *t_offd, *e_valu, *e_vect, *e_rwrk;
62: PetscBLASInt *e_iblk, *e_splt, *e_iwrk;
64: Mat Qmat, Mmat;
65: Vec r, z, p, d;
66: PC pc;
68: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
69: PetscReal alpha, beta, kappa, rz, rzm1;
70: PetscReal rr, r2, piv, step;
71: PetscReal vl, vu;
72: PetscReal coef1, coef2, coef3, root1, root2, obj1, obj2;
73: PetscReal norm_t, norm_w, pert;
75: PetscInt i, j, max_cg_its, max_lanczos_its, max_newton_its, sigma;
76: PetscBLASInt t_size = 0, l_size = 0, il, iu, info;
77: PetscBLASInt nrhs, nldb;
79: PetscBLASInt e_valus = 0, e_splts;
80: PetscBool diagonalscale;
82: /***************************************************************************/
83: /* Check the arguments and parameters. */
84: /***************************************************************************/
86: PCGetDiagonalScale(ksp->pc, &diagonalscale);
90: /***************************************************************************/
91: /* Get the workspace vectors and initialize variables */
92: /***************************************************************************/
94: r2 = cg->radius * cg->radius;
95: r = ksp->work[0];
96: z = ksp->work[1];
97: p = ksp->work[2];
98: d = ksp->vec_sol;
99: pc = ksp->pc;
101: PCGetOperators(pc, &Qmat, &Mmat);
103: VecGetSize(d, &max_cg_its);
104: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
105: max_lanczos_its = cg->max_lanczos_its;
106: max_newton_its = cg->max_newton_its;
107: ksp->its = 0;
109: /***************************************************************************/
110: /* Initialize objective function direction, and minimum eigenvalue. */
111: /***************************************************************************/
113: cg->o_fcn = 0.0;
115: VecSet(d, 0.0); /* d = 0 */
116: cg->norm_d = 0.0;
118: cg->e_min = 0.0;
119: cg->lambda = 0.0;
121: /***************************************************************************/
122: /* The first phase of GLTR performs a standard conjugate gradient method, */
123: /* but stores the values required for the Lanczos matrix. We switch to */
124: /* the Lanczos when the conjugate gradient method breaks down. Check the */
125: /* right-hand side for numerical problems. The check for not-a-number and */
126: /* infinite values need be performed only once. */
127: /***************************************************************************/
129: VecCopy(ksp->vec_rhs, r); /* r = -grad */
130: VecDot(r, r, &rr); /* rr = r^T r */
131: KSPCheckDot(ksp, rr);
133: /***************************************************************************/
134: /* Check the preconditioner for numerical problems and for positive */
135: /* definiteness. The check for not-a-number and infinite values need be */
136: /* performed only once. */
137: /***************************************************************************/
139: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
140: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
141: if (PetscIsInfOrNanScalar(rz)) {
142: /*************************************************************************/
143: /* The preconditioner contains not-a-number or an infinite value. */
144: /* Return the gradient direction intersected with the trust region. */
145: /*************************************************************************/
147: ksp->reason = KSP_DIVERGED_NANORINF;
148: PetscInfo(ksp, "KSPCGSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz);
150: if (cg->radius) {
151: if (r2 >= rr) {
152: alpha = 1.0;
153: cg->norm_d = PetscSqrtReal(rr);
154: } else {
155: alpha = PetscSqrtReal(r2 / rr);
156: cg->norm_d = cg->radius;
157: }
159: VecAXPY(d, alpha, r); /* d = d + alpha r */
161: /***********************************************************************/
162: /* Compute objective function. */
163: /***********************************************************************/
165: KSP_MatMult(ksp, Qmat, d, z);
166: VecAYPX(z, -0.5, ksp->vec_rhs);
167: VecDot(d, z, &cg->o_fcn);
168: cg->o_fcn = -cg->o_fcn;
169: ++ksp->its;
170: }
171: return 0;
172: }
174: if (rz < 0.0) {
175: /*************************************************************************/
176: /* The preconditioner is indefinite. Because this is the first */
177: /* and we do not have a direction yet, we use the gradient step. Note */
178: /* that we cannot use the preconditioned norm when computing the step */
179: /* because the matrix is indefinite. */
180: /*************************************************************************/
182: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
183: PetscInfo(ksp, "KSPCGSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz);
185: if (cg->radius) {
186: if (r2 >= rr) {
187: alpha = 1.0;
188: cg->norm_d = PetscSqrtReal(rr);
189: } else {
190: alpha = PetscSqrtReal(r2 / rr);
191: cg->norm_d = cg->radius;
192: }
194: VecAXPY(d, alpha, r); /* d = d + alpha r */
196: /***********************************************************************/
197: /* Compute objective function. */
198: /***********************************************************************/
200: KSP_MatMult(ksp, Qmat, d, z);
201: VecAYPX(z, -0.5, ksp->vec_rhs);
202: VecDot(d, z, &cg->o_fcn);
203: cg->o_fcn = -cg->o_fcn;
204: ++ksp->its;
205: }
206: return 0;
207: }
209: /***************************************************************************/
210: /* As far as we know, the preconditioner is positive semidefinite. */
211: /* Compute and log the residual. Check convergence because this */
212: /* initializes things, but do not terminate until at least one conjugate */
213: /* gradient iteration has been performed. */
214: /***************************************************************************/
216: cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
218: switch (ksp->normtype) {
219: case KSP_NORM_PRECONDITIONED:
220: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
221: break;
223: case KSP_NORM_UNPRECONDITIONED:
224: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
225: break;
227: case KSP_NORM_NATURAL:
228: norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
229: break;
231: default:
232: norm_r = 0.0;
233: break;
234: }
236: KSPLogResidualHistory(ksp, norm_r);
237: KSPMonitor(ksp, ksp->its, norm_r);
238: ksp->rnorm = norm_r;
240: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
242: /***************************************************************************/
243: /* Compute the first direction and update the iteration. */
244: /***************************************************************************/
246: VecCopy(z, p); /* p = z */
247: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
248: ++ksp->its;
250: /***************************************************************************/
251: /* Check the matrix for numerical problems. */
252: /***************************************************************************/
254: VecDot(p, z, &kappa); /* kappa = p^T Q p */
255: if (PetscIsInfOrNanScalar(kappa)) {
256: /*************************************************************************/
257: /* The matrix produced not-a-number or an infinite value. In this case, */
258: /* we must stop and use the gradient direction. This condition need */
259: /* only be checked once. */
260: /*************************************************************************/
262: ksp->reason = KSP_DIVERGED_NANORINF;
263: PetscInfo(ksp, "KSPCGSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa);
265: if (cg->radius) {
266: if (r2 >= rr) {
267: alpha = 1.0;
268: cg->norm_d = PetscSqrtReal(rr);
269: } else {
270: alpha = PetscSqrtReal(r2 / rr);
271: cg->norm_d = cg->radius;
272: }
274: VecAXPY(d, alpha, r); /* d = d + alpha r */
276: /***********************************************************************/
277: /* Compute objective function. */
278: /***********************************************************************/
280: KSP_MatMult(ksp, Qmat, d, z);
281: VecAYPX(z, -0.5, ksp->vec_rhs);
282: VecDot(d, z, &cg->o_fcn);
283: cg->o_fcn = -cg->o_fcn;
284: ++ksp->its;
285: }
286: return 0;
287: }
289: /***************************************************************************/
290: /* Initialize variables for calculating the norm of the direction and for */
291: /* the Lanczos tridiagonal matrix. Note that we track the diagonal value */
292: /* of the Cholesky factorization of the Lanczos matrix in order to */
293: /* determine when negative curvature is encountered. */
294: /***************************************************************************/
296: dMp = 0.0;
297: norm_d = 0.0;
298: switch (cg->dtype) {
299: case GLTR_PRECONDITIONED_DIRECTION:
300: norm_p = rz;
301: break;
303: default:
304: VecDot(p, p, &norm_p);
305: break;
306: }
308: cg->diag[t_size] = kappa / rz;
309: cg->offd[t_size] = 0.0;
310: ++t_size;
312: piv = 1.0;
314: /***************************************************************************/
315: /* Check for breakdown of the conjugate gradient method; this occurs when */
316: /* kappa is zero. */
317: /***************************************************************************/
319: if (PetscAbsReal(kappa) <= 0.0) {
320: /*************************************************************************/
321: /* The curvature is zero. In this case, we must stop and use follow */
322: /* the direction of negative curvature since the Lanczos matrix is zero. */
323: /*************************************************************************/
325: ksp->reason = KSP_DIVERGED_BREAKDOWN;
326: PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: kappa=%g\n", (double)kappa);
328: if (cg->radius && norm_p > 0.0) {
329: /***********************************************************************/
330: /* Follow direction of negative curvature to the boundary of the */
331: /* trust region. */
332: /***********************************************************************/
334: step = PetscSqrtReal(r2 / norm_p);
335: cg->norm_d = cg->radius;
337: VecAXPY(d, step, p); /* d = d + step p */
339: /***********************************************************************/
340: /* Update objective function. */
341: /***********************************************************************/
343: cg->o_fcn += step * (0.5 * step * kappa - rz);
344: } else if (cg->radius) {
345: /***********************************************************************/
346: /* The norm of the preconditioned direction is zero; use the gradient */
347: /* step. */
348: /***********************************************************************/
350: if (r2 >= rr) {
351: alpha = 1.0;
352: cg->norm_d = PetscSqrtReal(rr);
353: } else {
354: alpha = PetscSqrtReal(r2 / rr);
355: cg->norm_d = cg->radius;
356: }
358: VecAXPY(d, alpha, r); /* d = d + alpha r */
360: /***********************************************************************/
361: /* Compute objective function. */
362: /***********************************************************************/
364: KSP_MatMult(ksp, Qmat, d, z);
365: VecAYPX(z, -0.5, ksp->vec_rhs);
366: VecDot(d, z, &cg->o_fcn);
367: cg->o_fcn = -cg->o_fcn;
368: ++ksp->its;
369: }
370: return 0;
371: }
373: /***************************************************************************/
374: /* Begin the first part of the GLTR algorithm which runs the conjugate */
375: /* gradient method until either the problem is solved, we encounter the */
376: /* boundary of the trust region, or the conjugate gradient method breaks */
377: /* down. */
378: /***************************************************************************/
380: while (1) {
381: /*************************************************************************/
382: /* Know that kappa is nonzero, because we have not broken down, so we */
383: /* can compute the steplength. */
384: /*************************************************************************/
386: alpha = rz / kappa;
387: cg->alpha[l_size] = alpha;
389: /*************************************************************************/
390: /* Compute the diagonal value of the Cholesky factorization of the */
391: /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
392: /* This indicates a direction of negative curvature. */
393: /*************************************************************************/
395: piv = cg->diag[l_size] - cg->offd[l_size] * cg->offd[l_size] / piv;
396: if (piv <= 0.0) {
397: /***********************************************************************/
398: /* In this case, the matrix is indefinite and we have encountered */
399: /* a direction of negative curvature. Follow the direction to the */
400: /* boundary of the trust region. */
401: /***********************************************************************/
403: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
404: PetscInfo(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa);
406: if (cg->radius && norm_p > 0.0) {
407: /*********************************************************************/
408: /* Follow direction of negative curvature to boundary. */
409: /*********************************************************************/
411: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
412: cg->norm_d = cg->radius;
414: VecAXPY(d, step, p); /* d = d + step p */
416: /*********************************************************************/
417: /* Update objective function. */
418: /*********************************************************************/
420: cg->o_fcn += step * (0.5 * step * kappa - rz);
421: } else if (cg->radius) {
422: /*********************************************************************/
423: /* The norm of the direction is zero; there is nothing to follow. */
424: /*********************************************************************/
425: }
426: break;
427: }
429: /*************************************************************************/
430: /* Compute the steplength and check for intersection with the trust */
431: /* region. */
432: /*************************************************************************/
434: norm_dp1 = norm_d + alpha * (2.0 * dMp + alpha * norm_p);
435: if (cg->radius && norm_dp1 >= r2) {
436: /***********************************************************************/
437: /* In this case, the matrix is positive definite as far as we know. */
438: /* However, the full step goes beyond the trust region. */
439: /***********************************************************************/
441: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
442: PetscInfo(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius);
444: if (norm_p > 0.0) {
445: /*********************************************************************/
446: /* Follow the direction to the boundary of the trust region. */
447: /*********************************************************************/
449: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
450: cg->norm_d = cg->radius;
452: VecAXPY(d, step, p); /* d = d + step p */
454: /*********************************************************************/
455: /* Update objective function. */
456: /*********************************************************************/
458: cg->o_fcn += step * (0.5 * step * kappa - rz);
459: } else {
460: /*********************************************************************/
461: /* The norm of the direction is zero; there is nothing to follow. */
462: /*********************************************************************/
463: }
464: break;
465: }
467: /*************************************************************************/
468: /* Now we can update the direction and residual. */
469: /*************************************************************************/
471: VecAXPY(d, alpha, p); /* d = d + alpha p */
472: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
473: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
475: switch (cg->dtype) {
476: case GLTR_PRECONDITIONED_DIRECTION:
477: norm_d = norm_dp1;
478: break;
480: default:
481: VecDot(d, d, &norm_d);
482: break;
483: }
484: cg->norm_d = PetscSqrtReal(norm_d);
486: /*************************************************************************/
487: /* Update objective function. */
488: /*************************************************************************/
490: cg->o_fcn -= 0.5 * alpha * rz;
492: /*************************************************************************/
493: /* Check that the preconditioner appears positive semidefinite. */
494: /*************************************************************************/
496: rzm1 = rz;
497: VecDot(r, z, &rz); /* rz = r^T z */
498: if (rz < 0.0) {
499: /***********************************************************************/
500: /* The preconditioner is indefinite. */
501: /***********************************************************************/
503: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
504: PetscInfo(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz);
505: break;
506: }
508: /*************************************************************************/
509: /* As far as we know, the preconditioner is positive semidefinite. */
510: /* Compute the residual and check for convergence. */
511: /*************************************************************************/
513: cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
515: switch (ksp->normtype) {
516: case KSP_NORM_PRECONDITIONED:
517: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
518: break;
520: case KSP_NORM_UNPRECONDITIONED:
521: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
522: break;
524: case KSP_NORM_NATURAL:
525: norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M */
526: break;
528: default:
529: norm_r = 0.0;
530: break;
531: }
533: KSPLogResidualHistory(ksp, norm_r);
534: KSPMonitor(ksp, ksp->its, norm_r);
535: ksp->rnorm = norm_r;
537: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
538: if (ksp->reason) {
539: /***********************************************************************/
540: /* The method has converged. */
541: /***********************************************************************/
543: PetscInfo(ksp, "KSPCGSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
544: break;
545: }
547: /*************************************************************************/
548: /* We have not converged yet. Check for breakdown. */
549: /*************************************************************************/
551: beta = rz / rzm1;
552: if (PetscAbsReal(beta) <= 0.0) {
553: /***********************************************************************/
554: /* Conjugate gradients has broken down. */
555: /***********************************************************************/
557: ksp->reason = KSP_DIVERGED_BREAKDOWN;
558: PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta);
559: break;
560: }
562: /*************************************************************************/
563: /* Check iteration limit. */
564: /*************************************************************************/
566: if (ksp->its >= max_cg_its) {
567: ksp->reason = KSP_DIVERGED_ITS;
568: PetscInfo(ksp, "KSPCGSolve_GLTR: iterlim: its=%" PetscInt_FMT "\n", ksp->its);
569: break;
570: }
572: /*************************************************************************/
573: /* Update p and the norms. */
574: /*************************************************************************/
576: cg->beta[l_size] = beta;
577: VecAYPX(p, beta, z); /* p = z + beta p */
579: switch (cg->dtype) {
580: case GLTR_PRECONDITIONED_DIRECTION:
581: dMp = beta * (dMp + alpha * norm_p);
582: norm_p = beta * (rzm1 + beta * norm_p);
583: break;
585: default:
586: VecDot(d, p, &dMp);
587: VecDot(p, p, &norm_p);
588: break;
589: }
591: /*************************************************************************/
592: /* Compute the new direction and update the iteration. */
593: /*************************************************************************/
595: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
596: VecDot(p, z, &kappa); /* kappa = p^T Q p */
597: ++ksp->its;
599: /*************************************************************************/
600: /* Update the Lanczos tridiagonal matrix. */
601: /*************************************************************************/
603: ++l_size;
604: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
605: cg->diag[t_size] = kappa / rz + beta / alpha;
606: ++t_size;
608: /*************************************************************************/
609: /* Check for breakdown of the conjugate gradient method; this occurs */
610: /* when kappa is zero. */
611: /*************************************************************************/
613: if (PetscAbsReal(kappa) <= 0.0) {
614: /***********************************************************************/
615: /* The method breaks down; move along the direction as if the matrix */
616: /* were indefinite. */
617: /***********************************************************************/
619: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
620: PetscInfo(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa);
622: if (cg->radius && norm_p > 0.0) {
623: /*********************************************************************/
624: /* Follow direction to boundary. */
625: /*********************************************************************/
627: step = (PetscSqrtReal(dMp * dMp + norm_p * (r2 - norm_d)) - dMp) / norm_p;
628: cg->norm_d = cg->radius;
630: VecAXPY(d, step, p); /* d = d + step p */
632: /*********************************************************************/
633: /* Update objective function. */
634: /*********************************************************************/
636: cg->o_fcn += step * (0.5 * step * kappa - rz);
637: } else if (cg->radius) {
638: /*********************************************************************/
639: /* The norm of the direction is zero; there is nothing to follow. */
640: /*********************************************************************/
641: }
642: break;
643: }
644: }
646: /***************************************************************************/
647: /* Check to see if we need to continue with the Lanczos method. */
648: /***************************************************************************/
650: if (!cg->radius) {
651: /*************************************************************************/
652: /* There is no radius. Therefore, we cannot move along the boundary. */
653: /*************************************************************************/
654: return 0;
655: }
657: if (KSP_CONVERGED_CG_NEG_CURVE != ksp->reason) {
658: /*************************************************************************/
659: /* The method either converged to an interior point, hit the boundary of */
660: /* the trust-region without encountering a direction of negative */
661: /* curvature or the iteration limit was reached. */
662: /*************************************************************************/
663: return 0;
664: }
666: /***************************************************************************/
667: /* Switch to constructing the Lanczos basis by way of the conjugate */
668: /* directions. */
669: /***************************************************************************/
671: for (i = 0; i < max_lanczos_its; ++i) {
672: /*************************************************************************/
673: /* Check for breakdown of the conjugate gradient method; this occurs */
674: /* when kappa is zero. */
675: /*************************************************************************/
677: if (PetscAbsReal(kappa) <= 0.0) {
678: ksp->reason = KSP_DIVERGED_BREAKDOWN;
679: PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa);
680: break;
681: }
683: /*************************************************************************/
684: /* Update the direction and residual. */
685: /*************************************************************************/
687: alpha = rz / kappa;
688: cg->alpha[l_size] = alpha;
690: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
691: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
693: /*************************************************************************/
694: /* Check that the preconditioner appears positive semidefinite. */
695: /*************************************************************************/
697: rzm1 = rz;
698: VecDot(r, z, &rz); /* rz = r^T z */
699: if (rz < 0.0) {
700: /***********************************************************************/
701: /* The preconditioner is indefinite. */
702: /***********************************************************************/
704: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
705: PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz);
706: break;
707: }
709: /*************************************************************************/
710: /* As far as we know, the preconditioner is positive definite. Compute */
711: /* the residual. Do NOT check for convergence. */
712: /*************************************************************************/
714: cg->norm_r[l_size + 1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
716: switch (ksp->normtype) {
717: case KSP_NORM_PRECONDITIONED:
718: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
719: break;
721: case KSP_NORM_UNPRECONDITIONED:
722: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
723: break;
725: case KSP_NORM_NATURAL:
726: norm_r = cg->norm_r[l_size + 1]; /* norm_r = |r|_M */
727: break;
729: default:
730: norm_r = 0.0;
731: break;
732: }
734: KSPLogResidualHistory(ksp, norm_r);
735: KSPMonitor(ksp, ksp->its, norm_r);
736: ksp->rnorm = norm_r;
738: /*************************************************************************/
739: /* Check for breakdown. */
740: /*************************************************************************/
742: beta = rz / rzm1;
743: if (PetscAbsReal(beta) <= 0.0) {
744: /***********************************************************************/
745: /* Conjugate gradients has broken down. */
746: /***********************************************************************/
748: ksp->reason = KSP_DIVERGED_BREAKDOWN;
749: PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta);
750: break;
751: }
753: /*************************************************************************/
754: /* Update p and the norms. */
755: /*************************************************************************/
757: cg->beta[l_size] = beta;
758: VecAYPX(p, beta, z); /* p = z + beta p */
760: /*************************************************************************/
761: /* Compute the new direction and update the iteration. */
762: /*************************************************************************/
764: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
765: VecDot(p, z, &kappa); /* kappa = p^T Q p */
766: ++ksp->its;
768: /*************************************************************************/
769: /* Update the Lanczos tridiagonal matrix. */
770: /*************************************************************************/
772: ++l_size;
773: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
774: cg->diag[t_size] = kappa / rz + beta / alpha;
775: ++t_size;
776: }
778: /***************************************************************************/
779: /* We have the Lanczos basis, solve the tridiagonal trust-region problem */
780: /* to obtain the Lanczos direction. We know that the solution lies on */
781: /* the boundary of the trust region. We start by checking that the */
782: /* workspace allocated is large enough. */
783: /***************************************************************************/
784: /* Note that the current version only computes the solution by using the */
785: /* preconditioned direction. Need to think about how to do the */
786: /* unpreconditioned direction calculation. */
787: /***************************************************************************/
789: if (t_size > cg->alloced) {
790: if (cg->alloced) {
791: PetscFree(cg->rwork);
792: PetscFree(cg->iwork);
793: cg->alloced += cg->init_alloc;
794: } else {
795: cg->alloced = cg->init_alloc;
796: }
798: while (t_size > cg->alloced) cg->alloced += cg->init_alloc;
800: cg->alloced = PetscMin(cg->alloced, t_size);
801: PetscMalloc2(10 * cg->alloced, &cg->rwork, 5 * cg->alloced, &cg->iwork);
802: }
804: /***************************************************************************/
805: /* Set up the required vectors. */
806: /***************************************************************************/
808: t_soln = cg->rwork + 0 * t_size; /* Solution */
809: t_diag = cg->rwork + 1 * t_size; /* Diagonal of T */
810: t_offd = cg->rwork + 2 * t_size; /* Off-diagonal of T */
811: e_valu = cg->rwork + 3 * t_size; /* Eigenvalues of T */
812: e_vect = cg->rwork + 4 * t_size; /* Eigenvector of T */
813: e_rwrk = cg->rwork + 5 * t_size; /* Eigen workspace */
815: e_iblk = cg->iwork + 0 * t_size; /* Eigen blocks */
816: e_splt = cg->iwork + 1 * t_size; /* Eigen splits */
817: e_iwrk = cg->iwork + 2 * t_size; /* Eigen workspace */
819: /***************************************************************************/
820: /* Compute the minimum eigenvalue of T. */
821: /***************************************************************************/
823: vl = 0.0;
824: vu = 0.0;
825: il = 1;
826: iu = 1;
828: 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));
830: if ((0 != info) || (1 != e_valus)) {
831: /*************************************************************************/
832: /* Calculation of the minimum eigenvalue failed. Return the */
833: /* Steihaug-Toint direction. */
834: /*************************************************************************/
836: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n");
837: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
838: return 0;
839: }
841: cg->e_min = e_valu[0];
843: /***************************************************************************/
844: /* Compute the initial value of lambda to make (T + lambda I) positive */
845: /* definite. */
846: /***************************************************************************/
848: pert = cg->init_pert;
849: if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
851: while (1) {
852: for (i = 0; i < t_size; ++i) {
853: t_diag[i] = cg->diag[i] + cg->lambda;
854: t_offd[i] = cg->offd[i];
855: }
857: PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
859: if (0 == info) break;
861: pert += pert;
862: cg->lambda = cg->lambda * (1.0 + pert) + pert;
863: }
865: /***************************************************************************/
866: /* Compute the initial step and its norm. */
867: /***************************************************************************/
869: nrhs = 1;
870: nldb = t_size;
872: t_soln[0] = -cg->norm_r[0];
873: for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
875: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
877: if (0 != info) {
878: /*************************************************************************/
879: /* Calculation of the initial step failed; return the Steihaug-Toint */
880: /* direction. */
881: /*************************************************************************/
883: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
884: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
885: return 0;
886: }
888: norm_t = 0.;
889: for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
890: norm_t = PetscSqrtReal(norm_t);
892: /***************************************************************************/
893: /* Determine the case we are in. */
894: /***************************************************************************/
896: if (norm_t <= cg->radius) {
897: /*************************************************************************/
898: /* The step is within the trust region; check if we are in the hard case */
899: /* and need to move to the boundary by following a direction of negative */
900: /* curvature. */
901: /*************************************************************************/
903: if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
904: /***********************************************************************/
905: /* This is the hard case; compute the eigenvector associated with the */
906: /* minimum eigenvalue and move along this direction to the boundary. */
907: /***********************************************************************/
909: 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));
911: if (0 != info) {
912: /*********************************************************************/
913: /* Calculation of the minimum eigenvalue failed. Return the */
914: /* Steihaug-Toint direction. */
915: /*********************************************************************/
917: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n");
918: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
919: return 0;
920: }
922: coef1 = 0.0;
923: coef2 = 0.0;
924: coef3 = -cg->radius * cg->radius;
925: for (i = 0; i < t_size; ++i) {
926: coef1 += e_vect[i] * e_vect[i];
927: coef2 += e_vect[i] * t_soln[i];
928: coef3 += t_soln[i] * t_soln[i];
929: }
931: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
932: root1 = (-coef2 + coef3) / coef1;
933: root2 = (-coef2 - coef3) / coef1;
935: /***********************************************************************/
936: /* Compute objective value for (t_soln + root1 * e_vect) */
937: /***********************************************************************/
939: for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
941: obj1 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
942: 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]);
943: obj1 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
945: /***********************************************************************/
946: /* Compute objective value for (t_soln + root2 * e_vect) */
947: /***********************************************************************/
949: for (i = 0; i < t_size; ++i) e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
951: obj2 = e_rwrk[0] * (0.5 * (cg->diag[0] * e_rwrk[0] + cg->offd[1] * e_rwrk[1]) + cg->norm_r[0]);
952: 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]);
953: obj2 += 0.5 * e_rwrk[i] * (cg->offd[i] * e_rwrk[i - 1] + cg->diag[i] * e_rwrk[i]);
955: /***********************************************************************/
956: /* Choose the point with the best objective function value. */
957: /***********************************************************************/
959: if (obj1 <= obj2) {
960: for (i = 0; i < t_size; ++i) t_soln[i] += root1 * e_vect[i];
961: } else {
962: for (i = 0; i < t_size; ++i) t_soln[i] += root2 * e_vect[i];
963: }
964: } else {
965: /***********************************************************************/
966: /* The matrix is positive definite or there was no room to move; the */
967: /* solution is already contained in t_soln. */
968: /***********************************************************************/
969: }
970: } else {
971: /*************************************************************************/
972: /* The step is outside the trust-region. Compute the correct value for */
973: /* lambda by performing Newton's method. */
974: /*************************************************************************/
976: for (i = 0; i < max_newton_its; ++i) {
977: /***********************************************************************/
978: /* Check for convergence. */
979: /***********************************************************************/
981: if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
983: /***********************************************************************/
984: /* Compute the update. */
985: /***********************************************************************/
987: PetscArraycpy(e_rwrk, t_soln, t_size);
989: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
991: if (0 != info) {
992: /*********************************************************************/
993: /* Calculation of the step failed; return the Steihaug-Toint */
994: /* direction. */
995: /*********************************************************************/
997: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
998: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
999: return 0;
1000: }
1002: /***********************************************************************/
1003: /* Modify lambda. */
1004: /***********************************************************************/
1006: norm_w = 0.;
1007: for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
1009: cg->lambda += (norm_t - cg->radius) / cg->radius * (norm_t * norm_t) / norm_w;
1011: /***********************************************************************/
1012: /* Factor T + lambda I */
1013: /***********************************************************************/
1015: for (j = 0; j < t_size; ++j) {
1016: t_diag[j] = cg->diag[j] + cg->lambda;
1017: t_offd[j] = cg->offd[j];
1018: }
1020: PetscCallBLAS("LAPACKpttrf", LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1022: if (0 != info) {
1023: /*********************************************************************/
1024: /* Calculation of factorization failed; return the Steihaug-Toint */
1025: /* direction. */
1026: /*********************************************************************/
1028: PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n");
1029: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1030: return 0;
1031: }
1033: /***********************************************************************/
1034: /* Compute the new step and its norm. */
1035: /***********************************************************************/
1037: t_soln[0] = -cg->norm_r[0];
1038: for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
1040: PetscCallBLAS("LAPACKpttrs", LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1042: if (0 != info) {
1043: /*********************************************************************/
1044: /* Calculation of the step failed; return the Steihaug-Toint */
1045: /* direction. */
1046: /*********************************************************************/
1048: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
1049: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1050: return 0;
1051: }
1053: norm_t = 0.;
1054: for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1055: norm_t = PetscSqrtReal(norm_t);
1056: }
1058: /*************************************************************************/
1059: /* Check for convergence. */
1060: /*************************************************************************/
1062: if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1063: /***********************************************************************/
1064: /* Newton method failed to converge in iteration limit. */
1065: /***********************************************************************/
1067: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n");
1068: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1069: return 0;
1070: }
1071: }
1073: /***************************************************************************/
1074: /* Recover the norm of the direction and objective function value. */
1075: /***************************************************************************/
1077: cg->norm_d = norm_t;
1079: cg->o_fcn = t_soln[0] * (0.5 * (cg->diag[0] * t_soln[0] + cg->offd[1] * t_soln[1]) + cg->norm_r[0]);
1080: 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]);
1081: cg->o_fcn += 0.5 * t_soln[i] * (cg->offd[i] * t_soln[i - 1] + cg->diag[i] * t_soln[i]);
1083: /***************************************************************************/
1084: /* Recover the direction. */
1085: /***************************************************************************/
1087: sigma = -1;
1089: /***************************************************************************/
1090: /* Start conjugate gradient method from the beginning */
1091: /***************************************************************************/
1093: VecCopy(ksp->vec_rhs, r); /* r = -grad */
1094: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1096: /***************************************************************************/
1097: /* Accumulate Q * s */
1098: /***************************************************************************/
1100: VecCopy(z, d);
1101: VecScale(d, sigma * t_soln[0] / cg->norm_r[0]);
1103: /***************************************************************************/
1104: /* Compute the first direction. */
1105: /***************************************************************************/
1107: VecCopy(z, p); /* p = z */
1108: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1109: ++ksp->its;
1111: for (i = 0; i < l_size - 1; ++i) {
1112: /*************************************************************************/
1113: /* Update the residual and direction. */
1114: /*************************************************************************/
1116: alpha = cg->alpha[i];
1117: if (alpha >= 0.0) sigma = -sigma;
1119: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1120: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1122: /*************************************************************************/
1123: /* Accumulate Q * s */
1124: /*************************************************************************/
1126: VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z);
1128: /*************************************************************************/
1129: /* Update p. */
1130: /*************************************************************************/
1132: beta = cg->beta[i];
1133: VecAYPX(p, beta, z); /* p = z + beta p */
1134: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1135: ++ksp->its;
1136: }
1138: /***************************************************************************/
1139: /* Update the residual and direction. */
1140: /***************************************************************************/
1142: alpha = cg->alpha[i];
1143: if (alpha >= 0.0) sigma = -sigma;
1145: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1146: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1148: /***************************************************************************/
1149: /* Accumulate Q * s */
1150: /***************************************************************************/
1152: VecAXPY(d, sigma * t_soln[i + 1] / cg->norm_r[i + 1], z);
1154: /***************************************************************************/
1155: /* Set the termination reason. */
1156: /***************************************************************************/
1158: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1159: return 0;
1160: #endif
1161: }
1163: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
1164: {
1165: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1166: PetscInt max_its;
1168: /***************************************************************************/
1169: /* Determine the total maximum number of iterations. */
1170: /***************************************************************************/
1172: max_its = ksp->max_it + cg->max_lanczos_its + 1;
1174: /***************************************************************************/
1175: /* Set work vectors needed by conjugate gradient method and allocate */
1176: /* workspace for Lanczos matrix. */
1177: /***************************************************************************/
1179: KSPSetWorkVecs(ksp, 3);
1180: if (cg->diag) {
1181: PetscArrayzero(cg->diag, max_its);
1182: PetscArrayzero(cg->offd, max_its);
1183: PetscArrayzero(cg->alpha, max_its);
1184: PetscArrayzero(cg->beta, max_its);
1185: PetscArrayzero(cg->norm_r, max_its);
1186: } else {
1187: PetscCalloc5(max_its, &cg->diag, max_its, &cg->offd, max_its, &cg->alpha, max_its, &cg->beta, max_its, &cg->norm_r);
1188: }
1189: return 0;
1190: }
1192: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
1193: {
1194: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1196: /***************************************************************************/
1197: /* Free memory allocated for the data. */
1198: /***************************************************************************/
1200: PetscFree5(cg->diag, cg->offd, cg->alpha, cg->beta, cg->norm_r);
1201: if (cg->alloced) PetscFree2(cg->rwork, cg->iwork);
1203: /***************************************************************************/
1204: /* Clear composed functions */
1205: /***************************************************************************/
1207: PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", NULL);
1208: PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", NULL);
1209: PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", NULL);
1210: PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", NULL);
1211: PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", NULL);
1213: /***************************************************************************/
1214: /* Destroy KSP object. */
1215: /***************************************************************************/
1216: KSPDestroyDefault(ksp);
1217: return 0;
1218: }
1220: static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
1221: {
1222: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1224: cg->radius = radius;
1225: return 0;
1226: }
1228: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1229: {
1230: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1232: *norm_d = cg->norm_d;
1233: return 0;
1234: }
1236: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1237: {
1238: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1240: *o_fcn = cg->o_fcn;
1241: return 0;
1242: }
1244: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1245: {
1246: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1248: *e_min = cg->e_min;
1249: return 0;
1250: }
1252: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1253: {
1254: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1256: *lambda = cg->lambda;
1257: return 0;
1258: }
1260: static PetscErrorCode KSPCGSetFromOptions_GLTR(KSP ksp, PetscOptionItems *PetscOptionsObject)
1261: {
1262: KSPCG_GLTR *cg = (KSPCG_GLTR *)ksp->data;
1264: PetscOptionsHeadBegin(PetscOptionsObject, "KSP GLTR options");
1266: PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL);
1268: PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
1270: PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1271: PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1272: PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);
1274: PetscOptionsInt("-ksp_cg_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL);
1275: PetscOptionsInt("-ksp_cg_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL);
1277: PetscOptionsHeadEnd();
1278: return 0;
1279: }
1281: /*MC
1282: KSPGLTR - Code to run conjugate gradient method subject to a constraint
1283: on the solution norm.
1285: Options Database Keys:
1286: . -ksp_cg_radius <r> - Trust Region Radius
1288: Level: developer
1290: Notes:
1291: Uses preconditioned conjugate gradient to compute
1292: an approximate minimizer of the quadratic function
1294: q(s) = g^T * s + .5 * s^T * H * s
1296: subject to the trust region constraint
1298: || s || <= delta,
1300: where
1302: delta is the trust region radius,
1303: g is the gradient vector,
1304: H is the Hessian approximation,
1305: M is the positive definite preconditioner matrix.
1307: `KSPConvergedReason` may have the additional values
1308: .vb
1309: KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1310: KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step.
1311: .ve
1313: The operator and the preconditioner supplied must be symmetric and positive definite.
1315: This is rarely used directly, it is used in Trust Region methods for nonlinear equations, `SNESNEWTONTR`
1317: Reference:
1318: . * - Gould, N. and Lucidi, S. and Roma, M. and Toint, P., Solving the Trust-Region Subproblem using the Lanczos Method,
1319: SIAM Journal on Optimization, volume 9, number 2, 1999, 504-525
1321: .seealso: [](chapter_ksp), `KSPQCG`, `KSPNASH`, `KSPSTCG`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPCGSetRadius()`, `KSPCGGetNormD()`, `KSPCGGetObjFcn()`, `KSPGLTRGetMinEig()`, `KSPGLTRGetLambda()`, `KSPCG`
1322: M*/
1324: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1325: {
1326: KSPCG_GLTR *cg;
1328: PetscNew(&cg);
1329: cg->radius = 0.0;
1330: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1332: cg->init_pert = 1.0e-8;
1333: cg->eigen_tol = 1.0e-10;
1334: cg->newton_tol = 1.0e-6;
1336: cg->alloced = 0;
1337: cg->init_alloc = 1024;
1339: cg->max_lanczos_its = 20;
1340: cg->max_newton_its = 10;
1342: ksp->data = (void *)cg;
1343: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 3);
1344: KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
1345: KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2);
1346: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);
1348: /***************************************************************************/
1349: /* Sets the functions that are associated with this data structure */
1350: /* (in C++ this is the same as defining virtual functions). */
1351: /***************************************************************************/
1353: ksp->ops->setup = KSPCGSetUp_GLTR;
1354: ksp->ops->solve = KSPCGSolve_GLTR;
1355: ksp->ops->destroy = KSPCGDestroy_GLTR;
1356: ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1357: ksp->ops->buildsolution = KSPBuildSolutionDefault;
1358: ksp->ops->buildresidual = KSPBuildResidualDefault;
1359: ksp->ops->view = NULL;
1361: PetscObjectComposeFunction((PetscObject)ksp, "KSPCGSetRadius_C", KSPCGSetRadius_GLTR);
1362: PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetNormD_C", KSPCGGetNormD_GLTR);
1363: PetscObjectComposeFunction((PetscObject)ksp, "KSPCGGetObjFcn_C", KSPCGGetObjFcn_GLTR);
1364: PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetMinEig_C", KSPGLTRGetMinEig_GLTR);
1365: PetscObjectComposeFunction((PetscObject)ksp, "KSPGLTRGetLambda_C", KSPGLTRGetLambda_GLTR);
1366: return 0;
1367: }