Actual source code: gltr.c
petsc-3.7.7 2017-09-25
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: KSPCheckDot(ksp,rr);
226: /***************************************************************************/
227: /* Check the preconditioner for numerical problems and for positive */
228: /* definiteness. The check for not-a-number and infinite values need be */
229: /* performed only once. */
230: /***************************************************************************/
232: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
233: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
234: if (PetscIsInfOrNanScalar(rz)) {
235: /*************************************************************************/
236: /* The preconditioner contains not-a-number or an infinite value. */
237: /* Return the gradient direction intersected with the trust region. */
238: /*************************************************************************/
240: ksp->reason = KSP_DIVERGED_NANORINF;
241: PetscInfo1(ksp, "KSPSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz);
243: if (cg->radius) {
244: if (r2 >= rr) {
245: alpha = 1.0;
246: cg->norm_d = PetscSqrtReal(rr);
247: } else {
248: alpha = PetscSqrtReal(r2 / rr);
249: cg->norm_d = cg->radius;
250: }
252: VecAXPY(d, alpha, r); /* d = d + alpha r */
254: /***********************************************************************/
255: /* Compute objective function. */
256: /***********************************************************************/
258: KSP_MatMult(ksp, Qmat, d, z);
259: VecAYPX(z, -0.5, ksp->vec_rhs);
260: VecDot(d, z, &cg->o_fcn);
261: cg->o_fcn = -cg->o_fcn;
262: ++ksp->its;
263: }
264: return(0);
265: }
267: if (rz < 0.0) {
268: /*************************************************************************/
269: /* The preconditioner is indefinite. Because this is the first */
270: /* and we do not have a direction yet, we use the gradient step. Note */
271: /* that we cannot use the preconditioned norm when computing the step */
272: /* because the matrix is indefinite. */
273: /*************************************************************************/
275: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
276: PetscInfo1(ksp, "KSPSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz);
278: if (cg->radius) {
279: if (r2 >= rr) {
280: alpha = 1.0;
281: cg->norm_d = PetscSqrtReal(rr);
282: } else {
283: alpha = PetscSqrtReal(r2 / rr);
284: cg->norm_d = cg->radius;
285: }
287: VecAXPY(d, alpha, r); /* d = d + alpha r */
289: /***********************************************************************/
290: /* Compute objective function. */
291: /***********************************************************************/
293: KSP_MatMult(ksp, Qmat, d, z);
294: VecAYPX(z, -0.5, ksp->vec_rhs);
295: VecDot(d, z, &cg->o_fcn);
296: cg->o_fcn = -cg->o_fcn;
297: ++ksp->its;
298: }
299: return(0);
300: }
302: /***************************************************************************/
303: /* As far as we know, the preconditioner is positive semidefinite. */
304: /* Compute and log the residual. Check convergence because this */
305: /* initializes things, but do not terminate until at least one conjugate */
306: /* gradient iteration has been performed. */
307: /***************************************************************************/
309: cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
311: switch (ksp->normtype) {
312: case KSP_NORM_PRECONDITIONED:
313: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
314: break;
316: case KSP_NORM_UNPRECONDITIONED:
317: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
318: break;
320: case KSP_NORM_NATURAL:
321: norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
322: break;
324: default:
325: norm_r = 0.0;
326: break;
327: }
329: KSPLogResidualHistory(ksp, norm_r);
330: KSPMonitor(ksp, ksp->its, norm_r);
331: ksp->rnorm = norm_r;
333: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
335: /***************************************************************************/
336: /* Compute the first direction and update the iteration. */
337: /***************************************************************************/
339: VecCopy(z, p); /* p = z */
340: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
341: ++ksp->its;
343: /***************************************************************************/
344: /* Check the matrix for numerical problems. */
345: /***************************************************************************/
347: VecDot(p, z, &kappa); /* kappa = p^T Q p */
348: if (PetscIsInfOrNanScalar(kappa)) {
349: /*************************************************************************/
350: /* The matrix produced not-a-number or an infinite value. In this case, */
351: /* we must stop and use the gradient direction. This condition need */
352: /* only be checked once. */
353: /*************************************************************************/
355: ksp->reason = KSP_DIVERGED_NANORINF;
356: PetscInfo1(ksp, "KSPSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa);
358: if (cg->radius) {
359: if (r2 >= rr) {
360: alpha = 1.0;
361: cg->norm_d = PetscSqrtReal(rr);
362: } else {
363: alpha = PetscSqrtReal(r2 / rr);
364: cg->norm_d = cg->radius;
365: }
367: VecAXPY(d, alpha, r); /* d = d + alpha r */
369: /***********************************************************************/
370: /* Compute objective function. */
371: /***********************************************************************/
373: KSP_MatMult(ksp, Qmat, d, z);
374: VecAYPX(z, -0.5, ksp->vec_rhs);
375: VecDot(d, z, &cg->o_fcn);
376: cg->o_fcn = -cg->o_fcn;
377: ++ksp->its;
378: }
379: return(0);
380: }
382: /***************************************************************************/
383: /* Initialize variables for calculating the norm of the direction and for */
384: /* the Lanczos tridiagonal matrix. Note that we track the diagonal value */
385: /* of the Cholesky factorization of the Lanczos matrix in order to */
386: /* determine when negative curvature is encountered. */
387: /***************************************************************************/
389: dMp = 0.0;
390: norm_d = 0.0;
391: switch (cg->dtype) {
392: case GLTR_PRECONDITIONED_DIRECTION:
393: norm_p = rz;
394: break;
396: default:
397: VecDot(p, p, &norm_p);
398: break;
399: }
401: cg->diag[t_size] = kappa / rz;
402: cg->offd[t_size] = 0.0;
403: ++t_size;
405: piv = 1.0;
407: /***************************************************************************/
408: /* Check for breakdown of the conjugate gradient method; this occurs when */
409: /* kappa is zero. */
410: /***************************************************************************/
412: if (PetscAbsReal(kappa) <= 0.0) {
413: /*************************************************************************/
414: /* The curvature is zero. In this case, we must stop and use follow */
415: /* the direction of negative curvature since the Lanczos matrix is zero. */
416: /*************************************************************************/
418: ksp->reason = KSP_DIVERGED_BREAKDOWN;
419: PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: kappa=%g\n", (double)kappa);
421: if (cg->radius && norm_p > 0.0) {
422: /***********************************************************************/
423: /* Follow direction of negative curvature to the boundary of the */
424: /* trust region. */
425: /***********************************************************************/
427: step = PetscSqrtReal(r2 / norm_p);
428: cg->norm_d = cg->radius;
430: VecAXPY(d, step, p); /* d = d + step p */
432: /***********************************************************************/
433: /* Update objective function. */
434: /***********************************************************************/
436: cg->o_fcn += step * (0.5 * step * kappa - rz);
437: } else if (cg->radius) {
438: /***********************************************************************/
439: /* The norm of the preconditioned direction is zero; use the gradient */
440: /* step. */
441: /***********************************************************************/
443: if (r2 >= rr) {
444: alpha = 1.0;
445: cg->norm_d = PetscSqrtReal(rr);
446: } else {
447: alpha = PetscSqrtReal(r2 / rr);
448: cg->norm_d = cg->radius;
449: }
451: VecAXPY(d, alpha, r); /* d = d + alpha r */
453: /***********************************************************************/
454: /* Compute objective function. */
455: /***********************************************************************/
457: KSP_MatMult(ksp, Qmat, d, z);
458: VecAYPX(z, -0.5, ksp->vec_rhs);
459: VecDot(d, z, &cg->o_fcn);
460: cg->o_fcn = -cg->o_fcn;
461: ++ksp->its;
462: }
463: return(0);
464: }
466: /***************************************************************************/
467: /* Begin the first part of the GLTR algorithm which runs the conjugate */
468: /* gradient method until either the problem is solved, we encounter the */
469: /* boundary of the trust region, or the conjugate gradient method breaks */
470: /* down. */
471: /***************************************************************************/
473: while (1) {
474: /*************************************************************************/
475: /* Know that kappa is nonzero, because we have not broken down, so we */
476: /* can compute the steplength. */
477: /*************************************************************************/
479: alpha = rz / kappa;
480: cg->alpha[l_size] = alpha;
482: /*************************************************************************/
483: /* Compute the diagonal value of the Cholesky factorization of the */
484: /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
485: /* This indicates a direction of negative curvature. */
486: /*************************************************************************/
488: piv = cg->diag[l_size] - cg->offd[l_size]*cg->offd[l_size] / piv;
489: if (piv <= 0.0) {
490: /***********************************************************************/
491: /* In this case, the matrix is indefinite and we have encountered */
492: /* a direction of negative curvature. Follow the direction to the */
493: /* boundary of the trust region. */
494: /***********************************************************************/
496: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
497: PetscInfo1(ksp, "KSPSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa);
499: if (cg->radius && norm_p > 0.0) {
500: /*********************************************************************/
501: /* Follow direction of negative curvature to boundary. */
502: /*********************************************************************/
504: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
505: cg->norm_d = cg->radius;
507: VecAXPY(d, step, p); /* d = d + step p */
509: /*********************************************************************/
510: /* Update objective function. */
511: /*********************************************************************/
513: cg->o_fcn += step * (0.5 * step * kappa - rz);
514: } else if (cg->radius) {
515: /*********************************************************************/
516: /* The norm of the direction is zero; there is nothing to follow. */
517: /*********************************************************************/
518: }
519: break;
520: }
522: /*************************************************************************/
523: /* Compute the steplength and check for intersection with the trust */
524: /* region. */
525: /*************************************************************************/
527: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
528: if (cg->radius && norm_dp1 >= r2) {
529: /***********************************************************************/
530: /* In this case, the matrix is positive definite as far as we know. */
531: /* However, the full step goes beyond the trust region. */
532: /***********************************************************************/
534: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
535: PetscInfo1(ksp, "KSPSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius);
537: if (norm_p > 0.0) {
538: /*********************************************************************/
539: /* Follow the direction to the boundary of the trust region. */
540: /*********************************************************************/
542: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
543: cg->norm_d = cg->radius;
545: VecAXPY(d, step, p); /* d = d + step p */
547: /*********************************************************************/
548: /* Update objective function. */
549: /*********************************************************************/
551: cg->o_fcn += step * (0.5 * step * kappa - rz);
552: } else {
553: /*********************************************************************/
554: /* The norm of the direction is zero; there is nothing to follow. */
555: /*********************************************************************/
556: }
557: break;
558: }
560: /*************************************************************************/
561: /* Now we can update the direction and residual. */
562: /*************************************************************************/
564: VecAXPY(d, alpha, p); /* d = d + alpha p */
565: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
566: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
568: switch (cg->dtype) {
569: case GLTR_PRECONDITIONED_DIRECTION:
570: norm_d = norm_dp1;
571: break;
573: default:
574: VecDot(d, d, &norm_d);
575: break;
576: }
577: cg->norm_d = PetscSqrtReal(norm_d);
579: /*************************************************************************/
580: /* Update objective function. */
581: /*************************************************************************/
583: cg->o_fcn -= 0.5 * alpha * rz;
585: /*************************************************************************/
586: /* Check that the preconditioner appears positive semidefinite. */
587: /*************************************************************************/
589: rzm1 = rz;
590: VecDot(r, z, &rz); /* rz = r^T z */
591: if (rz < 0.0) {
592: /***********************************************************************/
593: /* The preconditioner is indefinite. */
594: /***********************************************************************/
596: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
597: PetscInfo1(ksp, "KSPSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz);
598: break;
599: }
601: /*************************************************************************/
602: /* As far as we know, the preconditioner is positive semidefinite. */
603: /* Compute the residual and check for convergence. */
604: /*************************************************************************/
606: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
608: switch (ksp->normtype) {
609: case KSP_NORM_PRECONDITIONED:
610: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
611: break;
613: case KSP_NORM_UNPRECONDITIONED:
614: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
615: break;
617: case KSP_NORM_NATURAL:
618: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
619: break;
621: default:
622: norm_r = 0.0;
623: break;
624: }
626: KSPLogResidualHistory(ksp, norm_r);
627: KSPMonitor(ksp, ksp->its, norm_r);
628: ksp->rnorm = norm_r;
630: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
631: if (ksp->reason) {
632: /***********************************************************************/
633: /* The method has converged. */
634: /***********************************************************************/
636: PetscInfo2(ksp, "KSPSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
637: break;
638: }
640: /*************************************************************************/
641: /* We have not converged yet. Check for breakdown. */
642: /*************************************************************************/
644: beta = rz / rzm1;
645: if (PetscAbsReal(beta) <= 0.0) {
646: /***********************************************************************/
647: /* Conjugate gradients has broken down. */
648: /***********************************************************************/
650: ksp->reason = KSP_DIVERGED_BREAKDOWN;
651: PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n", (double)beta);
652: break;
653: }
655: /*************************************************************************/
656: /* Check iteration limit. */
657: /*************************************************************************/
659: if (ksp->its >= max_cg_its) {
660: ksp->reason = KSP_DIVERGED_ITS;
661: PetscInfo1(ksp, "KSPSolve_GLTR: iterlim: its=%D\n", ksp->its);
662: break;
663: }
665: /*************************************************************************/
666: /* Update p and the norms. */
667: /*************************************************************************/
669: cg->beta[l_size] = beta;
670: VecAYPX(p, beta, z); /* p = z + beta p */
672: switch (cg->dtype) {
673: case GLTR_PRECONDITIONED_DIRECTION:
674: dMp = beta*(dMp + alpha*norm_p);
675: norm_p = beta*(rzm1 + beta*norm_p);
676: break;
678: default:
679: VecDot(d, p, &dMp);
680: VecDot(p, p, &norm_p);
681: break;
682: }
684: /*************************************************************************/
685: /* Compute the new direction and update the iteration. */
686: /*************************************************************************/
688: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
689: VecDot(p, z, &kappa); /* kappa = p^T Q p */
690: ++ksp->its;
692: /*************************************************************************/
693: /* Update the Lanczos tridiagonal matrix. */
694: /*************************************************************************/
696: ++l_size;
697: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
698: cg->diag[t_size] = kappa / rz + beta / alpha;
699: ++t_size;
701: /*************************************************************************/
702: /* Check for breakdown of the conjugate gradient method; this occurs */
703: /* when kappa is zero. */
704: /*************************************************************************/
706: if (PetscAbsReal(kappa) <= 0.0) {
707: /***********************************************************************/
708: /* The method breaks down; move along the direction as if the matrix */
709: /* were indefinite. */
710: /***********************************************************************/
712: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
713: PetscInfo1(ksp, "KSPSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa);
715: if (cg->radius && norm_p > 0.0) {
716: /*********************************************************************/
717: /* Follow direction to boundary. */
718: /*********************************************************************/
720: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
721: cg->norm_d = cg->radius;
723: VecAXPY(d, step, p); /* d = d + step p */
725: /*********************************************************************/
726: /* Update objective function. */
727: /*********************************************************************/
729: cg->o_fcn += step * (0.5 * step * kappa - rz);
730: } else if (cg->radius) {
731: /*********************************************************************/
732: /* The norm of the direction is zero; there is nothing to follow. */
733: /*********************************************************************/
734: }
735: break;
736: }
737: }
739: /***************************************************************************/
740: /* Check to see if we need to continue with the Lanczos method. */
741: /***************************************************************************/
743: if (!cg->radius) {
744: /*************************************************************************/
745: /* There is no radius. Therefore, we cannot move along the boundary. */
746: /*************************************************************************/
747: return(0);
748: }
750: if (KSP_CONVERGED_CG_NEG_CURVE != ksp->reason) {
751: /*************************************************************************/
752: /* The method either converged to an interior point, hit the boundary of */
753: /* the trust-region without encountering a direction of negative */
754: /* curvature or the iteration limit was reached. */
755: /*************************************************************************/
756: return(0);
757: }
759: /***************************************************************************/
760: /* Switch to contructing the Lanczos basis by way of the conjugate */
761: /* directions. */
762: /***************************************************************************/
764: for (i = 0; i < max_lanczos_its; ++i) {
765: /*************************************************************************/
766: /* Check for breakdown of the conjugate gradient method; this occurs */
767: /* when kappa is zero. */
768: /*************************************************************************/
770: if (PetscAbsReal(kappa) <= 0.0) {
771: ksp->reason = KSP_DIVERGED_BREAKDOWN;
772: PetscInfo1(ksp, "KSPSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa);
773: break;
774: }
776: /*************************************************************************/
777: /* Update the direction and residual. */
778: /*************************************************************************/
780: alpha = rz / kappa;
781: cg->alpha[l_size] = alpha;
783: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
784: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
786: /*************************************************************************/
787: /* Check that the preconditioner appears positive semidefinite. */
788: /*************************************************************************/
790: rzm1 = rz;
791: VecDot(r, z, &rz); /* rz = r^T z */
792: if (rz < 0.0) {
793: /***********************************************************************/
794: /* The preconditioner is indefinite. */
795: /***********************************************************************/
797: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
798: PetscInfo1(ksp, "KSPSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz);
799: break;
800: }
802: /*************************************************************************/
803: /* As far as we know, the preconditioner is positive definite. Compute */
804: /* the residual. Do NOT check for convergence. */
805: /*************************************************************************/
807: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
809: switch (ksp->normtype) {
810: case KSP_NORM_PRECONDITIONED:
811: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
812: break;
814: case KSP_NORM_UNPRECONDITIONED:
815: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
816: break;
818: case KSP_NORM_NATURAL:
819: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
820: break;
822: default:
823: norm_r = 0.0;
824: break;
825: }
827: KSPLogResidualHistory(ksp, norm_r);
828: KSPMonitor(ksp, ksp->its, norm_r);
829: ksp->rnorm = norm_r;
831: /*************************************************************************/
832: /* Check for breakdown. */
833: /*************************************************************************/
835: beta = rz / rzm1;
836: if (PetscAbsReal(beta) <= 0.0) {
837: /***********************************************************************/
838: /* Conjugate gradients has broken down. */
839: /***********************************************************************/
841: ksp->reason = KSP_DIVERGED_BREAKDOWN;
842: PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n",(double) beta);
843: break;
844: }
846: /*************************************************************************/
847: /* Update p and the norms. */
848: /*************************************************************************/
850: cg->beta[l_size] = beta;
851: VecAYPX(p, beta, z); /* p = z + beta p */
853: /*************************************************************************/
854: /* Compute the new direction and update the iteration. */
855: /*************************************************************************/
857: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
858: VecDot(p, z, &kappa); /* kappa = p^T Q p */
859: ++ksp->its;
861: /*************************************************************************/
862: /* Update the Lanczos tridiagonal matrix. */
863: /*************************************************************************/
865: ++l_size;
866: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
867: cg->diag[t_size] = kappa / rz + beta / alpha;
868: ++t_size;
869: }
871: /***************************************************************************/
872: /* We have the Lanczos basis, solve the tridiagonal trust-region problem */
873: /* to obtain the Lanczos direction. We know that the solution lies on */
874: /* the boundary of the trust region. We start by checking that the */
875: /* workspace allocated is large enough. */
876: /***************************************************************************/
877: /* Note that the current version only computes the solution by using the */
878: /* preconditioned direction. Need to think about how to do the */
879: /* unpreconditioned direction calculation. */
880: /***************************************************************************/
882: if (t_size > cg->alloced) {
883: if (cg->alloced) {
884: PetscFree(cg->rwork);
885: PetscFree(cg->iwork);
886: cg->alloced += cg->init_alloc;
887: } else {
888: cg->alloced = cg->init_alloc;
889: }
891: while (t_size > cg->alloced) {
892: cg->alloced += cg->init_alloc;
893: }
895: cg->alloced = PetscMin(cg->alloced, t_size);
896: PetscMalloc2(10*cg->alloced, &cg->rwork,5*cg->alloced, &cg->iwork);
897: }
899: /***************************************************************************/
900: /* Set up the required vectors. */
901: /***************************************************************************/
903: t_soln = cg->rwork + 0*t_size; /* Solution */
904: t_diag = cg->rwork + 1*t_size; /* Diagonal of T */
905: t_offd = cg->rwork + 2*t_size; /* Off-diagonal of T */
906: e_valu = cg->rwork + 3*t_size; /* Eigenvalues of T */
907: e_vect = cg->rwork + 4*t_size; /* Eigenvector of T */
908: e_rwrk = cg->rwork + 5*t_size; /* Eigen workspace */
910: e_iblk = cg->iwork + 0*t_size; /* Eigen blocks */
911: e_splt = cg->iwork + 1*t_size; /* Eigen splits */
912: e_iwrk = cg->iwork + 2*t_size; /* Eigen workspace */
914: /***************************************************************************/
915: /* Compute the minimum eigenvalue of T. */
916: /***************************************************************************/
918: vl = 0.0;
919: vu = 0.0;
920: il = 1;
921: iu = 1;
923: #if defined(PETSC_MISSING_LAPACK_STEBZ)
924: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEBZ - Lapack routine is unavailable.");
925: #else
926: 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));
928: if ((0 != info) || (1 != e_valus)) {
929: /*************************************************************************/
930: /* Calculation of the minimum eigenvalue failed. Return the */
931: /* Steihaug-Toint direction. */
932: /*************************************************************************/
934: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvalue.\n");
935: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
936: return(0);
937: }
939: cg->e_min = e_valu[0];
941: /***************************************************************************/
942: /* Compute the initial value of lambda to make (T + lamba I) positive */
943: /* definite. */
944: /***************************************************************************/
946: pert = cg->init_pert;
947: if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
948: #endif
950: while (1) {
951: for (i = 0; i < t_size; ++i) {
952: t_diag[i] = cg->diag[i] + cg->lambda;
953: t_offd[i] = cg->offd[i];
954: }
956: #if defined(PETSC_MISSING_LAPACK_PTTRF)
957: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
958: #else
959: PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
961: if (0 == info) break;
963: pert += pert;
964: cg->lambda = cg->lambda * (1.0 + pert) + pert;
965: #endif
966: }
968: /***************************************************************************/
969: /* Compute the initial step and its norm. */
970: /***************************************************************************/
972: nrhs = 1;
973: nldb = t_size;
975: t_soln[0] = -cg->norm_r[0];
976: for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
978: #if defined(PETSC_MISSING_LAPACK_PTTRS)
979: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
980: #else
981: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
982: #endif
984: if (0 != info) {
985: /*************************************************************************/
986: /* Calculation of the initial step failed; return the Steihaug-Toint */
987: /* direction. */
988: /*************************************************************************/
990: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
991: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
992: return(0);
993: }
995: norm_t = 0.;
996: for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
997: norm_t = PetscSqrtReal(norm_t);
999: /***************************************************************************/
1000: /* Determine the case we are in. */
1001: /***************************************************************************/
1003: if (norm_t <= cg->radius) {
1004: /*************************************************************************/
1005: /* The step is within the trust region; check if we are in the hard case */
1006: /* and need to move to the boundary by following a direction of negative */
1007: /* curvature. */
1008: /*************************************************************************/
1010: if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
1011: /***********************************************************************/
1012: /* This is the hard case; compute the eigenvector associated with the */
1013: /* minimum eigenvalue and move along this direction to the boundary. */
1014: /***********************************************************************/
1016: #if defined(PETSC_MISSING_LAPACK_STEIN)
1017: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEIN - Lapack routine is unavailable.");
1018: #else
1019: 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));
1020: #endif
1022: if (0 != info) {
1023: /*********************************************************************/
1024: /* Calculation of the minimum eigenvalue failed. Return the */
1025: /* Steihaug-Toint direction. */
1026: /*********************************************************************/
1028: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvector.\n");
1029: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1030: return(0);
1031: }
1033: coef1 = 0.0;
1034: coef2 = 0.0;
1035: coef3 = -cg->radius * cg->radius;
1036: for (i = 0; i < t_size; ++i) {
1037: coef1 += e_vect[i] * e_vect[i];
1038: coef2 += e_vect[i] * t_soln[i];
1039: coef3 += t_soln[i] * t_soln[i];
1040: }
1042: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
1043: root1 = (-coef2 + coef3) / coef1;
1044: root2 = (-coef2 - coef3) / coef1;
1046: /***********************************************************************/
1047: /* Compute objective value for (t_soln + root1 * e_vect) */
1048: /***********************************************************************/
1050: for (i = 0; i < t_size; ++i) {
1051: e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
1052: }
1054: obj1 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1055: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1056: for (i = 1; i < t_size - 1; ++i) {
1057: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1058: cg->diag[i]*e_rwrk[i]+
1059: cg->offd[i+1]*e_rwrk[i+1]);
1060: }
1061: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1062: cg->diag[i]*e_rwrk[i]);
1064: /***********************************************************************/
1065: /* Compute objective value for (t_soln + root2 * e_vect) */
1066: /***********************************************************************/
1068: for (i = 0; i < t_size; ++i) {
1069: e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
1070: }
1072: obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1073: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1074: for (i = 1; i < t_size - 1; ++i) {
1075: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1076: cg->diag[i]*e_rwrk[i]+
1077: cg->offd[i+1]*e_rwrk[i+1]);
1078: }
1079: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1080: cg->diag[i]*e_rwrk[i]);
1082: /***********************************************************************/
1083: /* Choose the point with the best objective function value. */
1084: /***********************************************************************/
1086: if (obj1 <= obj2) {
1087: for (i = 0; i < t_size; ++i) {
1088: t_soln[i] += root1 * e_vect[i];
1089: }
1090: }
1091: else {
1092: for (i = 0; i < t_size; ++i) {
1093: t_soln[i] += root2 * e_vect[i];
1094: }
1095: }
1096: } else {
1097: /***********************************************************************/
1098: /* The matrix is positive definite or there was no room to move; the */
1099: /* solution is already contained in t_soln. */
1100: /***********************************************************************/
1101: }
1102: } else {
1103: /*************************************************************************/
1104: /* The step is outside the trust-region. Compute the correct value for */
1105: /* lambda by performing Newton's method. */
1106: /*************************************************************************/
1108: for (i = 0; i < max_newton_its; ++i) {
1109: /***********************************************************************/
1110: /* Check for convergence. */
1111: /***********************************************************************/
1113: if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
1115: /***********************************************************************/
1116: /* Compute the update. */
1117: /***********************************************************************/
1119: PetscMemcpy(e_rwrk, t_soln, sizeof(PetscReal)*t_size);
1121: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1122: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1123: #else
1124: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
1125: #endif
1127: if (0 != info) {
1128: /*********************************************************************/
1129: /* Calculation of the step failed; return the Steihaug-Toint */
1130: /* direction. */
1131: /*********************************************************************/
1133: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1134: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1135: return(0);
1136: }
1138: /***********************************************************************/
1139: /* Modify lambda. */
1140: /***********************************************************************/
1142: norm_w = 0.;
1143: for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
1145: cg->lambda += (norm_t - cg->radius)/cg->radius * (norm_t * norm_t) / norm_w;
1147: /***********************************************************************/
1148: /* Factor T + lambda I */
1149: /***********************************************************************/
1151: for (j = 0; j < t_size; ++j) {
1152: t_diag[j] = cg->diag[j] + cg->lambda;
1153: t_offd[j] = cg->offd[j];
1154: }
1156: #if defined(PETSC_MISSING_LAPACK_PTTRF)
1157: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
1158: #else
1159: PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1160: #endif
1162: if (0 != info) {
1163: /*********************************************************************/
1164: /* Calculation of factorization failed; return the Steihaug-Toint */
1165: /* direction. */
1166: /*********************************************************************/
1168: PetscInfo(ksp, "KSPSolve_GLTR: factorization failed.\n");
1169: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1170: return(0);
1171: }
1173: /***********************************************************************/
1174: /* Compute the new step and its norm. */
1175: /***********************************************************************/
1177: t_soln[0] = -cg->norm_r[0];
1178: for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
1180: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1181: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1182: #else
1183: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1184: #endif
1186: if (0 != info) {
1187: /*********************************************************************/
1188: /* Calculation of the step failed; return the Steihaug-Toint */
1189: /* direction. */
1190: /*********************************************************************/
1192: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1193: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1194: return(0);
1195: }
1197: norm_t = 0.;
1198: for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1199: norm_t = PetscSqrtReal(norm_t);
1200: }
1202: /*************************************************************************/
1203: /* Check for convergence. */
1204: /*************************************************************************/
1206: if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1207: /***********************************************************************/
1208: /* Newton method failed to converge in iteration limit. */
1209: /***********************************************************************/
1211: PetscInfo(ksp, "KSPSolve_GLTR: failed to converge.\n");
1212: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1213: return(0);
1214: }
1215: }
1217: /***************************************************************************/
1218: /* Recover the norm of the direction and objective function value. */
1219: /***************************************************************************/
1221: cg->norm_d = norm_t;
1223: cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1224: for (i = 1; i < t_size - 1; ++i) {
1225: 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]);
1226: }
1227: cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);
1229: /***************************************************************************/
1230: /* Recover the direction. */
1231: /***************************************************************************/
1233: sigma = -1;
1235: /***************************************************************************/
1236: /* Start conjugate gradient method from the beginning */
1237: /***************************************************************************/
1239: VecCopy(ksp->vec_rhs, r); /* r = -grad */
1240: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1242: /***************************************************************************/
1243: /* Accumulate Q * s */
1244: /***************************************************************************/
1246: VecCopy(z, d);
1247: VecScale(d, sigma * t_soln[0] / cg->norm_r[0]);
1249: /***************************************************************************/
1250: /* Compute the first direction. */
1251: /***************************************************************************/
1253: VecCopy(z, p); /* p = z */
1254: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1255: ++ksp->its;
1257: for (i = 0; i < l_size - 1; ++i) {
1258: /*************************************************************************/
1259: /* Update the residual and direction. */
1260: /*************************************************************************/
1262: alpha = cg->alpha[i];
1263: if (alpha >= 0.0) sigma = -sigma;
1265: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1266: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1268: /*************************************************************************/
1269: /* Accumulate Q * s */
1270: /*************************************************************************/
1272: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1274: /*************************************************************************/
1275: /* Update p. */
1276: /*************************************************************************/
1278: beta = cg->beta[i];
1279: VecAYPX(p, beta, z); /* p = z + beta p */
1280: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1281: ++ksp->its;
1282: }
1284: /***************************************************************************/
1285: /* Update the residual and direction. */
1286: /***************************************************************************/
1288: alpha = cg->alpha[i];
1289: if (alpha >= 0.0) sigma = -sigma;
1291: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1292: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1294: /***************************************************************************/
1295: /* Accumulate Q * s */
1296: /***************************************************************************/
1298: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1300: /***************************************************************************/
1301: /* Set the termination reason. */
1302: /***************************************************************************/
1304: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1305: return(0);
1306: #endif
1307: }
1311: PetscErrorCode KSPSetUp_GLTR(KSP ksp)
1312: {
1313: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1314: PetscInt max_its;
1318: /***************************************************************************/
1319: /* Determine the total maximum number of iterations. */
1320: /***************************************************************************/
1322: max_its = ksp->max_it + cg->max_lanczos_its + 1;
1324: /***************************************************************************/
1325: /* Set work vectors needed by conjugate gradient method and allocate */
1326: /* workspace for Lanczos matrix. */
1327: /***************************************************************************/
1329: KSPSetWorkVecs(ksp, 3);
1331: PetscCalloc5(max_its,&cg->diag,max_its,&cg->offd,max_its,&cg->alpha,max_its,&cg->beta,max_its,&cg->norm_r);
1332: PetscLogObjectMemory((PetscObject)ksp, 5*max_its*sizeof(PetscReal));
1333: return(0);
1334: }
1338: PetscErrorCode KSPDestroy_GLTR(KSP ksp)
1339: {
1340: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1344: /***************************************************************************/
1345: /* Free memory allocated for the data. */
1346: /***************************************************************************/
1348: PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1349: if (cg->alloced) {
1350: PetscFree2(cg->rwork,cg->iwork);
1351: }
1353: /***************************************************************************/
1354: /* Clear composed functions */
1355: /***************************************************************************/
1357: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRSetRadius_C",NULL);
1358: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetNormD_C",NULL);
1359: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetObjFcn_C",NULL);
1360: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",NULL);
1361: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",NULL);
1363: /***************************************************************************/
1364: /* Destroy KSP object. */
1365: /***************************************************************************/
1366: KSPDestroyDefault(ksp);
1367: return(0);
1368: }
1372: static PetscErrorCode KSPGLTRSetRadius_GLTR(KSP ksp, PetscReal radius)
1373: {
1374: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1377: cg->radius = radius;
1378: return(0);
1379: }
1383: static PetscErrorCode KSPGLTRGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1384: {
1385: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1388: *norm_d = cg->norm_d;
1389: return(0);
1390: }
1394: static PetscErrorCode KSPGLTRGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1395: {
1396: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1399: *o_fcn = cg->o_fcn;
1400: return(0);
1401: }
1405: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1406: {
1407: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1410: *e_min = cg->e_min;
1411: return(0);
1412: }
1416: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1417: {
1418: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1421: *lambda = cg->lambda;
1422: return(0);
1423: }
1427: PetscErrorCode KSPSetFromOptions_GLTR(PetscOptionItems *PetscOptionsObject,KSP ksp)
1428: {
1430: KSP_GLTR *cg = (KSP_GLTR*)ksp->data;
1433: PetscOptionsHead(PetscOptionsObject,"KSP GLTR options");
1435: PetscOptionsReal("-ksp_gltr_radius", "Trust Region Radius", "KSPGLTRSetRadius", cg->radius, &cg->radius, NULL);
1437: PetscOptionsReal("-ksp_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1438: PetscOptionsReal("-ksp_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1439: PetscOptionsReal("-ksp_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);
1441: PetscOptionsInt("-ksp_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL);
1442: PetscOptionsInt("-ksp_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL);
1444: PetscOptionsEList("-ksp_gltr_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
1446: PetscOptionsTail();
1447: return(0);
1448: }
1450: /*MC
1451: KSPGLTR - Code to run conjugate gradient method subject to a constraint
1452: on the solution norm. This is used in Trust Region methods for
1453: nonlinear equations, SNESNEWTONTR
1455: Options Database Keys:
1456: . -ksp_gltr_radius <r> - Trust Region Radius
1458: Notes: This is rarely used directly
1460: Use preconditioned conjugate gradient to compute
1461: an approximate minimizer of the quadratic function
1463: q(s) = g^T * s + .5 * s^T * H * s
1465: subject to the trust region constraint
1467: || s || <= delta,
1469: where
1471: delta is the trust region radius,
1472: g is the gradient vector,
1473: H is the Hessian approximation,
1474: M is the positive definite preconditioner matrix.
1476: KSPConvergedReason may be
1477: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1478: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1479: $ other KSP converged/diverged reasons
1481: Notes:
1482: The preconditioner supplied should be symmetric and positive definite.
1484: Level: developer
1486: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGLTRSetRadius(), KSPGLTRGetNormD(), KSPGLTRGetObjFcn(), KSPGLTRGetMinEig(), KSPGLTRGetLambda()
1487: M*/
1491: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1492: {
1494: KSP_GLTR *cg;
1497: PetscNewLog(ksp,&cg);
1498: cg->radius = 0.0;
1499: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1501: cg->init_pert = 1.0e-8;
1502: cg->eigen_tol = 1.0e-10;
1503: cg->newton_tol = 1.0e-6;
1505: cg->alloced = 0;
1506: cg->init_alloc = 1024;
1508: cg->max_lanczos_its = 20;
1509: cg->max_newton_its = 10;
1511: ksp->data = (void*) cg;
1512: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
1513: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
1514: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
1516: /***************************************************************************/
1517: /* Sets the functions that are associated with this data structure */
1518: /* (in C++ this is the same as defining virtual functions). */
1519: /***************************************************************************/
1521: ksp->ops->setup = KSPSetUp_GLTR;
1522: ksp->ops->solve = KSPSolve_GLTR;
1523: ksp->ops->destroy = KSPDestroy_GLTR;
1524: ksp->ops->setfromoptions = KSPSetFromOptions_GLTR;
1525: ksp->ops->buildsolution = KSPBuildSolutionDefault;
1526: ksp->ops->buildresidual = KSPBuildResidualDefault;
1527: ksp->ops->view = 0;
1529: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRSetRadius_C",KSPGLTRSetRadius_GLTR);
1530: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetNormD_C", KSPGLTRGetNormD_GLTR);
1531: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetObjFcn_C",KSPGLTRGetObjFcn_GLTR);
1532: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",KSPGLTRGetMinEig_GLTR);
1533: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",KSPGLTRGetLambda_GLTR);
1534: return(0);
1535: }