Actual source code: gltr.c
petsc-3.5.4 2015-05-23
2: #include <petsc-private/kspimpl.h> /*I "petscksp.h" I*/
3: #include <../src/ksp/ksp/impls/cg/gltr/gltrimpl.h>
4: #include <petscblaslapack.h>
6: #define GLTR_PRECONDITIONED_DIRECTION 0
7: #define GLTR_UNPRECONDITIONED_DIRECTION 1
8: #define GLTR_DIRECTION_TYPES 2
10: static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};
14: /*@
15: KSPGLTRSetRadius - Sets the radius of the trust region.
17: Logically Collective on KSP
19: Input Parameters:
20: + ksp - the iterative context
21: - radius - the trust region radius (Infinity is the default)
23: Options Database Key:
24: . -ksp_gltr_radius <r>
26: Level: advanced
28: .keywords: KSP, GLTR, set, trust region radius
29: @*/
30: PetscErrorCode KSPGLTRSetRadius(KSP ksp, PetscReal radius)
31: {
36: if (radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Radius negative");
38: PetscTryMethod(ksp,"KSPGLTRSetRadius_C",(KSP,PetscReal),(ksp,radius));
39: return(0);
40: }
44: /*@
45: KSPGLTRGetNormD - Get norm of the direction.
47: Collective
49: Input Parameters:
50: + ksp - the iterative context
51: - norm_d - the norm of the direction
53: Level: advanced
55: .keywords: KSP, GLTR, get, norm direction
56: @*/
57: PetscErrorCode KSPGLTRGetNormD(KSP ksp, PetscReal *norm_d)
58: {
63: PetscUseMethod(ksp,"KSPGLTRGetNormD_C",(KSP,PetscReal*),(ksp,norm_d));
64: return(0);
65: }
69: /*@
70: KSPGLTRGetObjFcn - Get objective function value.
72: Note Collective
74: Input Parameters:
75: + ksp - the iterative context
76: - o_fcn - the objective function value
78: Level: advanced
80: .keywords: KSP, GLTR, get, objective function
81: @*/
82: PetscErrorCode KSPGLTRGetObjFcn(KSP ksp, PetscReal *o_fcn)
83: {
88: PetscUseMethod(ksp,"KSPGLTRGetObjFcn_C",(KSP,PetscReal*),(ksp,o_fcn));
89: return(0);
90: }
94: /*@
95: KSPGLTRGetMinEig - Get minimum eigenvalue.
97: Collective on KSP
99: Input Parameters:
100: + ksp - the iterative context
101: - e_min - the minimum eigenvalue
103: Level: advanced
105: .keywords: KSP, GLTR, get, minimum eigenvalue
106: @*/
107: PetscErrorCode KSPGLTRGetMinEig(KSP ksp, PetscReal *e_min)
108: {
113: PetscUseMethod(ksp,"KSPGLTRGetMinEig_C",(KSP,PetscReal*),(ksp,e_min));
114: return(0);
115: }
119: /*@
120: KSPGLTRGetLambda - Get multiplier on trust-region constraint.
122: Not Collective
124: Input Parameters:
125: + ksp - the iterative context
126: - lambda - the multiplier
128: Level: advanced
130: .keywords: KSP, GLTR, get, multiplier
131: @*/
132: PetscErrorCode KSPGLTRGetLambda(KSP ksp, PetscReal *lambda)
133: {
138: PetscUseMethod(ksp,"KSPGLTRGetLambda_C",(KSP,PetscReal*),(ksp,lambda));
139: return(0);
140: }
144: PetscErrorCode KSPSolve_GLTR(KSP ksp)
145: {
146: #if defined(PETSC_USE_COMPLEX)
147: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "GLTR is not available for complex systems");
148: #else
149: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
150: PetscReal *t_soln, *t_diag, *t_offd, *e_valu, *e_vect, *e_rwrk;
151: PetscBLASInt *e_iblk, *e_splt, *e_iwrk;
154: Mat Qmat, Mmat;
155: Vec r, z, p, d;
156: PC pc;
158: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
159: PetscReal alpha, beta, kappa, rz, rzm1;
160: PetscReal rr, r2, piv, step;
161: PetscReal vl, vu;
162: PetscReal coef1, coef2, coef3, root1, root2, obj1, obj2;
163: PetscReal norm_t, norm_w, pert;
165: PetscInt i, j, max_cg_its, max_lanczos_its, max_newton_its, sigma;
166: PetscBLASInt t_size = 0, l_size = 0, il, iu, info;
167: PetscBLASInt nrhs, nldb;
169: #if !defined(PETSC_MISSING_LAPACK_STEBZ)
170: PetscBLASInt e_valus, e_splts;
171: #endif
172: PetscBool diagonalscale;
175: /***************************************************************************/
176: /* Check the arguments and parameters. */
177: /***************************************************************************/
179: PCGetDiagonalScale(ksp->pc, &diagonalscale);
180: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
181: if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
183: /***************************************************************************/
184: /* Get the workspace vectors and initialize variables */
185: /***************************************************************************/
187: r2 = cg->radius * cg->radius;
188: r = ksp->work[0];
189: z = ksp->work[1];
190: p = ksp->work[2];
191: d = ksp->vec_sol;
192: pc = ksp->pc;
194: PCGetOperators(pc, &Qmat, &Mmat);
196: VecGetSize(d, &max_cg_its);
197: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
198: max_lanczos_its = cg->max_lanczos_its;
199: max_newton_its = cg->max_newton_its;
200: ksp->its = 0;
202: /***************************************************************************/
203: /* Initialize objective function direction, and minimum eigenvalue. */
204: /***************************************************************************/
206: cg->o_fcn = 0.0;
208: VecSet(d, 0.0); /* d = 0 */
209: cg->norm_d = 0.0;
211: cg->e_min = 0.0;
212: cg->lambda = 0.0;
214: /***************************************************************************/
215: /* The first phase of GLTR performs a standard conjugate gradient method, */
216: /* but stores the values required for the Lanczos matrix. We switch to */
217: /* the Lanczos when the conjugate gradient method breaks down. Check the */
218: /* right-hand side for numerical problems. The check for not-a-number and */
219: /* infinite values need be performed only once. */
220: /***************************************************************************/
222: VecCopy(ksp->vec_rhs, r); /* r = -grad */
223: VecDot(r, r, &rr); /* rr = r^T r */
224: if (PetscIsInfOrNanScalar(rr)) {
225: /*************************************************************************/
226: /* The right-hand side contains not-a-number or an infinite value. */
227: /* The gradient step does not work; return a zero value for the step. */
228: /*************************************************************************/
230: ksp->reason = KSP_DIVERGED_NANORINF;
231: PetscInfo1(ksp, "KSPSolve_GLTR: bad right-hand side: rr=%g\n", (double)rr);
232: return(0);
233: }
235: /***************************************************************************/
236: /* Check the preconditioner for numerical problems and for positive */
237: /* definiteness. The check for not-a-number and infinite values need be */
238: /* performed only once. */
239: /***************************************************************************/
241: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
242: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
243: if (PetscIsInfOrNanScalar(rz)) {
244: /*************************************************************************/
245: /* The preconditioner contains not-a-number or an infinite value. */
246: /* Return the gradient direction intersected with the trust region. */
247: /*************************************************************************/
249: ksp->reason = KSP_DIVERGED_NANORINF;
250: PetscInfo1(ksp, "KSPSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz);
252: if (cg->radius) {
253: if (r2 >= rr) {
254: alpha = 1.0;
255: cg->norm_d = PetscSqrtReal(rr);
256: } else {
257: alpha = PetscSqrtReal(r2 / rr);
258: cg->norm_d = cg->radius;
259: }
261: VecAXPY(d, alpha, r); /* d = d + alpha r */
263: /***********************************************************************/
264: /* Compute objective function. */
265: /***********************************************************************/
267: KSP_MatMult(ksp, Qmat, d, z);
268: VecAYPX(z, -0.5, ksp->vec_rhs);
269: VecDot(d, z, &cg->o_fcn);
270: cg->o_fcn = -cg->o_fcn;
271: ++ksp->its;
272: }
273: return(0);
274: }
276: if (rz < 0.0) {
277: /*************************************************************************/
278: /* The preconditioner is indefinite. Because this is the first */
279: /* and we do not have a direction yet, we use the gradient step. Note */
280: /* that we cannot use the preconditioned norm when computing the step */
281: /* because the matrix is indefinite. */
282: /*************************************************************************/
284: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
285: PetscInfo1(ksp, "KSPSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz);
287: if (cg->radius) {
288: if (r2 >= rr) {
289: alpha = 1.0;
290: cg->norm_d = PetscSqrtReal(rr);
291: } else {
292: alpha = PetscSqrtReal(r2 / rr);
293: cg->norm_d = cg->radius;
294: }
296: VecAXPY(d, alpha, r); /* d = d + alpha r */
298: /***********************************************************************/
299: /* Compute objective function. */
300: /***********************************************************************/
302: KSP_MatMult(ksp, Qmat, d, z);
303: VecAYPX(z, -0.5, ksp->vec_rhs);
304: VecDot(d, z, &cg->o_fcn);
305: cg->o_fcn = -cg->o_fcn;
306: ++ksp->its;
307: }
308: return(0);
309: }
311: /***************************************************************************/
312: /* As far as we know, the preconditioner is positive semidefinite. */
313: /* Compute and log the residual. Check convergence because this */
314: /* initializes things, but do not terminate until at least one conjugate */
315: /* gradient iteration has been performed. */
316: /***************************************************************************/
318: cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
320: switch (ksp->normtype) {
321: case KSP_NORM_PRECONDITIONED:
322: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
323: break;
325: case KSP_NORM_UNPRECONDITIONED:
326: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
327: break;
329: case KSP_NORM_NATURAL:
330: norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
331: break;
333: default:
334: norm_r = 0.0;
335: break;
336: }
338: KSPLogResidualHistory(ksp, norm_r);
339: KSPMonitor(ksp, ksp->its, norm_r);
340: ksp->rnorm = norm_r;
342: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
344: /***************************************************************************/
345: /* Compute the first direction and update the iteration. */
346: /***************************************************************************/
348: VecCopy(z, p); /* p = z */
349: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
350: ++ksp->its;
352: /***************************************************************************/
353: /* Check the matrix for numerical problems. */
354: /***************************************************************************/
356: VecDot(p, z, &kappa); /* kappa = p^T Q p */
357: if (PetscIsInfOrNanScalar(kappa)) {
358: /*************************************************************************/
359: /* The matrix produced not-a-number or an infinite value. In this case, */
360: /* we must stop and use the gradient direction. This condition need */
361: /* only be checked once. */
362: /*************************************************************************/
364: ksp->reason = KSP_DIVERGED_NANORINF;
365: PetscInfo1(ksp, "KSPSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa);
367: if (cg->radius) {
368: if (r2 >= rr) {
369: alpha = 1.0;
370: cg->norm_d = PetscSqrtReal(rr);
371: } else {
372: alpha = PetscSqrtReal(r2 / rr);
373: cg->norm_d = cg->radius;
374: }
376: VecAXPY(d, alpha, r); /* d = d + alpha r */
378: /***********************************************************************/
379: /* Compute objective function. */
380: /***********************************************************************/
382: KSP_MatMult(ksp, Qmat, d, z);
383: VecAYPX(z, -0.5, ksp->vec_rhs);
384: VecDot(d, z, &cg->o_fcn);
385: cg->o_fcn = -cg->o_fcn;
386: ++ksp->its;
387: }
388: return(0);
389: }
391: /***************************************************************************/
392: /* Initialize variables for calculating the norm of the direction and for */
393: /* the Lanczos tridiagonal matrix. Note that we track the diagonal value */
394: /* of the Cholesky factorization of the Lanczos matrix in order to */
395: /* determine when negative curvature is encountered. */
396: /***************************************************************************/
398: dMp = 0.0;
399: norm_d = 0.0;
400: switch (cg->dtype) {
401: case GLTR_PRECONDITIONED_DIRECTION:
402: norm_p = rz;
403: break;
405: default:
406: VecDot(p, p, &norm_p);
407: break;
408: }
410: cg->diag[t_size] = kappa / rz;
411: cg->offd[t_size] = 0.0;
412: ++t_size;
414: piv = 1.0;
416: /***************************************************************************/
417: /* Check for breakdown of the conjugate gradient method; this occurs when */
418: /* kappa is zero. */
419: /***************************************************************************/
421: if (PetscAbsReal(kappa) <= 0.0) {
422: /*************************************************************************/
423: /* The curvature is zero. In this case, we must stop and use follow */
424: /* the direction of negative curvature since the Lanczos matrix is zero. */
425: /*************************************************************************/
427: ksp->reason = KSP_DIVERGED_BREAKDOWN;
428: PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: kappa=%g\n", (double)kappa);
430: if (cg->radius && norm_p > 0.0) {
431: /***********************************************************************/
432: /* Follow direction of negative curvature to the boundary of the */
433: /* trust region. */
434: /***********************************************************************/
436: step = PetscSqrtReal(r2 / norm_p);
437: cg->norm_d = cg->radius;
439: VecAXPY(d, step, p); /* d = d + step p */
441: /***********************************************************************/
442: /* Update objective function. */
443: /***********************************************************************/
445: cg->o_fcn += step * (0.5 * step * kappa - rz);
446: } else if (cg->radius) {
447: /***********************************************************************/
448: /* The norm of the preconditioned direction is zero; use the gradient */
449: /* step. */
450: /***********************************************************************/
452: if (r2 >= rr) {
453: alpha = 1.0;
454: cg->norm_d = PetscSqrtReal(rr);
455: } else {
456: alpha = PetscSqrtReal(r2 / rr);
457: cg->norm_d = cg->radius;
458: }
460: VecAXPY(d, alpha, r); /* d = d + alpha r */
462: /***********************************************************************/
463: /* Compute objective function. */
464: /***********************************************************************/
466: KSP_MatMult(ksp, Qmat, d, z);
467: VecAYPX(z, -0.5, ksp->vec_rhs);
468: VecDot(d, z, &cg->o_fcn);
469: cg->o_fcn = -cg->o_fcn;
470: ++ksp->its;
471: }
472: return(0);
473: }
475: /***************************************************************************/
476: /* Begin the first part of the GLTR algorithm which runs the conjugate */
477: /* gradient method until either the problem is solved, we encounter the */
478: /* boundary of the trust region, or the conjugate gradient method breaks */
479: /* down. */
480: /***************************************************************************/
482: while (1) {
483: /*************************************************************************/
484: /* Know that kappa is nonzero, because we have not broken down, so we */
485: /* can compute the steplength. */
486: /*************************************************************************/
488: alpha = rz / kappa;
489: cg->alpha[l_size] = alpha;
491: /*************************************************************************/
492: /* Compute the diagonal value of the Cholesky factorization of the */
493: /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
494: /* This indicates a direction of negative curvature. */
495: /*************************************************************************/
497: piv = cg->diag[l_size] - cg->offd[l_size]*cg->offd[l_size] / piv;
498: if (piv <= 0.0) {
499: /***********************************************************************/
500: /* In this case, the matrix is indefinite and we have encountered */
501: /* a direction of negative curvature. Follow the direction to the */
502: /* boundary of the trust region. */
503: /***********************************************************************/
505: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
506: PetscInfo1(ksp, "KSPSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa);
508: if (cg->radius && norm_p > 0.0) {
509: /*********************************************************************/
510: /* Follow direction of negative curvature to boundary. */
511: /*********************************************************************/
513: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
514: cg->norm_d = cg->radius;
516: VecAXPY(d, step, p); /* d = d + step p */
518: /*********************************************************************/
519: /* Update objective function. */
520: /*********************************************************************/
522: cg->o_fcn += step * (0.5 * step * kappa - rz);
523: } else if (cg->radius) {
524: /*********************************************************************/
525: /* The norm of the direction is zero; there is nothing to follow. */
526: /*********************************************************************/
527: }
528: break;
529: }
531: /*************************************************************************/
532: /* Compute the steplength and check for intersection with the trust */
533: /* region. */
534: /*************************************************************************/
536: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
537: if (cg->radius && norm_dp1 >= r2) {
538: /***********************************************************************/
539: /* In this case, the matrix is positive definite as far as we know. */
540: /* However, the full step goes beyond the trust region. */
541: /***********************************************************************/
543: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
544: PetscInfo1(ksp, "KSPSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius);
546: if (norm_p > 0.0) {
547: /*********************************************************************/
548: /* Follow the direction to the boundary of the trust region. */
549: /*********************************************************************/
551: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
552: cg->norm_d = cg->radius;
554: VecAXPY(d, step, p); /* d = d + step p */
556: /*********************************************************************/
557: /* Update objective function. */
558: /*********************************************************************/
560: cg->o_fcn += step * (0.5 * step * kappa - rz);
561: } else {
562: /*********************************************************************/
563: /* The norm of the direction is zero; there is nothing to follow. */
564: /*********************************************************************/
565: }
566: break;
567: }
569: /*************************************************************************/
570: /* Now we can update the direction and residual. */
571: /*************************************************************************/
573: VecAXPY(d, alpha, p); /* d = d + alpha p */
574: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
575: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
577: switch (cg->dtype) {
578: case GLTR_PRECONDITIONED_DIRECTION:
579: norm_d = norm_dp1;
580: break;
582: default:
583: VecDot(d, d, &norm_d);
584: break;
585: }
586: cg->norm_d = PetscSqrtReal(norm_d);
588: /*************************************************************************/
589: /* Update objective function. */
590: /*************************************************************************/
592: cg->o_fcn -= 0.5 * alpha * rz;
594: /*************************************************************************/
595: /* Check that the preconditioner appears positive semidefinite. */
596: /*************************************************************************/
598: rzm1 = rz;
599: VecDot(r, z, &rz); /* rz = r^T z */
600: if (rz < 0.0) {
601: /***********************************************************************/
602: /* The preconditioner is indefinite. */
603: /***********************************************************************/
605: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
606: PetscInfo1(ksp, "KSPSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz);
607: break;
608: }
610: /*************************************************************************/
611: /* As far as we know, the preconditioner is positive semidefinite. */
612: /* Compute the residual and check for convergence. */
613: /*************************************************************************/
615: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
617: switch (ksp->normtype) {
618: case KSP_NORM_PRECONDITIONED:
619: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
620: break;
622: case KSP_NORM_UNPRECONDITIONED:
623: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
624: break;
626: case KSP_NORM_NATURAL:
627: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
628: break;
630: default:
631: norm_r = 0.0;
632: break;
633: }
635: KSPLogResidualHistory(ksp, norm_r);
636: KSPMonitor(ksp, ksp->its, norm_r);
637: ksp->rnorm = norm_r;
639: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
640: if (ksp->reason) {
641: /***********************************************************************/
642: /* The method has converged. */
643: /***********************************************************************/
645: PetscInfo2(ksp, "KSPSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
646: break;
647: }
649: /*************************************************************************/
650: /* We have not converged yet. Check for breakdown. */
651: /*************************************************************************/
653: beta = rz / rzm1;
654: if (PetscAbsReal(beta) <= 0.0) {
655: /***********************************************************************/
656: /* Conjugate gradients has broken down. */
657: /***********************************************************************/
659: ksp->reason = KSP_DIVERGED_BREAKDOWN;
660: PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n", (double)beta);
661: break;
662: }
664: /*************************************************************************/
665: /* Check iteration limit. */
666: /*************************************************************************/
668: if (ksp->its >= max_cg_its) {
669: ksp->reason = KSP_DIVERGED_ITS;
670: PetscInfo1(ksp, "KSPSolve_GLTR: iterlim: its=%D\n", ksp->its);
671: break;
672: }
674: /*************************************************************************/
675: /* Update p and the norms. */
676: /*************************************************************************/
678: cg->beta[l_size] = beta;
679: VecAYPX(p, beta, z); /* p = z + beta p */
681: switch (cg->dtype) {
682: case GLTR_PRECONDITIONED_DIRECTION:
683: dMp = beta*(dMp + alpha*norm_p);
684: norm_p = beta*(rzm1 + beta*norm_p);
685: break;
687: default:
688: VecDot(d, p, &dMp);
689: VecDot(p, p, &norm_p);
690: break;
691: }
693: /*************************************************************************/
694: /* Compute the new direction and update the iteration. */
695: /*************************************************************************/
697: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
698: VecDot(p, z, &kappa); /* kappa = p^T Q p */
699: ++ksp->its;
701: /*************************************************************************/
702: /* Update the Lanczos tridiagonal matrix. */
703: /*************************************************************************/
705: ++l_size;
706: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
707: cg->diag[t_size] = kappa / rz + beta / alpha;
708: ++t_size;
710: /*************************************************************************/
711: /* Check for breakdown of the conjugate gradient method; this occurs */
712: /* when kappa is zero. */
713: /*************************************************************************/
715: if (PetscAbsReal(kappa) <= 0.0) {
716: /***********************************************************************/
717: /* The method breaks down; move along the direction as if the matrix */
718: /* were indefinite. */
719: /***********************************************************************/
721: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
722: PetscInfo1(ksp, "KSPSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa);
724: if (cg->radius && norm_p > 0.0) {
725: /*********************************************************************/
726: /* Follow direction to boundary. */
727: /*********************************************************************/
729: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
730: cg->norm_d = cg->radius;
732: VecAXPY(d, step, p); /* d = d + step p */
734: /*********************************************************************/
735: /* Update objective function. */
736: /*********************************************************************/
738: cg->o_fcn += step * (0.5 * step * kappa - rz);
739: } else if (cg->radius) {
740: /*********************************************************************/
741: /* The norm of the direction is zero; there is nothing to follow. */
742: /*********************************************************************/
743: }
744: break;
745: }
746: }
748: /***************************************************************************/
749: /* Check to see if we need to continue with the Lanczos method. */
750: /***************************************************************************/
752: if (!cg->radius) {
753: /*************************************************************************/
754: /* There is no radius. Therefore, we cannot move along the boundary. */
755: /*************************************************************************/
756: return(0);
757: }
759: if (KSP_CONVERGED_CG_NEG_CURVE != ksp->reason) {
760: /*************************************************************************/
761: /* The method either converged to an interior point, hit the boundary of */
762: /* the trust-region without encountering a direction of negative */
763: /* curvature or the iteration limit was reached. */
764: /*************************************************************************/
765: return(0);
766: }
768: /***************************************************************************/
769: /* Switch to contructing the Lanczos basis by way of the conjugate */
770: /* directions. */
771: /***************************************************************************/
773: for (i = 0; i < max_lanczos_its; ++i) {
774: /*************************************************************************/
775: /* Check for breakdown of the conjugate gradient method; this occurs */
776: /* when kappa is zero. */
777: /*************************************************************************/
779: if (PetscAbsReal(kappa) <= 0.0) {
780: ksp->reason = KSP_DIVERGED_BREAKDOWN;
781: PetscInfo1(ksp, "KSPSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa);
782: break;
783: }
785: /*************************************************************************/
786: /* Update the direction and residual. */
787: /*************************************************************************/
789: alpha = rz / kappa;
790: cg->alpha[l_size] = alpha;
792: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
793: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
795: /*************************************************************************/
796: /* Check that the preconditioner appears positive semidefinite. */
797: /*************************************************************************/
799: rzm1 = rz;
800: VecDot(r, z, &rz); /* rz = r^T z */
801: if (rz < 0.0) {
802: /***********************************************************************/
803: /* The preconditioner is indefinite. */
804: /***********************************************************************/
806: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
807: PetscInfo1(ksp, "KSPSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz);
808: break;
809: }
811: /*************************************************************************/
812: /* As far as we know, the preconditioner is positive definite. Compute */
813: /* the residual. Do NOT check for convergence. */
814: /*************************************************************************/
816: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
818: switch (ksp->normtype) {
819: case KSP_NORM_PRECONDITIONED:
820: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
821: break;
823: case KSP_NORM_UNPRECONDITIONED:
824: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
825: break;
827: case KSP_NORM_NATURAL:
828: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
829: break;
831: default:
832: norm_r = 0.0;
833: break;
834: }
836: KSPLogResidualHistory(ksp, norm_r);
837: KSPMonitor(ksp, ksp->its, norm_r);
838: ksp->rnorm = norm_r;
840: /*************************************************************************/
841: /* Check for breakdown. */
842: /*************************************************************************/
844: beta = rz / rzm1;
845: if (PetscAbsReal(beta) <= 0.0) {
846: /***********************************************************************/
847: /* Conjugate gradients has broken down. */
848: /***********************************************************************/
850: ksp->reason = KSP_DIVERGED_BREAKDOWN;
851: PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n",(double) beta);
852: break;
853: }
855: /*************************************************************************/
856: /* Update p and the norms. */
857: /*************************************************************************/
859: cg->beta[l_size] = beta;
860: VecAYPX(p, beta, z); /* p = z + beta p */
862: /*************************************************************************/
863: /* Compute the new direction and update the iteration. */
864: /*************************************************************************/
866: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
867: VecDot(p, z, &kappa); /* kappa = p^T Q p */
868: ++ksp->its;
870: /*************************************************************************/
871: /* Update the Lanczos tridiagonal matrix. */
872: /*************************************************************************/
874: ++l_size;
875: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
876: cg->diag[t_size] = kappa / rz + beta / alpha;
877: ++t_size;
878: }
880: /***************************************************************************/
881: /* We have the Lanczos basis, solve the tridiagonal trust-region problem */
882: /* to obtain the Lanczos direction. We know that the solution lies on */
883: /* the boundary of the trust region. We start by checking that the */
884: /* workspace allocated is large enough. */
885: /***************************************************************************/
886: /* Note that the current version only computes the solution by using the */
887: /* preconditioned direction. Need to think about how to do the */
888: /* unpreconditioned direction calculation. */
889: /***************************************************************************/
891: if (t_size > cg->alloced) {
892: if (cg->alloced) {
893: PetscFree(cg->rwork);
894: PetscFree(cg->iwork);
895: cg->alloced += cg->init_alloc;
896: } else {
897: cg->alloced = cg->init_alloc;
898: }
900: while (t_size > cg->alloced) {
901: cg->alloced += cg->init_alloc;
902: }
904: cg->alloced = PetscMin(cg->alloced, t_size);
905: PetscMalloc2(10*cg->alloced, &cg->rwork,5*cg->alloced, &cg->iwork);
906: }
908: /***************************************************************************/
909: /* Set up the required vectors. */
910: /***************************************************************************/
912: t_soln = cg->rwork + 0*t_size; /* Solution */
913: t_diag = cg->rwork + 1*t_size; /* Diagonal of T */
914: t_offd = cg->rwork + 2*t_size; /* Off-diagonal of T */
915: e_valu = cg->rwork + 3*t_size; /* Eigenvalues of T */
916: e_vect = cg->rwork + 4*t_size; /* Eigenvector of T */
917: e_rwrk = cg->rwork + 5*t_size; /* Eigen workspace */
919: e_iblk = cg->iwork + 0*t_size; /* Eigen blocks */
920: e_splt = cg->iwork + 1*t_size; /* Eigen splits */
921: e_iwrk = cg->iwork + 2*t_size; /* Eigen workspace */
923: /***************************************************************************/
924: /* Compute the minimum eigenvalue of T. */
925: /***************************************************************************/
927: vl = 0.0;
928: vu = 0.0;
929: il = 1;
930: iu = 1;
932: #if defined(PETSC_MISSING_LAPACK_STEBZ)
933: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEBZ - Lapack routine is unavailable.");
934: #else
935: 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));
937: if ((0 != info) || (1 != e_valus)) {
938: /*************************************************************************/
939: /* Calculation of the minimum eigenvalue failed. Return the */
940: /* Steihaug-Toint direction. */
941: /*************************************************************************/
943: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvalue.\n");
944: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
945: return(0);
946: }
948: cg->e_min = e_valu[0];
950: /***************************************************************************/
951: /* Compute the initial value of lambda to make (T + lamba I) positive */
952: /* definite. */
953: /***************************************************************************/
955: pert = cg->init_pert;
956: if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
957: #endif
959: while (1) {
960: for (i = 0; i < t_size; ++i) {
961: t_diag[i] = cg->diag[i] + cg->lambda;
962: t_offd[i] = cg->offd[i];
963: }
965: #if defined(PETSC_MISSING_LAPACK_PTTRF)
966: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
967: #else
968: PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
970: if (0 == info) break;
972: pert += pert;
973: cg->lambda = cg->lambda * (1.0 + pert) + pert;
974: #endif
975: }
977: /***************************************************************************/
978: /* Compute the initial step and its norm. */
979: /***************************************************************************/
981: nrhs = 1;
982: nldb = t_size;
984: t_soln[0] = -cg->norm_r[0];
985: for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
987: #if defined(PETSC_MISSING_LAPACK_PTTRS)
988: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
989: #else
990: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
991: #endif
993: if (0 != info) {
994: /*************************************************************************/
995: /* Calculation of the initial step failed; return the Steihaug-Toint */
996: /* direction. */
997: /*************************************************************************/
999: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1000: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1001: return(0);
1002: }
1004: norm_t = 0.;
1005: for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
1006: norm_t = PetscSqrtReal(norm_t);
1008: /***************************************************************************/
1009: /* Determine the case we are in. */
1010: /***************************************************************************/
1012: if (norm_t <= cg->radius) {
1013: /*************************************************************************/
1014: /* The step is within the trust region; check if we are in the hard case */
1015: /* and need to move to the boundary by following a direction of negative */
1016: /* curvature. */
1017: /*************************************************************************/
1019: if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
1020: /***********************************************************************/
1021: /* This is the hard case; compute the eigenvector associated with the */
1022: /* minimum eigenvalue and move along this direction to the boundary. */
1023: /***********************************************************************/
1025: #if defined(PETSC_MISSING_LAPACK_STEIN)
1026: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEIN - Lapack routine is unavailable.");
1027: #else
1028: 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));
1029: #endif
1031: if (0 != info) {
1032: /*********************************************************************/
1033: /* Calculation of the minimum eigenvalue failed. Return the */
1034: /* Steihaug-Toint direction. */
1035: /*********************************************************************/
1037: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvector.\n");
1038: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1039: return(0);
1040: }
1042: coef1 = 0.0;
1043: coef2 = 0.0;
1044: coef3 = -cg->radius * cg->radius;
1045: for (i = 0; i < t_size; ++i) {
1046: coef1 += e_vect[i] * e_vect[i];
1047: coef2 += e_vect[i] * t_soln[i];
1048: coef3 += t_soln[i] * t_soln[i];
1049: }
1051: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
1052: root1 = (-coef2 + coef3) / coef1;
1053: root2 = (-coef2 - coef3) / coef1;
1055: /***********************************************************************/
1056: /* Compute objective value for (t_soln + root1 * e_vect) */
1057: /***********************************************************************/
1059: for (i = 0; i < t_size; ++i) {
1060: e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
1061: }
1063: obj1 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1064: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1065: for (i = 1; i < t_size - 1; ++i) {
1066: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1067: cg->diag[i]*e_rwrk[i]+
1068: cg->offd[i+1]*e_rwrk[i+1]);
1069: }
1070: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1071: cg->diag[i]*e_rwrk[i]);
1073: /***********************************************************************/
1074: /* Compute objective value for (t_soln + root2 * e_vect) */
1075: /***********************************************************************/
1077: for (i = 0; i < t_size; ++i) {
1078: e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
1079: }
1081: obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1082: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1083: for (i = 1; i < t_size - 1; ++i) {
1084: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1085: cg->diag[i]*e_rwrk[i]+
1086: cg->offd[i+1]*e_rwrk[i+1]);
1087: }
1088: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1089: cg->diag[i]*e_rwrk[i]);
1091: /***********************************************************************/
1092: /* Choose the point with the best objective function value. */
1093: /***********************************************************************/
1095: if (obj1 <= obj2) {
1096: for (i = 0; i < t_size; ++i) {
1097: t_soln[i] += root1 * e_vect[i];
1098: }
1099: }
1100: else {
1101: for (i = 0; i < t_size; ++i) {
1102: t_soln[i] += root2 * e_vect[i];
1103: }
1104: }
1105: } else {
1106: /***********************************************************************/
1107: /* The matrix is positive definite or there was no room to move; the */
1108: /* solution is already contained in t_soln. */
1109: /***********************************************************************/
1110: }
1111: } else {
1112: /*************************************************************************/
1113: /* The step is outside the trust-region. Compute the correct value for */
1114: /* lambda by performing Newton's method. */
1115: /*************************************************************************/
1117: for (i = 0; i < max_newton_its; ++i) {
1118: /***********************************************************************/
1119: /* Check for convergence. */
1120: /***********************************************************************/
1122: if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
1124: /***********************************************************************/
1125: /* Compute the update. */
1126: /***********************************************************************/
1128: PetscMemcpy(e_rwrk, t_soln, sizeof(PetscReal)*t_size);
1130: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1131: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1132: #else
1133: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
1134: #endif
1136: if (0 != info) {
1137: /*********************************************************************/
1138: /* Calculation of the step failed; return the Steihaug-Toint */
1139: /* direction. */
1140: /*********************************************************************/
1142: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1143: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1144: return(0);
1145: }
1147: /***********************************************************************/
1148: /* Modify lambda. */
1149: /***********************************************************************/
1151: norm_w = 0.;
1152: for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
1154: cg->lambda += (norm_t - cg->radius)/cg->radius * (norm_t * norm_t) / norm_w;
1156: /***********************************************************************/
1157: /* Factor T + lambda I */
1158: /***********************************************************************/
1160: for (j = 0; j < t_size; ++j) {
1161: t_diag[j] = cg->diag[j] + cg->lambda;
1162: t_offd[j] = cg->offd[j];
1163: }
1165: #if defined(PETSC_MISSING_LAPACK_PTTRF)
1166: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
1167: #else
1168: PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1169: #endif
1171: if (0 != info) {
1172: /*********************************************************************/
1173: /* Calculation of factorization failed; return the Steihaug-Toint */
1174: /* direction. */
1175: /*********************************************************************/
1177: PetscInfo(ksp, "KSPSolve_GLTR: factorization failed.\n");
1178: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1179: return(0);
1180: }
1182: /***********************************************************************/
1183: /* Compute the new step and its norm. */
1184: /***********************************************************************/
1186: t_soln[0] = -cg->norm_r[0];
1187: for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
1189: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1190: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1191: #else
1192: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1193: #endif
1195: if (0 != info) {
1196: /*********************************************************************/
1197: /* Calculation of the step failed; return the Steihaug-Toint */
1198: /* direction. */
1199: /*********************************************************************/
1201: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1202: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1203: return(0);
1204: }
1206: norm_t = 0.;
1207: for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1208: norm_t = PetscSqrtReal(norm_t);
1209: }
1211: /*************************************************************************/
1212: /* Check for convergence. */
1213: /*************************************************************************/
1215: if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1216: /***********************************************************************/
1217: /* Newton method failed to converge in iteration limit. */
1218: /***********************************************************************/
1220: PetscInfo(ksp, "KSPSolve_GLTR: failed to converge.\n");
1221: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1222: return(0);
1223: }
1224: }
1226: /***************************************************************************/
1227: /* Recover the norm of the direction and objective function value. */
1228: /***************************************************************************/
1230: cg->norm_d = norm_t;
1232: cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1233: for (i = 1; i < t_size - 1; ++i) {
1234: 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]);
1235: }
1236: cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);
1238: /***************************************************************************/
1239: /* Recover the direction. */
1240: /***************************************************************************/
1242: sigma = -1;
1244: /***************************************************************************/
1245: /* Start conjugate gradient method from the beginning */
1246: /***************************************************************************/
1248: VecCopy(ksp->vec_rhs, r); /* r = -grad */
1249: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1251: /***************************************************************************/
1252: /* Accumulate Q * s */
1253: /***************************************************************************/
1255: VecCopy(z, d);
1256: VecScale(d, sigma * t_soln[0] / cg->norm_r[0]);
1258: /***************************************************************************/
1259: /* Compute the first direction. */
1260: /***************************************************************************/
1262: VecCopy(z, p); /* p = z */
1263: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1264: ++ksp->its;
1266: for (i = 0; i < l_size - 1; ++i) {
1267: /*************************************************************************/
1268: /* Update the residual and direction. */
1269: /*************************************************************************/
1271: alpha = cg->alpha[i];
1272: if (alpha >= 0.0) sigma = -sigma;
1274: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1275: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1277: /*************************************************************************/
1278: /* Accumulate Q * s */
1279: /*************************************************************************/
1281: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1283: /*************************************************************************/
1284: /* Update p. */
1285: /*************************************************************************/
1287: beta = cg->beta[i];
1288: VecAYPX(p, beta, z); /* p = z + beta p */
1289: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1290: ++ksp->its;
1291: }
1293: /***************************************************************************/
1294: /* Update the residual and direction. */
1295: /***************************************************************************/
1297: alpha = cg->alpha[i];
1298: if (alpha >= 0.0) sigma = -sigma;
1300: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1301: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1303: /***************************************************************************/
1304: /* Accumulate Q * s */
1305: /***************************************************************************/
1307: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1309: /***************************************************************************/
1310: /* Set the termination reason. */
1311: /***************************************************************************/
1313: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1314: return(0);
1315: #endif
1316: }
1320: PetscErrorCode KSPSetUp_GLTR(KSP ksp)
1321: {
1322: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1323: PetscInt max_its;
1327: /***************************************************************************/
1328: /* Determine the total maximum number of iterations. */
1329: /***************************************************************************/
1331: max_its = ksp->max_it + cg->max_lanczos_its + 1;
1333: /***************************************************************************/
1334: /* Set work vectors needed by conjugate gradient method and allocate */
1335: /* workspace for Lanczos matrix. */
1336: /***************************************************************************/
1338: KSPSetWorkVecs(ksp, 3);
1340: PetscCalloc5(max_its,&cg->diag,max_its,&cg->offd,max_its,&cg->alpha,max_its,&cg->beta,max_its,&cg->norm_r);
1341: PetscLogObjectMemory((PetscObject)ksp, 5*max_its*sizeof(PetscReal));
1342: return(0);
1343: }
1347: PetscErrorCode KSPDestroy_GLTR(KSP ksp)
1348: {
1349: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1353: /***************************************************************************/
1354: /* Free memory allocated for the data. */
1355: /***************************************************************************/
1357: PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1358: if (cg->alloced) {
1359: PetscFree2(cg->rwork,cg->iwork);
1360: }
1362: /***************************************************************************/
1363: /* Clear composed functions */
1364: /***************************************************************************/
1366: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRSetRadius_C",NULL);
1367: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetNormD_C",NULL);
1368: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetObjFcn_C",NULL);
1369: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",NULL);
1370: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",NULL);
1372: /***************************************************************************/
1373: /* Destroy KSP object. */
1374: /***************************************************************************/
1375: KSPDestroyDefault(ksp);
1376: return(0);
1377: }
1381: static PetscErrorCode KSPGLTRSetRadius_GLTR(KSP ksp, PetscReal radius)
1382: {
1383: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1386: cg->radius = radius;
1387: return(0);
1388: }
1392: static PetscErrorCode KSPGLTRGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1393: {
1394: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1397: *norm_d = cg->norm_d;
1398: return(0);
1399: }
1403: static PetscErrorCode KSPGLTRGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1404: {
1405: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1408: *o_fcn = cg->o_fcn;
1409: return(0);
1410: }
1414: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1415: {
1416: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1419: *e_min = cg->e_min;
1420: return(0);
1421: }
1425: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1426: {
1427: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1430: *lambda = cg->lambda;
1431: return(0);
1432: }
1436: PetscErrorCode KSPSetFromOptions_GLTR(KSP ksp)
1437: {
1439: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1442: PetscOptionsHead("KSP GLTR options");
1444: PetscOptionsReal("-ksp_gltr_radius", "Trust Region Radius", "KSPGLTRSetRadius", cg->radius, &cg->radius, NULL);
1446: PetscOptionsReal("-ksp_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1447: PetscOptionsReal("-ksp_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1448: PetscOptionsReal("-ksp_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);
1450: PetscOptionsInt("-ksp_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL);
1451: PetscOptionsInt("-ksp_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL);
1453: PetscOptionsEList("-ksp_gltr_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
1455: PetscOptionsTail();
1456: return(0);
1457: }
1459: /*MC
1460: KSPGLTR - Code to run conjugate gradient method subject to a constraint
1461: on the solution norm. This is used in Trust Region methods for
1462: nonlinear equations, SNESNEWTONTR
1464: Options Database Keys:
1465: . -ksp_gltr_radius <r> - Trust Region Radius
1467: Notes: This is rarely used directly
1469: Use preconditioned conjugate gradient to compute
1470: an approximate minimizer of the quadratic function
1472: q(s) = g^T * s + .5 * s^T * H * s
1474: subject to the trust region constraint
1476: || s || <= delta,
1478: where
1480: delta is the trust region radius,
1481: g is the gradient vector,
1482: H is the Hessian approximation,
1483: M is the positive definite preconditioner matrix.
1485: KSPConvergedReason may be
1486: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1487: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1488: $ other KSP converged/diverged reasons
1490: Notes:
1491: The preconditioner supplied should be symmetric and positive definite.
1493: Level: developer
1495: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGLTRSetRadius(), KSPGLTRGetNormD(), KSPGLTRGetObjFcn(), KSPGLTRGetMinEig(), KSPGLTRGetLambda()
1496: M*/
1500: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1501: {
1503: KSP_GLTR *cg;
1506: PetscNewLog(ksp,&cg);
1507: cg->radius = 0.0;
1508: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1510: cg->init_pert = 1.0e-8;
1511: cg->eigen_tol = 1.0e-10;
1512: cg->newton_tol = 1.0e-6;
1514: cg->alloced = 0;
1515: cg->init_alloc = 1024;
1517: cg->max_lanczos_its = 20;
1518: cg->max_newton_its = 10;
1520: ksp->data = (void*) cg;
1521: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
1522: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
1523: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
1525: /***************************************************************************/
1526: /* Sets the functions that are associated with this data structure */
1527: /* (in C++ this is the same as defining virtual functions). */
1528: /***************************************************************************/
1530: ksp->ops->setup = KSPSetUp_GLTR;
1531: ksp->ops->solve = KSPSolve_GLTR;
1532: ksp->ops->destroy = KSPDestroy_GLTR;
1533: ksp->ops->setfromoptions = KSPSetFromOptions_GLTR;
1534: ksp->ops->buildsolution = KSPBuildSolutionDefault;
1535: ksp->ops->buildresidual = KSPBuildResidualDefault;
1536: ksp->ops->view = 0;
1538: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRSetRadius_C",KSPGLTRSetRadius_GLTR);
1539: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetNormD_C", KSPGLTRGetNormD_GLTR);
1540: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetObjFcn_C",KSPGLTRGetObjFcn_GLTR);
1541: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",KSPGLTRGetMinEig_GLTR);
1542: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",KSPGLTRGetLambda_GLTR);
1543: return(0);
1544: }