Actual source code: gltr.c
petsc-3.4.5 2014-06-29
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: MatStructure pflag;
155: Mat Qmat, Mmat;
156: Vec r, z, p, d;
157: PC pc;
159: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
160: PetscReal alpha, beta, kappa, rz, rzm1;
161: PetscReal rr, r2, piv, step;
162: PetscReal vl, vu;
163: PetscReal coef1, coef2, coef3, root1, root2, obj1, obj2;
164: PetscReal norm_t, norm_w, pert;
166: PetscInt i, j, max_cg_its, max_lanczos_its, max_newton_its, sigma;
167: PetscBLASInt t_size = 0, l_size = 0, il, iu, info;
168: PetscBLASInt nrhs, nldb;
170: #if !defined(PETSC_MISSING_LAPACK_STEBZ)
171: PetscBLASInt e_valus, e_splts;
172: #endif
173: PetscBool diagonalscale;
176: /***************************************************************************/
177: /* Check the arguments and parameters. */
178: /***************************************************************************/
180: PCGetDiagonalScale(ksp->pc, &diagonalscale);
181: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
182: if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
184: /***************************************************************************/
185: /* Get the workspace vectors and initialize variables */
186: /***************************************************************************/
188: r2 = cg->radius * cg->radius;
189: r = ksp->work[0];
190: z = ksp->work[1];
191: p = ksp->work[2];
192: d = ksp->vec_sol;
193: pc = ksp->pc;
195: PCGetOperators(pc, &Qmat, &Mmat, &pflag);
197: VecGetSize(d, &max_cg_its);
198: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
199: max_lanczos_its = cg->max_lanczos_its;
200: max_newton_its = cg->max_newton_its;
201: ksp->its = 0;
203: /***************************************************************************/
204: /* Initialize objective function direction, and minimum eigenvalue. */
205: /***************************************************************************/
207: cg->o_fcn = 0.0;
209: VecSet(d, 0.0); /* d = 0 */
210: cg->norm_d = 0.0;
212: cg->e_min = 0.0;
213: cg->lambda = 0.0;
215: /***************************************************************************/
216: /* The first phase of GLTR performs a standard conjugate gradient method, */
217: /* but stores the values required for the Lanczos matrix. We switch to */
218: /* the Lanczos when the conjugate gradient method breaks down. Check the */
219: /* right-hand side for numerical problems. The check for not-a-number and */
220: /* infinite values need be performed only once. */
221: /***************************************************************************/
223: VecCopy(ksp->vec_rhs, r); /* r = -grad */
224: VecDot(r, r, &rr); /* rr = r^T r */
225: if (PetscIsInfOrNanScalar(rr)) {
226: /*************************************************************************/
227: /* The right-hand side contains not-a-number or an infinite value. */
228: /* The gradient step does not work; return a zero value for the step. */
229: /*************************************************************************/
231: ksp->reason = KSP_DIVERGED_NANORINF;
232: PetscInfo1(ksp, "KSPSolve_GLTR: bad right-hand side: rr=%g\n", rr);
233: return(0);
234: }
236: /***************************************************************************/
237: /* Check the preconditioner for numerical problems and for positive */
238: /* definiteness. The check for not-a-number and infinite values need be */
239: /* performed only once. */
240: /***************************************************************************/
242: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
243: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
244: if (PetscIsInfOrNanScalar(rz)) {
245: /*************************************************************************/
246: /* The preconditioner contains not-a-number or an infinite value. */
247: /* Return the gradient direction intersected with the trust region. */
248: /*************************************************************************/
250: ksp->reason = KSP_DIVERGED_NANORINF;
251: PetscInfo1(ksp, "KSPSolve_GLTR: bad preconditioner: rz=%g\n", rz);
253: if (cg->radius) {
254: if (r2 >= rr) {
255: alpha = 1.0;
256: cg->norm_d = PetscSqrtReal(rr);
257: } else {
258: alpha = PetscSqrtReal(r2 / rr);
259: cg->norm_d = cg->radius;
260: }
262: VecAXPY(d, alpha, r); /* d = d + alpha r */
264: /***********************************************************************/
265: /* Compute objective function. */
266: /***********************************************************************/
268: KSP_MatMult(ksp, Qmat, d, z);
269: VecAYPX(z, -0.5, ksp->vec_rhs);
270: VecDot(d, z, &cg->o_fcn);
271: cg->o_fcn = -cg->o_fcn;
272: ++ksp->its;
273: }
274: return(0);
275: }
277: if (rz < 0.0) {
278: /*************************************************************************/
279: /* The preconditioner is indefinite. Because this is the first */
280: /* and we do not have a direction yet, we use the gradient step. Note */
281: /* that we cannot use the preconditioned norm when computing the step */
282: /* because the matrix is indefinite. */
283: /*************************************************************************/
285: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
286: PetscInfo1(ksp, "KSPSolve_GLTR: indefinite preconditioner: rz=%g\n", rz);
288: if (cg->radius) {
289: if (r2 >= rr) {
290: alpha = 1.0;
291: cg->norm_d = PetscSqrtReal(rr);
292: } else {
293: alpha = PetscSqrtReal(r2 / rr);
294: cg->norm_d = cg->radius;
295: }
297: VecAXPY(d, alpha, r); /* d = d + alpha r */
299: /***********************************************************************/
300: /* Compute objective function. */
301: /***********************************************************************/
303: KSP_MatMult(ksp, Qmat, d, z);
304: VecAYPX(z, -0.5, ksp->vec_rhs);
305: VecDot(d, z, &cg->o_fcn);
306: cg->o_fcn = -cg->o_fcn;
307: ++ksp->its;
308: }
309: return(0);
310: }
312: /***************************************************************************/
313: /* As far as we know, the preconditioner is positive semidefinite. */
314: /* Compute and log the residual. Check convergence because this */
315: /* initializes things, but do not terminate until at least one conjugate */
316: /* gradient iteration has been performed. */
317: /***************************************************************************/
319: cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
321: switch (ksp->normtype) {
322: case KSP_NORM_PRECONDITIONED:
323: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
324: break;
326: case KSP_NORM_UNPRECONDITIONED:
327: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
328: break;
330: case KSP_NORM_NATURAL:
331: norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
332: break;
334: default:
335: norm_r = 0.0;
336: break;
337: }
339: KSPLogResidualHistory(ksp, norm_r);
340: KSPMonitor(ksp, ksp->its, norm_r);
341: ksp->rnorm = norm_r;
343: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
345: /***************************************************************************/
346: /* Compute the first direction and update the iteration. */
347: /***************************************************************************/
349: VecCopy(z, p); /* p = z */
350: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
351: ++ksp->its;
353: /***************************************************************************/
354: /* Check the matrix for numerical problems. */
355: /***************************************************************************/
357: VecDot(p, z, &kappa); /* kappa = p^T Q p */
358: if (PetscIsInfOrNanScalar(kappa)) {
359: /*************************************************************************/
360: /* The matrix produced not-a-number or an infinite value. In this case, */
361: /* we must stop and use the gradient direction. This condition need */
362: /* only be checked once. */
363: /*************************************************************************/
365: ksp->reason = KSP_DIVERGED_NANORINF;
366: PetscInfo1(ksp, "KSPSolve_GLTR: bad matrix: kappa=%g\n", kappa);
368: if (cg->radius) {
369: if (r2 >= rr) {
370: alpha = 1.0;
371: cg->norm_d = PetscSqrtReal(rr);
372: } else {
373: alpha = PetscSqrtReal(r2 / rr);
374: cg->norm_d = cg->radius;
375: }
377: VecAXPY(d, alpha, r); /* d = d + alpha r */
379: /***********************************************************************/
380: /* Compute objective function. */
381: /***********************************************************************/
383: KSP_MatMult(ksp, Qmat, d, z);
384: VecAYPX(z, -0.5, ksp->vec_rhs);
385: VecDot(d, z, &cg->o_fcn);
386: cg->o_fcn = -cg->o_fcn;
387: ++ksp->its;
388: }
389: return(0);
390: }
392: /***************************************************************************/
393: /* Initialize variables for calculating the norm of the direction and for */
394: /* the Lanczos tridiagonal matrix. Note that we track the diagonal value */
395: /* of the Cholesky factorization of the Lanczos matrix in order to */
396: /* determine when negative curvature is encountered. */
397: /***************************************************************************/
399: dMp = 0.0;
400: norm_d = 0.0;
401: switch (cg->dtype) {
402: case GLTR_PRECONDITIONED_DIRECTION:
403: norm_p = rz;
404: break;
406: default:
407: VecDot(p, p, &norm_p);
408: break;
409: }
411: cg->diag[t_size] = kappa / rz;
412: cg->offd[t_size] = 0.0;
413: ++t_size;
415: piv = 1.0;
417: /***************************************************************************/
418: /* Check for breakdown of the conjugate gradient method; this occurs when */
419: /* kappa is zero. */
420: /***************************************************************************/
422: if (fabs(kappa) <= 0.0) {
423: /*************************************************************************/
424: /* The curvature is zero. In this case, we must stop and use follow */
425: /* the direction of negative curvature since the Lanczos matrix is zero. */
426: /*************************************************************************/
428: ksp->reason = KSP_DIVERGED_BREAKDOWN;
429: PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: kappa=%g\n", kappa);
431: if (cg->radius && norm_p > 0.0) {
432: /***********************************************************************/
433: /* Follow direction of negative curvature to the boundary of the */
434: /* trust region. */
435: /***********************************************************************/
437: step = PetscSqrtReal(r2 / norm_p);
438: cg->norm_d = cg->radius;
440: VecAXPY(d, step, p); /* d = d + step p */
442: /***********************************************************************/
443: /* Update objective function. */
444: /***********************************************************************/
446: cg->o_fcn += step * (0.5 * step * kappa - rz);
447: } else if (cg->radius) {
448: /***********************************************************************/
449: /* The norm of the preconditioned direction is zero; use the gradient */
450: /* step. */
451: /***********************************************************************/
453: if (r2 >= rr) {
454: alpha = 1.0;
455: cg->norm_d = PetscSqrtReal(rr);
456: } else {
457: alpha = PetscSqrtReal(r2 / rr);
458: cg->norm_d = cg->radius;
459: }
461: VecAXPY(d, alpha, r); /* d = d + alpha r */
463: /***********************************************************************/
464: /* Compute objective function. */
465: /***********************************************************************/
467: KSP_MatMult(ksp, Qmat, d, z);
468: VecAYPX(z, -0.5, ksp->vec_rhs);
469: VecDot(d, z, &cg->o_fcn);
470: cg->o_fcn = -cg->o_fcn;
471: ++ksp->its;
472: }
473: return(0);
474: }
476: /***************************************************************************/
477: /* Begin the first part of the GLTR algorithm which runs the conjugate */
478: /* gradient method until either the problem is solved, we encounter the */
479: /* boundary of the trust region, or the conjugate gradient method breaks */
480: /* down. */
481: /***************************************************************************/
483: while (1) {
484: /*************************************************************************/
485: /* Know that kappa is nonzero, because we have not broken down, so we */
486: /* can compute the steplength. */
487: /*************************************************************************/
489: alpha = rz / kappa;
490: cg->alpha[l_size] = alpha;
492: /*************************************************************************/
493: /* Compute the diagonal value of the Cholesky factorization of the */
494: /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
495: /* This indicates a direction of negative curvature. */
496: /*************************************************************************/
498: piv = cg->diag[l_size] - cg->offd[l_size]*cg->offd[l_size] / piv;
499: if (piv <= 0.0) {
500: /***********************************************************************/
501: /* In this case, the matrix is indefinite and we have encountered */
502: /* a direction of negative curvature. Follow the direction to the */
503: /* boundary of the trust region. */
504: /***********************************************************************/
506: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
507: PetscInfo1(ksp, "KSPSolve_GLTR: negative curvature: kappa=%g\n", kappa);
509: if (cg->radius && norm_p > 0.0) {
510: /*********************************************************************/
511: /* Follow direction of negative curvature to boundary. */
512: /*********************************************************************/
514: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
515: cg->norm_d = cg->radius;
517: VecAXPY(d, step, p); /* d = d + step p */
519: /*********************************************************************/
520: /* Update objective function. */
521: /*********************************************************************/
523: cg->o_fcn += step * (0.5 * step * kappa - rz);
524: } else if (cg->radius) {
525: /*********************************************************************/
526: /* The norm of the direction is zero; there is nothing to follow. */
527: /*********************************************************************/
528: }
529: break;
530: }
532: /*************************************************************************/
533: /* Compute the steplength and check for intersection with the trust */
534: /* region. */
535: /*************************************************************************/
537: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
538: if (cg->radius && norm_dp1 >= r2) {
539: /***********************************************************************/
540: /* In this case, the matrix is positive definite as far as we know. */
541: /* However, the full step goes beyond the trust region. */
542: /***********************************************************************/
544: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
545: PetscInfo1(ksp, "KSPSolve_GLTR: constrained step: radius=%g\n", cg->radius);
547: if (norm_p > 0.0) {
548: /*********************************************************************/
549: /* Follow the direction to the boundary of the trust region. */
550: /*********************************************************************/
552: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
553: cg->norm_d = cg->radius;
555: VecAXPY(d, step, p); /* d = d + step p */
557: /*********************************************************************/
558: /* Update objective function. */
559: /*********************************************************************/
561: cg->o_fcn += step * (0.5 * step * kappa - rz);
562: } else {
563: /*********************************************************************/
564: /* The norm of the direction is zero; there is nothing to follow. */
565: /*********************************************************************/
566: }
567: break;
568: }
570: /*************************************************************************/
571: /* Now we can update the direction and residual. */
572: /*************************************************************************/
574: VecAXPY(d, alpha, p); /* d = d + alpha p */
575: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
576: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
578: switch (cg->dtype) {
579: case GLTR_PRECONDITIONED_DIRECTION:
580: norm_d = norm_dp1;
581: break;
583: default:
584: VecDot(d, d, &norm_d);
585: break;
586: }
587: cg->norm_d = PetscSqrtReal(norm_d);
589: /*************************************************************************/
590: /* Update objective function. */
591: /*************************************************************************/
593: cg->o_fcn -= 0.5 * alpha * rz;
595: /*************************************************************************/
596: /* Check that the preconditioner appears positive semidefinite. */
597: /*************************************************************************/
599: rzm1 = rz;
600: VecDot(r, z, &rz); /* rz = r^T z */
601: if (rz < 0.0) {
602: /***********************************************************************/
603: /* The preconditioner is indefinite. */
604: /***********************************************************************/
606: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
607: PetscInfo1(ksp, "KSPSolve_GLTR: cg indefinite preconditioner: rz=%g\n", rz);
608: break;
609: }
611: /*************************************************************************/
612: /* As far as we know, the preconditioner is positive semidefinite. */
613: /* Compute the residual and check for convergence. */
614: /*************************************************************************/
616: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
618: switch (ksp->normtype) {
619: case KSP_NORM_PRECONDITIONED:
620: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
621: break;
623: case KSP_NORM_UNPRECONDITIONED:
624: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
625: break;
627: case KSP_NORM_NATURAL:
628: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
629: break;
631: default:
632: norm_r = 0.0;
633: break;
634: }
636: KSPLogResidualHistory(ksp, norm_r);
637: KSPMonitor(ksp, ksp->its, norm_r);
638: ksp->rnorm = norm_r;
640: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
641: if (ksp->reason) {
642: /***********************************************************************/
643: /* The method has converged. */
644: /***********************************************************************/
646: PetscInfo2(ksp, "KSPSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", norm_r, cg->radius);
647: break;
648: }
650: /*************************************************************************/
651: /* We have not converged yet. Check for breakdown. */
652: /*************************************************************************/
654: beta = rz / rzm1;
655: if (fabs(beta) <= 0.0) {
656: /***********************************************************************/
657: /* Conjugate gradients has broken down. */
658: /***********************************************************************/
660: ksp->reason = KSP_DIVERGED_BREAKDOWN;
661: PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n", beta);
662: break;
663: }
665: /*************************************************************************/
666: /* Check iteration limit. */
667: /*************************************************************************/
669: if (ksp->its >= max_cg_its) {
670: ksp->reason = KSP_DIVERGED_ITS;
671: PetscInfo1(ksp, "KSPSolve_GLTR: iterlim: its=%d\n", ksp->its);
672: break;
673: }
675: /*************************************************************************/
676: /* Update p and the norms. */
677: /*************************************************************************/
679: cg->beta[l_size] = beta;
680: VecAYPX(p, beta, z); /* p = z + beta p */
682: switch (cg->dtype) {
683: case GLTR_PRECONDITIONED_DIRECTION:
684: dMp = beta*(dMp + alpha*norm_p);
685: norm_p = beta*(rzm1 + beta*norm_p);
686: break;
688: default:
689: VecDot(d, p, &dMp);
690: VecDot(p, p, &norm_p);
691: break;
692: }
694: /*************************************************************************/
695: /* Compute the new direction and update the iteration. */
696: /*************************************************************************/
698: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
699: VecDot(p, z, &kappa); /* kappa = p^T Q p */
700: ++ksp->its;
702: /*************************************************************************/
703: /* Update the the Lanczos tridiagonal matrix. */
704: /*************************************************************************/
706: ++l_size;
707: cg->offd[t_size] = PetscSqrtReal(beta) / fabs(alpha);
708: cg->diag[t_size] = kappa / rz + beta / alpha;
709: ++t_size;
711: /*************************************************************************/
712: /* Check for breakdown of the conjugate gradient method; this occurs */
713: /* when kappa is zero. */
714: /*************************************************************************/
716: if (fabs(kappa) <= 0.0) {
717: /***********************************************************************/
718: /* The method breaks down; move along the direction as if the matrix */
719: /* were indefinite. */
720: /***********************************************************************/
722: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
723: PetscInfo1(ksp, "KSPSolve_GLTR: cg breakdown: kappa=%g\n", kappa);
725: if (cg->radius && norm_p > 0.0) {
726: /*********************************************************************/
727: /* Follow direction to boundary. */
728: /*********************************************************************/
730: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
731: cg->norm_d = cg->radius;
733: VecAXPY(d, step, p); /* d = d + step p */
735: /*********************************************************************/
736: /* Update objective function. */
737: /*********************************************************************/
739: cg->o_fcn += step * (0.5 * step * kappa - rz);
740: } else if (cg->radius) {
741: /*********************************************************************/
742: /* The norm of the direction is zero; there is nothing to follow. */
743: /*********************************************************************/
744: }
745: break;
746: }
747: }
749: /***************************************************************************/
750: /* Check to see if we need to continue with the Lanczos method. */
751: /***************************************************************************/
753: if (!cg->radius) {
754: /*************************************************************************/
755: /* There is no radius. Therefore, we cannot move along the boundary. */
756: /*************************************************************************/
757: return(0);
758: }
760: if (KSP_CONVERGED_CG_NEG_CURVE != ksp->reason) {
761: /*************************************************************************/
762: /* The method either converged to an interior point, hit the boundary of */
763: /* the trust-region without encountering a direction of negative */
764: /* curvature or the iteration limit was reached. */
765: /*************************************************************************/
766: return(0);
767: }
769: /***************************************************************************/
770: /* Switch to contructing the Lanczos basis by way of the conjugate */
771: /* directions. */
772: /***************************************************************************/
774: for (i = 0; i < max_lanczos_its; ++i) {
775: /*************************************************************************/
776: /* Check for breakdown of the conjugate gradient method; this occurs */
777: /* when kappa is zero. */
778: /*************************************************************************/
780: if (fabs(kappa) <= 0.0) {
781: ksp->reason = KSP_DIVERGED_BREAKDOWN;
782: PetscInfo1(ksp, "KSPSolve_GLTR: lanczos breakdown: kappa=%g\n", kappa);
783: break;
784: }
786: /*************************************************************************/
787: /* Update the direction and residual. */
788: /*************************************************************************/
790: alpha = rz / kappa;
791: cg->alpha[l_size] = alpha;
793: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
794: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
796: /*************************************************************************/
797: /* Check that the preconditioner appears positive semidefinite. */
798: /*************************************************************************/
800: rzm1 = rz;
801: VecDot(r, z, &rz); /* rz = r^T z */
802: if (rz < 0.0) {
803: /***********************************************************************/
804: /* The preconditioner is indefinite. */
805: /***********************************************************************/
807: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
808: PetscInfo1(ksp, "KSPSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", rz);
809: break;
810: }
812: /*************************************************************************/
813: /* As far as we know, the preconditioner is positive definite. Compute */
814: /* the residual. Do NOT check for convergence. */
815: /*************************************************************************/
817: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
819: switch (ksp->normtype) {
820: case KSP_NORM_PRECONDITIONED:
821: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
822: break;
824: case KSP_NORM_UNPRECONDITIONED:
825: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
826: break;
828: case KSP_NORM_NATURAL:
829: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
830: break;
832: default:
833: norm_r = 0.0;
834: break;
835: }
837: KSPLogResidualHistory(ksp, norm_r);
838: KSPMonitor(ksp, ksp->its, norm_r);
839: ksp->rnorm = norm_r;
841: /*************************************************************************/
842: /* Check for breakdown. */
843: /*************************************************************************/
845: beta = rz / rzm1;
846: if (fabs(beta) <= 0.0) {
847: /***********************************************************************/
848: /* Conjugate gradients has broken down. */
849: /***********************************************************************/
851: ksp->reason = KSP_DIVERGED_BREAKDOWN;
852: PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n", beta);
853: break;
854: }
856: /*************************************************************************/
857: /* Update p and the norms. */
858: /*************************************************************************/
860: cg->beta[l_size] = beta;
861: VecAYPX(p, beta, z); /* p = z + beta p */
863: /*************************************************************************/
864: /* Compute the new direction and update the iteration. */
865: /*************************************************************************/
867: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
868: VecDot(p, z, &kappa); /* kappa = p^T Q p */
869: ++ksp->its;
871: /*************************************************************************/
872: /* Update the Lanczos tridiagonal matrix. */
873: /*************************************************************************/
875: ++l_size;
876: cg->offd[t_size] = PetscSqrtReal(beta) / fabs(alpha);
877: cg->diag[t_size] = kappa / rz + beta / alpha;
878: ++t_size;
879: }
881: /***************************************************************************/
882: /* We have the Lanczos basis, solve the tridiagonal trust-region problem */
883: /* to obtain the Lanczos direction. We know that the solution lies on */
884: /* the boundary of the trust region. We start by checking that the */
885: /* workspace allocated is large enough. */
886: /***************************************************************************/
887: /* Note that the current version only computes the solution by using the */
888: /* preconditioned direction. Need to think about how to do the */
889: /* unpreconditioned direction calculation. */
890: /***************************************************************************/
892: if (t_size > cg->alloced) {
893: if (cg->alloced) {
894: PetscFree(cg->rwork);
895: PetscFree(cg->iwork);
896: cg->alloced += cg->init_alloc;
897: } else {
898: cg->alloced = cg->init_alloc;
899: }
901: while (t_size > cg->alloced) {
902: cg->alloced += cg->init_alloc;
903: }
905: cg->alloced = PetscMin(cg->alloced, t_size);
906: PetscMalloc2(10*cg->alloced,PetscReal, &cg->rwork,5*cg->alloced,PetscBLASInt, &cg->iwork);
907: }
909: /***************************************************************************/
910: /* Set up the required vectors. */
911: /***************************************************************************/
913: t_soln = cg->rwork + 0*t_size; /* Solution */
914: t_diag = cg->rwork + 1*t_size; /* Diagonal of T */
915: t_offd = cg->rwork + 2*t_size; /* Off-diagonal of T */
916: e_valu = cg->rwork + 3*t_size; /* Eigenvalues of T */
917: e_vect = cg->rwork + 4*t_size; /* Eigenvector of T */
918: e_rwrk = cg->rwork + 5*t_size; /* Eigen workspace */
920: e_iblk = cg->iwork + 0*t_size; /* Eigen blocks */
921: e_splt = cg->iwork + 1*t_size; /* Eigen splits */
922: e_iwrk = cg->iwork + 2*t_size; /* Eigen workspace */
924: /***************************************************************************/
925: /* Compute the minimum eigenvalue of T. */
926: /***************************************************************************/
928: vl = 0.0;
929: vu = 0.0;
930: il = 1;
931: iu = 1;
933: #if defined(PETSC_MISSING_LAPACK_STEBZ)
934: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEBZ - Lapack routine is unavailable.");
935: #else
936: 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));
938: if ((0 != info) || (1 != e_valus)) {
939: /*************************************************************************/
940: /* Calculation of the minimum eigenvalue failed. Return the */
941: /* Steihaug-Toint direction. */
942: /*************************************************************************/
944: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvalue.\n");
945: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
946: return(0);
947: }
949: cg->e_min = e_valu[0];
951: /***************************************************************************/
952: /* Compute the initial value of lambda to make (T + lamba I) positive */
953: /* definite. */
954: /***************************************************************************/
956: pert = cg->init_pert;
957: if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
958: #endif
960: while (1) {
961: for (i = 0; i < t_size; ++i) {
962: t_diag[i] = cg->diag[i] + cg->lambda;
963: t_offd[i] = cg->offd[i];
964: }
966: #if defined(PETSC_MISSING_LAPACK_PTTRF)
967: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
968: #else
969: PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
971: if (0 == info) break;
973: pert += pert;
974: cg->lambda = cg->lambda * (1.0 + pert) + pert;
975: #endif
976: }
978: /***************************************************************************/
979: /* Compute the initial step and its norm. */
980: /***************************************************************************/
982: nrhs = 1;
983: nldb = t_size;
985: t_soln[0] = -cg->norm_r[0];
986: for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
988: #if defined(PETSC_MISSING_LAPACK_PTTRS)
989: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
990: #else
991: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
992: #endif
994: if (0 != info) {
995: /*************************************************************************/
996: /* Calculation of the initial step failed; return the Steihaug-Toint */
997: /* direction. */
998: /*************************************************************************/
1000: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1001: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1002: return(0);
1003: }
1005: norm_t = 0.;
1006: for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
1007: norm_t = PetscSqrtReal(norm_t);
1009: /***************************************************************************/
1010: /* Determine the case we are in. */
1011: /***************************************************************************/
1013: if (norm_t <= cg->radius) {
1014: /*************************************************************************/
1015: /* The step is within the trust region; check if we are in the hard case */
1016: /* and need to move to the boundary by following a direction of negative */
1017: /* curvature. */
1018: /*************************************************************************/
1020: if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
1021: /***********************************************************************/
1022: /* This is the hard case; compute the eigenvector associated with the */
1023: /* minimum eigenvalue and move along this direction to the boundary. */
1024: /***********************************************************************/
1026: #if defined(PETSC_MISSING_LAPACK_STEIN)
1027: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEIN - Lapack routine is unavailable.");
1028: #else
1029: 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));
1030: #endif
1032: if (0 != info) {
1033: /*********************************************************************/
1034: /* Calculation of the minimum eigenvalue failed. Return the */
1035: /* Steihaug-Toint direction. */
1036: /*********************************************************************/
1038: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvector.\n");
1039: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1040: return(0);
1041: }
1043: coef1 = 0.0;
1044: coef2 = 0.0;
1045: coef3 = -cg->radius * cg->radius;
1046: for (i = 0; i < t_size; ++i) {
1047: coef1 += e_vect[i] * e_vect[i];
1048: coef2 += e_vect[i] * t_soln[i];
1049: coef3 += t_soln[i] * t_soln[i];
1050: }
1052: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
1053: root1 = (-coef2 + coef3) / coef1;
1054: root2 = (-coef2 - coef3) / coef1;
1056: /***********************************************************************/
1057: /* Compute objective value for (t_soln + root1 * e_vect) */
1058: /***********************************************************************/
1060: for (i = 0; i < t_size; ++i) {
1061: e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
1062: }
1064: obj1 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1065: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1066: for (i = 1; i < t_size - 1; ++i) {
1067: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1068: cg->diag[i]*e_rwrk[i]+
1069: cg->offd[i+1]*e_rwrk[i+1]);
1070: }
1071: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1072: cg->diag[i]*e_rwrk[i]);
1074: /***********************************************************************/
1075: /* Compute objective value for (t_soln + root2 * e_vect) */
1076: /***********************************************************************/
1078: for (i = 0; i < t_size; ++i) {
1079: e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
1080: }
1082: obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1083: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1084: for (i = 1; i < t_size - 1; ++i) {
1085: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1086: cg->diag[i]*e_rwrk[i]+
1087: cg->offd[i+1]*e_rwrk[i+1]);
1088: }
1089: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1090: cg->diag[i]*e_rwrk[i]);
1092: /***********************************************************************/
1093: /* Choose the point with the best objective function value. */
1094: /***********************************************************************/
1096: if (obj1 <= obj2) {
1097: for (i = 0; i < t_size; ++i) {
1098: t_soln[i] += root1 * e_vect[i];
1099: }
1100: }
1101: else {
1102: for (i = 0; i < t_size; ++i) {
1103: t_soln[i] += root2 * e_vect[i];
1104: }
1105: }
1106: } else {
1107: /***********************************************************************/
1108: /* The matrix is positive definite or there was no room to move; the */
1109: /* solution is already contained in t_soln. */
1110: /***********************************************************************/
1111: }
1112: } else {
1113: /*************************************************************************/
1114: /* The step is outside the trust-region. Compute the correct value for */
1115: /* lambda by performing Newton's method. */
1116: /*************************************************************************/
1118: for (i = 0; i < max_newton_its; ++i) {
1119: /***********************************************************************/
1120: /* Check for convergence. */
1121: /***********************************************************************/
1123: if (fabs(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
1125: /***********************************************************************/
1126: /* Compute the update. */
1127: /***********************************************************************/
1129: PetscMemcpy(e_rwrk, t_soln, sizeof(PetscReal)*t_size);
1131: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1132: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1133: #else
1134: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
1135: #endif
1137: if (0 != info) {
1138: /*********************************************************************/
1139: /* Calculation of the step failed; return the Steihaug-Toint */
1140: /* direction. */
1141: /*********************************************************************/
1143: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1144: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1145: return(0);
1146: }
1148: /***********************************************************************/
1149: /* Modify lambda. */
1150: /***********************************************************************/
1152: norm_w = 0.;
1153: for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
1155: cg->lambda += (norm_t - cg->radius)/cg->radius * (norm_t * norm_t) / norm_w;
1157: /***********************************************************************/
1158: /* Factor T + lambda I */
1159: /***********************************************************************/
1161: for (j = 0; j < t_size; ++j) {
1162: t_diag[j] = cg->diag[j] + cg->lambda;
1163: t_offd[j] = cg->offd[j];
1164: }
1166: #if defined(PETSC_MISSING_LAPACK_PTTRF)
1167: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
1168: #else
1169: PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1170: #endif
1172: if (0 != info) {
1173: /*********************************************************************/
1174: /* Calculation of factorization failed; return the Steihaug-Toint */
1175: /* direction. */
1176: /*********************************************************************/
1178: PetscInfo(ksp, "KSPSolve_GLTR: factorization failed.\n");
1179: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1180: return(0);
1181: }
1183: /***********************************************************************/
1184: /* Compute the new step and its norm. */
1185: /***********************************************************************/
1187: t_soln[0] = -cg->norm_r[0];
1188: for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
1190: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1191: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1192: #else
1193: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1194: #endif
1196: if (0 != info) {
1197: /*********************************************************************/
1198: /* Calculation of the step failed; return the Steihaug-Toint */
1199: /* direction. */
1200: /*********************************************************************/
1202: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1203: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1204: return(0);
1205: }
1207: norm_t = 0.;
1208: for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1209: norm_t = PetscSqrtReal(norm_t);
1210: }
1212: /*************************************************************************/
1213: /* Check for convergence. */
1214: /*************************************************************************/
1216: if (fabs(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1217: /***********************************************************************/
1218: /* Newton method failed to converge in iteration limit. */
1219: /***********************************************************************/
1221: PetscInfo(ksp, "KSPSolve_GLTR: failed to converge.\n");
1222: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1223: return(0);
1224: }
1225: }
1227: /***************************************************************************/
1228: /* Recover the norm of the direction and objective function value. */
1229: /***************************************************************************/
1231: cg->norm_d = norm_t;
1233: cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1234: for (i = 1; i < t_size - 1; ++i) {
1235: 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]);
1236: }
1237: cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);
1239: /***************************************************************************/
1240: /* Recover the direction. */
1241: /***************************************************************************/
1243: sigma = -1;
1245: /***************************************************************************/
1246: /* Start conjugate gradient method from the beginning */
1247: /***************************************************************************/
1249: VecCopy(ksp->vec_rhs, r); /* r = -grad */
1250: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1252: /***************************************************************************/
1253: /* Accumulate Q * s */
1254: /***************************************************************************/
1256: VecCopy(z, d);
1257: VecScale(d, sigma * t_soln[0] / cg->norm_r[0]);
1259: /***************************************************************************/
1260: /* Compute the first direction. */
1261: /***************************************************************************/
1263: VecCopy(z, p); /* p = z */
1264: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1265: ++ksp->its;
1267: for (i = 0; i < l_size - 1; ++i) {
1268: /*************************************************************************/
1269: /* Update the residual and direction. */
1270: /*************************************************************************/
1272: alpha = cg->alpha[i];
1273: if (alpha >= 0.0) sigma = -sigma;
1275: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1276: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1278: /*************************************************************************/
1279: /* Accumulate Q * s */
1280: /*************************************************************************/
1282: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1284: /*************************************************************************/
1285: /* Update p. */
1286: /*************************************************************************/
1288: beta = cg->beta[i];
1289: VecAYPX(p, beta, z); /* p = z + beta p */
1290: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1291: ++ksp->its;
1292: }
1294: /***************************************************************************/
1295: /* Update the residual and direction. */
1296: /***************************************************************************/
1298: alpha = cg->alpha[i];
1299: if (alpha >= 0.0) sigma = -sigma;
1301: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1302: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1304: /***************************************************************************/
1305: /* Accumulate Q * s */
1306: /***************************************************************************/
1308: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1310: /***************************************************************************/
1311: /* Set the termination reason. */
1312: /***************************************************************************/
1314: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1315: return(0);
1316: #endif
1317: }
1321: PetscErrorCode KSPSetUp_GLTR(KSP ksp)
1322: {
1323: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1324: PetscInt max_its;
1328: /***************************************************************************/
1329: /* Determine the total maximum number of iterations. */
1330: /***************************************************************************/
1332: max_its = ksp->max_it + cg->max_lanczos_its + 1;
1334: /***************************************************************************/
1335: /* Set work vectors needed by conjugate gradient method and allocate */
1336: /* workspace for Lanczos matrix. */
1337: /***************************************************************************/
1339: KSPSetWorkVecs(ksp, 3);
1341: PetscMalloc5(max_its,PetscReal,&cg->diag,max_its,PetscReal,&cg->offd,max_its,PetscReal,&cg->alpha,max_its,PetscReal,&cg->beta,max_its,PetscReal,&cg->norm_r);
1342: PetscMemzero(cg->diag, max_its*sizeof(PetscReal));
1343: PetscMemzero(cg->offd, max_its*sizeof(PetscReal));
1344: PetscMemzero(cg->alpha, max_its*sizeof(PetscReal));
1345: PetscMemzero(cg->beta, max_its*sizeof(PetscReal));
1346: PetscMemzero(cg->norm_r, max_its*sizeof(PetscReal));
1347: PetscLogObjectMemory(ksp, 5*max_its*sizeof(PetscReal));
1348: return(0);
1349: }
1353: PetscErrorCode KSPDestroy_GLTR(KSP ksp)
1354: {
1355: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1359: /***************************************************************************/
1360: /* Free memory allocated for the data. */
1361: /***************************************************************************/
1363: PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1364: if (cg->alloced) {
1365: PetscFree2(cg->rwork,cg->iwork);
1366: }
1368: /***************************************************************************/
1369: /* Clear composed functions */
1370: /***************************************************************************/
1372: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRSetRadius_C",NULL);
1373: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetNormD_C",NULL);
1374: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetObjFcn_C",NULL);
1375: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",NULL);
1376: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",NULL);
1378: /***************************************************************************/
1379: /* Destroy KSP object. */
1380: /***************************************************************************/
1381: KSPDestroyDefault(ksp);
1382: return(0);
1383: }
1387: static PetscErrorCode KSPGLTRSetRadius_GLTR(KSP ksp, PetscReal radius)
1388: {
1389: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1392: cg->radius = radius;
1393: return(0);
1394: }
1398: static PetscErrorCode KSPGLTRGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1399: {
1400: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1403: *norm_d = cg->norm_d;
1404: return(0);
1405: }
1409: static PetscErrorCode KSPGLTRGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1410: {
1411: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1414: *o_fcn = cg->o_fcn;
1415: return(0);
1416: }
1420: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1421: {
1422: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1425: *e_min = cg->e_min;
1426: return(0);
1427: }
1431: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1432: {
1433: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1436: *lambda = cg->lambda;
1437: return(0);
1438: }
1442: PetscErrorCode KSPSetFromOptions_GLTR(KSP ksp)
1443: {
1445: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1448: PetscOptionsHead("KSP GLTR options");
1450: PetscOptionsReal("-ksp_gltr_radius", "Trust Region Radius", "KSPGLTRSetRadius", cg->radius, &cg->radius, NULL);
1452: PetscOptionsReal("-ksp_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1453: PetscOptionsReal("-ksp_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1454: PetscOptionsReal("-ksp_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);
1456: PetscOptionsInt("-ksp_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL);
1457: PetscOptionsInt("-ksp_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL);
1459: PetscOptionsEList("-ksp_gltr_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
1461: PetscOptionsTail();
1462: return(0);
1463: }
1465: /*MC
1466: KSPGLTR - Code to run conjugate gradient method subject to a constraint
1467: on the solution norm. This is used in Trust Region methods for
1468: nonlinear equations, SNESNEWTONTR
1470: Options Database Keys:
1471: . -ksp_gltr_radius <r> - Trust Region Radius
1473: Notes: This is rarely used directly
1475: Use preconditioned conjugate gradient to compute
1476: an approximate minimizer of the quadratic function
1478: q(s) = g^T * s + .5 * s^T * H * s
1480: subject to the trust region constraint
1482: || s || <= delta,
1484: where
1486: delta is the trust region radius,
1487: g is the gradient vector,
1488: H is the Hessian approximation,
1489: M is the positive definite preconditioner matrix.
1491: KSPConvergedReason may be
1492: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1493: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1494: $ other KSP converged/diverged reasons
1496: Notes:
1497: The preconditioner supplied should be symmetric and positive definite.
1499: Level: developer
1501: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGLTRSetRadius(), KSPGLTRGetNormD(), KSPGLTRGetObjFcn(), KSPGLTRGetMinEig(), KSPGLTRGetLambda()
1502: M*/
1506: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1507: {
1509: KSP_GLTR *cg;
1512: PetscNewLog(ksp, KSP_GLTR, &cg);
1513: cg->radius = 0.0;
1514: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1516: cg->init_pert = 1.0e-8;
1517: cg->eigen_tol = 1.0e-10;
1518: cg->newton_tol = 1.0e-6;
1520: cg->alloced = 0;
1521: cg->init_alloc = 1024;
1523: cg->max_lanczos_its = 20;
1524: cg->max_newton_its = 10;
1526: ksp->data = (void*) cg;
1527: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
1528: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,1);
1529: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);
1531: /***************************************************************************/
1532: /* Sets the functions that are associated with this data structure */
1533: /* (in C++ this is the same as defining virtual functions). */
1534: /***************************************************************************/
1536: ksp->ops->setup = KSPSetUp_GLTR;
1537: ksp->ops->solve = KSPSolve_GLTR;
1538: ksp->ops->destroy = KSPDestroy_GLTR;
1539: ksp->ops->setfromoptions = KSPSetFromOptions_GLTR;
1540: ksp->ops->buildsolution = KSPBuildSolutionDefault;
1541: ksp->ops->buildresidual = KSPBuildResidualDefault;
1542: ksp->ops->view = 0;
1544: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRSetRadius_C",KSPGLTRSetRadius_GLTR);
1545: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetNormD_C", KSPGLTRGetNormD_GLTR);
1546: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetObjFcn_C",KSPGLTRGetObjFcn_GLTR);
1547: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",KSPGLTRGetMinEig_GLTR);
1548: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",KSPGLTRGetLambda_GLTR);
1549: return(0);
1550: }