Actual source code: gltr.c
petsc-3.10.5 2019-03-28
2: #include <../src/ksp/ksp/impls/cg/gltr/gltrimpl.h>
3: #include <petscblaslapack.h>
5: #define GLTR_PRECONDITIONED_DIRECTION 0
6: #define GLTR_UNPRECONDITIONED_DIRECTION 1
7: #define GLTR_DIRECTION_TYPES 2
9: static const char *DType_Table[64] = {"preconditioned", "unpreconditioned"};
11: /*@
12: KSPCGGLTRGetMinEig - Get minimum eigenvalue.
14: Collective on KSP
16: Input Parameters:
17: + ksp - the iterative context
18: - e_min - the minimum eigenvalue
20: Level: advanced
22: .keywords: KSP, GLTR, get, minimum eigenvalue
23: @*/
24: PetscErrorCode KSPCGGLTRGetMinEig(KSP ksp, PetscReal *e_min)
25: {
30: PetscUseMethod(ksp,"KSPCGGLTRGetMinEig_C",(KSP,PetscReal*),(ksp,e_min));
31: return(0);
32: }
34: /*@
35: KSPCGGLTRGetLambda - Get multiplier on trust-region constraint.
37: Not Collective
39: Input Parameters:
40: + ksp - the iterative context
41: - lambda - the multiplier
43: Level: advanced
45: .keywords: KSP, GLTR, get, multiplier
46: @*/
47: PetscErrorCode KSPCGGLTRGetLambda(KSP ksp, PetscReal *lambda)
48: {
53: PetscUseMethod(ksp,"KSPCGGLTRGetLambda_C",(KSP,PetscReal*),(ksp,lambda));
54: return(0);
55: }
57: static PetscErrorCode KSPCGSolve_GLTR(KSP ksp)
58: {
59: #if defined(PETSC_USE_COMPLEX)
60: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "GLTR is not available for complex systems");
61: #else
62: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
63: PetscReal *t_soln, *t_diag, *t_offd, *e_valu, *e_vect, *e_rwrk;
64: PetscBLASInt *e_iblk, *e_splt, *e_iwrk;
67: Mat Qmat, Mmat;
68: Vec r, z, p, d;
69: PC pc;
71: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
72: PetscReal alpha, beta, kappa, rz, rzm1;
73: PetscReal rr, r2, piv, step;
74: PetscReal vl, vu;
75: PetscReal coef1, coef2, coef3, root1, root2, obj1, obj2;
76: PetscReal norm_t, norm_w, pert;
78: PetscInt i, j, max_cg_its, max_lanczos_its, max_newton_its, sigma;
79: PetscBLASInt t_size = 0, l_size = 0, il, iu, info;
80: PetscBLASInt nrhs, nldb;
82: #if !defined(PETSC_MISSING_LAPACK_STEBZ)
83: PetscBLASInt e_valus, e_splts;
84: #endif
85: PetscBool diagonalscale;
88: /***************************************************************************/
89: /* Check the arguments and parameters. */
90: /***************************************************************************/
92: PCGetDiagonalScale(ksp->pc, &diagonalscale);
93: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
94: if (cg->radius < 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE, "Input error: radius < 0");
96: /***************************************************************************/
97: /* Get the workspace vectors and initialize variables */
98: /***************************************************************************/
100: r2 = cg->radius * cg->radius;
101: r = ksp->work[0];
102: z = ksp->work[1];
103: p = ksp->work[2];
104: d = ksp->vec_sol;
105: pc = ksp->pc;
107: PCGetOperators(pc, &Qmat, &Mmat);
109: VecGetSize(d, &max_cg_its);
110: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
111: max_lanczos_its = cg->max_lanczos_its;
112: max_newton_its = cg->max_newton_its;
113: ksp->its = 0;
115: /***************************************************************************/
116: /* Initialize objective function direction, and minimum eigenvalue. */
117: /***************************************************************************/
119: cg->o_fcn = 0.0;
121: VecSet(d, 0.0); /* d = 0 */
122: cg->norm_d = 0.0;
124: cg->e_min = 0.0;
125: cg->lambda = 0.0;
127: /***************************************************************************/
128: /* The first phase of GLTR performs a standard conjugate gradient method, */
129: /* but stores the values required for the Lanczos matrix. We switch to */
130: /* the Lanczos when the conjugate gradient method breaks down. Check the */
131: /* right-hand side for numerical problems. The check for not-a-number and */
132: /* infinite values need be performed only once. */
133: /***************************************************************************/
135: VecCopy(ksp->vec_rhs, r); /* r = -grad */
136: VecDot(r, r, &rr); /* rr = r^T r */
137: KSPCheckDot(ksp,rr);
139: /***************************************************************************/
140: /* Check the preconditioner for numerical problems and for positive */
141: /* definiteness. The check for not-a-number and infinite values need be */
142: /* performed only once. */
143: /***************************************************************************/
145: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
146: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
147: if (PetscIsInfOrNanScalar(rz)) {
148: /*************************************************************************/
149: /* The preconditioner contains not-a-number or an infinite value. */
150: /* Return the gradient direction intersected with the trust region. */
151: /*************************************************************************/
153: ksp->reason = KSP_DIVERGED_NANORINF;
154: PetscInfo1(ksp, "KSPCGSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz);
156: if (cg->radius) {
157: if (r2 >= rr) {
158: alpha = 1.0;
159: cg->norm_d = PetscSqrtReal(rr);
160: } else {
161: alpha = PetscSqrtReal(r2 / rr);
162: cg->norm_d = cg->radius;
163: }
165: VecAXPY(d, alpha, r); /* d = d + alpha r */
167: /***********************************************************************/
168: /* Compute objective function. */
169: /***********************************************************************/
171: KSP_MatMult(ksp, Qmat, d, z);
172: VecAYPX(z, -0.5, ksp->vec_rhs);
173: VecDot(d, z, &cg->o_fcn);
174: cg->o_fcn = -cg->o_fcn;
175: ++ksp->its;
176: }
177: return(0);
178: }
180: if (rz < 0.0) {
181: /*************************************************************************/
182: /* The preconditioner is indefinite. Because this is the first */
183: /* and we do not have a direction yet, we use the gradient step. Note */
184: /* that we cannot use the preconditioned norm when computing the step */
185: /* because the matrix is indefinite. */
186: /*************************************************************************/
188: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
189: PetscInfo1(ksp, "KSPCGSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz);
191: if (cg->radius) {
192: if (r2 >= rr) {
193: alpha = 1.0;
194: cg->norm_d = PetscSqrtReal(rr);
195: } else {
196: alpha = PetscSqrtReal(r2 / rr);
197: cg->norm_d = cg->radius;
198: }
200: VecAXPY(d, alpha, r); /* d = d + alpha r */
202: /***********************************************************************/
203: /* Compute objective function. */
204: /***********************************************************************/
206: KSP_MatMult(ksp, Qmat, d, z);
207: VecAYPX(z, -0.5, ksp->vec_rhs);
208: VecDot(d, z, &cg->o_fcn);
209: cg->o_fcn = -cg->o_fcn;
210: ++ksp->its;
211: }
212: return(0);
213: }
215: /***************************************************************************/
216: /* As far as we know, the preconditioner is positive semidefinite. */
217: /* Compute and log the residual. Check convergence because this */
218: /* initializes things, but do not terminate until at least one conjugate */
219: /* gradient iteration has been performed. */
220: /***************************************************************************/
222: cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
224: switch (ksp->normtype) {
225: case KSP_NORM_PRECONDITIONED:
226: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
227: break;
229: case KSP_NORM_UNPRECONDITIONED:
230: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
231: break;
233: case KSP_NORM_NATURAL:
234: norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
235: break;
237: default:
238: norm_r = 0.0;
239: break;
240: }
242: KSPLogResidualHistory(ksp, norm_r);
243: KSPMonitor(ksp, ksp->its, norm_r);
244: ksp->rnorm = norm_r;
246: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
248: /***************************************************************************/
249: /* Compute the first direction and update the iteration. */
250: /***************************************************************************/
252: VecCopy(z, p); /* p = z */
253: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
254: ++ksp->its;
256: /***************************************************************************/
257: /* Check the matrix for numerical problems. */
258: /***************************************************************************/
260: VecDot(p, z, &kappa); /* kappa = p^T Q p */
261: if (PetscIsInfOrNanScalar(kappa)) {
262: /*************************************************************************/
263: /* The matrix produced not-a-number or an infinite value. In this case, */
264: /* we must stop and use the gradient direction. This condition need */
265: /* only be checked once. */
266: /*************************************************************************/
268: ksp->reason = KSP_DIVERGED_NANORINF;
269: PetscInfo1(ksp, "KSPCGSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa);
271: if (cg->radius) {
272: if (r2 >= rr) {
273: alpha = 1.0;
274: cg->norm_d = PetscSqrtReal(rr);
275: } else {
276: alpha = PetscSqrtReal(r2 / rr);
277: cg->norm_d = cg->radius;
278: }
280: VecAXPY(d, alpha, r); /* d = d + alpha r */
282: /***********************************************************************/
283: /* Compute objective function. */
284: /***********************************************************************/
286: KSP_MatMult(ksp, Qmat, d, z);
287: VecAYPX(z, -0.5, ksp->vec_rhs);
288: VecDot(d, z, &cg->o_fcn);
289: cg->o_fcn = -cg->o_fcn;
290: ++ksp->its;
291: }
292: return(0);
293: }
295: /***************************************************************************/
296: /* Initialize variables for calculating the norm of the direction and for */
297: /* the Lanczos tridiagonal matrix. Note that we track the diagonal value */
298: /* of the Cholesky factorization of the Lanczos matrix in order to */
299: /* determine when negative curvature is encountered. */
300: /***************************************************************************/
302: dMp = 0.0;
303: norm_d = 0.0;
304: switch (cg->dtype) {
305: case GLTR_PRECONDITIONED_DIRECTION:
306: norm_p = rz;
307: break;
309: default:
310: VecDot(p, p, &norm_p);
311: break;
312: }
314: cg->diag[t_size] = kappa / rz;
315: cg->offd[t_size] = 0.0;
316: ++t_size;
318: piv = 1.0;
320: /***************************************************************************/
321: /* Check for breakdown of the conjugate gradient method; this occurs when */
322: /* kappa is zero. */
323: /***************************************************************************/
325: if (PetscAbsReal(kappa) <= 0.0) {
326: /*************************************************************************/
327: /* The curvature is zero. In this case, we must stop and use follow */
328: /* the direction of negative curvature since the Lanczos matrix is zero. */
329: /*************************************************************************/
331: ksp->reason = KSP_DIVERGED_BREAKDOWN;
332: PetscInfo1(ksp, "KSPCGSolve_GLTR: breakdown: kappa=%g\n", (double)kappa);
334: if (cg->radius && norm_p > 0.0) {
335: /***********************************************************************/
336: /* Follow direction of negative curvature to the boundary of the */
337: /* trust region. */
338: /***********************************************************************/
340: step = PetscSqrtReal(r2 / norm_p);
341: cg->norm_d = cg->radius;
343: VecAXPY(d, step, p); /* d = d + step p */
345: /***********************************************************************/
346: /* Update objective function. */
347: /***********************************************************************/
349: cg->o_fcn += step * (0.5 * step * kappa - rz);
350: } else if (cg->radius) {
351: /***********************************************************************/
352: /* The norm of the preconditioned direction is zero; use the gradient */
353: /* step. */
354: /***********************************************************************/
356: if (r2 >= rr) {
357: alpha = 1.0;
358: cg->norm_d = PetscSqrtReal(rr);
359: } else {
360: alpha = PetscSqrtReal(r2 / rr);
361: cg->norm_d = cg->radius;
362: }
364: VecAXPY(d, alpha, r); /* d = d + alpha r */
366: /***********************************************************************/
367: /* Compute objective function. */
368: /***********************************************************************/
370: KSP_MatMult(ksp, Qmat, d, z);
371: VecAYPX(z, -0.5, ksp->vec_rhs);
372: VecDot(d, z, &cg->o_fcn);
373: cg->o_fcn = -cg->o_fcn;
374: ++ksp->its;
375: }
376: return(0);
377: }
379: /***************************************************************************/
380: /* Begin the first part of the GLTR algorithm which runs the conjugate */
381: /* gradient method until either the problem is solved, we encounter the */
382: /* boundary of the trust region, or the conjugate gradient method breaks */
383: /* down. */
384: /***************************************************************************/
386: while (1) {
387: /*************************************************************************/
388: /* Know that kappa is nonzero, because we have not broken down, so we */
389: /* can compute the steplength. */
390: /*************************************************************************/
392: alpha = rz / kappa;
393: cg->alpha[l_size] = alpha;
395: /*************************************************************************/
396: /* Compute the diagonal value of the Cholesky factorization of the */
397: /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
398: /* This indicates a direction of negative curvature. */
399: /*************************************************************************/
401: piv = cg->diag[l_size] - cg->offd[l_size]*cg->offd[l_size] / piv;
402: if (piv <= 0.0) {
403: /***********************************************************************/
404: /* In this case, the matrix is indefinite and we have encountered */
405: /* a direction of negative curvature. Follow the direction to the */
406: /* boundary of the trust region. */
407: /***********************************************************************/
409: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
410: PetscInfo1(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa);
412: if (cg->radius && norm_p > 0.0) {
413: /*********************************************************************/
414: /* Follow direction of negative curvature to boundary. */
415: /*********************************************************************/
417: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
418: cg->norm_d = cg->radius;
420: VecAXPY(d, step, p); /* d = d + step p */
422: /*********************************************************************/
423: /* Update objective function. */
424: /*********************************************************************/
426: cg->o_fcn += step * (0.5 * step * kappa - rz);
427: } else if (cg->radius) {
428: /*********************************************************************/
429: /* The norm of the direction is zero; there is nothing to follow. */
430: /*********************************************************************/
431: }
432: break;
433: }
435: /*************************************************************************/
436: /* Compute the steplength and check for intersection with the trust */
437: /* region. */
438: /*************************************************************************/
440: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
441: if (cg->radius && norm_dp1 >= r2) {
442: /***********************************************************************/
443: /* In this case, the matrix is positive definite as far as we know. */
444: /* However, the full step goes beyond the trust region. */
445: /***********************************************************************/
447: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
448: PetscInfo1(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius);
450: if (norm_p > 0.0) {
451: /*********************************************************************/
452: /* Follow the direction to the boundary of the trust region. */
453: /*********************************************************************/
455: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
456: cg->norm_d = cg->radius;
458: VecAXPY(d, step, p); /* d = d + step p */
460: /*********************************************************************/
461: /* Update objective function. */
462: /*********************************************************************/
464: cg->o_fcn += step * (0.5 * step * kappa - rz);
465: } else {
466: /*********************************************************************/
467: /* The norm of the direction is zero; there is nothing to follow. */
468: /*********************************************************************/
469: }
470: break;
471: }
473: /*************************************************************************/
474: /* Now we can update the direction and residual. */
475: /*************************************************************************/
477: VecAXPY(d, alpha, p); /* d = d + alpha p */
478: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
479: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
481: switch (cg->dtype) {
482: case GLTR_PRECONDITIONED_DIRECTION:
483: norm_d = norm_dp1;
484: break;
486: default:
487: VecDot(d, d, &norm_d);
488: break;
489: }
490: cg->norm_d = PetscSqrtReal(norm_d);
492: /*************************************************************************/
493: /* Update objective function. */
494: /*************************************************************************/
496: cg->o_fcn -= 0.5 * alpha * rz;
498: /*************************************************************************/
499: /* Check that the preconditioner appears positive semidefinite. */
500: /*************************************************************************/
502: rzm1 = rz;
503: VecDot(r, z, &rz); /* rz = r^T z */
504: if (rz < 0.0) {
505: /***********************************************************************/
506: /* The preconditioner is indefinite. */
507: /***********************************************************************/
509: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
510: PetscInfo1(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz);
511: break;
512: }
514: /*************************************************************************/
515: /* As far as we know, the preconditioner is positive semidefinite. */
516: /* Compute the residual and check for convergence. */
517: /*************************************************************************/
519: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
521: switch (ksp->normtype) {
522: case KSP_NORM_PRECONDITIONED:
523: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
524: break;
526: case KSP_NORM_UNPRECONDITIONED:
527: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
528: break;
530: case KSP_NORM_NATURAL:
531: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
532: break;
534: default:
535: norm_r = 0.0;
536: break;
537: }
539: KSPLogResidualHistory(ksp, norm_r);
540: KSPMonitor(ksp, ksp->its, norm_r);
541: ksp->rnorm = norm_r;
543: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
544: if (ksp->reason) {
545: /***********************************************************************/
546: /* The method has converged. */
547: /***********************************************************************/
549: PetscInfo2(ksp, "KSPCGSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
550: break;
551: }
553: /*************************************************************************/
554: /* We have not converged yet. Check for breakdown. */
555: /*************************************************************************/
557: beta = rz / rzm1;
558: if (PetscAbsReal(beta) <= 0.0) {
559: /***********************************************************************/
560: /* Conjugate gradients has broken down. */
561: /***********************************************************************/
563: ksp->reason = KSP_DIVERGED_BREAKDOWN;
564: PetscInfo1(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta);
565: break;
566: }
568: /*************************************************************************/
569: /* Check iteration limit. */
570: /*************************************************************************/
572: if (ksp->its >= max_cg_its) {
573: ksp->reason = KSP_DIVERGED_ITS;
574: PetscInfo1(ksp, "KSPCGSolve_GLTR: iterlim: its=%D\n", ksp->its);
575: break;
576: }
578: /*************************************************************************/
579: /* Update p and the norms. */
580: /*************************************************************************/
582: cg->beta[l_size] = beta;
583: VecAYPX(p, beta, z); /* p = z + beta p */
585: switch (cg->dtype) {
586: case GLTR_PRECONDITIONED_DIRECTION:
587: dMp = beta*(dMp + alpha*norm_p);
588: norm_p = beta*(rzm1 + beta*norm_p);
589: break;
591: default:
592: VecDot(d, p, &dMp);
593: VecDot(p, p, &norm_p);
594: break;
595: }
597: /*************************************************************************/
598: /* Compute the new direction and update the iteration. */
599: /*************************************************************************/
601: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
602: VecDot(p, z, &kappa); /* kappa = p^T Q p */
603: ++ksp->its;
605: /*************************************************************************/
606: /* Update the Lanczos tridiagonal matrix. */
607: /*************************************************************************/
609: ++l_size;
610: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
611: cg->diag[t_size] = kappa / rz + beta / alpha;
612: ++t_size;
614: /*************************************************************************/
615: /* Check for breakdown of the conjugate gradient method; this occurs */
616: /* when kappa is zero. */
617: /*************************************************************************/
619: if (PetscAbsReal(kappa) <= 0.0) {
620: /***********************************************************************/
621: /* The method breaks down; move along the direction as if the matrix */
622: /* were indefinite. */
623: /***********************************************************************/
625: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
626: PetscInfo1(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa);
628: if (cg->radius && norm_p > 0.0) {
629: /*********************************************************************/
630: /* Follow direction to boundary. */
631: /*********************************************************************/
633: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
634: cg->norm_d = cg->radius;
636: VecAXPY(d, step, p); /* d = d + step p */
638: /*********************************************************************/
639: /* Update objective function. */
640: /*********************************************************************/
642: cg->o_fcn += step * (0.5 * step * kappa - rz);
643: } else if (cg->radius) {
644: /*********************************************************************/
645: /* The norm of the direction is zero; there is nothing to follow. */
646: /*********************************************************************/
647: }
648: break;
649: }
650: }
652: /***************************************************************************/
653: /* Check to see if we need to continue with the Lanczos method. */
654: /***************************************************************************/
656: if (!cg->radius) {
657: /*************************************************************************/
658: /* There is no radius. Therefore, we cannot move along the boundary. */
659: /*************************************************************************/
660: return(0);
661: }
663: if (KSP_CONVERGED_CG_NEG_CURVE != ksp->reason) {
664: /*************************************************************************/
665: /* The method either converged to an interior point, hit the boundary of */
666: /* the trust-region without encountering a direction of negative */
667: /* curvature or the iteration limit was reached. */
668: /*************************************************************************/
669: return(0);
670: }
672: /***************************************************************************/
673: /* Switch to contructing the Lanczos basis by way of the conjugate */
674: /* directions. */
675: /***************************************************************************/
677: for (i = 0; i < max_lanczos_its; ++i) {
678: /*************************************************************************/
679: /* Check for breakdown of the conjugate gradient method; this occurs */
680: /* when kappa is zero. */
681: /*************************************************************************/
683: if (PetscAbsReal(kappa) <= 0.0) {
684: ksp->reason = KSP_DIVERGED_BREAKDOWN;
685: PetscInfo1(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa);
686: break;
687: }
689: /*************************************************************************/
690: /* Update the direction and residual. */
691: /*************************************************************************/
693: alpha = rz / kappa;
694: cg->alpha[l_size] = alpha;
696: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
697: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
699: /*************************************************************************/
700: /* Check that the preconditioner appears positive semidefinite. */
701: /*************************************************************************/
703: rzm1 = rz;
704: VecDot(r, z, &rz); /* rz = r^T z */
705: if (rz < 0.0) {
706: /***********************************************************************/
707: /* The preconditioner is indefinite. */
708: /***********************************************************************/
710: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
711: PetscInfo1(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz);
712: break;
713: }
715: /*************************************************************************/
716: /* As far as we know, the preconditioner is positive definite. Compute */
717: /* the residual. Do NOT check for convergence. */
718: /*************************************************************************/
720: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
722: switch (ksp->normtype) {
723: case KSP_NORM_PRECONDITIONED:
724: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
725: break;
727: case KSP_NORM_UNPRECONDITIONED:
728: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
729: break;
731: case KSP_NORM_NATURAL:
732: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
733: break;
735: default:
736: norm_r = 0.0;
737: break;
738: }
740: KSPLogResidualHistory(ksp, norm_r);
741: KSPMonitor(ksp, ksp->its, norm_r);
742: ksp->rnorm = norm_r;
744: /*************************************************************************/
745: /* Check for breakdown. */
746: /*************************************************************************/
748: beta = rz / rzm1;
749: if (PetscAbsReal(beta) <= 0.0) {
750: /***********************************************************************/
751: /* Conjugate gradients has broken down. */
752: /***********************************************************************/
754: ksp->reason = KSP_DIVERGED_BREAKDOWN;
755: PetscInfo1(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n",(double) beta);
756: break;
757: }
759: /*************************************************************************/
760: /* Update p and the norms. */
761: /*************************************************************************/
763: cg->beta[l_size] = beta;
764: VecAYPX(p, beta, z); /* p = z + beta p */
766: /*************************************************************************/
767: /* Compute the new direction and update the iteration. */
768: /*************************************************************************/
770: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
771: VecDot(p, z, &kappa); /* kappa = p^T Q p */
772: ++ksp->its;
774: /*************************************************************************/
775: /* Update the Lanczos tridiagonal matrix. */
776: /*************************************************************************/
778: ++l_size;
779: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
780: cg->diag[t_size] = kappa / rz + beta / alpha;
781: ++t_size;
782: }
784: /***************************************************************************/
785: /* We have the Lanczos basis, solve the tridiagonal trust-region problem */
786: /* to obtain the Lanczos direction. We know that the solution lies on */
787: /* the boundary of the trust region. We start by checking that the */
788: /* workspace allocated is large enough. */
789: /***************************************************************************/
790: /* Note that the current version only computes the solution by using the */
791: /* preconditioned direction. Need to think about how to do the */
792: /* unpreconditioned direction calculation. */
793: /***************************************************************************/
795: if (t_size > cg->alloced) {
796: if (cg->alloced) {
797: PetscFree(cg->rwork);
798: PetscFree(cg->iwork);
799: cg->alloced += cg->init_alloc;
800: } else {
801: cg->alloced = cg->init_alloc;
802: }
804: while (t_size > cg->alloced) {
805: cg->alloced += cg->init_alloc;
806: }
808: cg->alloced = PetscMin(cg->alloced, t_size);
809: PetscMalloc2(10*cg->alloced, &cg->rwork,5*cg->alloced, &cg->iwork);
810: }
812: /***************************************************************************/
813: /* Set up the required vectors. */
814: /***************************************************************************/
816: t_soln = cg->rwork + 0*t_size; /* Solution */
817: t_diag = cg->rwork + 1*t_size; /* Diagonal of T */
818: t_offd = cg->rwork + 2*t_size; /* Off-diagonal of T */
819: e_valu = cg->rwork + 3*t_size; /* Eigenvalues of T */
820: e_vect = cg->rwork + 4*t_size; /* Eigenvector of T */
821: e_rwrk = cg->rwork + 5*t_size; /* Eigen workspace */
823: e_iblk = cg->iwork + 0*t_size; /* Eigen blocks */
824: e_splt = cg->iwork + 1*t_size; /* Eigen splits */
825: e_iwrk = cg->iwork + 2*t_size; /* Eigen workspace */
827: /***************************************************************************/
828: /* Compute the minimum eigenvalue of T. */
829: /***************************************************************************/
831: vl = 0.0;
832: vu = 0.0;
833: il = 1;
834: iu = 1;
836: #if defined(PETSC_MISSING_LAPACK_STEBZ)
837: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEBZ - Lapack routine is unavailable.");
838: #else
839: 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));
841: if ((0 != info) || (1 != e_valus)) {
842: /*************************************************************************/
843: /* Calculation of the minimum eigenvalue failed. Return the */
844: /* Steihaug-Toint direction. */
845: /*************************************************************************/
847: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n");
848: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
849: return(0);
850: }
852: cg->e_min = e_valu[0];
854: /***************************************************************************/
855: /* Compute the initial value of lambda to make (T + lamba I) positive */
856: /* definite. */
857: /***************************************************************************/
859: pert = cg->init_pert;
860: if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
861: #endif
863: while (1) {
864: for (i = 0; i < t_size; ++i) {
865: t_diag[i] = cg->diag[i] + cg->lambda;
866: t_offd[i] = cg->offd[i];
867: }
869: #if defined(PETSC_MISSING_LAPACK_PTTRF)
870: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
871: #else
872: PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
874: if (0 == info) break;
876: pert += pert;
877: cg->lambda = cg->lambda * (1.0 + pert) + pert;
878: #endif
879: }
881: /***************************************************************************/
882: /* Compute the initial step and its norm. */
883: /***************************************************************************/
885: nrhs = 1;
886: nldb = t_size;
888: t_soln[0] = -cg->norm_r[0];
889: for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
891: #if defined(PETSC_MISSING_LAPACK_PTTRS)
892: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
893: #else
894: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
895: #endif
897: if (0 != info) {
898: /*************************************************************************/
899: /* Calculation of the initial step failed; return the Steihaug-Toint */
900: /* direction. */
901: /*************************************************************************/
903: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
904: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
905: return(0);
906: }
908: norm_t = 0.;
909: for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
910: norm_t = PetscSqrtReal(norm_t);
912: /***************************************************************************/
913: /* Determine the case we are in. */
914: /***************************************************************************/
916: if (norm_t <= cg->radius) {
917: /*************************************************************************/
918: /* The step is within the trust region; check if we are in the hard case */
919: /* and need to move to the boundary by following a direction of negative */
920: /* curvature. */
921: /*************************************************************************/
923: if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
924: /***********************************************************************/
925: /* This is the hard case; compute the eigenvector associated with the */
926: /* minimum eigenvalue and move along this direction to the boundary. */
927: /***********************************************************************/
929: #if defined(PETSC_MISSING_LAPACK_STEIN)
930: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"STEIN - Lapack routine is unavailable.");
931: #else
932: 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));
933: #endif
935: if (0 != info) {
936: /*********************************************************************/
937: /* Calculation of the minimum eigenvalue failed. Return the */
938: /* Steihaug-Toint direction. */
939: /*********************************************************************/
941: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n");
942: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
943: return(0);
944: }
946: coef1 = 0.0;
947: coef2 = 0.0;
948: coef3 = -cg->radius * cg->radius;
949: for (i = 0; i < t_size; ++i) {
950: coef1 += e_vect[i] * e_vect[i];
951: coef2 += e_vect[i] * t_soln[i];
952: coef3 += t_soln[i] * t_soln[i];
953: }
955: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
956: root1 = (-coef2 + coef3) / coef1;
957: root2 = (-coef2 - coef3) / coef1;
959: /***********************************************************************/
960: /* Compute objective value for (t_soln + root1 * e_vect) */
961: /***********************************************************************/
963: for (i = 0; i < t_size; ++i) {
964: e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
965: }
967: obj1 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
968: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
969: for (i = 1; i < t_size - 1; ++i) {
970: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
971: cg->diag[i]*e_rwrk[i]+
972: cg->offd[i+1]*e_rwrk[i+1]);
973: }
974: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
975: cg->diag[i]*e_rwrk[i]);
977: /***********************************************************************/
978: /* Compute objective value for (t_soln + root2 * e_vect) */
979: /***********************************************************************/
981: for (i = 0; i < t_size; ++i) {
982: e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
983: }
985: obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
986: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
987: for (i = 1; i < t_size - 1; ++i) {
988: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
989: cg->diag[i]*e_rwrk[i]+
990: cg->offd[i+1]*e_rwrk[i+1]);
991: }
992: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
993: cg->diag[i]*e_rwrk[i]);
995: /***********************************************************************/
996: /* Choose the point with the best objective function value. */
997: /***********************************************************************/
999: if (obj1 <= obj2) {
1000: for (i = 0; i < t_size; ++i) {
1001: t_soln[i] += root1 * e_vect[i];
1002: }
1003: }
1004: else {
1005: for (i = 0; i < t_size; ++i) {
1006: t_soln[i] += root2 * e_vect[i];
1007: }
1008: }
1009: } else {
1010: /***********************************************************************/
1011: /* The matrix is positive definite or there was no room to move; the */
1012: /* solution is already contained in t_soln. */
1013: /***********************************************************************/
1014: }
1015: } else {
1016: /*************************************************************************/
1017: /* The step is outside the trust-region. Compute the correct value for */
1018: /* lambda by performing Newton's method. */
1019: /*************************************************************************/
1021: for (i = 0; i < max_newton_its; ++i) {
1022: /***********************************************************************/
1023: /* Check for convergence. */
1024: /***********************************************************************/
1026: if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
1028: /***********************************************************************/
1029: /* Compute the update. */
1030: /***********************************************************************/
1032: PetscMemcpy(e_rwrk, t_soln, sizeof(PetscReal)*t_size);
1034: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1035: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1036: #else
1037: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
1038: #endif
1040: if (0 != info) {
1041: /*********************************************************************/
1042: /* Calculation of the step failed; return the Steihaug-Toint */
1043: /* direction. */
1044: /*********************************************************************/
1046: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
1047: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1048: return(0);
1049: }
1051: /***********************************************************************/
1052: /* Modify lambda. */
1053: /***********************************************************************/
1055: norm_w = 0.;
1056: for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
1058: cg->lambda += (norm_t - cg->radius)/cg->radius * (norm_t * norm_t) / norm_w;
1060: /***********************************************************************/
1061: /* Factor T + lambda I */
1062: /***********************************************************************/
1064: for (j = 0; j < t_size; ++j) {
1065: t_diag[j] = cg->diag[j] + cg->lambda;
1066: t_offd[j] = cg->offd[j];
1067: }
1069: #if defined(PETSC_MISSING_LAPACK_PTTRF)
1070: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRF - Lapack routine is unavailable.");
1071: #else
1072: PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1073: #endif
1075: if (0 != info) {
1076: /*********************************************************************/
1077: /* Calculation of factorization failed; return the Steihaug-Toint */
1078: /* direction. */
1079: /*********************************************************************/
1081: PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n");
1082: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1083: return(0);
1084: }
1086: /***********************************************************************/
1087: /* Compute the new step and its norm. */
1088: /***********************************************************************/
1090: t_soln[0] = -cg->norm_r[0];
1091: for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
1093: #if defined(PETSC_MISSING_LAPACK_PTTRS)
1094: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"PTTRS - Lapack routine is unavailable.");
1095: #else
1096: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1097: #endif
1099: if (0 != info) {
1100: /*********************************************************************/
1101: /* Calculation of the step failed; return the Steihaug-Toint */
1102: /* direction. */
1103: /*********************************************************************/
1105: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
1106: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1107: return(0);
1108: }
1110: norm_t = 0.;
1111: for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1112: norm_t = PetscSqrtReal(norm_t);
1113: }
1115: /*************************************************************************/
1116: /* Check for convergence. */
1117: /*************************************************************************/
1119: if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1120: /***********************************************************************/
1121: /* Newton method failed to converge in iteration limit. */
1122: /***********************************************************************/
1124: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n");
1125: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1126: return(0);
1127: }
1128: }
1130: /***************************************************************************/
1131: /* Recover the norm of the direction and objective function value. */
1132: /***************************************************************************/
1134: cg->norm_d = norm_t;
1136: cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1137: for (i = 1; i < t_size - 1; ++i) {
1138: 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]);
1139: }
1140: cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);
1142: /***************************************************************************/
1143: /* Recover the direction. */
1144: /***************************************************************************/
1146: sigma = -1;
1148: /***************************************************************************/
1149: /* Start conjugate gradient method from the beginning */
1150: /***************************************************************************/
1152: VecCopy(ksp->vec_rhs, r); /* r = -grad */
1153: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1155: /***************************************************************************/
1156: /* Accumulate Q * s */
1157: /***************************************************************************/
1159: VecCopy(z, d);
1160: VecScale(d, sigma * t_soln[0] / cg->norm_r[0]);
1162: /***************************************************************************/
1163: /* Compute the first direction. */
1164: /***************************************************************************/
1166: VecCopy(z, p); /* p = z */
1167: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1168: ++ksp->its;
1170: for (i = 0; i < l_size - 1; ++i) {
1171: /*************************************************************************/
1172: /* Update the residual and direction. */
1173: /*************************************************************************/
1175: alpha = cg->alpha[i];
1176: if (alpha >= 0.0) sigma = -sigma;
1178: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1179: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1181: /*************************************************************************/
1182: /* Accumulate Q * s */
1183: /*************************************************************************/
1185: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1187: /*************************************************************************/
1188: /* Update p. */
1189: /*************************************************************************/
1191: beta = cg->beta[i];
1192: VecAYPX(p, beta, z); /* p = z + beta p */
1193: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1194: ++ksp->its;
1195: }
1197: /***************************************************************************/
1198: /* Update the residual and direction. */
1199: /***************************************************************************/
1201: alpha = cg->alpha[i];
1202: if (alpha >= 0.0) sigma = -sigma;
1204: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1205: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1207: /***************************************************************************/
1208: /* Accumulate Q * s */
1209: /***************************************************************************/
1211: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1213: /***************************************************************************/
1214: /* Set the termination reason. */
1215: /***************************************************************************/
1217: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1218: return(0);
1219: #endif
1220: }
1222: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
1223: {
1224: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1225: PetscInt max_its;
1229: /***************************************************************************/
1230: /* Determine the total maximum number of iterations. */
1231: /***************************************************************************/
1233: max_its = ksp->max_it + cg->max_lanczos_its + 1;
1235: /***************************************************************************/
1236: /* Set work vectors needed by conjugate gradient method and allocate */
1237: /* workspace for Lanczos matrix. */
1238: /***************************************************************************/
1240: KSPSetWorkVecs(ksp, 3);
1241: if (cg->diag) {
1242: PetscMemzero(cg->diag, max_its*sizeof(PetscReal));
1243: PetscMemzero(cg->offd, max_its*sizeof(PetscReal));
1244: PetscMemzero(cg->alpha, max_its*sizeof(PetscReal));
1245: PetscMemzero(cg->beta, max_its*sizeof(PetscReal));
1246: PetscMemzero(cg->norm_r, max_its*sizeof(PetscReal));
1247: } else {
1248: PetscCalloc5(max_its,&cg->diag,max_its,&cg->offd,max_its,&cg->alpha,max_its,&cg->beta,max_its,&cg->norm_r);
1249: PetscLogObjectMemory((PetscObject)ksp, 5*max_its*sizeof(PetscReal));
1250: }
1251: return(0);
1252: }
1254: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
1255: {
1256: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1260: /***************************************************************************/
1261: /* Free memory allocated for the data. */
1262: /***************************************************************************/
1264: PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1265: if (cg->alloced) {
1266: PetscFree2(cg->rwork,cg->iwork);
1267: }
1269: /***************************************************************************/
1270: /* Clear composed functions */
1271: /***************************************************************************/
1273: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",NULL);
1274: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",NULL);
1275: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",NULL);
1276: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGLTRGetMinEig_C",NULL);
1277: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGLTRGetLambda_C",NULL);
1279: /***************************************************************************/
1280: /* Destroy KSP object. */
1281: /***************************************************************************/
1282: KSPDestroyDefault(ksp);
1283: return(0);
1284: }
1286: static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
1287: {
1288: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1291: cg->radius = radius;
1292: return(0);
1293: }
1295: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1296: {
1297: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1300: *norm_d = cg->norm_d;
1301: return(0);
1302: }
1304: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1305: {
1306: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1309: *o_fcn = cg->o_fcn;
1310: return(0);
1311: }
1313: static PetscErrorCode KSPCGGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1314: {
1315: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1318: *e_min = cg->e_min;
1319: return(0);
1320: }
1322: static PetscErrorCode KSPCGGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1323: {
1324: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1327: *lambda = cg->lambda;
1328: return(0);
1329: }
1331: static PetscErrorCode KSPCGSetFromOptions_GLTR(PetscOptionItems *PetscOptionsObject,KSP ksp)
1332: {
1334: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1337: PetscOptionsHead(PetscOptionsObject,"KSP GLTR options");
1339: PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL);
1341: PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
1343: PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1344: PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1345: PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);
1347: PetscOptionsInt("-ksp_cg_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL);
1348: PetscOptionsInt("-ksp_cg_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL);
1350: PetscOptionsTail();
1351: return(0);
1352: }
1354: /*MC
1355: KSPCGGLTR - Code to run conjugate gradient method subject to a constraint
1356: on the solution norm. This is used in Trust Region methods for
1357: nonlinear equations, SNESNEWTONTR
1359: Options Database Keys:
1360: . -ksp_cg_radius <r> - Trust Region Radius
1362: Notes:
1363: This is rarely used directly
1365: Use preconditioned conjugate gradient to compute
1366: an approximate minimizer of the quadratic function
1368: q(s) = g^T * s + .5 * s^T * H * s
1370: subject to the trust region constraint
1372: || s || <= delta,
1374: where
1376: delta is the trust region radius,
1377: g is the gradient vector,
1378: H is the Hessian approximation,
1379: M is the positive definite preconditioner matrix.
1381: KSPConvergedReason may be
1382: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1383: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1384: $ other KSP converged/diverged reasons
1386: Notes:
1387: The preconditioner supplied should be symmetric and positive definite.
1389: Level: developer
1391: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCGSetRadius(), KSPCGGetNormD(), KSPCGGetObjFcn(), KSPCGGLTRGetMinEig(), KSPCGGLTRGetLambda()
1392: M*/
1394: PETSC_EXTERN PetscErrorCode KSPCreate_CGGLTR(KSP ksp)
1395: {
1397: KSPCG_GLTR *cg;
1400: PetscNewLog(ksp,&cg);
1401: cg->radius = 0.0;
1402: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1404: cg->init_pert = 1.0e-8;
1405: cg->eigen_tol = 1.0e-10;
1406: cg->newton_tol = 1.0e-6;
1408: cg->alloced = 0;
1409: cg->init_alloc = 1024;
1411: cg->max_lanczos_its = 20;
1412: cg->max_newton_its = 10;
1414: ksp->data = (void*) cg;
1415: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
1416: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
1417: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
1419: /***************************************************************************/
1420: /* Sets the functions that are associated with this data structure */
1421: /* (in C++ this is the same as defining virtual functions). */
1422: /***************************************************************************/
1424: ksp->ops->setup = KSPCGSetUp_GLTR;
1425: ksp->ops->solve = KSPCGSolve_GLTR;
1426: ksp->ops->destroy = KSPCGDestroy_GLTR;
1427: ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1428: ksp->ops->buildsolution = KSPBuildSolutionDefault;
1429: ksp->ops->buildresidual = KSPBuildResidualDefault;
1430: ksp->ops->view = 0;
1432: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",KSPCGSetRadius_GLTR);
1433: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C", KSPCGGetNormD_GLTR);
1434: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",KSPCGGetObjFcn_GLTR);
1435: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGLTRGetMinEig_C",KSPCGGLTRGetMinEig_GLTR);
1436: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGLTRGetLambda_C",KSPCGGLTRGetLambda_GLTR);
1437: return(0);
1438: }