Actual source code: gltr.c
petsc-3.14.6 2021-03-30
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.
14: Collective on ksp
16: Input Parameters:
17: + ksp - the iterative context
18: - e_min - the minimum eigenvalue
20: Level: advanced
22: @*/
23: PetscErrorCode KSPGLTRGetMinEig(KSP ksp, PetscReal *e_min)
24: {
29: PetscUseMethod(ksp,"KSPGLTRGetMinEig_C",(KSP,PetscReal*),(ksp,e_min));
30: return(0);
31: }
33: /*@
34: KSPGLTRGetLambda - Get multiplier on trust-region constraint.
36: Not Collective
38: Input Parameters:
39: + ksp - the iterative context
40: - lambda - the multiplier
42: Level: advanced
44: @*/
45: PetscErrorCode KSPGLTRGetLambda(KSP ksp, PetscReal *lambda)
46: {
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;
65: Mat Qmat, Mmat;
66: Vec r, z, p, d;
67: PC pc;
69: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
70: PetscReal alpha, beta, kappa, rz, rzm1;
71: PetscReal rr, r2, piv, step;
72: PetscReal vl, vu;
73: PetscReal coef1, coef2, coef3, root1, root2, obj1, obj2;
74: PetscReal norm_t, norm_w, pert;
76: PetscInt i, j, max_cg_its, max_lanczos_its, max_newton_its, sigma;
77: PetscBLASInt t_size = 0, l_size = 0, il, iu, info;
78: PetscBLASInt nrhs, nldb;
80: PetscBLASInt e_valus=0, e_splts;
81: PetscBool diagonalscale;
84: /***************************************************************************/
85: /* Check the arguments and parameters. */
86: /***************************************************************************/
88: PCGetDiagonalScale(ksp->pc, &diagonalscale);
89: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
90: if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
92: /***************************************************************************/
93: /* Get the workspace vectors and initialize variables */
94: /***************************************************************************/
96: r2 = cg->radius * cg->radius;
97: r = ksp->work[0];
98: z = ksp->work[1];
99: p = ksp->work[2];
100: d = ksp->vec_sol;
101: pc = ksp->pc;
103: PCGetOperators(pc, &Qmat, &Mmat);
105: VecGetSize(d, &max_cg_its);
106: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
107: max_lanczos_its = cg->max_lanczos_its;
108: max_newton_its = cg->max_newton_its;
109: ksp->its = 0;
111: /***************************************************************************/
112: /* Initialize objective function direction, and minimum eigenvalue. */
113: /***************************************************************************/
115: cg->o_fcn = 0.0;
117: VecSet(d, 0.0); /* d = 0 */
118: cg->norm_d = 0.0;
120: cg->e_min = 0.0;
121: cg->lambda = 0.0;
123: /***************************************************************************/
124: /* The first phase of GLTR performs a standard conjugate gradient method, */
125: /* but stores the values required for the Lanczos matrix. We switch to */
126: /* the Lanczos when the conjugate gradient method breaks down. Check the */
127: /* right-hand side for numerical problems. The check for not-a-number and */
128: /* infinite values need be performed only once. */
129: /***************************************************************************/
131: VecCopy(ksp->vec_rhs, r); /* r = -grad */
132: VecDot(r, r, &rr); /* rr = r^T r */
133: KSPCheckDot(ksp,rr);
135: /***************************************************************************/
136: /* Check the preconditioner for numerical problems and for positive */
137: /* definiteness. The check for not-a-number and infinite values need be */
138: /* performed only once. */
139: /***************************************************************************/
141: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
142: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
143: if (PetscIsInfOrNanScalar(rz)) {
144: /*************************************************************************/
145: /* The preconditioner contains not-a-number or an infinite value. */
146: /* Return the gradient direction intersected with the trust region. */
147: /*************************************************************************/
149: ksp->reason = KSP_DIVERGED_NANORINF;
150: PetscInfo1(ksp, "KSPCGSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz);
152: if (cg->radius) {
153: if (r2 >= rr) {
154: alpha = 1.0;
155: cg->norm_d = PetscSqrtReal(rr);
156: } else {
157: alpha = PetscSqrtReal(r2 / rr);
158: cg->norm_d = cg->radius;
159: }
161: VecAXPY(d, alpha, r); /* d = d + alpha r */
163: /***********************************************************************/
164: /* Compute objective function. */
165: /***********************************************************************/
167: KSP_MatMult(ksp, Qmat, d, z);
168: VecAYPX(z, -0.5, ksp->vec_rhs);
169: VecDot(d, z, &cg->o_fcn);
170: cg->o_fcn = -cg->o_fcn;
171: ++ksp->its;
172: }
173: return(0);
174: }
176: if (rz < 0.0) {
177: /*************************************************************************/
178: /* The preconditioner is indefinite. Because this is the first */
179: /* and we do not have a direction yet, we use the gradient step. Note */
180: /* that we cannot use the preconditioned norm when computing the step */
181: /* because the matrix is indefinite. */
182: /*************************************************************************/
184: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
185: PetscInfo1(ksp, "KSPCGSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz);
187: if (cg->radius) {
188: if (r2 >= rr) {
189: alpha = 1.0;
190: cg->norm_d = PetscSqrtReal(rr);
191: } else {
192: alpha = PetscSqrtReal(r2 / rr);
193: cg->norm_d = cg->radius;
194: }
196: VecAXPY(d, alpha, r); /* d = d + alpha r */
198: /***********************************************************************/
199: /* Compute objective function. */
200: /***********************************************************************/
202: KSP_MatMult(ksp, Qmat, d, z);
203: VecAYPX(z, -0.5, ksp->vec_rhs);
204: VecDot(d, z, &cg->o_fcn);
205: cg->o_fcn = -cg->o_fcn;
206: ++ksp->its;
207: }
208: return(0);
209: }
211: /***************************************************************************/
212: /* As far as we know, the preconditioner is positive semidefinite. */
213: /* Compute and log the residual. Check convergence because this */
214: /* initializes things, but do not terminate until at least one conjugate */
215: /* gradient iteration has been performed. */
216: /***************************************************************************/
218: cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
220: switch (ksp->normtype) {
221: case KSP_NORM_PRECONDITIONED:
222: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
223: break;
225: case KSP_NORM_UNPRECONDITIONED:
226: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
227: break;
229: case KSP_NORM_NATURAL:
230: norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
231: break;
233: default:
234: norm_r = 0.0;
235: break;
236: }
238: KSPLogResidualHistory(ksp, norm_r);
239: KSPMonitor(ksp, ksp->its, norm_r);
240: ksp->rnorm = norm_r;
242: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
244: /***************************************************************************/
245: /* Compute the first direction and update the iteration. */
246: /***************************************************************************/
248: VecCopy(z, p); /* p = z */
249: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
250: ++ksp->its;
252: /***************************************************************************/
253: /* Check the matrix for numerical problems. */
254: /***************************************************************************/
256: VecDot(p, z, &kappa); /* kappa = p^T Q p */
257: if (PetscIsInfOrNanScalar(kappa)) {
258: /*************************************************************************/
259: /* The matrix produced not-a-number or an infinite value. In this case, */
260: /* we must stop and use the gradient direction. This condition need */
261: /* only be checked once. */
262: /*************************************************************************/
264: ksp->reason = KSP_DIVERGED_NANORINF;
265: PetscInfo1(ksp, "KSPCGSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa);
267: if (cg->radius) {
268: if (r2 >= rr) {
269: alpha = 1.0;
270: cg->norm_d = PetscSqrtReal(rr);
271: } else {
272: alpha = PetscSqrtReal(r2 / rr);
273: cg->norm_d = cg->radius;
274: }
276: VecAXPY(d, alpha, r); /* d = d + alpha r */
278: /***********************************************************************/
279: /* Compute objective function. */
280: /***********************************************************************/
282: KSP_MatMult(ksp, Qmat, d, z);
283: VecAYPX(z, -0.5, ksp->vec_rhs);
284: VecDot(d, z, &cg->o_fcn);
285: cg->o_fcn = -cg->o_fcn;
286: ++ksp->its;
287: }
288: return(0);
289: }
291: /***************************************************************************/
292: /* Initialize variables for calculating the norm of the direction and for */
293: /* the Lanczos tridiagonal matrix. Note that we track the diagonal value */
294: /* of the Cholesky factorization of the Lanczos matrix in order to */
295: /* determine when negative curvature is encountered. */
296: /***************************************************************************/
298: dMp = 0.0;
299: norm_d = 0.0;
300: switch (cg->dtype) {
301: case GLTR_PRECONDITIONED_DIRECTION:
302: norm_p = rz;
303: break;
305: default:
306: VecDot(p, p, &norm_p);
307: break;
308: }
310: cg->diag[t_size] = kappa / rz;
311: cg->offd[t_size] = 0.0;
312: ++t_size;
314: piv = 1.0;
316: /***************************************************************************/
317: /* Check for breakdown of the conjugate gradient method; this occurs when */
318: /* kappa is zero. */
319: /***************************************************************************/
321: if (PetscAbsReal(kappa) <= 0.0) {
322: /*************************************************************************/
323: /* The curvature is zero. In this case, we must stop and use follow */
324: /* the direction of negative curvature since the Lanczos matrix is zero. */
325: /*************************************************************************/
327: ksp->reason = KSP_DIVERGED_BREAKDOWN;
328: PetscInfo1(ksp, "KSPCGSolve_GLTR: breakdown: kappa=%g\n", (double)kappa);
330: if (cg->radius && norm_p > 0.0) {
331: /***********************************************************************/
332: /* Follow direction of negative curvature to the boundary of the */
333: /* trust region. */
334: /***********************************************************************/
336: step = PetscSqrtReal(r2 / norm_p);
337: cg->norm_d = cg->radius;
339: VecAXPY(d, step, p); /* d = d + step p */
341: /***********************************************************************/
342: /* Update objective function. */
343: /***********************************************************************/
345: cg->o_fcn += step * (0.5 * step * kappa - rz);
346: } else if (cg->radius) {
347: /***********************************************************************/
348: /* The norm of the preconditioned direction is zero; use the gradient */
349: /* step. */
350: /***********************************************************************/
352: if (r2 >= rr) {
353: alpha = 1.0;
354: cg->norm_d = PetscSqrtReal(rr);
355: } else {
356: alpha = PetscSqrtReal(r2 / rr);
357: cg->norm_d = cg->radius;
358: }
360: VecAXPY(d, alpha, r); /* d = d + alpha r */
362: /***********************************************************************/
363: /* Compute objective function. */
364: /***********************************************************************/
366: KSP_MatMult(ksp, Qmat, d, z);
367: VecAYPX(z, -0.5, ksp->vec_rhs);
368: VecDot(d, z, &cg->o_fcn);
369: cg->o_fcn = -cg->o_fcn;
370: ++ksp->its;
371: }
372: return(0);
373: }
375: /***************************************************************************/
376: /* Begin the first part of the GLTR algorithm which runs the conjugate */
377: /* gradient method until either the problem is solved, we encounter the */
378: /* boundary of the trust region, or the conjugate gradient method breaks */
379: /* down. */
380: /***************************************************************************/
382: while (1) {
383: /*************************************************************************/
384: /* Know that kappa is nonzero, because we have not broken down, so we */
385: /* can compute the steplength. */
386: /*************************************************************************/
388: alpha = rz / kappa;
389: cg->alpha[l_size] = alpha;
391: /*************************************************************************/
392: /* Compute the diagonal value of the Cholesky factorization of the */
393: /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
394: /* This indicates a direction of negative curvature. */
395: /*************************************************************************/
397: piv = cg->diag[l_size] - cg->offd[l_size]*cg->offd[l_size] / piv;
398: if (piv <= 0.0) {
399: /***********************************************************************/
400: /* In this case, the matrix is indefinite and we have encountered */
401: /* a direction of negative curvature. Follow the direction to the */
402: /* boundary of the trust region. */
403: /***********************************************************************/
405: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
406: PetscInfo1(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa);
408: if (cg->radius && norm_p > 0.0) {
409: /*********************************************************************/
410: /* Follow direction of negative curvature to boundary. */
411: /*********************************************************************/
413: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
414: cg->norm_d = cg->radius;
416: VecAXPY(d, step, p); /* d = d + step p */
418: /*********************************************************************/
419: /* Update objective function. */
420: /*********************************************************************/
422: cg->o_fcn += step * (0.5 * step * kappa - rz);
423: } else if (cg->radius) {
424: /*********************************************************************/
425: /* The norm of the direction is zero; there is nothing to follow. */
426: /*********************************************************************/
427: }
428: break;
429: }
431: /*************************************************************************/
432: /* Compute the steplength and check for intersection with the trust */
433: /* region. */
434: /*************************************************************************/
436: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
437: if (cg->radius && norm_dp1 >= r2) {
438: /***********************************************************************/
439: /* In this case, the matrix is positive definite as far as we know. */
440: /* However, the full step goes beyond the trust region. */
441: /***********************************************************************/
443: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
444: PetscInfo1(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius);
446: if (norm_p > 0.0) {
447: /*********************************************************************/
448: /* Follow the direction to the boundary of the trust region. */
449: /*********************************************************************/
451: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
452: cg->norm_d = cg->radius;
454: VecAXPY(d, step, p); /* d = d + step p */
456: /*********************************************************************/
457: /* Update objective function. */
458: /*********************************************************************/
460: cg->o_fcn += step * (0.5 * step * kappa - rz);
461: } else {
462: /*********************************************************************/
463: /* The norm of the direction is zero; there is nothing to follow. */
464: /*********************************************************************/
465: }
466: break;
467: }
469: /*************************************************************************/
470: /* Now we can update the direction and residual. */
471: /*************************************************************************/
473: VecAXPY(d, alpha, p); /* d = d + alpha p */
474: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
475: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
477: switch (cg->dtype) {
478: case GLTR_PRECONDITIONED_DIRECTION:
479: norm_d = norm_dp1;
480: break;
482: default:
483: VecDot(d, d, &norm_d);
484: break;
485: }
486: cg->norm_d = PetscSqrtReal(norm_d);
488: /*************************************************************************/
489: /* Update objective function. */
490: /*************************************************************************/
492: cg->o_fcn -= 0.5 * alpha * rz;
494: /*************************************************************************/
495: /* Check that the preconditioner appears positive semidefinite. */
496: /*************************************************************************/
498: rzm1 = rz;
499: VecDot(r, z, &rz); /* rz = r^T z */
500: if (rz < 0.0) {
501: /***********************************************************************/
502: /* The preconditioner is indefinite. */
503: /***********************************************************************/
505: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
506: PetscInfo1(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz);
507: break;
508: }
510: /*************************************************************************/
511: /* As far as we know, the preconditioner is positive semidefinite. */
512: /* Compute the residual and check for convergence. */
513: /*************************************************************************/
515: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
517: switch (ksp->normtype) {
518: case KSP_NORM_PRECONDITIONED:
519: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
520: break;
522: case KSP_NORM_UNPRECONDITIONED:
523: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
524: break;
526: case KSP_NORM_NATURAL:
527: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
528: break;
530: default:
531: norm_r = 0.0;
532: break;
533: }
535: KSPLogResidualHistory(ksp, norm_r);
536: KSPMonitor(ksp, ksp->its, norm_r);
537: ksp->rnorm = norm_r;
539: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
540: if (ksp->reason) {
541: /***********************************************************************/
542: /* The method has converged. */
543: /***********************************************************************/
545: PetscInfo2(ksp, "KSPCGSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
546: break;
547: }
549: /*************************************************************************/
550: /* We have not converged yet. Check for breakdown. */
551: /*************************************************************************/
553: beta = rz / rzm1;
554: if (PetscAbsReal(beta) <= 0.0) {
555: /***********************************************************************/
556: /* Conjugate gradients has broken down. */
557: /***********************************************************************/
559: ksp->reason = KSP_DIVERGED_BREAKDOWN;
560: PetscInfo1(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta);
561: break;
562: }
564: /*************************************************************************/
565: /* Check iteration limit. */
566: /*************************************************************************/
568: if (ksp->its >= max_cg_its) {
569: ksp->reason = KSP_DIVERGED_ITS;
570: PetscInfo1(ksp, "KSPCGSolve_GLTR: iterlim: its=%D\n", ksp->its);
571: break;
572: }
574: /*************************************************************************/
575: /* Update p and the norms. */
576: /*************************************************************************/
578: cg->beta[l_size] = beta;
579: VecAYPX(p, beta, z); /* p = z + beta p */
581: switch (cg->dtype) {
582: case GLTR_PRECONDITIONED_DIRECTION:
583: dMp = beta*(dMp + alpha*norm_p);
584: norm_p = beta*(rzm1 + beta*norm_p);
585: break;
587: default:
588: VecDot(d, p, &dMp);
589: VecDot(p, p, &norm_p);
590: break;
591: }
593: /*************************************************************************/
594: /* Compute the new direction and update the iteration. */
595: /*************************************************************************/
597: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
598: VecDot(p, z, &kappa); /* kappa = p^T Q p */
599: ++ksp->its;
601: /*************************************************************************/
602: /* Update the Lanczos tridiagonal matrix. */
603: /*************************************************************************/
605: ++l_size;
606: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
607: cg->diag[t_size] = kappa / rz + beta / alpha;
608: ++t_size;
610: /*************************************************************************/
611: /* Check for breakdown of the conjugate gradient method; this occurs */
612: /* when kappa is zero. */
613: /*************************************************************************/
615: if (PetscAbsReal(kappa) <= 0.0) {
616: /***********************************************************************/
617: /* The method breaks down; move along the direction as if the matrix */
618: /* were indefinite. */
619: /***********************************************************************/
621: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
622: PetscInfo1(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa);
624: if (cg->radius && norm_p > 0.0) {
625: /*********************************************************************/
626: /* Follow direction to boundary. */
627: /*********************************************************************/
629: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
630: cg->norm_d = cg->radius;
632: VecAXPY(d, step, p); /* d = d + step p */
634: /*********************************************************************/
635: /* Update objective function. */
636: /*********************************************************************/
638: cg->o_fcn += step * (0.5 * step * kappa - rz);
639: } else if (cg->radius) {
640: /*********************************************************************/
641: /* The norm of the direction is zero; there is nothing to follow. */
642: /*********************************************************************/
643: }
644: break;
645: }
646: }
648: /***************************************************************************/
649: /* Check to see if we need to continue with the Lanczos method. */
650: /***************************************************************************/
652: if (!cg->radius) {
653: /*************************************************************************/
654: /* There is no radius. Therefore, we cannot move along the boundary. */
655: /*************************************************************************/
656: return(0);
657: }
659: if (KSP_CONVERGED_CG_NEG_CURVE != ksp->reason) {
660: /*************************************************************************/
661: /* The method either converged to an interior point, hit the boundary of */
662: /* the trust-region without encountering a direction of negative */
663: /* curvature or the iteration limit was reached. */
664: /*************************************************************************/
665: return(0);
666: }
668: /***************************************************************************/
669: /* Switch to contructing the Lanczos basis by way of the conjugate */
670: /* directions. */
671: /***************************************************************************/
673: for (i = 0; i < max_lanczos_its; ++i) {
674: /*************************************************************************/
675: /* Check for breakdown of the conjugate gradient method; this occurs */
676: /* when kappa is zero. */
677: /*************************************************************************/
679: if (PetscAbsReal(kappa) <= 0.0) {
680: ksp->reason = KSP_DIVERGED_BREAKDOWN;
681: PetscInfo1(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa);
682: break;
683: }
685: /*************************************************************************/
686: /* Update the direction and residual. */
687: /*************************************************************************/
689: alpha = rz / kappa;
690: cg->alpha[l_size] = alpha;
692: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
693: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
695: /*************************************************************************/
696: /* Check that the preconditioner appears positive semidefinite. */
697: /*************************************************************************/
699: rzm1 = rz;
700: VecDot(r, z, &rz); /* rz = r^T z */
701: if (rz < 0.0) {
702: /***********************************************************************/
703: /* The preconditioner is indefinite. */
704: /***********************************************************************/
706: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
707: PetscInfo1(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz);
708: break;
709: }
711: /*************************************************************************/
712: /* As far as we know, the preconditioner is positive definite. Compute */
713: /* the residual. Do NOT check for convergence. */
714: /*************************************************************************/
716: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
718: switch (ksp->normtype) {
719: case KSP_NORM_PRECONDITIONED:
720: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
721: break;
723: case KSP_NORM_UNPRECONDITIONED:
724: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
725: break;
727: case KSP_NORM_NATURAL:
728: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
729: break;
731: default:
732: norm_r = 0.0;
733: break;
734: }
736: KSPLogResidualHistory(ksp, norm_r);
737: KSPMonitor(ksp, ksp->its, norm_r);
738: ksp->rnorm = norm_r;
740: /*************************************************************************/
741: /* Check for breakdown. */
742: /*************************************************************************/
744: beta = rz / rzm1;
745: if (PetscAbsReal(beta) <= 0.0) {
746: /***********************************************************************/
747: /* Conjugate gradients has broken down. */
748: /***********************************************************************/
750: ksp->reason = KSP_DIVERGED_BREAKDOWN;
751: PetscInfo1(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n",(double) beta);
752: break;
753: }
755: /*************************************************************************/
756: /* Update p and the norms. */
757: /*************************************************************************/
759: cg->beta[l_size] = beta;
760: VecAYPX(p, beta, z); /* p = z + beta p */
762: /*************************************************************************/
763: /* Compute the new direction and update the iteration. */
764: /*************************************************************************/
766: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
767: VecDot(p, z, &kappa); /* kappa = p^T Q p */
768: ++ksp->its;
770: /*************************************************************************/
771: /* Update the Lanczos tridiagonal matrix. */
772: /*************************************************************************/
774: ++l_size;
775: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
776: cg->diag[t_size] = kappa / rz + beta / alpha;
777: ++t_size;
778: }
780: /***************************************************************************/
781: /* We have the Lanczos basis, solve the tridiagonal trust-region problem */
782: /* to obtain the Lanczos direction. We know that the solution lies on */
783: /* the boundary of the trust region. We start by checking that the */
784: /* workspace allocated is large enough. */
785: /***************************************************************************/
786: /* Note that the current version only computes the solution by using the */
787: /* preconditioned direction. Need to think about how to do the */
788: /* unpreconditioned direction calculation. */
789: /***************************************************************************/
791: if (t_size > cg->alloced) {
792: if (cg->alloced) {
793: PetscFree(cg->rwork);
794: PetscFree(cg->iwork);
795: cg->alloced += cg->init_alloc;
796: } else {
797: cg->alloced = cg->init_alloc;
798: }
800: while (t_size > cg->alloced) {
801: cg->alloced += cg->init_alloc;
802: }
804: cg->alloced = PetscMin(cg->alloced, t_size);
805: PetscMalloc2(10*cg->alloced, &cg->rwork,5*cg->alloced, &cg->iwork);
806: }
808: /***************************************************************************/
809: /* Set up the required vectors. */
810: /***************************************************************************/
812: t_soln = cg->rwork + 0*t_size; /* Solution */
813: t_diag = cg->rwork + 1*t_size; /* Diagonal of T */
814: t_offd = cg->rwork + 2*t_size; /* Off-diagonal of T */
815: e_valu = cg->rwork + 3*t_size; /* Eigenvalues of T */
816: e_vect = cg->rwork + 4*t_size; /* Eigenvector of T */
817: e_rwrk = cg->rwork + 5*t_size; /* Eigen workspace */
819: e_iblk = cg->iwork + 0*t_size; /* Eigen blocks */
820: e_splt = cg->iwork + 1*t_size; /* Eigen splits */
821: e_iwrk = cg->iwork + 2*t_size; /* Eigen workspace */
823: /***************************************************************************/
824: /* Compute the minimum eigenvalue of T. */
825: /***************************************************************************/
827: vl = 0.0;
828: vu = 0.0;
829: il = 1;
830: iu = 1;
832: PetscStackCallBLAS("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));
834: if ((0 != info) || (1 != e_valus)) {
835: /*************************************************************************/
836: /* Calculation of the minimum eigenvalue failed. Return the */
837: /* Steihaug-Toint direction. */
838: /*************************************************************************/
840: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n");
841: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
842: return(0);
843: }
845: cg->e_min = e_valu[0];
847: /***************************************************************************/
848: /* Compute the initial value of lambda to make (T + lamba I) positive */
849: /* definite. */
850: /***************************************************************************/
852: pert = cg->init_pert;
853: if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
855: while (1) {
856: for (i = 0; i < t_size; ++i) {
857: t_diag[i] = cg->diag[i] + cg->lambda;
858: t_offd[i] = cg->offd[i];
859: }
861: PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
863: if (0 == info) break;
865: pert += pert;
866: cg->lambda = cg->lambda * (1.0 + pert) + pert;
867: }
869: /***************************************************************************/
870: /* Compute the initial step and its norm. */
871: /***************************************************************************/
873: nrhs = 1;
874: nldb = t_size;
876: t_soln[0] = -cg->norm_r[0];
877: for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
879: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
881: if (0 != info) {
882: /*************************************************************************/
883: /* Calculation of the initial step failed; return the Steihaug-Toint */
884: /* direction. */
885: /*************************************************************************/
887: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
888: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
889: return(0);
890: }
892: norm_t = 0.;
893: for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
894: norm_t = PetscSqrtReal(norm_t);
896: /***************************************************************************/
897: /* Determine the case we are in. */
898: /***************************************************************************/
900: if (norm_t <= cg->radius) {
901: /*************************************************************************/
902: /* The step is within the trust region; check if we are in the hard case */
903: /* and need to move to the boundary by following a direction of negative */
904: /* curvature. */
905: /*************************************************************************/
907: if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
908: /***********************************************************************/
909: /* This is the hard case; compute the eigenvector associated with the */
910: /* minimum eigenvalue and move along this direction to the boundary. */
911: /***********************************************************************/
913: PetscStackCallBLAS("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));
915: if (0 != info) {
916: /*********************************************************************/
917: /* Calculation of the minimum eigenvalue failed. Return the */
918: /* Steihaug-Toint direction. */
919: /*********************************************************************/
921: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n");
922: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
923: return(0);
924: }
926: coef1 = 0.0;
927: coef2 = 0.0;
928: coef3 = -cg->radius * cg->radius;
929: for (i = 0; i < t_size; ++i) {
930: coef1 += e_vect[i] * e_vect[i];
931: coef2 += e_vect[i] * t_soln[i];
932: coef3 += t_soln[i] * t_soln[i];
933: }
935: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
936: root1 = (-coef2 + coef3) / coef1;
937: root2 = (-coef2 - coef3) / coef1;
939: /***********************************************************************/
940: /* Compute objective value for (t_soln + root1 * e_vect) */
941: /***********************************************************************/
943: for (i = 0; i < t_size; ++i) {
944: e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
945: }
947: obj1 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
948: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
949: for (i = 1; i < t_size - 1; ++i) {
950: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
951: cg->diag[i]*e_rwrk[i]+
952: cg->offd[i+1]*e_rwrk[i+1]);
953: }
954: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
955: cg->diag[i]*e_rwrk[i]);
957: /***********************************************************************/
958: /* Compute objective value for (t_soln + root2 * e_vect) */
959: /***********************************************************************/
961: for (i = 0; i < t_size; ++i) {
962: e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
963: }
965: obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
966: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
967: for (i = 1; i < t_size - 1; ++i) {
968: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
969: cg->diag[i]*e_rwrk[i]+
970: cg->offd[i+1]*e_rwrk[i+1]);
971: }
972: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
973: cg->diag[i]*e_rwrk[i]);
975: /***********************************************************************/
976: /* Choose the point with the best objective function value. */
977: /***********************************************************************/
979: if (obj1 <= obj2) {
980: for (i = 0; i < t_size; ++i) {
981: t_soln[i] += root1 * e_vect[i];
982: }
983: }
984: else {
985: for (i = 0; i < t_size; ++i) {
986: t_soln[i] += root2 * e_vect[i];
987: }
988: }
989: } else {
990: /***********************************************************************/
991: /* The matrix is positive definite or there was no room to move; the */
992: /* solution is already contained in t_soln. */
993: /***********************************************************************/
994: }
995: } else {
996: /*************************************************************************/
997: /* The step is outside the trust-region. Compute the correct value for */
998: /* lambda by performing Newton's method. */
999: /*************************************************************************/
1001: for (i = 0; i < max_newton_its; ++i) {
1002: /***********************************************************************/
1003: /* Check for convergence. */
1004: /***********************************************************************/
1006: if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
1008: /***********************************************************************/
1009: /* Compute the update. */
1010: /***********************************************************************/
1012: PetscArraycpy(e_rwrk, t_soln, t_size);
1014: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
1016: if (0 != info) {
1017: /*********************************************************************/
1018: /* Calculation of the step failed; return the Steihaug-Toint */
1019: /* direction. */
1020: /*********************************************************************/
1022: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
1023: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1024: return(0);
1025: }
1027: /***********************************************************************/
1028: /* Modify lambda. */
1029: /***********************************************************************/
1031: norm_w = 0.;
1032: for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
1034: cg->lambda += (norm_t - cg->radius)/cg->radius * (norm_t * norm_t) / norm_w;
1036: /***********************************************************************/
1037: /* Factor T + lambda I */
1038: /***********************************************************************/
1040: for (j = 0; j < t_size; ++j) {
1041: t_diag[j] = cg->diag[j] + cg->lambda;
1042: t_offd[j] = cg->offd[j];
1043: }
1045: PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1047: if (0 != info) {
1048: /*********************************************************************/
1049: /* Calculation of factorization failed; return the Steihaug-Toint */
1050: /* direction. */
1051: /*********************************************************************/
1053: PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n");
1054: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1055: return(0);
1056: }
1058: /***********************************************************************/
1059: /* Compute the new step and its norm. */
1060: /***********************************************************************/
1062: t_soln[0] = -cg->norm_r[0];
1063: for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
1065: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1067: if (0 != info) {
1068: /*********************************************************************/
1069: /* Calculation of the step failed; return the Steihaug-Toint */
1070: /* direction. */
1071: /*********************************************************************/
1073: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
1074: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1075: return(0);
1076: }
1078: norm_t = 0.;
1079: for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1080: norm_t = PetscSqrtReal(norm_t);
1081: }
1083: /*************************************************************************/
1084: /* Check for convergence. */
1085: /*************************************************************************/
1087: if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1088: /***********************************************************************/
1089: /* Newton method failed to converge in iteration limit. */
1090: /***********************************************************************/
1092: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n");
1093: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1094: return(0);
1095: }
1096: }
1098: /***************************************************************************/
1099: /* Recover the norm of the direction and objective function value. */
1100: /***************************************************************************/
1102: cg->norm_d = norm_t;
1104: cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1105: for (i = 1; i < t_size - 1; ++i) {
1106: 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]);
1107: }
1108: cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);
1110: /***************************************************************************/
1111: /* Recover the direction. */
1112: /***************************************************************************/
1114: sigma = -1;
1116: /***************************************************************************/
1117: /* Start conjugate gradient method from the beginning */
1118: /***************************************************************************/
1120: VecCopy(ksp->vec_rhs, r); /* r = -grad */
1121: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1123: /***************************************************************************/
1124: /* Accumulate Q * s */
1125: /***************************************************************************/
1127: VecCopy(z, d);
1128: VecScale(d, sigma * t_soln[0] / cg->norm_r[0]);
1130: /***************************************************************************/
1131: /* Compute the first direction. */
1132: /***************************************************************************/
1134: VecCopy(z, p); /* p = z */
1135: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1136: ++ksp->its;
1138: for (i = 0; i < l_size - 1; ++i) {
1139: /*************************************************************************/
1140: /* Update the residual and direction. */
1141: /*************************************************************************/
1143: alpha = cg->alpha[i];
1144: if (alpha >= 0.0) sigma = -sigma;
1146: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1147: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1149: /*************************************************************************/
1150: /* Accumulate Q * s */
1151: /*************************************************************************/
1153: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1155: /*************************************************************************/
1156: /* Update p. */
1157: /*************************************************************************/
1159: beta = cg->beta[i];
1160: VecAYPX(p, beta, z); /* p = z + beta p */
1161: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1162: ++ksp->its;
1163: }
1165: /***************************************************************************/
1166: /* Update the residual and direction. */
1167: /***************************************************************************/
1169: alpha = cg->alpha[i];
1170: if (alpha >= 0.0) sigma = -sigma;
1172: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1173: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1175: /***************************************************************************/
1176: /* Accumulate Q * s */
1177: /***************************************************************************/
1179: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1181: /***************************************************************************/
1182: /* Set the termination reason. */
1183: /***************************************************************************/
1185: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1186: return(0);
1187: #endif
1188: }
1190: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
1191: {
1192: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1193: PetscInt max_its;
1197: /***************************************************************************/
1198: /* Determine the total maximum number of iterations. */
1199: /***************************************************************************/
1201: max_its = ksp->max_it + cg->max_lanczos_its + 1;
1203: /***************************************************************************/
1204: /* Set work vectors needed by conjugate gradient method and allocate */
1205: /* workspace for Lanczos matrix. */
1206: /***************************************************************************/
1208: KSPSetWorkVecs(ksp, 3);
1209: if (cg->diag) {
1210: PetscArrayzero(cg->diag, max_its);
1211: PetscArrayzero(cg->offd, max_its);
1212: PetscArrayzero(cg->alpha, max_its);
1213: PetscArrayzero(cg->beta, max_its);
1214: PetscArrayzero(cg->norm_r, max_its);
1215: } else {
1216: PetscCalloc5(max_its,&cg->diag,max_its,&cg->offd,max_its,&cg->alpha,max_its,&cg->beta,max_its,&cg->norm_r);
1217: PetscLogObjectMemory((PetscObject)ksp, 5*max_its*sizeof(PetscReal));
1218: }
1219: return(0);
1220: }
1222: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
1223: {
1224: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1228: /***************************************************************************/
1229: /* Free memory allocated for the data. */
1230: /***************************************************************************/
1232: PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1233: if (cg->alloced) {
1234: PetscFree2(cg->rwork,cg->iwork);
1235: }
1237: /***************************************************************************/
1238: /* Clear composed functions */
1239: /***************************************************************************/
1241: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",NULL);
1242: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",NULL);
1243: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",NULL);
1244: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",NULL);
1245: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",NULL);
1247: /***************************************************************************/
1248: /* Destroy KSP object. */
1249: /***************************************************************************/
1250: KSPDestroyDefault(ksp);
1251: return(0);
1252: }
1254: static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
1255: {
1256: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1259: cg->radius = radius;
1260: return(0);
1261: }
1263: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1264: {
1265: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1268: *norm_d = cg->norm_d;
1269: return(0);
1270: }
1272: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1273: {
1274: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1277: *o_fcn = cg->o_fcn;
1278: return(0);
1279: }
1281: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1282: {
1283: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1286: *e_min = cg->e_min;
1287: return(0);
1288: }
1290: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1291: {
1292: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1295: *lambda = cg->lambda;
1296: return(0);
1297: }
1299: static PetscErrorCode KSPCGSetFromOptions_GLTR(PetscOptionItems *PetscOptionsObject,KSP ksp)
1300: {
1302: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1305: PetscOptionsHead(PetscOptionsObject,"KSP GLTR options");
1307: PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL);
1309: PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
1311: PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1312: PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1313: PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);
1315: PetscOptionsInt("-ksp_cg_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL);
1316: PetscOptionsInt("-ksp_cg_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL);
1318: PetscOptionsTail();
1319: return(0);
1320: }
1322: /*MC
1323: KSPGLTR - Code to run conjugate gradient method subject to a constraint
1324: on the solution norm. This is used in Trust Region methods for
1325: nonlinear equations, SNESNEWTONTR
1327: Options Database Keys:
1328: . -ksp_cg_radius <r> - Trust Region Radius
1330: Notes:
1331: This is rarely used directly
1333: Use preconditioned conjugate gradient to compute
1334: an approximate minimizer of the quadratic function
1336: q(s) = g^T * s + .5 * s^T * H * s
1338: subject to the trust region constraint
1340: || s || <= delta,
1342: where
1344: delta is the trust region radius,
1345: g is the gradient vector,
1346: H is the Hessian approximation,
1347: M is the positive definite preconditioner matrix.
1349: KSPConvergedReason may be
1350: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1351: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1352: $ other KSP converged/diverged reasons
1354: Notes:
1355: The preconditioner supplied should be symmetric and positive definite.
1357: Reference:
1358: Gould, N. and Lucidi, S. and Roma, M. and Toint, P., Solving the Trust-Region Subproblem using the Lanczos Method,
1359: SIAM Journal on Optimization, volume 9, number 2, 1999, 504-525
1361: Level: developer
1363: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCGSetRadius(), KSPCGGetNormD(), KSPCGGetObjFcn(), KSPGLTRGetMinEig(), KSPGLTRGetLambda()
1364: M*/
1366: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1367: {
1369: KSPCG_GLTR *cg;
1372: PetscNewLog(ksp,&cg);
1373: cg->radius = 0.0;
1374: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1376: cg->init_pert = 1.0e-8;
1377: cg->eigen_tol = 1.0e-10;
1378: cg->newton_tol = 1.0e-6;
1380: cg->alloced = 0;
1381: cg->init_alloc = 1024;
1383: cg->max_lanczos_its = 20;
1384: cg->max_newton_its = 10;
1386: ksp->data = (void*) cg;
1387: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
1388: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
1389: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
1390: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
1392: /***************************************************************************/
1393: /* Sets the functions that are associated with this data structure */
1394: /* (in C++ this is the same as defining virtual functions). */
1395: /***************************************************************************/
1397: ksp->ops->setup = KSPCGSetUp_GLTR;
1398: ksp->ops->solve = KSPCGSolve_GLTR;
1399: ksp->ops->destroy = KSPCGDestroy_GLTR;
1400: ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1401: ksp->ops->buildsolution = KSPBuildSolutionDefault;
1402: ksp->ops->buildresidual = KSPBuildResidualDefault;
1403: ksp->ops->view = NULL;
1405: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",KSPCGSetRadius_GLTR);
1406: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C", KSPCGGetNormD_GLTR);
1407: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",KSPCGGetObjFcn_GLTR);
1408: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",KSPGLTRGetMinEig_GLTR);
1409: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",KSPGLTRGetLambda_GLTR);
1410: return(0);
1411: }