Actual source code: gpcglinesearch.c
1: #include <petsc/private/taolinesearchimpl.h>
2: #include <../src/tao/linesearch/impls/gpcglinesearch/gpcglinesearch.h>
4: /* ---------------------------------------------------------- */
6: static PetscErrorCode TaoLineSearchDestroy_GPCG(TaoLineSearch ls)
7: {
8: TaoLineSearch_GPCG *ctx = (TaoLineSearch_GPCG *)ls->data;
10: VecDestroy(&ctx->W1);
11: VecDestroy(&ctx->W2);
12: VecDestroy(&ctx->Gold);
13: VecDestroy(&ctx->x);
14: PetscFree(ls->data);
15: return 0;
16: }
18: /*------------------------------------------------------------*/
19: static PetscErrorCode TaoLineSearchView_GPCG(TaoLineSearch ls, PetscViewer viewer)
20: {
21: PetscBool isascii;
23: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
24: if (isascii) PetscViewerASCIIPrintf(viewer, " GPCG Line search");
25: return 0;
26: }
28: /*------------------------------------------------------------*/
29: static PetscErrorCode TaoLineSearchApply_GPCG(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s)
30: {
31: TaoLineSearch_GPCG *neP = (TaoLineSearch_GPCG *)ls->data;
32: PetscInt i;
33: PetscBool g_computed = PETSC_FALSE; /* to prevent extra gradient computation */
34: PetscReal d1, finit, actred, prered, rho, gdx;
36: /* ls->stepmin - lower bound for step */
37: /* ls->stepmax - upper bound for step */
38: /* ls->rtol - relative tolerance for an acceptable step */
39: /* ls->ftol - tolerance for sufficient decrease condition */
40: /* ls->gtol - tolerance for curvature condition */
41: /* ls->nfeval - number of function evaluations */
42: /* ls->nfeval - number of function/gradient evaluations */
43: /* ls->max_funcs - maximum number of function evaluations */
45: TaoLineSearchMonitor(ls, 0, *f, 0.0);
47: ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
48: ls->step = ls->initstep;
49: if (!neP->W2) {
50: VecDuplicate(x, &neP->W2);
51: VecDuplicate(x, &neP->W1);
52: VecDuplicate(x, &neP->Gold);
53: neP->x = x;
54: PetscObjectReference((PetscObject)neP->x);
55: } else if (x != neP->x) {
56: VecDestroy(&neP->x);
57: VecDestroy(&neP->W1);
58: VecDestroy(&neP->W2);
59: VecDestroy(&neP->Gold);
60: VecDuplicate(x, &neP->W1);
61: VecDuplicate(x, &neP->W2);
62: VecDuplicate(x, &neP->Gold);
63: PetscObjectDereference((PetscObject)neP->x);
64: neP->x = x;
65: PetscObjectReference((PetscObject)neP->x);
66: }
68: VecDot(g, s, &gdx);
69: if (gdx > 0) {
70: PetscInfo(ls, "Line search error: search direction is not descent direction. dot(g,s) = %g\n", (double)gdx);
71: ls->reason = TAOLINESEARCH_FAILED_ASCENT;
72: return 0;
73: }
74: VecCopy(x, neP->W2);
75: VecCopy(g, neP->Gold);
76: if (ls->bounded) {
77: /* Compute the smallest steplength that will make one nonbinding variable equal the bound */
78: VecStepBoundInfo(x, s, ls->lower, ls->upper, &rho, &actred, &d1);
79: ls->step = PetscMin(ls->step, d1);
80: }
81: rho = 0;
82: actred = 0;
84: if (ls->step < 0) {
85: PetscInfo(ls, "Line search error: initial step parameter %g< 0\n", (double)ls->step);
86: ls->reason = TAOLINESEARCH_HALTED_OTHER;
87: return 0;
88: }
90: /* Initialization */
91: finit = *f;
92: for (i = 0; i < ls->max_funcs; i++) {
93: /* Force the step to be within the bounds */
94: ls->step = PetscMax(ls->step, ls->stepmin);
95: ls->step = PetscMin(ls->step, ls->stepmax);
97: VecWAXPY(neP->W2, ls->step, s, x);
98: if (ls->bounded) {
99: /* Make sure new vector is numerically within bounds */
100: VecMedian(neP->W2, ls->lower, ls->upper, neP->W2);
101: }
103: /* Gradient is not needed here. Unless there is a separate
104: gradient routine, compute it here anyway to prevent recomputing at
105: the end of the line search */
106: if (ls->hasobjective) {
107: TaoLineSearchComputeObjective(ls, neP->W2, f);
108: g_computed = PETSC_FALSE;
109: } else if (ls->usegts) {
110: TaoLineSearchComputeObjectiveAndGTS(ls, neP->W2, f, &gdx);
111: g_computed = PETSC_FALSE;
112: } else {
113: TaoLineSearchComputeObjectiveAndGradient(ls, neP->W2, f, g);
114: g_computed = PETSC_TRUE;
115: }
117: TaoLineSearchMonitor(ls, i + 1, *f, ls->step);
119: if (0 == i) ls->f_fullstep = *f;
121: actred = *f - finit;
122: VecWAXPY(neP->W1, -1.0, x, neP->W2); /* W1 = W2 - X */
123: VecDot(neP->W1, neP->Gold, &prered);
125: if (PetscAbsReal(prered) < 1.0e-100) prered = 1.0e-12;
126: rho = actred / prered;
128: /*
129: If sufficient progress has been obtained, accept the
130: point. Otherwise, backtrack.
131: */
133: if (actred > 0) {
134: PetscInfo(ls, "Step resulted in ascent, rejecting.\n");
135: ls->step = (ls->step) / 2;
136: } else if (rho > ls->ftol) {
137: break;
138: } else {
139: ls->step = (ls->step) / 2;
140: }
142: /* Convergence testing */
144: if (ls->step <= ls->stepmin || ls->step >= ls->stepmax) {
145: ls->reason = TAOLINESEARCH_HALTED_OTHER;
146: PetscInfo(ls, "Rounding errors may prevent further progress. May not be a step satisfying\n");
147: PetscInfo(ls, "sufficient decrease and curvature conditions. Tolerances may be too small.\n");
148: break;
149: }
150: if (ls->step == ls->stepmax) {
151: PetscInfo(ls, "Step is at the upper bound, stepmax (%g)\n", (double)ls->stepmax);
152: ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND;
153: break;
154: }
155: if (ls->step == ls->stepmin) {
156: PetscInfo(ls, "Step is at the lower bound, stepmin (%g)\n", (double)ls->stepmin);
157: ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
158: break;
159: }
160: if ((ls->nfeval + ls->nfgeval) >= ls->max_funcs) {
161: PetscInfo(ls, "Number of line search function evals (%" PetscInt_FMT ") > maximum (%" PetscInt_FMT ")\n", ls->nfeval + ls->nfgeval, ls->max_funcs);
162: ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
163: break;
164: }
165: if ((neP->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) {
166: PetscInfo(ls, "Relative width of interval of uncertainty is at most rtol (%g)\n", (double)ls->rtol);
167: ls->reason = TAOLINESEARCH_HALTED_RTOL;
168: break;
169: }
170: }
171: PetscInfo(ls, "%" PetscInt_FMT " function evals in line search, step = %g\n", ls->nfeval + ls->nfgeval, (double)ls->step);
172: /* set new solution vector and compute gradient if necessary */
173: VecCopy(neP->W2, x);
174: if (ls->reason == TAOLINESEARCH_CONTINUE_ITERATING) ls->reason = TAOLINESEARCH_SUCCESS;
175: if (!g_computed) TaoLineSearchComputeGradient(ls, x, g);
176: return 0;
177: }
179: /* ---------------------------------------------------------- */
181: /*MC
182: TAOLINESEARCHGPCG - Special line-search method for the Gradient-Projected Conjugate Gradient (TAOGPCG) algorithm.
183: Should not be used with any other algorithm.
185: Level: developer
187: .keywords: Tao, linesearch
188: M*/
189: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_GPCG(TaoLineSearch ls)
190: {
191: TaoLineSearch_GPCG *neP;
193: ls->ftol = 0.05;
194: ls->rtol = 0.0;
195: ls->gtol = 0.0;
196: ls->stepmin = 1.0e-20;
197: ls->stepmax = 1.0e+20;
198: ls->nfeval = 0;
199: ls->max_funcs = 30;
200: ls->step = 1.0;
202: PetscNew(&neP);
203: neP->bracket = 0;
204: neP->infoc = 1;
205: ls->data = (void *)neP;
207: ls->ops->setup = NULL;
208: ls->ops->reset = NULL;
209: ls->ops->apply = TaoLineSearchApply_GPCG;
210: ls->ops->view = TaoLineSearchView_GPCG;
211: ls->ops->destroy = TaoLineSearchDestroy_GPCG;
212: ls->ops->setfromoptions = NULL;
213: ls->ops->monitor = NULL;
214: return 0;
215: }