1 //------------------------------------------------------------------------
2 // Copyright (C) 1993:
3 // J.C. Meza
4 // Sandia National Laboratories
5 // meza@california.sandia.gov
6 //------------------------------------------------------------------------
7
8 #include "NLF.h"
9 #include "TOLS.h"
10 #include "cblas.h"
11
12 using namespace std;
13 using Teuchos::SerialDenseVector;
14 using Teuchos::SerialDenseMatrix;
15 using Teuchos::SerialSymDenseMatrix;
16
17 namespace OPTPP {
18
reset()19 void NLF1::reset() // Reset parameter values
20 {
21 init_flag = false;
22 nfevals = ngevals = 0;
23 #ifdef OPTPP_HAVE_MPI
24 SpecFlag = Spec1;
25 #else
26 SpecFlag = NoSpec;
27 #endif
28 application.reset();
29 }
30
initFcn()31 void NLF1::initFcn() // Initialize Function
32 {
33 if (init_flag == false) {
34 init_fcn(dim, mem_xc);
35 init_flag = true;
36 }
37 else {
38 cerr << "NLF1:initFcn: Warning - initialization called twice\n";
39 init_fcn(dim, mem_xc);
40 }
41 }
42
evalF()43 double NLF1::evalF() // Evaluate Function
44 {
45 int result = 0;
46 SerialDenseVector<int,double> gtmp(dim);
47
48 double time0 = get_wall_clock_time();
49 // *** CHANGE *** //
50 if (!application.getF(mem_xc,fvalue)) {
51 fcn_v(NLPFunction, dim, mem_xc, fvalue, gtmp, result, vptr);
52 application.update(result,dim,mem_xc,fvalue,gtmp);
53 nfevals++;
54 }
55 // *** CHANGE *** //
56 function_time = get_wall_clock_time() - time0;
57
58 if (debug_)
59 cout << "NLF1::evalF()\n"
60 << "nfevals = " << nfevals << "\n"
61 << "fvalue = " << fvalue << "\n"
62 << "function time = " << function_time << "\n";
63 return fvalue;
64 }
65
evalF(const SerialDenseVector<int,double> & x)66 double NLF1::evalF(const SerialDenseVector<int,double>& x) // Evaluate Function at x
67 {
68 int result = 0;
69 double fx;
70 SerialDenseVector<int,double> gtmp(dim);
71
72 double time0 = get_wall_clock_time();
73 // *** CHANGE *** //
74 if (!application.getF(x,fx)) {
75 fcn_v(NLPFunction, dim, x, fx, gtmp, result, vptr);
76 application.update(result,dim,x,fx,gtmp);
77 nfevals++;
78 }
79 // *** CHANGE *** //
80 function_time = get_wall_clock_time() - time0;
81
82 if (debug_)
83 cout << "NLF1::evalF(x)\n"
84 << "nfevals = " << nfevals << "\n"
85 << "fvalue = " << fx << "\n"
86 << "function time = " << function_time << "\n";
87 return fx;
88 }
89
evalG()90 SerialDenseVector<int,double> NLF1::evalG() // Evaluate the gradient
91 {
92 int result = 0;
93 double fx;
94
95 // *** CHANGE *** //
96 if (!application.getGrad(mem_xc,mem_grad)) {
97 fcn_v(NLPGradient, dim, mem_xc, fx, mem_grad, result, vptr);
98 application.update(result,dim,mem_xc,fx,mem_grad);
99 ngevals++;
100 }
101 // *** CHANGE *** //
102 return mem_grad;
103 }
104
evalG(const SerialDenseVector<int,double> & x)105 SerialDenseVector<int,double> NLF1::evalG(const SerialDenseVector<int,double>& x) // Evaluate the gradient at x
106 {
107 int result = 0 ;
108 double fx;
109 SerialDenseVector<int,double> gx(dim);
110
111 // *** CHANGE *** //
112 if (!application.getGrad(x,gx)) {
113 fcn_v(NLPGradient, dim, x, fx, gx, result, vptr);
114 application.update(result,dim,x,fx,gx);
115 ngevals++;
116 }
117 // *** CHANGE *** //
118 return gx;
119 }
120
evalH()121 SerialSymDenseMatrix<int,double> NLF1::evalH() // Evaluate the Hessian
122 {
123 SerialDenseVector<int,double> sx(dim);
124 SerialSymDenseMatrix<int,double> Hessian(dim);
125
126 sx = 1.0;
127 Hessian = FDHessian(sx);
128 return Hessian;
129 }
130
evalH(SerialDenseVector<int,double> & x)131 SerialSymDenseMatrix<int,double> NLF1::evalH(SerialDenseVector<int,double>& x) // Evaluate the Hessian at x
132 {
133 SerialSymDenseMatrix<int,double> Hessian(dim);
134
135 Hessian = FDHessian(x);
136 return Hessian;
137 }
138
eval()139 void NLF1::eval() // Evaluate Function and Gradient
140 {
141 int mode = NLPFunction | NLPGradient, result = 0;
142
143 double time0 = get_wall_clock_time();
144 // *** CHANGE *** //
145 if (!application.getF(mem_xc,fvalue) || !application.getGrad(mem_xc,mem_grad)) {
146 fcn_v(mode, dim, mem_xc, fvalue, mem_grad, result, vptr);
147 application.update(result,dim,mem_xc,fvalue,mem_grad);
148 nfevals++; ngevals++;
149 }
150 // *** CHANGE *** //
151
152 function_time = get_wall_clock_time() - time0;
153
154 if (debug_)
155 cout << "NLF1::eval()\n"
156 << "mode = " << mode << "\n"
157 << "nfevals = " << nfevals << "\n"
158 << "fvalue = " << fvalue << "\n"
159 << "function time = " << function_time << "\n";
160 }
161
evalLagrangian(const SerialDenseVector<int,double> & xc,SerialDenseVector<int,double> & multiplier,const SerialDenseVector<int,double> & type)162 double NLF1::evalLagrangian(const SerialDenseVector<int,double>& xc ,
163 SerialDenseVector<int,double>& multiplier,
164 const SerialDenseVector<int,double>& type)
165 {
166 double result = evalF(xc);
167 if( hasConstraints()){
168 SerialDenseVector<int,double> resid(constraint_->evalResidual(xc));
169 // SerialDenseVector<int,double> resid(constraint_->getNumOfCons());
170 // resid = constraint_->evalResidual(xc);
171 result -= resid.dot(multiplier);
172 }
173 return result;
174 }
175
evalLagrangianGradient(const SerialDenseVector<int,double> & xc,const SerialDenseVector<int,double> & multiplier,const SerialDenseVector<int,double> & type)176 SerialDenseVector<int,double> NLF1::evalLagrangianGradient(const SerialDenseVector<int,double>& xc,
177 const SerialDenseVector<int,double>& multiplier,
178 const SerialDenseVector<int,double>& type)
179 {
180 SerialDenseVector<int,double> grad = evalG(xc);
181 if(hasConstraints()){
182 SerialDenseVector<int,double> tmult(multiplier);
183 // tmult = multiplier;
184 int alpha = grad.length();
185 SerialDenseVector<int,double> tmult2(alpha);
186 //tmult *= -1;
187 // grad += constraint_->evalGradient(xc)*tmult;
188 tmult2.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS,-1.0,constraint_->evalGradient(xc),tmult,0.0);
189 grad+=tmult2;
190 }
191 return grad;
192
193 }
194
195
evalCF(const SerialDenseVector<int,double> & x)196 SerialDenseVector<int,double> NLF1::evalCF(const SerialDenseVector<int,double>& x) // Evaluate Function at x
197 {
198 int result = 0;
199 SerialDenseVector<int,double> cfx(ncnln);
200 SerialDenseMatrix<int,double> gtmp(dim,ncnln);
201
202 double time0 = get_wall_clock_time();
203 // *** CHANGE *** //
204 if (!application.getCF(x,cfx)) {
205 confcn(NLPFunction, dim, x, cfx, gtmp, result);
206 application.constraint_update(result,dim,ncnln,x,cfx,gtmp);
207 }
208 // *** CHANGE *** //
209 function_time = get_wall_clock_time() - time0;
210
211 if (debug_)
212 cout << "NLF1::evalCF(x)\n"
213 << "nfevals = " << nfevals << "\n"
214 // << "fvalue = " << cfx << "\n"
215 << "function time = " << function_time << "\n";
216 return cfx;
217 }
218
evalCG(const SerialDenseVector<int,double> & x)219 SerialDenseMatrix<int,double> NLF1::evalCG(const SerialDenseVector<int,double>& x) // Evaluate the gradient at x
220 {
221 int result = 0 ;
222 SerialDenseVector<int,double> cfx(ncnln);
223 SerialDenseMatrix<int,double> cgx(dim,ncnln);
224
225 // *** CHANGE *** //
226 if (!application.getCGrad(x,cgx)) {
227 confcn(NLPGradient, dim, x, cfx, cgx, result);
228 application.constraint_update(result,dim,ncnln,x,cfx,cgx);
229 }
230 // *** CHANGE *** //
231 return cgx;
232 }
233
evalCH(SerialDenseVector<int,double> & x)234 SerialSymDenseMatrix<int,double> NLF1::evalCH(SerialDenseVector<int,double>& x) // Evaluate the Hessian at x
235 {
236 SerialSymDenseMatrix<int,double> Hessian(dim);
237
238 // PJW This is a dummy routine. NIPS is the only algorithm which supports
239 // nonlinear constraints and currently this routine is not being accessed.
240 Hessian = 0.0;
241 return Hessian;
242 }
243
244
evalCH(SerialDenseVector<int,double> & x,int darg)245 OptppArray<SerialSymDenseMatrix<int,double> > NLF1::evalCH(SerialDenseVector<int,double>& x, int darg) // Evaluate the Hessian at x
246 {
247 OptppArray<SerialSymDenseMatrix<int,double> > Hessian(ncnln);
248 Hessian = CONFDHessian(x);
249 return Hessian;
250 }
251
evalC(const SerialDenseVector<int,double> & x)252 void NLF1::evalC(const SerialDenseVector<int,double>& x)
253 {
254 int mode = NLPFunction | NLPGradient, result = 0;
255 SerialDenseVector<int,double> cfx(ncnln);
256 SerialDenseMatrix<int,double> cgx(dim, ncnln);
257
258 double time0 = get_wall_clock_time();
259
260 if (!application.getCF(x, cfx) || !application.getCGrad(x, cgx)) {
261 confcn(mode, dim, x, cfx, cgx, result);
262 application.constraint_update(result, dim, ncnln, x, cfx, cgx);
263 }
264
265 function_time = get_wall_clock_time() - time0;
266 }
267
268 } // namespace OPTPP
269