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