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