Actual source code: gltr.c
petsc-3.3-p7 2013-05-11
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(((PetscObject)ksp)->comm,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: #ifdef PETSC_USE_COMPLEX
147: SETERRQ(((PetscObject)ksp)->comm,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(((PetscObject)ksp)->comm,PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
182: if (cg->radius < 0.0) SETERRQ(((PetscObject)ksp)->comm,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_NAN;
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_NAN;
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: }
258: else {
259: alpha = PetscSqrtReal(r2 / rr);
260: cg->norm_d = cg->radius;
261: }
263: VecAXPY(d, alpha, r); /* d = d + alpha r */
264:
265: /***********************************************************************/
266: /* Compute objective function. */
267: /***********************************************************************/
269: KSP_MatMult(ksp, Qmat, d, z);
270: VecAYPX(z, -0.5, ksp->vec_rhs);
271: VecDot(d, z, &cg->o_fcn);
272: cg->o_fcn = -cg->o_fcn;
273: ++ksp->its;
274: }
275: return(0);
276: }
278: if (rz < 0.0) {
279: /*************************************************************************/
280: /* The preconditioner is indefinite. Because this is the first */
281: /* and we do not have a direction yet, we use the gradient step. Note */
282: /* that we cannot use the preconditioned norm when computing the step */
283: /* because the matrix is indefinite. */
284: /*************************************************************************/
286: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
287: PetscInfo1(ksp, "KSPSolve_GLTR: indefinite preconditioner: rz=%g\n", rz);
289: if (cg->radius) {
290: if (r2 >= rr) {
291: alpha = 1.0;
292: cg->norm_d = PetscSqrtReal(rr);
293: }
294: else {
295: alpha = PetscSqrtReal(r2 / rr);
296: cg->norm_d = cg->radius;
297: }
299: VecAXPY(d, alpha, r); /* d = d + alpha r */
301: /***********************************************************************/
302: /* Compute objective function. */
303: /***********************************************************************/
305: KSP_MatMult(ksp, Qmat, d, z);
306: VecAYPX(z, -0.5, ksp->vec_rhs);
307: VecDot(d, z, &cg->o_fcn);
308: cg->o_fcn = -cg->o_fcn;
309: ++ksp->its;
310: }
311: return(0);
312: }
314: /***************************************************************************/
315: /* As far as we know, the preconditioner is positive semidefinite. */
316: /* Compute and log the residual. Check convergence because this */
317: /* initializes things, but do not terminate until at least one conjugate */
318: /* gradient iteration has been performed. */
319: /***************************************************************************/
321: cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
323: switch(ksp->normtype) {
324: case KSP_NORM_PRECONDITIONED:
325: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
326: break;
328: case KSP_NORM_UNPRECONDITIONED:
329: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
330: break;
332: case KSP_NORM_NATURAL:
333: norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
334: break;
336: default:
337: norm_r = 0.0;
338: break;
339: }
341: KSPLogResidualHistory(ksp, norm_r);
342: KSPMonitor(ksp, ksp->its, norm_r);
343: ksp->rnorm = norm_r;
345: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
347: /***************************************************************************/
348: /* Compute the first direction and update the iteration. */
349: /***************************************************************************/
351: VecCopy(z, p); /* p = z */
352: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
353: ++ksp->its;
355: /***************************************************************************/
356: /* Check the matrix for numerical problems. */
357: /***************************************************************************/
359: VecDot(p, z, &kappa); /* kappa = p^T Q p */
360: if (PetscIsInfOrNanScalar(kappa)) {
361: /*************************************************************************/
362: /* The matrix produced not-a-number or an infinite value. In this case, */
363: /* we must stop and use the gradient direction. This condition need */
364: /* only be checked once. */
365: /*************************************************************************/
367: ksp->reason = KSP_DIVERGED_NAN;
368: PetscInfo1(ksp, "KSPSolve_GLTR: bad matrix: kappa=%g\n", kappa);
370: if (cg->radius) {
371: if (r2 >= rr) {
372: alpha = 1.0;
373: cg->norm_d = PetscSqrtReal(rr);
374: }
375: else {
376: alpha = PetscSqrtReal(r2 / rr);
377: cg->norm_d = cg->radius;
378: }
380: VecAXPY(d, alpha, r); /* d = d + alpha r */
382: /***********************************************************************/
383: /* Compute objective function. */
384: /***********************************************************************/
386: KSP_MatMult(ksp, Qmat, d, z);
387: VecAYPX(z, -0.5, ksp->vec_rhs);
388: VecDot(d, z, &cg->o_fcn);
389: cg->o_fcn = -cg->o_fcn;
390: ++ksp->its;
391: }
392: return(0);
393: }
395: /***************************************************************************/
396: /* Initialize variables for calculating the norm of the direction and for */
397: /* the Lanczos tridiagonal matrix. Note that we track the diagonal value */
398: /* of the Cholesky factorization of the Lanczos matrix in order to */
399: /* determine when negative curvature is encountered. */
400: /***************************************************************************/
402: dMp = 0.0;
403: norm_d = 0.0;
404: switch(cg->dtype) {
405: case GLTR_PRECONDITIONED_DIRECTION:
406: norm_p = rz;
407: break;
409: default:
410: VecDot(p, p, &norm_p);
411: break;
412: }
414: cg->diag[t_size] = kappa / rz;
415: cg->offd[t_size] = 0.0;
416: ++t_size;
418: piv = 1.0;
420: /***************************************************************************/
421: /* Check for breakdown of the conjugate gradient method; this occurs when */
422: /* kappa is zero. */
423: /***************************************************************************/
425: if (fabs(kappa) <= 0.0) {
426: /*************************************************************************/
427: /* The curvature is zero. In this case, we must stop and use follow */
428: /* the direction of negative curvature since the Lanczos matrix is zero. */
429: /*************************************************************************/
431: ksp->reason = KSP_DIVERGED_BREAKDOWN;
432: PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: kappa=%g\n", kappa);
434: if (cg->radius && norm_p > 0.0) {
435: /***********************************************************************/
436: /* Follow direction of negative curvature to the boundary of the */
437: /* trust region. */
438: /***********************************************************************/
440: step = PetscSqrtReal(r2 / norm_p);
441: cg->norm_d = cg->radius;
443: VecAXPY(d, step, p); /* d = d + step p */
445: /***********************************************************************/
446: /* Update objective function. */
447: /***********************************************************************/
449: cg->o_fcn += step * (0.5 * step * kappa - rz);
450: }
451: else if (cg->radius) {
452: /***********************************************************************/
453: /* The norm of the preconditioned direction is zero; use the gradient */
454: /* step. */
455: /***********************************************************************/
457: if (r2 >= rr) {
458: alpha = 1.0;
459: cg->norm_d = PetscSqrtReal(rr);
460: }
461: else {
462: alpha = PetscSqrtReal(r2 / rr);
463: cg->norm_d = cg->radius;
464: }
466: VecAXPY(d, alpha, r); /* d = d + alpha r */
468: /***********************************************************************/
469: /* Compute objective function. */
470: /***********************************************************************/
472: KSP_MatMult(ksp, Qmat, d, z);
473: VecAYPX(z, -0.5, ksp->vec_rhs);
474: VecDot(d, z, &cg->o_fcn);
475: cg->o_fcn = -cg->o_fcn;
476: ++ksp->its;
477: }
478: return(0);
479: }
481: /***************************************************************************/
482: /* Begin the first part of the GLTR algorithm which runs the conjugate */
483: /* gradient method until either the problem is solved, we encounter the */
484: /* boundary of the trust region, or the conjugate gradient method breaks */
485: /* down. */
486: /***************************************************************************/
488: while(1) {
489: /*************************************************************************/
490: /* Know that kappa is nonzero, because we have not broken down, so we */
491: /* can compute the steplength. */
492: /*************************************************************************/
494: alpha = rz / kappa;
495: cg->alpha[l_size] = alpha;
497: /*************************************************************************/
498: /* Compute the diagonal value of the Cholesky factorization of the */
499: /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
500: /* This indicates a direction of negative curvature. */
501: /*************************************************************************/
503: piv = cg->diag[l_size] - cg->offd[l_size]*cg->offd[l_size] / piv;
504: if (piv <= 0.0) {
505: /***********************************************************************/
506: /* In this case, the matrix is indefinite and we have encountered */
507: /* a direction of negative curvature. Follow the direction to the */
508: /* boundary of the trust region. */
509: /***********************************************************************/
511: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
512: PetscInfo1(ksp, "KSPSolve_GLTR: negative curvature: kappa=%g\n", kappa);
514: if (cg->radius && norm_p > 0.0) {
515: /*********************************************************************/
516: /* Follow direction of negative curvature to boundary. */
517: /*********************************************************************/
519: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
520: cg->norm_d = cg->radius;
522: VecAXPY(d, step, p); /* d = d + step p */
524: /*********************************************************************/
525: /* Update objective function. */
526: /*********************************************************************/
528: cg->o_fcn += step * (0.5 * step * kappa - rz);
529: }
530: else if (cg->radius) {
531: /*********************************************************************/
532: /* The norm of the direction is zero; there is nothing to follow. */
533: /*********************************************************************/
534: }
535: break;
536: }
538: /*************************************************************************/
539: /* Compute the steplength and check for intersection with the trust */
540: /* region. */
541: /*************************************************************************/
543: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
544: if (cg->radius && norm_dp1 >= r2) {
545: /***********************************************************************/
546: /* In this case, the matrix is positive definite as far as we know. */
547: /* However, the full step goes beyond the trust region. */
548: /***********************************************************************/
550: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
551: PetscInfo1(ksp, "KSPSolve_GLTR: constrained step: radius=%g\n", cg->radius);
553: if (norm_p > 0.0) {
554: /*********************************************************************/
555: /* Follow the direction to the boundary of the trust region. */
556: /*********************************************************************/
558: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
559: cg->norm_d = cg->radius;
561: VecAXPY(d, step, p); /* d = d + step p */
563: /*********************************************************************/
564: /* Update objective function. */
565: /*********************************************************************/
567: cg->o_fcn += step * (0.5 * step * kappa - rz);
568: }
569: else {
570: /*********************************************************************/
571: /* The norm of the direction is zero; there is nothing to follow. */
572: /*********************************************************************/
573: }
574: break;
575: }
577: /*************************************************************************/
578: /* Now we can update the direction and residual. */
579: /*************************************************************************/
581: VecAXPY(d, alpha, p); /* d = d + alpha p */
582: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
583: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
585: switch(cg->dtype) {
586: case GLTR_PRECONDITIONED_DIRECTION:
587: norm_d = norm_dp1;
588: break;
590: default:
591: VecDot(d, d, &norm_d);
592: break;
593: }
594: cg->norm_d = PetscSqrtReal(norm_d);
596: /*************************************************************************/
597: /* Update objective function. */
598: /*************************************************************************/
600: cg->o_fcn -= 0.5 * alpha * rz;
602: /*************************************************************************/
603: /* Check that the preconditioner appears positive semidefinite. */
604: /*************************************************************************/
606: rzm1 = rz;
607: VecDot(r, z, &rz); /* rz = r^T z */
608: if (rz < 0.0) {
609: /***********************************************************************/
610: /* The preconditioner is indefinite. */
611: /***********************************************************************/
613: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
614: PetscInfo1(ksp, "KSPSolve_GLTR: cg indefinite preconditioner: rz=%g\n", rz);
615: break;
616: }
618: /*************************************************************************/
619: /* As far as we know, the preconditioner is positive semidefinite. */
620: /* Compute the residual and check for convergence. */
621: /*************************************************************************/
623: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
625: switch(ksp->normtype) {
626: case KSP_NORM_PRECONDITIONED:
627: VecNorm(z, NORM_2, &norm_r);/* norm_r = |z| */
628: break;
630: case KSP_NORM_UNPRECONDITIONED:
631: VecNorm(r, NORM_2, &norm_r);/* norm_r = |r| */
632: break;
634: case KSP_NORM_NATURAL:
635: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
636: break;
638: default:
639: norm_r = 0.0;
640: break;
641: }
643: KSPLogResidualHistory(ksp, norm_r);
644: KSPMonitor(ksp, ksp->its, norm_r);
645: ksp->rnorm = norm_r;
646:
647: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
648: if (ksp->reason) {
649: /***********************************************************************/
650: /* The method has converged. */
651: /***********************************************************************/
653: PetscInfo2(ksp, "KSPSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", norm_r, cg->radius);
654: break;
655: }
657: /*************************************************************************/
658: /* We have not converged yet. Check for breakdown. */
659: /*************************************************************************/
661: beta = rz / rzm1;
662: if (fabs(beta) <= 0.0) {
663: /***********************************************************************/
664: /* Conjugate gradients has broken down. */
665: /***********************************************************************/
667: ksp->reason = KSP_DIVERGED_BREAKDOWN;
668: PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n", beta);
669: break;
670: }
672: /*************************************************************************/
673: /* Check iteration limit. */
674: /*************************************************************************/
676: if (ksp->its >= max_cg_its) {
677: ksp->reason = KSP_DIVERGED_ITS;
678: PetscInfo1(ksp, "KSPSolve_GLTR: iterlim: its=%d\n", ksp->its);
679: break;
680: }
682: /*************************************************************************/
683: /* Update p and the norms. */
684: /*************************************************************************/
686: cg->beta[l_size] = beta;
687: VecAYPX(p, beta, z); /* p = z + beta p */
689: switch(cg->dtype) {
690: case GLTR_PRECONDITIONED_DIRECTION:
691: dMp = beta*(dMp + alpha*norm_p);
692: norm_p = beta*(rzm1 + beta*norm_p);
693: break;
695: default:
696: VecDot(d, p, &dMp);
697: VecDot(p, p, &norm_p);
698: break;
699: }
701: /*************************************************************************/
702: /* Compute the new direction and update the iteration. */
703: /*************************************************************************/
705: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
706: VecDot(p, z, &kappa); /* kappa = p^T Q p */
707: ++ksp->its;
709: /*************************************************************************/
710: /* Update the the Lanczos tridiagonal matrix. */
711: /*************************************************************************/
713: ++l_size;
714: cg->offd[t_size] = PetscSqrtReal(beta) / fabs(alpha);
715: cg->diag[t_size] = kappa / rz + beta / alpha;
716: ++t_size;
718: /*************************************************************************/
719: /* Check for breakdown of the conjugate gradient method; this occurs */
720: /* when kappa is zero. */
721: /*************************************************************************/
723: if (fabs(kappa) <= 0.0) {
724: /***********************************************************************/
725: /* The method breaks down; move along the direction as if the matrix */
726: /* were indefinite. */
727: /***********************************************************************/
729: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
730: PetscInfo1(ksp, "KSPSolve_GLTR: cg breakdown: kappa=%g\n", kappa);
732: if (cg->radius && norm_p > 0.0) {
733: /*********************************************************************/
734: /* Follow direction to boundary. */
735: /*********************************************************************/
737: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
738: cg->norm_d = cg->radius;
740: VecAXPY(d, step, p); /* d = d + step p */
742: /*********************************************************************/
743: /* Update objective function. */
744: /*********************************************************************/
746: cg->o_fcn += step * (0.5 * step * kappa - rz);
747: }
748: else if (cg->radius) {
749: /*********************************************************************/
750: /* The norm of the direction is zero; there is nothing to follow. */
751: /*********************************************************************/
752: }
753: break;
754: }
755: }
757: /***************************************************************************/
758: /* Check to see if we need to continue with the Lanczos method. */
759: /***************************************************************************/
761: if (!cg->radius) {
762: /*************************************************************************/
763: /* There is no radius. Therefore, we cannot move along the boundary. */
764: /*************************************************************************/
766: return(0);
767: }
769: if (KSP_CONVERGED_CG_NEG_CURVE != ksp->reason) {
770: /*************************************************************************/
771: /* The method either converged to an interior point, hit the boundary of */
772: /* the trust-region without encountering a direction of negative */
773: /* curvature or the iteration limit was reached. */
774: /*************************************************************************/
776: return(0);
777: }
779: /***************************************************************************/
780: /* Switch to contructing the Lanczos basis by way of the conjugate */
781: /* directions. */
782: /***************************************************************************/
784: for (i = 0; i < max_lanczos_its; ++i) {
785: /*************************************************************************/
786: /* Check for breakdown of the conjugate gradient method; this occurs */
787: /* when kappa is zero. */
788: /*************************************************************************/
790: if (fabs(kappa) <= 0.0) {
791: ksp->reason = KSP_DIVERGED_BREAKDOWN;
792: PetscInfo1(ksp, "KSPSolve_GLTR: lanczos breakdown: kappa=%g\n", kappa);
793: break;
794: }
796: /*************************************************************************/
797: /* Update the direction and residual. */
798: /*************************************************************************/
799:
800: alpha = rz / kappa;
801: cg->alpha[l_size] = alpha;
803: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
804: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
806: /*************************************************************************/
807: /* Check that the preconditioner appears positive semidefinite. */
808: /*************************************************************************/
810: rzm1 = rz;
811: VecDot(r, z, &rz); /* rz = r^T z */
812: if (rz < 0.0) {
813: /***********************************************************************/
814: /* The preconditioner is indefinite. */
815: /***********************************************************************/
817: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
818: PetscInfo1(ksp, "KSPSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", rz);
819: break;
820: }
822: /*************************************************************************/
823: /* As far as we know, the preconditioner is positive definite. Compute */
824: /* the residual. Do NOT check for convergence. */
825: /*************************************************************************/
827: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
829: switch(ksp->normtype) {
830: case KSP_NORM_PRECONDITIONED:
831: VecNorm(z, NORM_2, &norm_r);/* norm_r = |z| */
832: break;
834: case KSP_NORM_UNPRECONDITIONED:
835: VecNorm(r, NORM_2, &norm_r);/* norm_r = |r| */
836: break;
838: case KSP_NORM_NATURAL:
839: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
840: break;
842: default:
843: norm_r = 0.0;
844: break;
845: }
847: KSPLogResidualHistory(ksp, norm_r);
848: KSPMonitor(ksp, ksp->its, norm_r);
849: ksp->rnorm = norm_r;
850:
851: /*************************************************************************/
852: /* Check for breakdown. */
853: /*************************************************************************/
855: beta = rz / rzm1;
856: if (fabs(beta) <= 0.0) {
857: /***********************************************************************/
858: /* Conjugate gradients has broken down. */
859: /***********************************************************************/
861: ksp->reason = KSP_DIVERGED_BREAKDOWN;
862: PetscInfo1(ksp, "KSPSolve_GLTR: breakdown: beta=%g\n", beta);
863: break;
864: }
866: /*************************************************************************/
867: /* Update p and the norms. */
868: /*************************************************************************/
870: cg->beta[l_size] = beta;
871: VecAYPX(p, beta, z); /* p = z + beta p */
873: /*************************************************************************/
874: /* Compute the new direction and update the iteration. */
875: /*************************************************************************/
877: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
878: VecDot(p, z, &kappa); /* kappa = p^T Q p */
879: ++ksp->its;
881: /*************************************************************************/
882: /* Update the Lanczos tridiagonal matrix. */
883: /*************************************************************************/
885: ++l_size;
886: cg->offd[t_size] = PetscSqrtReal(beta) / fabs(alpha);
887: cg->diag[t_size] = kappa / rz + beta / alpha;
888: ++t_size;
889: }
891: /***************************************************************************/
892: /* We have the Lanczos basis, solve the tridiagonal trust-region problem */
893: /* to obtain the Lanczos direction. We know that the solution lies on */
894: /* the boundary of the trust region. We start by checking that the */
895: /* workspace allocated is large enough. */
896: /***************************************************************************/
897: /* Note that the current version only computes the solution by using the */
898: /* preconditioned direction. Need to think about how to do the */
899: /* unpreconditioned direction calculation. */
900: /***************************************************************************/
902: if (t_size > cg->alloced) {
903: if (cg->alloced) {
904: PetscFree(cg->rwork);
905: PetscFree(cg->iwork);
906: cg->alloced += cg->init_alloc;
907: } else {
908: cg->alloced = cg->init_alloc;
909: }
911: while (t_size > cg->alloced) {
912: cg->alloced += cg->init_alloc;
913: }
914:
915: cg->alloced = PetscMin(cg->alloced, t_size);
916: PetscMalloc2(10*cg->alloced,PetscReal, &cg->rwork,5*cg->alloced,PetscBLASInt, &cg->iwork);
917: }
919: /***************************************************************************/
920: /* Set up the required vectors. */
921: /***************************************************************************/
923: t_soln = cg->rwork + 0*t_size; /* Solution */
924: t_diag = cg->rwork + 1*t_size; /* Diagonal of T */
925: t_offd = cg->rwork + 2*t_size; /* Off-diagonal of T */
926: e_valu = cg->rwork + 3*t_size; /* Eigenvalues of T */
927: e_vect = cg->rwork + 4*t_size; /* Eigenvector of T */
928: e_rwrk = cg->rwork + 5*t_size; /* Eigen workspace */
929:
930: e_iblk = cg->iwork + 0*t_size; /* Eigen blocks */
931: e_splt = cg->iwork + 1*t_size; /* Eigen splits */
932: e_iwrk = cg->iwork + 2*t_size; /* Eigen workspace */
934: /***************************************************************************/
935: /* Compute the minimum eigenvalue of T. */
936: /***************************************************************************/
938: vl = 0.0;
939: vu = 0.0;
940: il = 1;
941: iu = 1;
943: #if defined(PETSC_MISSING_LAPACK_STEBZ)
944: SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"STEBZ - Lapack routine is unavailable.");
945: #else
946: LAPACKstebz_("I", "E", &t_size, &vl, &vu, &il, &iu, &cg->eigen_tol,
947: cg->diag, cg->offd + 1, &e_valus, &e_splts, e_valu,
948: e_iblk, e_splt, e_rwrk, e_iwrk, &info);
950: if ((0 != info) || (1 != e_valus)) {
951: /*************************************************************************/
952: /* Calculation of the minimum eigenvalue failed. Return the */
953: /* Steihaug-Toint direction. */
954: /*************************************************************************/
956: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvalue.\n");
957: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
958: return(0);
959: }
961: cg->e_min = e_valu[0];
963: /***************************************************************************/
964: /* Compute the initial value of lambda to make (T + lamba I) positive */
965: /* definite. */
966: /***************************************************************************/
968: pert = cg->init_pert;
969: if (e_valu[0] < 0.0) {
970: cg->lambda = pert - e_valu[0];
971: }
972: #endif
974: while(1) {
975: for (i = 0; i < t_size; ++i) {
976: t_diag[i] = cg->diag[i] + cg->lambda;
977: t_offd[i] = cg->offd[i];
978: }
980: #if defined(PETSC_MISSING_LAPACK_PTTRF)
981: SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
982: #else
983: LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info);
985: if (0 == info) {
986: break;
987: }
989: pert += pert;
990: cg->lambda = cg->lambda * (1.0 + pert) + pert;
991: #endif
992: }
994: /***************************************************************************/
995: /* Compute the initial step and its norm. */
996: /***************************************************************************/
998: nrhs = 1;
999: nldb = t_size;
1001: t_soln[0] = -cg->norm_r[0];
1002: for (i = 1; i < t_size; ++i) {
1003: t_soln[i] = 0.0;
1004: }
1006: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1007: SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1008: #else
1009: LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info);
1010: #endif
1012: if (0 != info) {
1013: /*************************************************************************/
1014: /* Calculation of the initial step failed; return the Steihaug-Toint */
1015: /* direction. */
1016: /*************************************************************************/
1018: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1019: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1020: return(0);
1021: }
1023: norm_t = 0.;
1024: for (i = 0; i < t_size; ++i) {
1025: norm_t += t_soln[i] * t_soln[i];
1026: }
1027: norm_t = PetscSqrtReal(norm_t);
1029: /***************************************************************************/
1030: /* Determine the case we are in. */
1031: /***************************************************************************/
1033: if (norm_t <= cg->radius) {
1034: /*************************************************************************/
1035: /* The step is within the trust region; check if we are in the hard case */
1036: /* and need to move to the boundary by following a direction of negative */
1037: /* curvature. */
1038: /*************************************************************************/
1039:
1040: if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
1041: /***********************************************************************/
1042: /* This is the hard case; compute the eigenvector associated with the */
1043: /* minimum eigenvalue and move along this direction to the boundary. */
1044: /***********************************************************************/
1046: #if defined(PETSC_MISSING_LAPACK_STEIN)
1047: SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"STEIN - Lapack routine is unavailable.");
1048: #else
1049: LAPACKstein_(&t_size, cg->diag, cg->offd + 1, &e_valus, e_valu,
1050: e_iblk, e_splt, e_vect, &nldb,
1051: e_rwrk, e_iwrk, e_iwrk + t_size, &info);
1052: #endif
1054: if (0 != info) {
1055: /*********************************************************************/
1056: /* Calculation of the minimum eigenvalue failed. Return the */
1057: /* Steihaug-Toint direction. */
1058: /*********************************************************************/
1059:
1060: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute eigenvector.\n");
1061: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1062: return(0);
1063: }
1064:
1065: coef1 = 0.0;
1066: coef2 = 0.0;
1067: coef3 = -cg->radius * cg->radius;
1068: for (i = 0; i < t_size; ++i) {
1069: coef1 += e_vect[i] * e_vect[i];
1070: coef2 += e_vect[i] * t_soln[i];
1071: coef3 += t_soln[i] * t_soln[i];
1072: }
1073:
1074: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
1075: root1 = (-coef2 + coef3) / coef1;
1076: root2 = (-coef2 - coef3) / coef1;
1078: /***********************************************************************/
1079: /* Compute objective value for (t_soln + root1 * e_vect) */
1080: /***********************************************************************/
1082: for (i = 0; i < t_size; ++i) {
1083: e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
1084: }
1085:
1086: obj1 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1087: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1088: for (i = 1; i < t_size - 1; ++i) {
1089: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1090: cg->diag[i]*e_rwrk[i]+
1091: cg->offd[i+1]*e_rwrk[i+1]);
1092: }
1093: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1094: cg->diag[i]*e_rwrk[i]);
1096: /***********************************************************************/
1097: /* Compute objective value for (t_soln + root2 * e_vect) */
1098: /***********************************************************************/
1100: for (i = 0; i < t_size; ++i) {
1101: e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
1102: }
1103:
1104: obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
1105: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
1106: for (i = 1; i < t_size - 1; ++i) {
1107: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1108: cg->diag[i]*e_rwrk[i]+
1109: cg->offd[i+1]*e_rwrk[i+1]);
1110: }
1111: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
1112: cg->diag[i]*e_rwrk[i]);
1113:
1114: /***********************************************************************/
1115: /* Choose the point with the best objective function value. */
1116: /***********************************************************************/
1118: if (obj1 <= obj2) {
1119: for (i = 0; i < t_size; ++i) {
1120: t_soln[i] += root1 * e_vect[i];
1121: }
1122: }
1123: else {
1124: for (i = 0; i < t_size; ++i) {
1125: t_soln[i] += root2 * e_vect[i];
1126: }
1127: }
1128: }
1129: else {
1130: /***********************************************************************/
1131: /* The matrix is positive definite or there was no room to move; the */
1132: /* solution is already contained in t_soln. */
1133: /***********************************************************************/
1134: }
1135: }
1136: else {
1137: /*************************************************************************/
1138: /* The step is outside the trust-region. Compute the correct value for */
1139: /* lambda by performing Newton's method. */
1140: /*************************************************************************/
1142: for (i = 0; i < max_newton_its; ++i) {
1143: /***********************************************************************/
1144: /* Check for convergence. */
1145: /***********************************************************************/
1147: if (fabs(norm_t - cg->radius) <= cg->newton_tol * cg->radius) {
1148: break;
1149: }
1151: /***********************************************************************/
1152: /* Compute the update. */
1153: /***********************************************************************/
1155: PetscMemcpy(e_rwrk, t_soln, sizeof(PetscReal)*t_size);
1157: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1158: SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1159: #else
1160: LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info);
1161: #endif
1163: if (0 != info) {
1164: /*********************************************************************/
1165: /* Calculation of the step failed; return the Steihaug-Toint */
1166: /* direction. */
1167: /*********************************************************************/
1169: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1170: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1171: return(0);
1172: }
1174: /***********************************************************************/
1175: /* Modify lambda. */
1176: /***********************************************************************/
1178: norm_w = 0.;
1179: for (j = 0; j < t_size; ++j) {
1180: norm_w += t_soln[j] * e_rwrk[j];
1181: }
1182:
1183: cg->lambda += (norm_t - cg->radius)/cg->radius * (norm_t * norm_t) / norm_w;
1185: /***********************************************************************/
1186: /* Factor T + lambda I */
1187: /***********************************************************************/
1188:
1189: for (j = 0; j < t_size; ++j) {
1190: t_diag[j] = cg->diag[j] + cg->lambda;
1191: t_offd[j] = cg->offd[j];
1192: }
1194: #if defined(PETSC_MISSING_LAPACK_PTTRF)
1195: SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
1196: #else
1197: LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info);
1198: #endif
1200: if (0 != info) {
1201: /*********************************************************************/
1202: /* Calculation of factorization failed; return the Steihaug-Toint */
1203: /* direction. */
1204: /*********************************************************************/
1206: PetscInfo(ksp, "KSPSolve_GLTR: factorization failed.\n");
1207: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1208: return(0);
1209: }
1211: /***********************************************************************/
1212: /* Compute the new step and its norm. */
1213: /***********************************************************************/
1215: t_soln[0] = -cg->norm_r[0];
1216: for (j = 1; j < t_size; ++j) {
1217: t_soln[j] = 0.0;
1218: }
1220: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1221: SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1222: #else
1223: LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info);
1224: #endif
1226: if (0 != info) {
1227: /*********************************************************************/
1228: /* Calculation of the step failed; return the Steihaug-Toint */
1229: /* direction. */
1230: /*********************************************************************/
1231:
1232: PetscInfo(ksp, "KSPSolve_GLTR: failed to compute step.\n");
1233: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1234: return(0);
1235: }
1237: norm_t = 0.;
1238: for (j = 0; j < t_size; ++j) {
1239: norm_t += t_soln[j] * t_soln[j];
1240: }
1241: norm_t = PetscSqrtReal(norm_t);
1242: }
1244: /*************************************************************************/
1245: /* Check for convergence. */
1246: /*************************************************************************/
1248: if (fabs(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1249: /***********************************************************************/
1250: /* Newton method failed to converge in iteration limit. */
1251: /***********************************************************************/
1253: PetscInfo(ksp, "KSPSolve_GLTR: failed to converge.\n");
1254: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1255: return(0);
1256: }
1257: }
1259: /***************************************************************************/
1260: /* Recover the norm of the direction and objective function value. */
1261: /***************************************************************************/
1263: cg->norm_d = norm_t;
1265: cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1266: for (i = 1; i < t_size - 1; ++i) {
1267: 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]);
1268: }
1269: cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);
1271: /***************************************************************************/
1272: /* Recover the direction. */
1273: /***************************************************************************/
1275: sigma = -1;
1277: /***************************************************************************/
1278: /* Start conjugate gradient method from the beginning */
1279: /***************************************************************************/
1281: VecCopy(ksp->vec_rhs, r); /* r = -grad */
1282: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1284: /***************************************************************************/
1285: /* Accumulate Q * s */
1286: /***************************************************************************/
1288: VecCopy(z, d);
1289: VecScale(d, sigma * t_soln[0] / cg->norm_r[0]);
1290:
1291: /***************************************************************************/
1292: /* Compute the first direction. */
1293: /***************************************************************************/
1295: VecCopy(z, p); /* p = z */
1296: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1297: ++ksp->its;
1299: for (i = 0; i < l_size - 1; ++i) {
1300: /*************************************************************************/
1301: /* Update the residual and direction. */
1302: /*************************************************************************/
1304: alpha = cg->alpha[i];
1305: if (alpha >= 0.0) {
1306: sigma = -sigma;
1307: }
1309: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1310: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1312: /*************************************************************************/
1313: /* Accumulate Q * s */
1314: /*************************************************************************/
1316: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1318: /*************************************************************************/
1319: /* Update p. */
1320: /*************************************************************************/
1322: beta = cg->beta[i];
1323: VecAYPX(p, beta, z); /* p = z + beta p */
1324: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1325: ++ksp->its;
1326: }
1328: /***************************************************************************/
1329: /* Update the residual and direction. */
1330: /***************************************************************************/
1332: alpha = cg->alpha[i];
1333: if (alpha >= 0.0) {
1334: sigma = -sigma;
1335: }
1337: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1338: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1340: /***************************************************************************/
1341: /* Accumulate Q * s */
1342: /***************************************************************************/
1344: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1346: /***************************************************************************/
1347: /* Set the termination reason. */
1348: /***************************************************************************/
1350: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1351: return(0);
1352: #endif
1353: }
1357: PetscErrorCode KSPSetUp_GLTR(KSP ksp)
1358: {
1359: KSP_GLTR *cg = (KSP_GLTR *)ksp->data;
1360: PetscInt max_its;
1364: /***************************************************************************/
1365: /* Determine the total maximum number of iterations. */
1366: /***************************************************************************/
1368: max_its = ksp->max_it + cg->max_lanczos_its + 1;
1370: /***************************************************************************/
1371: /* Set work vectors needed by conjugate gradient method and allocate */
1372: /* workspace for Lanczos matrix. */
1373: /***************************************************************************/
1375: KSPDefaultGetWork(ksp, 3);
1377: 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);
1378: PetscMemzero(cg->diag, max_its*sizeof(PetscReal));
1379: PetscMemzero(cg->offd, max_its*sizeof(PetscReal));
1380: PetscMemzero(cg->alpha, max_its*sizeof(PetscReal));
1381: PetscMemzero(cg->beta, max_its*sizeof(PetscReal));
1382: PetscMemzero(cg->norm_r, max_its*sizeof(PetscReal));
1383: PetscLogObjectMemory(ksp, 5*max_its*sizeof(PetscReal));
1384: return(0);
1385: }
1389: PetscErrorCode KSPDestroy_GLTR(KSP ksp)
1390: {
1391: KSP_GLTR *cg = (KSP_GLTR *)ksp->data;
1395: /***************************************************************************/
1396: /* Free memory allocated for the data. */
1397: /***************************************************************************/
1399: PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1400: if (cg->alloced) {
1401: PetscFree2(cg->rwork,cg->iwork);
1402: }
1404: /***************************************************************************/
1405: /* Clear composed functions */
1406: /***************************************************************************/
1408: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGLTRSetRadius_C","",PETSC_NULL);
1409: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGLTRGetNormD_C","",PETSC_NULL);
1410: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGLTRGetObjFcn_C","",PETSC_NULL);
1411: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGLTRGetMinEig_C","",PETSC_NULL);
1412: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPGLTRGetLambda_C","",PETSC_NULL);
1414: /***************************************************************************/
1415: /* Destroy KSP object. */
1416: /***************************************************************************/
1417: KSPDefaultDestroy(ksp);
1418: return(0);
1419: }
1421: EXTERN_C_BEGIN
1424: PetscErrorCode KSPGLTRSetRadius_GLTR(KSP ksp, PetscReal radius)
1425: {
1426: KSP_GLTR *cg = (KSP_GLTR *)ksp->data;
1429: cg->radius = radius;
1430: return(0);
1431: }
1435: PetscErrorCode KSPGLTRGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1436: {
1437: KSP_GLTR *cg = (KSP_GLTR *)ksp->data;
1440: *norm_d = cg->norm_d;
1441: return(0);
1442: }
1446: PetscErrorCode KSPGLTRGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1447: {
1448: KSP_GLTR *cg = (KSP_GLTR *)ksp->data;
1451: *o_fcn = cg->o_fcn;
1452: return(0);
1453: }
1457: PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1458: {
1459: KSP_GLTR *cg = (KSP_GLTR *)ksp->data;
1462: *e_min = cg->e_min;
1463: return(0);
1464: }
1468: PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1469: {
1470: KSP_GLTR *cg = (KSP_GLTR *)ksp->data;
1473: *lambda = cg->lambda;
1474: return(0);
1475: }
1476: EXTERN_C_END
1480: PetscErrorCode KSPSetFromOptions_GLTR(KSP ksp)
1481: {
1483: KSP_GLTR *cg = (KSP_GLTR *)ksp->data;
1486: PetscOptionsHead("KSP GLTR options");
1488: PetscOptionsReal("-ksp_gltr_radius", "Trust Region Radius", "KSPGLTRSetRadius", cg->radius, &cg->radius, PETSC_NULL);
1490: PetscOptionsReal("-ksp_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, PETSC_NULL);
1491: PetscOptionsReal("-ksp_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, PETSC_NULL);
1492: PetscOptionsReal("-ksp_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, PETSC_NULL);
1494: PetscOptionsInt("-ksp_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, PETSC_NULL);
1495: PetscOptionsInt("-ksp_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, PETSC_NULL);
1497: PetscOptionsEList("-ksp_gltr_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, PETSC_NULL);
1499: PetscOptionsTail();
1500: return(0);
1501: }
1503: /*MC
1504: KSPGLTR - Code to run conjugate gradient method subject to a constraint
1505: on the solution norm. This is used in Trust Region methods for
1506: nonlinear equations, SNESTR
1508: Options Database Keys:
1509: . -ksp_gltr_radius <r> - Trust Region Radius
1511: Notes: This is rarely used directly
1513: Use preconditioned conjugate gradient to compute
1514: an approximate minimizer of the quadratic function
1516: q(s) = g^T * s + .5 * s^T * H * s
1518: subject to the trust region constraint
1520: || s || <= delta,
1522: where
1524: delta is the trust region radius,
1525: g is the gradient vector,
1526: H is the Hessian approximation,
1527: M is the positive definite preconditioner matrix.
1529: KSPConvergedReason may be
1530: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1531: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1532: $ other KSP converged/diverged reasons
1534: Notes:
1535: The preconditioner supplied should be symmetric and positive definite.
1537: Level: developer
1539: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGLTRSetRadius(), KSPGLTRGetNormD(), KSPGLTRGetObjFcn(), KSPGLTRGetMinEig(), KSPGLTRGetLambda()
1540: M*/
1542: EXTERN_C_BEGIN
1545: PetscErrorCode KSPCreate_GLTR(KSP ksp)
1546: {
1548: KSP_GLTR *cg;
1551: PetscNewLog(ksp, KSP_GLTR, &cg);
1552: cg->radius = 0.0;
1553: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1555: cg->init_pert = 1.0e-8;
1556: cg->eigen_tol = 1.0e-10;
1557: cg->newton_tol = 1.0e-6;
1559: cg->alloced = 0;
1560: cg->init_alloc = 1024;
1562: cg->max_lanczos_its = 20;
1563: cg->max_newton_its = 10;
1565: ksp->data = (void *) cg;
1566: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
1567: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,1);
1568: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);
1570: /***************************************************************************/
1571: /* Sets the functions that are associated with this data structure */
1572: /* (in C++ this is the same as defining virtual functions). */
1573: /***************************************************************************/
1575: ksp->ops->setup = KSPSetUp_GLTR;
1576: ksp->ops->solve = KSPSolve_GLTR;
1577: ksp->ops->destroy = KSPDestroy_GLTR;
1578: ksp->ops->setfromoptions = KSPSetFromOptions_GLTR;
1579: ksp->ops->buildsolution = KSPDefaultBuildSolution;
1580: ksp->ops->buildresidual = KSPDefaultBuildResidual;
1581: ksp->ops->view = 0;
1583: PetscObjectComposeFunctionDynamic((PetscObject)ksp,
1584: "KSPGLTRSetRadius_C",
1585: "KSPGLTRSetRadius_GLTR",
1586: KSPGLTRSetRadius_GLTR);
1587: PetscObjectComposeFunctionDynamic((PetscObject)ksp,
1588: "KSPGLTRGetNormD_C",
1589: "KSPGLTRGetNormD_GLTR",
1590: KSPGLTRGetNormD_GLTR);
1591: PetscObjectComposeFunctionDynamic((PetscObject)ksp,
1592: "KSPGLTRGetObjFcn_C",
1593: "KSPGLTRGetObjFcn_GLTR",
1594: KSPGLTRGetObjFcn_GLTR);
1595: PetscObjectComposeFunctionDynamic((PetscObject)ksp,
1596: "KSPGLTRGetMinEig_C",
1597: "KSPGLTRGetMinEig_GLTR",
1598: KSPGLTRGetMinEig_GLTR);
1599: PetscObjectComposeFunctionDynamic((PetscObject)ksp,
1600: "KSPGLTRGetLambda_C",
1601: "KSPGLTRGetLambda_GLTR",
1602: KSPGLTRGetLambda_GLTR);
1603: return(0);
1604: }
1605: EXTERN_C_END