Actual source code: gltr.c
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: KSPGLTRGetMinEig - 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: @*/
23: PetscErrorCode KSPGLTRGetMinEig(KSP ksp, PetscReal *e_min)
24: {
26: PetscUseMethod(ksp,"KSPGLTRGetMinEig_C",(KSP,PetscReal*),(ksp,e_min));
27: return 0;
28: }
30: /*@
31: KSPGLTRGetLambda - Get multiplier on trust-region constraint.
33: Not Collective
35: Input Parameters:
36: + ksp - the iterative context
37: - lambda - the multiplier
39: Level: advanced
41: @*/
42: PetscErrorCode KSPGLTRGetLambda(KSP ksp, PetscReal *lambda)
43: {
45: PetscUseMethod(ksp,"KSPGLTRGetLambda_C",(KSP,PetscReal*),(ksp,lambda));
46: return 0;
47: }
49: static PetscErrorCode KSPCGSolve_GLTR(KSP ksp)
50: {
51: #if defined(PETSC_USE_COMPLEX)
52: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP, "GLTR is not available for complex systems");
53: #else
54: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
55: PetscReal *t_soln, *t_diag, *t_offd, *e_valu, *e_vect, *e_rwrk;
56: PetscBLASInt *e_iblk, *e_splt, *e_iwrk;
58: Mat Qmat, Mmat;
59: Vec r, z, p, d;
60: PC pc;
62: PetscReal norm_r, norm_d, norm_dp1, norm_p, dMp;
63: PetscReal alpha, beta, kappa, rz, rzm1;
64: PetscReal rr, r2, piv, step;
65: PetscReal vl, vu;
66: PetscReal coef1, coef2, coef3, root1, root2, obj1, obj2;
67: PetscReal norm_t, norm_w, pert;
69: PetscInt i, j, max_cg_its, max_lanczos_its, max_newton_its, sigma;
70: PetscBLASInt t_size = 0, l_size = 0, il, iu, info;
71: PetscBLASInt nrhs, nldb;
73: PetscBLASInt e_valus=0, e_splts;
74: PetscBool diagonalscale;
76: /***************************************************************************/
77: /* Check the arguments and parameters. */
78: /***************************************************************************/
80: PCGetDiagonalScale(ksp->pc, &diagonalscale);
84: /***************************************************************************/
85: /* Get the workspace vectors and initialize variables */
86: /***************************************************************************/
88: r2 = cg->radius * cg->radius;
89: r = ksp->work[0];
90: z = ksp->work[1];
91: p = ksp->work[2];
92: d = ksp->vec_sol;
93: pc = ksp->pc;
95: PCGetOperators(pc, &Qmat, &Mmat);
97: VecGetSize(d, &max_cg_its);
98: max_cg_its = PetscMin(max_cg_its, ksp->max_it);
99: max_lanczos_its = cg->max_lanczos_its;
100: max_newton_its = cg->max_newton_its;
101: ksp->its = 0;
103: /***************************************************************************/
104: /* Initialize objective function direction, and minimum eigenvalue. */
105: /***************************************************************************/
107: cg->o_fcn = 0.0;
109: VecSet(d, 0.0); /* d = 0 */
110: cg->norm_d = 0.0;
112: cg->e_min = 0.0;
113: cg->lambda = 0.0;
115: /***************************************************************************/
116: /* The first phase of GLTR performs a standard conjugate gradient method, */
117: /* but stores the values required for the Lanczos matrix. We switch to */
118: /* the Lanczos when the conjugate gradient method breaks down. Check the */
119: /* right-hand side for numerical problems. The check for not-a-number and */
120: /* infinite values need be performed only once. */
121: /***************************************************************************/
123: VecCopy(ksp->vec_rhs, r); /* r = -grad */
124: VecDot(r, r, &rr); /* rr = r^T r */
125: KSPCheckDot(ksp,rr);
127: /***************************************************************************/
128: /* Check the preconditioner for numerical problems and for positive */
129: /* definiteness. The check for not-a-number and infinite values need be */
130: /* performed only once. */
131: /***************************************************************************/
133: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
134: VecDot(r, z, &rz); /* rz = r^T inv(M) r */
135: if (PetscIsInfOrNanScalar(rz)) {
136: /*************************************************************************/
137: /* The preconditioner contains not-a-number or an infinite value. */
138: /* Return the gradient direction intersected with the trust region. */
139: /*************************************************************************/
141: ksp->reason = KSP_DIVERGED_NANORINF;
142: PetscInfo(ksp, "KSPCGSolve_GLTR: bad preconditioner: rz=%g\n", (double)rz);
144: if (cg->radius) {
145: if (r2 >= rr) {
146: alpha = 1.0;
147: cg->norm_d = PetscSqrtReal(rr);
148: } else {
149: alpha = PetscSqrtReal(r2 / rr);
150: cg->norm_d = cg->radius;
151: }
153: VecAXPY(d, alpha, r); /* d = d + alpha r */
155: /***********************************************************************/
156: /* Compute objective function. */
157: /***********************************************************************/
159: KSP_MatMult(ksp, Qmat, d, z);
160: VecAYPX(z, -0.5, ksp->vec_rhs);
161: VecDot(d, z, &cg->o_fcn);
162: cg->o_fcn = -cg->o_fcn;
163: ++ksp->its;
164: }
165: return 0;
166: }
168: if (rz < 0.0) {
169: /*************************************************************************/
170: /* The preconditioner is indefinite. Because this is the first */
171: /* and we do not have a direction yet, we use the gradient step. Note */
172: /* that we cannot use the preconditioned norm when computing the step */
173: /* because the matrix is indefinite. */
174: /*************************************************************************/
176: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
177: PetscInfo(ksp, "KSPCGSolve_GLTR: indefinite preconditioner: rz=%g\n", (double)rz);
179: if (cg->radius) {
180: if (r2 >= rr) {
181: alpha = 1.0;
182: cg->norm_d = PetscSqrtReal(rr);
183: } else {
184: alpha = PetscSqrtReal(r2 / rr);
185: cg->norm_d = cg->radius;
186: }
188: VecAXPY(d, alpha, r); /* d = d + alpha r */
190: /***********************************************************************/
191: /* Compute objective function. */
192: /***********************************************************************/
194: KSP_MatMult(ksp, Qmat, d, z);
195: VecAYPX(z, -0.5, ksp->vec_rhs);
196: VecDot(d, z, &cg->o_fcn);
197: cg->o_fcn = -cg->o_fcn;
198: ++ksp->its;
199: }
200: return 0;
201: }
203: /***************************************************************************/
204: /* As far as we know, the preconditioner is positive semidefinite. */
205: /* Compute and log the residual. Check convergence because this */
206: /* initializes things, but do not terminate until at least one conjugate */
207: /* gradient iteration has been performed. */
208: /***************************************************************************/
210: cg->norm_r[0] = PetscSqrtReal(rz); /* norm_r = |r|_M */
212: switch (ksp->normtype) {
213: case KSP_NORM_PRECONDITIONED:
214: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
215: break;
217: case KSP_NORM_UNPRECONDITIONED:
218: norm_r = PetscSqrtReal(rr); /* norm_r = |r| */
219: break;
221: case KSP_NORM_NATURAL:
222: norm_r = cg->norm_r[0]; /* norm_r = |r|_M */
223: break;
225: default:
226: norm_r = 0.0;
227: break;
228: }
230: KSPLogResidualHistory(ksp, norm_r);
231: KSPMonitor(ksp, ksp->its, norm_r);
232: ksp->rnorm = norm_r;
234: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
236: /***************************************************************************/
237: /* Compute the first direction and update the iteration. */
238: /***************************************************************************/
240: VecCopy(z, p); /* p = z */
241: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
242: ++ksp->its;
244: /***************************************************************************/
245: /* Check the matrix for numerical problems. */
246: /***************************************************************************/
248: VecDot(p, z, &kappa); /* kappa = p^T Q p */
249: if (PetscIsInfOrNanScalar(kappa)) {
250: /*************************************************************************/
251: /* The matrix produced not-a-number or an infinite value. In this case, */
252: /* we must stop and use the gradient direction. This condition need */
253: /* only be checked once. */
254: /*************************************************************************/
256: ksp->reason = KSP_DIVERGED_NANORINF;
257: PetscInfo(ksp, "KSPCGSolve_GLTR: bad matrix: kappa=%g\n", (double)kappa);
259: if (cg->radius) {
260: if (r2 >= rr) {
261: alpha = 1.0;
262: cg->norm_d = PetscSqrtReal(rr);
263: } else {
264: alpha = PetscSqrtReal(r2 / rr);
265: cg->norm_d = cg->radius;
266: }
268: VecAXPY(d, alpha, r); /* d = d + alpha r */
270: /***********************************************************************/
271: /* Compute objective function. */
272: /***********************************************************************/
274: KSP_MatMult(ksp, Qmat, d, z);
275: VecAYPX(z, -0.5, ksp->vec_rhs);
276: VecDot(d, z, &cg->o_fcn);
277: cg->o_fcn = -cg->o_fcn;
278: ++ksp->its;
279: }
280: return 0;
281: }
283: /***************************************************************************/
284: /* Initialize variables for calculating the norm of the direction and for */
285: /* the Lanczos tridiagonal matrix. Note that we track the diagonal value */
286: /* of the Cholesky factorization of the Lanczos matrix in order to */
287: /* determine when negative curvature is encountered. */
288: /***************************************************************************/
290: dMp = 0.0;
291: norm_d = 0.0;
292: switch (cg->dtype) {
293: case GLTR_PRECONDITIONED_DIRECTION:
294: norm_p = rz;
295: break;
297: default:
298: VecDot(p, p, &norm_p);
299: break;
300: }
302: cg->diag[t_size] = kappa / rz;
303: cg->offd[t_size] = 0.0;
304: ++t_size;
306: piv = 1.0;
308: /***************************************************************************/
309: /* Check for breakdown of the conjugate gradient method; this occurs when */
310: /* kappa is zero. */
311: /***************************************************************************/
313: if (PetscAbsReal(kappa) <= 0.0) {
314: /*************************************************************************/
315: /* The curvature is zero. In this case, we must stop and use follow */
316: /* the direction of negative curvature since the Lanczos matrix is zero. */
317: /*************************************************************************/
319: ksp->reason = KSP_DIVERGED_BREAKDOWN;
320: PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: kappa=%g\n", (double)kappa);
322: if (cg->radius && norm_p > 0.0) {
323: /***********************************************************************/
324: /* Follow direction of negative curvature to the boundary of the */
325: /* trust region. */
326: /***********************************************************************/
328: step = PetscSqrtReal(r2 / norm_p);
329: cg->norm_d = cg->radius;
331: VecAXPY(d, step, p); /* d = d + step p */
333: /***********************************************************************/
334: /* Update objective function. */
335: /***********************************************************************/
337: cg->o_fcn += step * (0.5 * step * kappa - rz);
338: } else if (cg->radius) {
339: /***********************************************************************/
340: /* The norm of the preconditioned direction is zero; use the gradient */
341: /* step. */
342: /***********************************************************************/
344: if (r2 >= rr) {
345: alpha = 1.0;
346: cg->norm_d = PetscSqrtReal(rr);
347: } else {
348: alpha = PetscSqrtReal(r2 / rr);
349: cg->norm_d = cg->radius;
350: }
352: VecAXPY(d, alpha, r); /* d = d + alpha r */
354: /***********************************************************************/
355: /* Compute objective function. */
356: /***********************************************************************/
358: KSP_MatMult(ksp, Qmat, d, z);
359: VecAYPX(z, -0.5, ksp->vec_rhs);
360: VecDot(d, z, &cg->o_fcn);
361: cg->o_fcn = -cg->o_fcn;
362: ++ksp->its;
363: }
364: return 0;
365: }
367: /***************************************************************************/
368: /* Begin the first part of the GLTR algorithm which runs the conjugate */
369: /* gradient method until either the problem is solved, we encounter the */
370: /* boundary of the trust region, or the conjugate gradient method breaks */
371: /* down. */
372: /***************************************************************************/
374: while (1) {
375: /*************************************************************************/
376: /* Know that kappa is nonzero, because we have not broken down, so we */
377: /* can compute the steplength. */
378: /*************************************************************************/
380: alpha = rz / kappa;
381: cg->alpha[l_size] = alpha;
383: /*************************************************************************/
384: /* Compute the diagonal value of the Cholesky factorization of the */
385: /* Lanczos matrix and check to see if the Lanczos matrix is indefinite. */
386: /* This indicates a direction of negative curvature. */
387: /*************************************************************************/
389: piv = cg->diag[l_size] - cg->offd[l_size]*cg->offd[l_size] / piv;
390: if (piv <= 0.0) {
391: /***********************************************************************/
392: /* In this case, the matrix is indefinite and we have encountered */
393: /* a direction of negative curvature. Follow the direction to the */
394: /* boundary of the trust region. */
395: /***********************************************************************/
397: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
398: PetscInfo(ksp, "KSPCGSolve_GLTR: negative curvature: kappa=%g\n", (double)kappa);
400: if (cg->radius && norm_p > 0.0) {
401: /*********************************************************************/
402: /* Follow direction of negative curvature to boundary. */
403: /*********************************************************************/
405: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
406: cg->norm_d = cg->radius;
408: VecAXPY(d, step, p); /* d = d + step p */
410: /*********************************************************************/
411: /* Update objective function. */
412: /*********************************************************************/
414: cg->o_fcn += step * (0.5 * step * kappa - rz);
415: } else if (cg->radius) {
416: /*********************************************************************/
417: /* The norm of the direction is zero; there is nothing to follow. */
418: /*********************************************************************/
419: }
420: break;
421: }
423: /*************************************************************************/
424: /* Compute the steplength and check for intersection with the trust */
425: /* region. */
426: /*************************************************************************/
428: norm_dp1 = norm_d + alpha*(2.0*dMp + alpha*norm_p);
429: if (cg->radius && norm_dp1 >= r2) {
430: /***********************************************************************/
431: /* In this case, the matrix is positive definite as far as we know. */
432: /* However, the full step goes beyond the trust region. */
433: /***********************************************************************/
435: ksp->reason = KSP_CONVERGED_CG_CONSTRAINED;
436: PetscInfo(ksp, "KSPCGSolve_GLTR: constrained step: radius=%g\n", (double)cg->radius);
438: if (norm_p > 0.0) {
439: /*********************************************************************/
440: /* Follow the direction to the boundary of the trust region. */
441: /*********************************************************************/
443: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
444: cg->norm_d = cg->radius;
446: VecAXPY(d, step, p); /* d = d + step p */
448: /*********************************************************************/
449: /* Update objective function. */
450: /*********************************************************************/
452: cg->o_fcn += step * (0.5 * step * kappa - rz);
453: } else {
454: /*********************************************************************/
455: /* The norm of the direction is zero; there is nothing to follow. */
456: /*********************************************************************/
457: }
458: break;
459: }
461: /*************************************************************************/
462: /* Now we can update the direction and residual. */
463: /*************************************************************************/
465: VecAXPY(d, alpha, p); /* d = d + alpha p */
466: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
467: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
469: switch (cg->dtype) {
470: case GLTR_PRECONDITIONED_DIRECTION:
471: norm_d = norm_dp1;
472: break;
474: default:
475: VecDot(d, d, &norm_d);
476: break;
477: }
478: cg->norm_d = PetscSqrtReal(norm_d);
480: /*************************************************************************/
481: /* Update objective function. */
482: /*************************************************************************/
484: cg->o_fcn -= 0.5 * alpha * rz;
486: /*************************************************************************/
487: /* Check that the preconditioner appears positive semidefinite. */
488: /*************************************************************************/
490: rzm1 = rz;
491: VecDot(r, z, &rz); /* rz = r^T z */
492: if (rz < 0.0) {
493: /***********************************************************************/
494: /* The preconditioner is indefinite. */
495: /***********************************************************************/
497: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
498: PetscInfo(ksp, "KSPCGSolve_GLTR: cg indefinite preconditioner: rz=%g\n", (double)rz);
499: break;
500: }
502: /*************************************************************************/
503: /* As far as we know, the preconditioner is positive semidefinite. */
504: /* Compute the residual and check for convergence. */
505: /*************************************************************************/
507: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
509: switch (ksp->normtype) {
510: case KSP_NORM_PRECONDITIONED:
511: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
512: break;
514: case KSP_NORM_UNPRECONDITIONED:
515: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
516: break;
518: case KSP_NORM_NATURAL:
519: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
520: break;
522: default:
523: norm_r = 0.0;
524: break;
525: }
527: KSPLogResidualHistory(ksp, norm_r);
528: KSPMonitor(ksp, ksp->its, norm_r);
529: ksp->rnorm = norm_r;
531: (*ksp->converged)(ksp, ksp->its, norm_r, &ksp->reason, ksp->cnvP);
532: if (ksp->reason) {
533: /***********************************************************************/
534: /* The method has converged. */
535: /***********************************************************************/
537: PetscInfo(ksp, "KSPCGSolve_GLTR: cg truncated step: rnorm=%g, radius=%g\n", (double)norm_r, (double)cg->radius);
538: break;
539: }
541: /*************************************************************************/
542: /* We have not converged yet. Check for breakdown. */
543: /*************************************************************************/
545: beta = rz / rzm1;
546: if (PetscAbsReal(beta) <= 0.0) {
547: /***********************************************************************/
548: /* Conjugate gradients has broken down. */
549: /***********************************************************************/
551: ksp->reason = KSP_DIVERGED_BREAKDOWN;
552: PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n", (double)beta);
553: break;
554: }
556: /*************************************************************************/
557: /* Check iteration limit. */
558: /*************************************************************************/
560: if (ksp->its >= max_cg_its) {
561: ksp->reason = KSP_DIVERGED_ITS;
562: PetscInfo(ksp, "KSPCGSolve_GLTR: iterlim: its=%D\n", ksp->its);
563: break;
564: }
566: /*************************************************************************/
567: /* Update p and the norms. */
568: /*************************************************************************/
570: cg->beta[l_size] = beta;
571: VecAYPX(p, beta, z); /* p = z + beta p */
573: switch (cg->dtype) {
574: case GLTR_PRECONDITIONED_DIRECTION:
575: dMp = beta*(dMp + alpha*norm_p);
576: norm_p = beta*(rzm1 + beta*norm_p);
577: break;
579: default:
580: VecDot(d, p, &dMp);
581: VecDot(p, p, &norm_p);
582: break;
583: }
585: /*************************************************************************/
586: /* Compute the new direction and update the iteration. */
587: /*************************************************************************/
589: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
590: VecDot(p, z, &kappa); /* kappa = p^T Q p */
591: ++ksp->its;
593: /*************************************************************************/
594: /* Update the Lanczos tridiagonal matrix. */
595: /*************************************************************************/
597: ++l_size;
598: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
599: cg->diag[t_size] = kappa / rz + beta / alpha;
600: ++t_size;
602: /*************************************************************************/
603: /* Check for breakdown of the conjugate gradient method; this occurs */
604: /* when kappa is zero. */
605: /*************************************************************************/
607: if (PetscAbsReal(kappa) <= 0.0) {
608: /***********************************************************************/
609: /* The method breaks down; move along the direction as if the matrix */
610: /* were indefinite. */
611: /***********************************************************************/
613: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
614: PetscInfo(ksp, "KSPCGSolve_GLTR: cg breakdown: kappa=%g\n", (double)kappa);
616: if (cg->radius && norm_p > 0.0) {
617: /*********************************************************************/
618: /* Follow direction to boundary. */
619: /*********************************************************************/
621: step = (PetscSqrtReal(dMp*dMp+norm_p*(r2-norm_d))-dMp)/norm_p;
622: cg->norm_d = cg->radius;
624: VecAXPY(d, step, p); /* d = d + step p */
626: /*********************************************************************/
627: /* Update objective function. */
628: /*********************************************************************/
630: cg->o_fcn += step * (0.5 * step * kappa - rz);
631: } else if (cg->radius) {
632: /*********************************************************************/
633: /* The norm of the direction is zero; there is nothing to follow. */
634: /*********************************************************************/
635: }
636: break;
637: }
638: }
640: /***************************************************************************/
641: /* Check to see if we need to continue with the Lanczos method. */
642: /***************************************************************************/
644: if (!cg->radius) {
645: /*************************************************************************/
646: /* There is no radius. Therefore, we cannot move along the boundary. */
647: /*************************************************************************/
648: return 0;
649: }
651: if (KSP_CONVERGED_CG_NEG_CURVE != ksp->reason) {
652: /*************************************************************************/
653: /* The method either converged to an interior point, hit the boundary of */
654: /* the trust-region without encountering a direction of negative */
655: /* curvature or the iteration limit was reached. */
656: /*************************************************************************/
657: return 0;
658: }
660: /***************************************************************************/
661: /* Switch to contructing the Lanczos basis by way of the conjugate */
662: /* directions. */
663: /***************************************************************************/
665: for (i = 0; i < max_lanczos_its; ++i) {
666: /*************************************************************************/
667: /* Check for breakdown of the conjugate gradient method; this occurs */
668: /* when kappa is zero. */
669: /*************************************************************************/
671: if (PetscAbsReal(kappa) <= 0.0) {
672: ksp->reason = KSP_DIVERGED_BREAKDOWN;
673: PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos breakdown: kappa=%g\n", (double)kappa);
674: break;
675: }
677: /*************************************************************************/
678: /* Update the direction and residual. */
679: /*************************************************************************/
681: alpha = rz / kappa;
682: cg->alpha[l_size] = alpha;
684: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
685: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
687: /*************************************************************************/
688: /* Check that the preconditioner appears positive semidefinite. */
689: /*************************************************************************/
691: rzm1 = rz;
692: VecDot(r, z, &rz); /* rz = r^T z */
693: if (rz < 0.0) {
694: /***********************************************************************/
695: /* The preconditioner is indefinite. */
696: /***********************************************************************/
698: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
699: PetscInfo(ksp, "KSPCGSolve_GLTR: lanczos indefinite preconditioner: rz=%g\n", (double)rz);
700: break;
701: }
703: /*************************************************************************/
704: /* As far as we know, the preconditioner is positive definite. Compute */
705: /* the residual. Do NOT check for convergence. */
706: /*************************************************************************/
708: cg->norm_r[l_size+1] = PetscSqrtReal(rz); /* norm_r = |r|_M */
710: switch (ksp->normtype) {
711: case KSP_NORM_PRECONDITIONED:
712: VecNorm(z, NORM_2, &norm_r); /* norm_r = |z| */
713: break;
715: case KSP_NORM_UNPRECONDITIONED:
716: VecNorm(r, NORM_2, &norm_r); /* norm_r = |r| */
717: break;
719: case KSP_NORM_NATURAL:
720: norm_r = cg->norm_r[l_size+1]; /* norm_r = |r|_M */
721: break;
723: default:
724: norm_r = 0.0;
725: break;
726: }
728: KSPLogResidualHistory(ksp, norm_r);
729: KSPMonitor(ksp, ksp->its, norm_r);
730: ksp->rnorm = norm_r;
732: /*************************************************************************/
733: /* Check for breakdown. */
734: /*************************************************************************/
736: beta = rz / rzm1;
737: if (PetscAbsReal(beta) <= 0.0) {
738: /***********************************************************************/
739: /* Conjugate gradients has broken down. */
740: /***********************************************************************/
742: ksp->reason = KSP_DIVERGED_BREAKDOWN;
743: PetscInfo(ksp, "KSPCGSolve_GLTR: breakdown: beta=%g\n",(double) beta);
744: break;
745: }
747: /*************************************************************************/
748: /* Update p and the norms. */
749: /*************************************************************************/
751: cg->beta[l_size] = beta;
752: VecAYPX(p, beta, z); /* p = z + beta p */
754: /*************************************************************************/
755: /* Compute the new direction and update the iteration. */
756: /*************************************************************************/
758: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
759: VecDot(p, z, &kappa); /* kappa = p^T Q p */
760: ++ksp->its;
762: /*************************************************************************/
763: /* Update the Lanczos tridiagonal matrix. */
764: /*************************************************************************/
766: ++l_size;
767: cg->offd[t_size] = PetscSqrtReal(beta) / PetscAbsReal(alpha);
768: cg->diag[t_size] = kappa / rz + beta / alpha;
769: ++t_size;
770: }
772: /***************************************************************************/
773: /* We have the Lanczos basis, solve the tridiagonal trust-region problem */
774: /* to obtain the Lanczos direction. We know that the solution lies on */
775: /* the boundary of the trust region. We start by checking that the */
776: /* workspace allocated is large enough. */
777: /***************************************************************************/
778: /* Note that the current version only computes the solution by using the */
779: /* preconditioned direction. Need to think about how to do the */
780: /* unpreconditioned direction calculation. */
781: /***************************************************************************/
783: if (t_size > cg->alloced) {
784: if (cg->alloced) {
785: PetscFree(cg->rwork);
786: PetscFree(cg->iwork);
787: cg->alloced += cg->init_alloc;
788: } else {
789: cg->alloced = cg->init_alloc;
790: }
792: while (t_size > cg->alloced) {
793: cg->alloced += cg->init_alloc;
794: }
796: cg->alloced = PetscMin(cg->alloced, t_size);
797: PetscMalloc2(10*cg->alloced, &cg->rwork,5*cg->alloced, &cg->iwork);
798: }
800: /***************************************************************************/
801: /* Set up the required vectors. */
802: /***************************************************************************/
804: t_soln = cg->rwork + 0*t_size; /* Solution */
805: t_diag = cg->rwork + 1*t_size; /* Diagonal of T */
806: t_offd = cg->rwork + 2*t_size; /* Off-diagonal of T */
807: e_valu = cg->rwork + 3*t_size; /* Eigenvalues of T */
808: e_vect = cg->rwork + 4*t_size; /* Eigenvector of T */
809: e_rwrk = cg->rwork + 5*t_size; /* Eigen workspace */
811: e_iblk = cg->iwork + 0*t_size; /* Eigen blocks */
812: e_splt = cg->iwork + 1*t_size; /* Eigen splits */
813: e_iwrk = cg->iwork + 2*t_size; /* Eigen workspace */
815: /***************************************************************************/
816: /* Compute the minimum eigenvalue of T. */
817: /***************************************************************************/
819: vl = 0.0;
820: vu = 0.0;
821: il = 1;
822: iu = 1;
824: 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));
826: if ((0 != info) || (1 != e_valus)) {
827: /*************************************************************************/
828: /* Calculation of the minimum eigenvalue failed. Return the */
829: /* Steihaug-Toint direction. */
830: /*************************************************************************/
832: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvalue.\n");
833: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
834: return 0;
835: }
837: cg->e_min = e_valu[0];
839: /***************************************************************************/
840: /* Compute the initial value of lambda to make (T + lamba I) positive */
841: /* definite. */
842: /***************************************************************************/
844: pert = cg->init_pert;
845: if (e_valu[0] < 0.0) cg->lambda = pert - e_valu[0];
847: while (1) {
848: for (i = 0; i < t_size; ++i) {
849: t_diag[i] = cg->diag[i] + cg->lambda;
850: t_offd[i] = cg->offd[i];
851: }
853: PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
855: if (0 == info) break;
857: pert += pert;
858: cg->lambda = cg->lambda * (1.0 + pert) + pert;
859: }
861: /***************************************************************************/
862: /* Compute the initial step and its norm. */
863: /***************************************************************************/
865: nrhs = 1;
866: nldb = t_size;
868: t_soln[0] = -cg->norm_r[0];
869: for (i = 1; i < t_size; ++i) t_soln[i] = 0.0;
871: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
873: if (0 != info) {
874: /*************************************************************************/
875: /* Calculation of the initial step failed; return the Steihaug-Toint */
876: /* direction. */
877: /*************************************************************************/
879: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
880: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
881: return 0;
882: }
884: norm_t = 0.;
885: for (i = 0; i < t_size; ++i) norm_t += t_soln[i] * t_soln[i];
886: norm_t = PetscSqrtReal(norm_t);
888: /***************************************************************************/
889: /* Determine the case we are in. */
890: /***************************************************************************/
892: if (norm_t <= cg->radius) {
893: /*************************************************************************/
894: /* The step is within the trust region; check if we are in the hard case */
895: /* and need to move to the boundary by following a direction of negative */
896: /* curvature. */
897: /*************************************************************************/
899: if ((e_valu[0] <= 0.0) && (norm_t < cg->radius)) {
900: /***********************************************************************/
901: /* This is the hard case; compute the eigenvector associated with the */
902: /* minimum eigenvalue and move along this direction to the boundary. */
903: /***********************************************************************/
905: 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));
907: if (0 != info) {
908: /*********************************************************************/
909: /* Calculation of the minimum eigenvalue failed. Return the */
910: /* Steihaug-Toint direction. */
911: /*********************************************************************/
913: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute eigenvector.\n");
914: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
915: return 0;
916: }
918: coef1 = 0.0;
919: coef2 = 0.0;
920: coef3 = -cg->radius * cg->radius;
921: for (i = 0; i < t_size; ++i) {
922: coef1 += e_vect[i] * e_vect[i];
923: coef2 += e_vect[i] * t_soln[i];
924: coef3 += t_soln[i] * t_soln[i];
925: }
927: coef3 = PetscSqrtReal(coef2 * coef2 - coef1 * coef3);
928: root1 = (-coef2 + coef3) / coef1;
929: root2 = (-coef2 - coef3) / coef1;
931: /***********************************************************************/
932: /* Compute objective value for (t_soln + root1 * e_vect) */
933: /***********************************************************************/
935: for (i = 0; i < t_size; ++i) {
936: e_rwrk[i] = t_soln[i] + root1 * e_vect[i];
937: }
939: obj1 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
940: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
941: for (i = 1; i < t_size - 1; ++i) {
942: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
943: cg->diag[i]*e_rwrk[i]+
944: cg->offd[i+1]*e_rwrk[i+1]);
945: }
946: obj1 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
947: cg->diag[i]*e_rwrk[i]);
949: /***********************************************************************/
950: /* Compute objective value for (t_soln + root2 * e_vect) */
951: /***********************************************************************/
953: for (i = 0; i < t_size; ++i) {
954: e_rwrk[i] = t_soln[i] + root2 * e_vect[i];
955: }
957: obj2 = e_rwrk[0]*(0.5*(cg->diag[0]*e_rwrk[0]+
958: cg->offd[1]*e_rwrk[1])+cg->norm_r[0]);
959: for (i = 1; i < t_size - 1; ++i) {
960: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
961: cg->diag[i]*e_rwrk[i]+
962: cg->offd[i+1]*e_rwrk[i+1]);
963: }
964: obj2 += 0.5*e_rwrk[i]*(cg->offd[i]*e_rwrk[i-1]+
965: cg->diag[i]*e_rwrk[i]);
967: /***********************************************************************/
968: /* Choose the point with the best objective function value. */
969: /***********************************************************************/
971: if (obj1 <= obj2) {
972: for (i = 0; i < t_size; ++i) {
973: t_soln[i] += root1 * e_vect[i];
974: }
975: }
976: else {
977: for (i = 0; i < t_size; ++i) {
978: t_soln[i] += root2 * e_vect[i];
979: }
980: }
981: } else {
982: /***********************************************************************/
983: /* The matrix is positive definite or there was no room to move; the */
984: /* solution is already contained in t_soln. */
985: /***********************************************************************/
986: }
987: } else {
988: /*************************************************************************/
989: /* The step is outside the trust-region. Compute the correct value for */
990: /* lambda by performing Newton's method. */
991: /*************************************************************************/
993: for (i = 0; i < max_newton_its; ++i) {
994: /***********************************************************************/
995: /* Check for convergence. */
996: /***********************************************************************/
998: if (PetscAbsReal(norm_t - cg->radius) <= cg->newton_tol * cg->radius) break;
1000: /***********************************************************************/
1001: /* Compute the update. */
1002: /***********************************************************************/
1004: PetscArraycpy(e_rwrk, t_soln, t_size);
1006: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, e_rwrk, &nldb, &info));
1008: if (0 != info) {
1009: /*********************************************************************/
1010: /* Calculation of the step failed; return the Steihaug-Toint */
1011: /* direction. */
1012: /*********************************************************************/
1014: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
1015: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1016: return 0;
1017: }
1019: /***********************************************************************/
1020: /* Modify lambda. */
1021: /***********************************************************************/
1023: norm_w = 0.;
1024: for (j = 0; j < t_size; ++j) norm_w += t_soln[j] * e_rwrk[j];
1026: cg->lambda += (norm_t - cg->radius)/cg->radius * (norm_t * norm_t) / norm_w;
1028: /***********************************************************************/
1029: /* Factor T + lambda I */
1030: /***********************************************************************/
1032: for (j = 0; j < t_size; ++j) {
1033: t_diag[j] = cg->diag[j] + cg->lambda;
1034: t_offd[j] = cg->offd[j];
1035: }
1037: PetscStackCallBLAS("LAPACKpttrf",LAPACKpttrf_(&t_size, t_diag, t_offd + 1, &info));
1039: if (0 != info) {
1040: /*********************************************************************/
1041: /* Calculation of factorization failed; return the Steihaug-Toint */
1042: /* direction. */
1043: /*********************************************************************/
1045: PetscInfo(ksp, "KSPCGSolve_GLTR: factorization failed.\n");
1046: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1047: return 0;
1048: }
1050: /***********************************************************************/
1051: /* Compute the new step and its norm. */
1052: /***********************************************************************/
1054: t_soln[0] = -cg->norm_r[0];
1055: for (j = 1; j < t_size; ++j) t_soln[j] = 0.0;
1057: PetscStackCallBLAS("LAPACKpttrs",LAPACKpttrs_(&t_size, &nrhs, t_diag, t_offd + 1, t_soln, &nldb, &info));
1059: if (0 != info) {
1060: /*********************************************************************/
1061: /* Calculation of the step failed; return the Steihaug-Toint */
1062: /* direction. */
1063: /*********************************************************************/
1065: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to compute step.\n");
1066: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1067: return 0;
1068: }
1070: norm_t = 0.;
1071: for (j = 0; j < t_size; ++j) norm_t += t_soln[j] * t_soln[j];
1072: norm_t = PetscSqrtReal(norm_t);
1073: }
1075: /*************************************************************************/
1076: /* Check for convergence. */
1077: /*************************************************************************/
1079: if (PetscAbsReal(norm_t - cg->radius) > cg->newton_tol * cg->radius) {
1080: /***********************************************************************/
1081: /* Newton method failed to converge in iteration limit. */
1082: /***********************************************************************/
1084: PetscInfo(ksp, "KSPCGSolve_GLTR: failed to converge.\n");
1085: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1086: return 0;
1087: }
1088: }
1090: /***************************************************************************/
1091: /* Recover the norm of the direction and objective function value. */
1092: /***************************************************************************/
1094: cg->norm_d = norm_t;
1096: cg->o_fcn = t_soln[0]*(0.5*(cg->diag[0]*t_soln[0]+cg->offd[1]*t_soln[1])+cg->norm_r[0]);
1097: for (i = 1; i < t_size - 1; ++i) {
1098: 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]);
1099: }
1100: cg->o_fcn += 0.5*t_soln[i]*(cg->offd[i]*t_soln[i-1]+cg->diag[i]*t_soln[i]);
1102: /***************************************************************************/
1103: /* Recover the direction. */
1104: /***************************************************************************/
1106: sigma = -1;
1108: /***************************************************************************/
1109: /* Start conjugate gradient method from the beginning */
1110: /***************************************************************************/
1112: VecCopy(ksp->vec_rhs, r); /* r = -grad */
1113: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1115: /***************************************************************************/
1116: /* Accumulate Q * s */
1117: /***************************************************************************/
1119: VecCopy(z, d);
1120: VecScale(d, sigma * t_soln[0] / cg->norm_r[0]);
1122: /***************************************************************************/
1123: /* Compute the first direction. */
1124: /***************************************************************************/
1126: VecCopy(z, p); /* p = z */
1127: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1128: ++ksp->its;
1130: for (i = 0; i < l_size - 1; ++i) {
1131: /*************************************************************************/
1132: /* Update the residual and direction. */
1133: /*************************************************************************/
1135: alpha = cg->alpha[i];
1136: if (alpha >= 0.0) sigma = -sigma;
1138: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1139: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1141: /*************************************************************************/
1142: /* Accumulate Q * s */
1143: /*************************************************************************/
1145: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1147: /*************************************************************************/
1148: /* Update p. */
1149: /*************************************************************************/
1151: beta = cg->beta[i];
1152: VecAYPX(p, beta, z); /* p = z + beta p */
1153: KSP_MatMult(ksp, Qmat, p, z); /* z = Q * p */
1154: ++ksp->its;
1155: }
1157: /***************************************************************************/
1158: /* Update the residual and direction. */
1159: /***************************************************************************/
1161: alpha = cg->alpha[i];
1162: if (alpha >= 0.0) sigma = -sigma;
1164: VecAXPY(r, -alpha, z); /* r = r - alpha Q p */
1165: KSP_PCApply(ksp, r, z); /* z = inv(M) r */
1167: /***************************************************************************/
1168: /* Accumulate Q * s */
1169: /***************************************************************************/
1171: VecAXPY(d, sigma * t_soln[i+1] / cg->norm_r[i+1], z);
1173: /***************************************************************************/
1174: /* Set the termination reason. */
1175: /***************************************************************************/
1177: ksp->reason = KSP_CONVERGED_CG_NEG_CURVE;
1178: return 0;
1179: #endif
1180: }
1182: static PetscErrorCode KSPCGSetUp_GLTR(KSP ksp)
1183: {
1184: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1185: PetscInt max_its;
1187: /***************************************************************************/
1188: /* Determine the total maximum number of iterations. */
1189: /***************************************************************************/
1191: max_its = ksp->max_it + cg->max_lanczos_its + 1;
1193: /***************************************************************************/
1194: /* Set work vectors needed by conjugate gradient method and allocate */
1195: /* workspace for Lanczos matrix. */
1196: /***************************************************************************/
1198: KSPSetWorkVecs(ksp, 3);
1199: if (cg->diag) {
1200: PetscArrayzero(cg->diag, max_its);
1201: PetscArrayzero(cg->offd, max_its);
1202: PetscArrayzero(cg->alpha, max_its);
1203: PetscArrayzero(cg->beta, max_its);
1204: PetscArrayzero(cg->norm_r, max_its);
1205: } else {
1206: PetscCalloc5(max_its,&cg->diag,max_its,&cg->offd,max_its,&cg->alpha,max_its,&cg->beta,max_its,&cg->norm_r);
1207: PetscLogObjectMemory((PetscObject)ksp, 5*max_its*sizeof(PetscReal));
1208: }
1209: return 0;
1210: }
1212: static PetscErrorCode KSPCGDestroy_GLTR(KSP ksp)
1213: {
1214: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1216: /***************************************************************************/
1217: /* Free memory allocated for the data. */
1218: /***************************************************************************/
1220: PetscFree5(cg->diag,cg->offd,cg->alpha,cg->beta,cg->norm_r);
1221: if (cg->alloced) {
1222: PetscFree2(cg->rwork,cg->iwork);
1223: }
1225: /***************************************************************************/
1226: /* Clear composed functions */
1227: /***************************************************************************/
1229: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",NULL);
1230: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C",NULL);
1231: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",NULL);
1232: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",NULL);
1233: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",NULL);
1235: /***************************************************************************/
1236: /* Destroy KSP object. */
1237: /***************************************************************************/
1238: KSPDestroyDefault(ksp);
1239: return 0;
1240: }
1242: static PetscErrorCode KSPCGSetRadius_GLTR(KSP ksp, PetscReal radius)
1243: {
1244: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1246: cg->radius = radius;
1247: return 0;
1248: }
1250: static PetscErrorCode KSPCGGetNormD_GLTR(KSP ksp, PetscReal *norm_d)
1251: {
1252: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1254: *norm_d = cg->norm_d;
1255: return 0;
1256: }
1258: static PetscErrorCode KSPCGGetObjFcn_GLTR(KSP ksp, PetscReal *o_fcn)
1259: {
1260: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1262: *o_fcn = cg->o_fcn;
1263: return 0;
1264: }
1266: static PetscErrorCode KSPGLTRGetMinEig_GLTR(KSP ksp, PetscReal *e_min)
1267: {
1268: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1270: *e_min = cg->e_min;
1271: return 0;
1272: }
1274: static PetscErrorCode KSPGLTRGetLambda_GLTR(KSP ksp, PetscReal *lambda)
1275: {
1276: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1278: *lambda = cg->lambda;
1279: return 0;
1280: }
1282: static PetscErrorCode KSPCGSetFromOptions_GLTR(PetscOptionItems *PetscOptionsObject,KSP ksp)
1283: {
1284: KSPCG_GLTR *cg = (KSPCG_GLTR*)ksp->data;
1286: PetscOptionsHead(PetscOptionsObject,"KSP GLTR options");
1288: PetscOptionsReal("-ksp_cg_radius", "Trust Region Radius", "KSPCGSetRadius", cg->radius, &cg->radius, NULL);
1290: PetscOptionsEList("-ksp_cg_dtype", "Norm used for direction", "", DType_Table, GLTR_DIRECTION_TYPES, DType_Table[cg->dtype], &cg->dtype, NULL);
1292: PetscOptionsReal("-ksp_cg_gltr_init_pert", "Initial perturbation", "", cg->init_pert, &cg->init_pert, NULL);
1293: PetscOptionsReal("-ksp_cg_gltr_eigen_tol", "Eigenvalue tolerance", "", cg->eigen_tol, &cg->eigen_tol, NULL);
1294: PetscOptionsReal("-ksp_cg_gltr_newton_tol", "Newton tolerance", "", cg->newton_tol, &cg->newton_tol, NULL);
1296: PetscOptionsInt("-ksp_cg_gltr_max_lanczos_its", "Maximum Lanczos Iters", "", cg->max_lanczos_its, &cg->max_lanczos_its, NULL);
1297: PetscOptionsInt("-ksp_cg_gltr_max_newton_its", "Maximum Newton Iters", "", cg->max_newton_its, &cg->max_newton_its, NULL);
1299: PetscOptionsTail();
1300: return 0;
1301: }
1303: /*MC
1304: KSPGLTR - Code to run conjugate gradient method subject to a constraint
1305: on the solution norm. This is used in Trust Region methods for
1306: nonlinear equations, SNESNEWTONTR
1308: Options Database Keys:
1309: . -ksp_cg_radius <r> - Trust Region Radius
1311: Notes:
1312: This is rarely used directly
1314: Use preconditioned conjugate gradient to compute
1315: an approximate minimizer of the quadratic function
1317: q(s) = g^T * s + .5 * s^T * H * s
1319: subject to the trust region constraint
1321: || s || <= delta,
1323: where
1325: delta is the trust region radius,
1326: g is the gradient vector,
1327: H is the Hessian approximation,
1328: M is the positive definite preconditioner matrix.
1330: KSPConvergedReason may be
1331: $ KSP_CONVERGED_CG_NEG_CURVE if convergence is reached along a negative curvature direction,
1332: $ KSP_CONVERGED_CG_CONSTRAINED if convergence is reached along a constrained step,
1333: $ other KSP converged/diverged reasons
1335: Notes:
1336: The preconditioner supplied should be symmetric and positive definite.
1338: Reference:
1339: Gould, N. and Lucidi, S. and Roma, M. and Toint, P., Solving the Trust-Region Subproblem using the Lanczos Method,
1340: SIAM Journal on Optimization, volume 9, number 2, 1999, 504-525
1342: Level: developer
1344: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCGSetRadius(), KSPCGGetNormD(), KSPCGGetObjFcn(), KSPGLTRGetMinEig(), KSPGLTRGetLambda()
1345: M*/
1347: PETSC_EXTERN PetscErrorCode KSPCreate_GLTR(KSP ksp)
1348: {
1349: KSPCG_GLTR *cg;
1351: PetscNewLog(ksp,&cg);
1352: cg->radius = 0.0;
1353: cg->dtype = GLTR_UNPRECONDITIONED_DIRECTION;
1355: cg->init_pert = 1.0e-8;
1356: cg->eigen_tol = 1.0e-10;
1357: cg->newton_tol = 1.0e-6;
1359: cg->alloced = 0;
1360: cg->init_alloc = 1024;
1362: cg->max_lanczos_its = 20;
1363: cg->max_newton_its = 10;
1365: ksp->data = (void*) cg;
1366: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,3);
1367: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
1368: KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,2);
1369: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
1371: /***************************************************************************/
1372: /* Sets the functions that are associated with this data structure */
1373: /* (in C++ this is the same as defining virtual functions). */
1374: /***************************************************************************/
1376: ksp->ops->setup = KSPCGSetUp_GLTR;
1377: ksp->ops->solve = KSPCGSolve_GLTR;
1378: ksp->ops->destroy = KSPCGDestroy_GLTR;
1379: ksp->ops->setfromoptions = KSPCGSetFromOptions_GLTR;
1380: ksp->ops->buildsolution = KSPBuildSolutionDefault;
1381: ksp->ops->buildresidual = KSPBuildResidualDefault;
1382: ksp->ops->view = NULL;
1384: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGSetRadius_C",KSPCGSetRadius_GLTR);
1385: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetNormD_C", KSPCGGetNormD_GLTR);
1386: PetscObjectComposeFunction((PetscObject)ksp,"KSPCGGetObjFcn_C",KSPCGGetObjFcn_GLTR);
1387: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetMinEig_C",KSPGLTRGetMinEig_GLTR);
1388: PetscObjectComposeFunction((PetscObject)ksp,"KSPGLTRGetLambda_C",KSPGLTRGetLambda_GLTR);
1389: return 0;
1390: }