1 //Possible Print Statement Issues
2
3 //JWG
4
5 //------------------------------------------------------------------------
6 // Copyright (C) 1993:
7 // J.C. Meza
8 // Sandia National Laboratories
9 // meza@california.sandia.gov
10 //------------------------------------------------------------------------
11
12 #include "Opt.h"
13 #include "ioformat.h"
14 #include "Teuchos_SerialDenseMatrix.hpp"
15 #include "Teuchos_SerialDenseVector.hpp"
16 #include "Teuchos_SerialSymDenseMatrix.hpp"
17
18 using Teuchos::SerialDenseMatrix;
19 using Teuchos::SerialSymDenseMatrix;
20
21 namespace OPTPP {
22
trustregion(NLP1 * nlp,std::ostream * fout,SerialSymDenseMatrix<int,double> & H,SerialDenseVector<int,double> & search_dir,SerialDenseVector<int,double> & sx,double & TR_size,double & step_length,double stpmax,double stpmin)23 int trustregion(NLP1* nlp, std::ostream *fout,
24 SerialSymDenseMatrix<int,double>& H, SerialDenseVector<int,double>& search_dir,
25 SerialDenseVector<int,double>& sx, double& TR_size, double& step_length,
26 double stpmax, double stpmin)
27 {
28 /****************************************************************************
29 * subroutine trustregion
30 *
31 * Purpose
32 * find a step which satisfies the Goldstein-Armijo line search conditions
33 *
34 * Compute the dogleg step
35 * Compute the predicted reduction, pred, of the quadratic model
36 * Compute the actual reduction, ared
37 * IF ared/pred > eta
38 * THEN x_vec = x_vec + d_vec
39 * TR_size >= TR_size
40 * Compute the gradient g_vec at the new point
41 * ELSE TR_size < ||d_vec||
42 *
43 * Parameters
44 * nlp --> pointer to nonlinear problem object
45 *
46 * search_dir --> Vector of length n which specifies the
47 * newton direction on input. On output it will
48 * contain the step
49 *
50 * step_length <-- is a nonnegative variable.
51 * On output step_length contains the step size taken
52 *
53 * ftol --> default Value = 1.e-4
54 * ftol should be smaller than 5.e-1
55 * suggested value = 1.e-4 for newton methods
56 * = 1.e-1 for more exact line searches
57 * xtol --> default Value = 2.2e-16
58 * gtol --> default Value = 0.9
59 * gtol should be greater than 1.e-4
60 *
61 * termination occurs when the sufficient decrease
62 * condition and the directional derivative condition are
63 * satisfied.
64 *
65 * stpmin and TR_size are nonnegative input variables which
66 * specify lower and upper bounds for the step.
67 * stpmin Default Value = 1.e-9
68 *
69 * Initial version Juan Meza November 1994
70 *
71 *
72 *****************************************************************************/
73
74 // Local variables
75
76 int n = nlp->getDim();
77 bool debug = nlp->getDebug();
78 bool modeOverride = nlp->getModeOverride();
79
80 SerialDenseVector<int,double> tgrad(n), newton_dir(n), tvec(n), xc(n), xtrial(n);
81 double fvalue, fplus, dnorm;
82 double eta1 = .001;
83 double eta2 = .1;
84 double eta3 = .75;
85 double rho_k;
86 int iter = 0;
87 int iter_max = 100;
88 double dd1, dd2;
89 double ared, pred;
90 int dog_step;
91 double TR_MAX = stpmax;
92 static const char *steps[] = {"C", "D", "N", "B"};
93 static bool accept;
94
95 //
96 // Initialize variables
97 //
98 fvalue = nlp->getF();
99 xc = nlp->getXc();
100 step_length = 1.0;
101 tgrad = nlp->getGrad();
102 newton_dir = search_dir;
103
104 if (debug) {
105 *fout << "\n***************************************";
106 *fout << "***************************************\n";
107 *fout << "\nComputeStep using trustregion\n";
108 *fout << "\tStep ||step|| ared pred TR_size \n";
109 }
110
111 while (iter < iter_max) {
112 iter++;
113 //
114 // Compute the dogleg step
115 //
116 search_dir = newton_dir;
117
118 dog_step = dogleg(nlp, fout, H, tgrad, search_dir, sx, dnorm, TR_size, stpmax);
119 step_length = dnorm;
120 //
121 // Compute pred = -g'd - 1/2 d'Hd
122 //
123 dd1 = tgrad.dot(search_dir);
124 tvec.multiply(Teuchos::LEFT_SIDE, 1.0, H, search_dir, 0.0);
125 dd2 = search_dir.dot(tvec);
126 pred = -dd1 - dd2/2.0;
127
128 //
129 // Compute objective function at trial point
130 //
131
132 xtrial = xc;
133 xtrial += search_dir;
134 if (modeOverride) {
135 nlp->setX(xtrial);
136 nlp->eval();
137 fplus = nlp->getF();
138 }
139 else
140 fplus = nlp->evalF(xtrial);
141 ared = fvalue - fplus;
142
143
144
145 //
146 // Should we take this step ?
147 //
148
149 rho_k = ared/pred;
150
151 accept = false;
152 if (rho_k >= eta1) accept = true;
153
154 //
155 // Update the trust region
156 //
157
158 if (accept) {
159
160 //
161 // Do the standard updating
162 //
163
164 if (rho_k <= eta2) {
165
166 // Model is just sort of bad
167 // New trust region will be TR_size/2
168
169 if (debug) {
170 *fout << "trustregion: rho_k = " << e(rho_k,14,4)
171 << " eta2 = " << e(eta2,14,4) << "\n";
172 }
173 TR_size = step_length / 2.0;
174
175 if (TR_size < stpmin) {
176 *fout << "***** Trust region too small to continue.\n";
177 nlp->setX(xc);
178 nlp->setF(fvalue);
179 nlp->setGrad(tgrad);
180 return(-1);
181 }
182 }
183
184 else if ((eta3 <= rho_k) && (rho_k <= (2.0 - eta3))) {
185
186 // Model is PRETTY good
187 // Double trust region
188
189 if (debug) {
190 *fout << "trustregion: rho_k = " << e(rho_k,14,4)
191 << " eta3 = " << e(eta3,14,4) << "\n";
192 }
193 TR_size = min(2.0*TR_size, TR_MAX);
194 }
195
196 else {
197
198 //
199 // All other cases
200 //
201
202 TR_size = max(2.0*step_length,TR_size);
203 TR_size = min(TR_size, TR_MAX);
204 if (debug) {
205 *fout << "trustregion: rho_k = " << e(rho_k,14,4) << "\n";
206 }
207 }
208 }
209
210 else {
211
212 // Model is REALLY bad
213 //
214
215 TR_size = step_length/10.0;
216
217 if (debug) {
218 *fout << "trustregion: rho_k = " << e(rho_k,14,4) << "n"
219 << " eta1 = " << e(eta1,14,4) << "\n";
220 }
221 if (TR_size < stpmin) {
222 *fout << "***** Trust region too small to continue.\n";
223 nlp->setX(xc);
224 nlp->setF(fvalue);
225 nlp->setGrad(tgrad);
226 return(-1);
227 }
228 }
229
230 //
231 // Accept/Reject Step
232 //
233
234 if (accept) {
235 //
236 // Update x, f, and grad
237 //
238
239 if (!modeOverride) {
240 nlp->setX(xtrial);
241 nlp->setF(fplus);
242 nlp->evalG();
243 }
244
245 if (debug) {
246 *fout << "\t Step ||step|| ared pred TR_size \n";
247 *fout << "Accept " << steps[dog_step] << e(step_length,14,4)
248 << e(ared,14,4) << e(pred,14,4) << e(TR_size,14,4) << "\n";
249 }
250 return(dog_step);
251 }
252
253 else {
254 //
255 // Reject step
256 //
257
258 if (debug) {
259 *fout << "\t Step ||step|| ared pred TR_size \n";
260 *fout << "Reject " << steps[dog_step] << e(step_length,14,4)
261 << e(ared,14,4) << e(pred,14,4) << e(TR_size,14,4) << "\n";
262 }
263 }
264 }
265 nlp->setX(xc);
266 nlp->setF(fvalue);
267 nlp->setGrad(tgrad);
268 return(-1); // Too many iterations
269 }
270
271 } // namespace OPTPP
272