1 //------------------------------------------------------------------------
2 // Copyright (C) 1993, 1994:
3 // J.C. Meza
4 // Sandia National Laboratories
5 // meza@california.sandia.gov
6 //------------------------------------------------------------------------
7 
8 #ifdef HAVE_CONFIG_H
9 #include "OPT++_config.h"
10 #endif
11 
12 #ifdef HAVE_STD
13 #include <cstdio>
14 #include <cstdlib>
15 #include <cmath>
16 #else
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <math.h>
20 #endif
21 
22 #if (defined(__sgi) || defined(__xlc__) || defined (__xlC__))
23 #define WANT_MATH
24 #else
25 #define WANT_STREAM
26 #define WANT_MATH
27 #endif
28 
29 #include "NLP.h"
30 #include "NLF.h"
31 
32 using std::cerr;
33 
34 using Teuchos::SerialDenseVector;
35 using Teuchos::SerialSymDenseMatrix;
36 
37 using namespace OPTPP;
38 
square(double x)39 static double square(double x) { return x*x; }
40 
init_rosen(int ndim,SerialDenseVector<int,double> & x)41 void init_rosen (int ndim, SerialDenseVector<int,double>& x)
42 {
43   if (ndim != 2)
44   {
45     exit (1);
46   }
47   x(0) = -1.2;
48   x(1) =  1.0;
49 }
rosen0(int n,const SerialDenseVector<int,double> & x,double & fx,int & result)50 void rosen0(int n, const SerialDenseVector<int,double>& x, double& fx, int& result)
51 { // Rosenbrock's function
52   double f1, f2, x1, x2;
53 
54   if (n != 2) return;
55 
56   x1 = x(0);
57   x2 = x(1);
58   f1 = (x2 - x1 * x1);
59   f2 = 1. - x1;
60 
61   fx  = 100.* f1*f1 + f2*f2;
62   result = NLPFunction;
63 }
64 
rosen(int mode,int n,const SerialDenseVector<int,double> & x,double & fx,SerialDenseVector<int,double> & g,int & result)65 void rosen(int mode, int n, const SerialDenseVector<int,double>& x, double& fx, SerialDenseVector<int,double>& g, int& result)
66 { // Rosenbrock's function
67   double f1, f2, x1, x2;
68 
69   if (n != 2) return;
70 
71   x1 = x(0);
72   x2 = x(1);
73   f1 = (x2 - x1 * x1);
74   f2 = 1. - x1;
75 
76   if (mode & NLPFunction) {
77     fx  = 100.* f1*f1 + f2*f2;
78     result = NLPFunction;
79   }
80   if (mode & NLPGradient) {
81     g(0) = -400.*f1*x1 - 2.*f2;
82     g(1) = 200.*f1;
83     result = NLPGradient;
84   }
85 }
86 
rosen2(int mode,int n,const SerialDenseVector<int,double> & x,double & fx,SerialDenseVector<int,double> & g,SerialSymDenseMatrix<int,double> & H,int & result)87 void rosen2(int mode, int n, const SerialDenseVector<int,double>& x, double& fx,
88 	SerialDenseVector<int,double>& g, SerialSymDenseMatrix<int,double>& H, int &result)
89 // Rosenbrock's function, n = 2 with first and second order derivatives
90 
91 {
92   //fcn_count++;
93 
94   double f1, f2, x1, x2;
95 
96 
97   // cout << "\nrosen2: mode = " << mode
98   //      << " count = " << fcn_count << "\n";
99   // for(int i=1; i<=n; i++)
100   //     cout << "x(" << i << ") = " << x(i) << "\n";
101 
102   if (n != 2) return;
103 
104   x1 = x(0);
105   x2 = x(1);
106   f1 = (x2 - x1 * x1);
107   f2 = 1. - x1;
108 
109   if (mode & NLPFunction) {
110     fx  = 100.* f1*f1 + f2*f2;
111     result = NLPFunction;
112   }
113   if (mode & NLPGradient) {
114     g(0) = -400.*f1*x1 - 2.*f2;
115     g(1) = 200.*f1;
116     result = NLPGradient;
117   }
118 
119   if (mode & NLPHessian) {
120 
121     f1 = (x2 - 3.0*x1*x1);
122     H(0,0) = -400.0*f1 + 2.0;
123     H(1,0) = -400.0*x1;
124     H(1,1) = 200.0;
125     result = NLPHessian;
126   }
127 }
128 
erosen1(int mode,int ndim,const SerialDenseVector<int,double> & x,double & fx,SerialDenseVector<int,double> & g,int & result)129 void erosen1(int mode, int ndim, const SerialDenseVector<int,double>& x, double& fx,
130 	SerialDenseVector<int,double>& g, int &result)
131 // Extended Rosenbrock's function, with analytic  derivatives
132 
133 {
134   static int i;
135   static double f1, f2, x1, x2;
136 
137   fx = 0.;
138 
139   for (i=0; i<ndim/2; ++i) {
140     x1 = x(2*i);
141     x2 = x(2*i+1);
142 
143     f1 = (x2 - x1 * x1);
144     f2 = 1. - x1;
145 
146     fx += 100.*square(f1) + square(f2);
147 
148     if (mode & NLPGradient) {
149       g(2*i) = -400.*f1*x1 - 2.*f2;
150       g(2*i+1) = 200.*f1;
151     }
152   }
153   result = NLPFunction;
154   if (mode & NLPGradient) result = NLPGradient;
155 
156 }
init_epowell(int ndim,SerialDenseVector<int,double> & x)157 void init_epowell (int ndim, SerialDenseVector<int,double>& x)
158 {
159   if (ndim % 4 != 0)
160   {
161     printf ("init_epowell: ndim must be multiple of 4, ndim = %d\n", ndim);
162     exit (1);
163   }
164 
165   for (int i=0; i<ndim/4; i++)
166   {
167     x(4*i    ) =  3;
168     x(4*i + 1) = -1;
169     x(4*i + 2) =  0;
170     x(4*i + 3) =  1;
171   }
172 }
epowell(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)173 void epowell(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
174 
175 {
176   // Extended Powell Singular Function, n variable, multiple of 4
177   // Standard initial guess  x0   = (psi1, psi2, psi3, psi4, ...)
178   //                         psi1 =  3
179   //                         psi2 = -1
180   //                         psi3 =  0
181   //                         psi4 =  1
182   // Solution                fmin = 0 at origin
183 
184   double f1, f2, f3, f4;
185 
186   if (ndim % 4 != 0) {
187     cerr << "epowell: ndim must be a multiple of 4, ndim = \n" << ndim;
188     return;
189   }
190 
191   fx = 0.0;
192   for (int i=0; i<ndim/4; i++) {
193     f1 = x(4*i) + 10*x(4*i+1);
194     f2 = sqrt(5.0) * ( x(4*i+2) - x(4*i+3));
195     f3 = square(x(4*i+1) - 2*x(4*i+2));
196     f4 = sqrt(10.0) * (x(4*i) - x(4*i+3));
197     fx += square(f1) + square(f2) + square(f3) + square(f4);
198   }
199   result = NLPFunction;
200 }
201 
init_trig(int ndim,SerialDenseVector<int,double> & x)202 void init_trig (int ndim, SerialDenseVector<int,double>& x)
203 {
204   for (int i=0; i<ndim; i++)
205   {
206     x(i) =  1/ double(ndim);
207   }
208 }
209 
trig(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)210 void trig(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
211 
212 {
213   // Trigonometric Function, ndim variable, m = ndim
214   // Standard initial guess  x0   = (1/ndim, ..., 1/ndim)
215   // Solution                fmin = 0
216 
217   double fi;
218 
219   fx = 0.0;
220 
221   for (int i=0; i<ndim; ++i) {
222     double sum = 0.0;
223     for (int j=0; j<ndim; ++j) {
224       sum += cos(x(j)) + double(i+1)*(1.0-cos(x(i)) - sin(x(i)));
225     }
226     fi = ndim - sum;
227     fx += fi*fi;
228   }
229 
230   result = NLPFunction;
231 }
init_erosen(int ndim,SerialDenseVector<int,double> & x)232 void init_erosen (int ndim, SerialDenseVector<int,double>& x)
233 {
234 
235   if (ndim % 2 != 0)
236   {
237     printf ("init_erosen: ndim must be even, ndim = %d\n", ndim);
238     exit (1);
239   }
240 
241   for (int i=0; i<ndim/2; ++i)
242   {
243     x(2*i)     = -1.2;
244     x(2*i + 1) =  1.0;
245   }
246 
247 }
erosen(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)248 void erosen(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
249 {
250 //  erosen
251 //  Extended rosenbrock function for any even value of n.
252 //                          2
253 //     f    (x) = 10(x   - x    )
254 //      2i-1          2i    2i-1
255 //
256 //     f  (x) = 1 - x
257 //      2i           2i-1
258 //                                           2         2
259 //     f(x) = sum from i = 1 to n/2 (f    (x)  + f  (x) ).
260 //                                    2i-1        2i
261 //
262 //     x  = (epsilon ) where epsilon    = -1.2, epsilon  = 1
263 //      0           i               2i-1               2i
264 //
265 //     f(x ) = 0 at x  = (1,...,1)
266 //        *          *
267 //
268 //     parameters
269 //
270 //     input
271 //        n             dimension of the problem to be solved must be even
272 //        x             point at which the function is to be evaluated
273 //     output
274 //        f             function value at x
275 
276 
277   static int i;
278   static double f1, f2, x1, x2;
279 
280   fx = 0.;
281 
282   for (i=0; i<ndim/2; ++i) {
283     x1 = x(2*i);
284     x2 = x(2*i+1);
285 
286     f1 = (x2 - x1 * x1);
287     f2 = 1. - x1;
288 
289     fx += 100.*square(f1) + square(f2);
290   }
291   result = NLPFunction;
292 }
293 
init_penalty1(int ndim,SerialDenseVector<int,double> & x)294 void init_penalty1 (int ndim, SerialDenseVector<int,double>& x)
295 {
296 
297   for (int i=0; i<ndim; i++) x(i) = double(i+1);
298 
299 }
penalty1(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)300 void penalty1(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
301 {
302 //  Penalty 1
303 //     parameters
304 //
305 //     input
306 //        ndim          dimension of the problem to be solved must be even
307 //        x             point at which the function is to be evaluated
308 //     output
309 //        f             function value at x
310 
311 
312   int i;
313   double fi;
314   double k1  = 1.e-5;
315   double sum = 0.0;
316 
317   fx = 0.;
318 
319   for (i=0; i<ndim; i++) sum += square(x(i));
320   fx = (sum - .25)*(sum - .25);
321 
322 
323   for (i=0; i<ndim; i++) {
324     fi = (x(i) - 1.0);
325     fx += k1*square(fi);
326   }
327   result = NLPFunction;
328 }
init_penalty2(int ndim,SerialDenseVector<int,double> & x)329 void init_penalty2 (int ndim, SerialDenseVector<int,double>& x)
330 {
331 
332   for (int i=0; i<ndim; i++) x(i) = 0.5;
333 
334 }
penalty2(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)335 void penalty2(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
336 {
337 //  Penalty 2
338 //     parameters
339 //
340 //     input
341 //        ndim          dimension of the problem to be solved must be even
342 //        x             point at which the function is to be evaluated
343 //     output
344 //        f             function value at x
345 
346 
347   int i;
348   double f1, fi, yi;
349 
350   double k1  = 1.e-5;
351   double sum = 0.0;
352 
353   fx = 0.0;
354 
355   // Take care of the two end cases first
356   //
357 
358   f1  = x(0) - 0.2;
359   fx  = square(f1);
360   for (i=0; i<ndim; i++) sum += double(ndim - i)*square(x(i));
361   sum = sum - 1.0;
362   fx += square(sum);
363 
364   //
365   // Now take care of all of the cases in between
366   //
367   for (i=1; i<ndim; i++) {
368     yi = exp(double(i+1)/10.0) + exp(double(i)/10.0);
369     fi = exp(x(i)/10.0) + exp(x(i-1)/10.0) - yi;
370     fx += k1*square(fi);
371   }
372   for (i=ndim; i<2*ndim-1; i++) {
373     fi = exp(x(i-ndim+1)/10.0) - exp(double(-1)/10.0);
374     fx += k1*square(fi);
375   }
376   result = NLPFunction;
377 }
init_vardim(int ndim,SerialDenseVector<int,double> & x)378 void init_vardim (int ndim, SerialDenseVector<int,double>& x)
379 {
380 
381   for (int i=0; i<ndim; i++) x(i) = 1.0 - (double(i+1)/double(ndim));
382 
383 }
vardim(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)384 void vardim(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
385 {
386 //  Variable dimensioned
387 //     parameters
388 //
389 //     input
390 //        ndim          dimension of the problem to be solved must be even
391 //        x             point at which the function is to be evaluated
392 //     output
393 //        f             function value at x
394 
395 
396   int i;
397   double fn_plus_1, fn_plus_2, fi;
398 
399   double sum = 0.0;
400 
401   fx = 0.;
402 
403   // Take care of the two end cases first
404   //
405   for (i=0; i<ndim; i++) sum += double(i+1)*(x(i) - 1.0);
406 
407   fn_plus_1  = sum;
408   fn_plus_2  = square(sum);
409   fx        += square(fn_plus_1) + square(fn_plus_2);
410 
411   //
412   // Now take care of all of the cases in between
413   //
414   for (i=0; i<ndim; i++) {
415     fi = x(i) - 1.0;
416     fx += square(fi);
417   }
418   result = NLPFunction;
419 }
init_broyden_tridiag(int ndim,SerialDenseVector<int,double> & x)420 void init_broyden_tridiag (int ndim, SerialDenseVector<int,double>& x)
421 {
422 
423   for (int i=0; i<ndim; i++) x(i)   = -1.0;
424 
425 }
broyden_tridiag(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)426 void broyden_tridiag(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
427 {
428 //  Broyden Tridiagonal function
429 //     parameters
430 //
431 //     input
432 //        ndim          dimension of the problem to be solved must be even
433 //        x             point at which the function is to be evaluated
434 //     output
435 //        f             function value at x
436 
437 
438   int i;
439   double f1, fi, fn;
440   double x0 = 0.0, xn_plus_1 = 0.0;
441 
442   fx = 0.;
443 
444   f1 = (3.0 - 2.0*x(0))*x(0) - x0 - 2.0*x(1) + 1.0;
445   fn = (3.0 - 2.0*x(ndim-1))*x(ndim-1) - x(ndim-2) - 2.0*xn_plus_1 + 1.0;
446 
447   fx += square(f1) + square(fn);
448   //
449   // Now take care of all of the cases in between
450   //
451   for (i=1; i<ndim-1; i++) {
452     fi  = (3.0 - 2.0*x(i))*x(i) - x(i-1) - 2.0*x(i+1) + 1.0;
453     fx += square(fi);
454   }
455   result = NLPFunction;
456 }
457 //-----------------------------------------------------------------------------
458 // Problems from Conn, Gould, Toint
459 //
460 //-----------------------------------------------------------------------------
init_gen_brown(int ndim,SerialDenseVector<int,double> & x)461 void init_gen_brown (int ndim, SerialDenseVector<int,double>& x)
462 {
463 
464   for (int i=0; i<ndim-1; i=i+2){
465     x(i)   = -1.0;
466     x(i+1) =  1.0;
467   }
468 
469 }
gen_brown(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)470 void gen_brown(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
471 {
472 //  Generalization of a function due to A. Brown
473 //  Conn, Gould, Toint, Test function 21, pp. 427
474 //     parameters
475 //
476 //     input
477 //        ndim          dimension of the problem to be solved must be even
478 //        x             point at which the function is to be evaluated
479 //     output
480 //        f             function value at x
481 
482 
483   int i;
484   double x0, x1, y1, y2;
485 
486   fx = 0.;
487 
488   //
489   for (i=0; i<ndim-1; i++) {
490     x0 = x(i);
491     x1 = x(i+1);
492     y1 = square(x1) + 1;
493     y2 = square(x0) + 1;
494     fx += pow(square(x0),y1) + pow(square(x1),y2);
495   }
496   result = NLPFunction;
497 }
init_chain_singular(int ndim,SerialDenseVector<int,double> & x)498 void init_chain_singular (int ndim, SerialDenseVector<int,double>& x)
499 {
500 
501   if (ndim % 4 != 0)
502   {
503     printf ("init_chain_singular: ndim must be multiple of 4, ndim = %d\n", ndim);
504     exit (1);
505   }
506 
507   for (int i=0; i<ndim/4; i++)
508   {
509     x(4*i    ) =  3;
510     x(4*i + 1) = -1;
511     x(4*i + 2) =  0;
512     x(4*i + 3) =  1;
513   }
514 
515 }
chain_singular(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)516 void chain_singular(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
517 {
518 //  Chained Singular Function
519 //  Conn, Gould, Toint, Test function 21, pp. 422
520 //     parameters
521 //
522 //     input
523 //        ndim          dimension of the problem to be solved must be even
524 //        x             point at which the function is to be evaluated
525 //     output
526 //        f             function value at x
527 
528 
529   int i;
530   double x0, x1, x2, x3;
531   double t1, t2, t3, t4;
532 
533   fx = 0.;
534 
535   if (ndim % 4 != 0)
536   {
537     printf ("chain_singular: ndim must be multiple of 4, ndim = %d\n", ndim);
538     exit (1);
539   }
540 
541   //
542   for (i=0; i<ndim-4; i=i+2) {
543     x0 = x(i);
544     x1 = x(i+1);
545     x2 = x(i+2);
546     x3 = x(i+3);
547     t1 = square(x0 + 10.0*x1);
548     t2 = 5.0*square(x2 - x3);
549     t3 = square(x1 - 2.0*x2)*square(x1 - 2.0*x2);
550     t4 = 10.0*square(x0 - 10.0*x3)*square(x0 - 10.0*x3);
551     fx += t1 + t2 + t3 + t4;
552   }
553   result = NLPFunction;
554 }
init_gen_wood(int ndim,SerialDenseVector<int,double> & x)555 void init_gen_wood (int ndim, SerialDenseVector<int,double>& x)
556 {
557 
558   if (ndim % 4 != 0)
559   {
560     printf ("init_gen_wood: ndim must be multiple of 4, ndim = %d\n", ndim);
561     exit (1);
562   }
563 
564   x(0) =  -3;
565   x(1) =  -1;
566   x(2) =  -3;
567   x(3) =  -1;
568 
569   for (int i=4; i<ndim; i=i+2)
570   {
571     x(i)   = -2;
572     x(i+1) =  0;
573   }
574   //  for (int i=1; i<=ndim; i++) x(i) = 1.1;
575 }
gen_wood(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)576 void gen_wood(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
577 {
578 //  Chained Singular Function
579 //  Conn, Gould, Toint, Test function 21, pp. 422
580 //     parameters
581 //
582 //     input
583 //        ndim          dimension of the problem to be solved must be even
584 //        x             point at which the function is to be evaluated
585 //     output
586 //        f             function value at x
587 
588 
589   int i;
590   double x0, x1, x2, x3;
591   double t1, t2, t3, t4, t5, t6;
592 
593   fx = 0.0;
594   if (ndim % 4 != 0)
595   {
596     printf ("gen_wood: ndim must be multiple of 4, ndim = %d\n", ndim);
597     exit (1);
598   }
599 
600   //
601   fx = 1.0;
602   for (i=0; i<ndim-4; i=i+4) {
603     x0 = x(i);
604     x1 = x(i+1);
605     x2 = x(i+2);
606     x3 = x(i+3);
607 
608     t1 = 100.0*square(x1 - square(x0));
609     t2 = square(1.0 - x0);
610     t3 = 90.0*square(x3 - square(x2));
611     t4 = square(1.0 - x2);
612     t5 = 10.0*square(x1 + x3 - 2.0);
613     t6 =  0.1*square(x1 - x3);
614     fx += t1 + t2 + t3 + t4 + t5 + t6;
615   }
616   result = NLPFunction;
617 }
init_chain_wood(int ndim,SerialDenseVector<int,double> & x)618 void init_chain_wood (int ndim, SerialDenseVector<int,double>& x)
619 {
620 
621   if (ndim % 4 != 0)
622   {
623     printf ("init_chain_wood: ndim must be multiple of 4, ndim = %d\n", ndim);
624     exit (1);
625   }
626 
627   x(0) =  -3;
628   x(1) =  -1;
629   x(2) =  -3;
630   x(3) =  -1;
631 
632   for (int i=4; i<ndim; i=i+2)
633   {
634     x(i)   = -2;
635     x(i+1) =  0;
636   }
637 
638 }
chain_wood(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)639 void chain_wood(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
640 {
641 //  Chained Singular Function
642 //  Conn, Gould, Toint, Test function 21, pp. 422
643 //     parameters
644 //
645 //     input
646 //        ndim          dimension of the problem to be solved must be even
647 //        x             point at which the function is to be evaluated
648 //     output
649 //        f             function value at x
650 
651 
652   int i;
653   double x0, x1, x2, x3;
654   double t1, t2, t3, t4, t5, t6;
655 
656   fx = 0.0;
657   if (ndim % 4 != 0)
658   {
659     printf ("chain_wood: ndim must be multiple of 4, ndim = %d\n", ndim);
660     exit (1);
661   }
662 
663   //
664   fx = 1.0;
665   for (i=0; i<ndim-4; i=i+2) {
666     x0 = x(i);
667     x1 = x(i+1);
668     x2 = x(i+2);
669     x3 = x(i+3);
670 
671     t1 = 100.0*square(x1 - square(x0));
672     t2 = square(1.0 - x0);
673     t3 = 90.0*square(x3 - square(x2));
674     t4 = square(1.0 - x2);
675     t5 = 10.0*square(x1 + x3 - 2.0);
676     t6 =  0.1*square(x1 - x3);
677     fx += t1 + t2 + t3 + t4 + t5 + t6;
678   }
679   result = NLPFunction;
680 }
init_broyden(int ndim,SerialDenseVector<int,double> & x)681 void init_broyden (int ndim, SerialDenseVector<int,double>& x)
682 {
683 
684   for (int i=0; i<ndim; i++) x(i)   = -1.0;
685 
686 }
broyden1a(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)687 void broyden1a(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
688 {
689 //  Broyden Tridiagonal function
690 //     parameters
691 //
692 //     input
693 //        ndim          dimension of the problem to be solved must be even
694 //        x             point at which the function is to be evaluated
695 //     output
696 //        f             function value at x
697 
698 
699   int i;
700   double f1, fi, fn;
701   double x0 = 0.0, xn_plus_1 = 0.0;
702   double p  = double(7)/double(3);
703 
704   fx = 1.0;
705 
706   // First the end cases
707   //
708   f1 = fabs((3.0 - 2.0*x(0))*x(0) - x0 - 2.0*x(1) + 1.0);
709   fn = fabs((3.0 - 2.0*x(ndim-1))*x(ndim-1) - x(ndim-2) - 2.0*xn_plus_1 + 1.0);
710 
711   fx += pow(f1,p) + pow(fn,p);
712   //
713   // Now take care of all of the cases in between
714   //
715   for (i=1; i<ndim-1; i++) {
716     fi  = fabs((3.0 - 2.0*x(i))*x(i) - x(i-1) - 2.0*x(i+1) + 1.0);
717     fx += pow(fi,p);
718   }
719   result = NLPFunction;
720 }
broyden1b(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)721 void broyden1b(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
722 {
723 //  Broyden Tridiagonal function
724 //     parameters
725 //
726 //     input
727 //        ndim          dimension of the problem to be solved must be even
728 //        x             point at which the function is to be evaluated
729 //     output
730 //        f             function value at x
731 
732 
733   int i;
734   double f1, fi, fn;
735   double x0  = 0.0, xn_plus_1 = 0.0;
736   double p   = 2.0;
737 
738   fx = 1.0;
739 
740   // First the end cases
741   //
742   f1 = fabs((3.0 - 2.0*x(0))*x(0) - x0 - 2.0*x(1) + 1.0);
743   fn = fabs((3.0 - 2.0*x(ndim-1))*x(ndim-1) - x(ndim-2) - 2.0*xn_plus_1 + 1.0);
744 
745   fx += pow(f1,p) + pow(fn,p);
746   //
747   // Now take care of all of the cases in between
748   //
749   for (i=1; i<ndim-1; i++) {
750     fi  = fabs((3.0 - 2.0*x(i))*x(i) - x(i-1) - 2.0*x(i+1) + 1.0);
751     fx += pow(fi,p);
752   }
753   result = NLPFunction;
754 }
broyden2a(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)755 void broyden2a(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
756 {
757 //  Generalization of Broyden Banded function, p = 7/3
758 //     parameters
759 //
760 //     input
761 //        ndim          dimension of the problem to be solved
762 //        x             point at which the function is to be evaluated
763 //     output
764 //        f             function value at x
765 
766 
767   int i, j;
768   double fi;
769   double sum;
770   double p  = double(7)/double(3);
771   int lower_bound, upper_bound;
772 
773   fx = 1.0;
774 
775   int ml = 5;
776   int mu = 1;
777 
778   for (i=0; i<ndim-1; i++) {
779 
780     lower_bound = max(1, i+1-ml);
781     upper_bound = min(ndim, i+1+mu);
782 
783     sum = 0.0;
784     for (j=0; j<ndim-1; j++) {
785       if ( (j+1 >= lower_bound) & (j+1 <= upper_bound))
786 	sum += x(j)*(1+x(j));
787     }
788 
789     fi  = fabs((2.0 + 5.0*square(x(i)))*x(i) + 1.0 - sum);
790     fx += pow(fi,p);
791   }
792   result = NLPFunction;
793 }
broyden2b(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)794 void broyden2b(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
795 {
796 //  Generalization of Broyden Banded function, p = 2
797 //     parameters
798 //
799 //     input
800 //        ndim          dimension of the problem to be solved must be even
801 //        x             point at which the function is to be evaluated
802 //     output
803 //        f             function value at x
804 
805 
806   int i, j;
807   double fi;
808   double sum;
809   double p  = 2.0;
810   int lower_bound, upper_bound;
811 
812   fx = 1.0;
813 
814   int ml = 5;
815   int mu = 1;
816 
817   for (i=0; i<ndim-1; i++) {
818 
819     lower_bound = max(1, i+1-ml);
820     upper_bound = min(ndim, i+1+mu);
821 
822     sum = 0.0;
823     for (j=0; j<ndim-1; j++) {
824       if ( (j+1 >= lower_bound) & (j+1 <= upper_bound))
825 	sum += x(j)*(1+x(j));
826     }
827 
828     fi  = fabs((2.0 + 5.0*square(x(i)))*x(i) + 1.0 - sum);
829     fx += pow(fi,p);
830   }
831   result = NLPFunction;
832 }
tointbroy(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)833 void tointbroy(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
834 {
835 //  Broyden Tridiagonal function
836 //     parameters
837 //
838 //     input
839 //        ndim          dimension of the problem to be solved must be even
840 //        x             point at which the function is to be evaluated
841 //     output
842 //        f             function value at x
843 
844 
845   int i;
846   double f1, fi, fn;
847   double x0 = 0.0, xn_plus_1 = 0.0;
848   double p  = double(7)/double(3);
849 
850   if (ndim % 2 != 0)
851   {
852     printf ("tointbroy: ndim must be multiple of 2, ndim = %d\n", ndim);
853     exit (1);
854   }
855 
856   fx = 1.0;
857 
858   // First the end cases
859   //
860   f1 = fabs((3.0 - 2.0*x(0))*x(0) - x0 - x(1) + 1.0);
861   fn = fabs((3.0 - 2.0*x(ndim-1))*x(ndim-1) - x(ndim-2) - xn_plus_1 + 1.0);
862 
863   fx += pow(f1,p) + pow(fn,p);
864   //
865   // Now take care of all of the cases in between
866   //
867   for (i=1; i<ndim-1; i++) {
868     fi  = fabs((3.0 - 2.0*x(i))*x(i) - x(i-1) - x(i+1) + 1.0);
869     fx += pow(fi,p);
870   }
871   //
872   // Add extra term due to Toint
873   //
874   int mid = ndim/2;
875   for (i=0; i<ndim/2-1; i++) {
876     fi  = fabs(x(i) + x(i+mid));
877     fx += pow(fi,p);
878   }
879   result = NLPFunction;
880 }
init_gen_cragg_levy(int ndim,SerialDenseVector<int,double> & x)881 void init_gen_cragg_levy (int ndim, SerialDenseVector<int,double>& x)
882 {
883 
884   if (ndim % 4 != 0)
885   {
886     printf ("init_gen_wood: ndim must be multiple of 4, ndim = %d\n", ndim);
887     exit (1);
888   }
889 
890   x(0) =  1;
891   for (int i=1; i<ndim; i++) x(i)   = 2;
892 }
gen_cragg_levy(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)893 void gen_cragg_levy(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
894 {
895 //  Generalized Cragg and Levy
896 //  Conn, Gould, Toint, Test function 17, pp. 422
897 //     parameters
898 //
899 //     input
900 //        ndim          dimension of the problem to be solved must be even
901 //        x             point at which the function is to be evaluated
902 //     output
903 //        f             function value at x
904 
905 
906   int i;
907   double x0, x1, x2, x3;
908   double t1, t2, t3, t4, t5;
909 
910   fx = 0.0;
911   if (ndim % 4 != 0)
912   {
913     printf ("gen_cragg_levy: ndim must be multiple of 4, ndim = %d\n", ndim);
914     exit (1);
915   }
916 
917   //
918   for (i=0; i<ndim-4; i=i+4) {
919     x0 = x(i);
920     x1 = x(i+1);
921     x2 = x(i+2);
922     x3 = x(i+3);
923 
924     t1 = pow((exp(x0) - x1),4.0);
925     t2 = 100.0*pow((x1 - x2),6.0);
926     t3 = pow((tan(x2 - x3)),4.0);
927     t4 = pow(x0,8.0);
928     t5 = square(x3 - 1.0);
929     fx += t1 + t2 + t3 + t4 + t5;
930   }
931   result = NLPFunction;
932 }
init_toint_trig(int ndim,SerialDenseVector<int,double> & x)933 void init_toint_trig (int ndim, SerialDenseVector<int,double>& x)
934 {
935 
936   for (int i=0; i<ndim; i++) x(i)   = 1.0;
937 
938 }void toint_trig(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
939 {
940 //  Toint Trig
941 //     parameters
942 //
943 //     input
944 //        ndim          dimension of the problem to be solved
945 //        x             point at which the function is to be evaluated
946 //     output
947 //        f             function value at x
948 
949 
950   int i, j;
951   SerialDenseVector<int,double> beta(ndim);
952   double alpha_ij, c_ij, t;
953   fx = 0.0;
954 
955   for (i=0; i<ndim; i++) {
956     beta(i) = 1.0 + double(i+1)/10.0;
957   }
958 
959   for (i=0; i<ndim; i++) {
960 
961     for (j=0; j<ndim; j++) {
962 
963       if ((abs(i-j))%4 == 0) {
964 	alpha_ij = 5.0*(1.0 + ((i+1)%5) + ((j+1)%5));
965 	c_ij     = double(i+j+2)/10.0;
966 	t        = beta(i)*x(i) + beta(j)*x(j) + c_ij;
967 	fx      += alpha_ij*sin(t);
968       }
969     }
970   }
971   result = NLPFunction;
972 }
init_chebyquad(int ndim,SerialDenseVector<int,double> & x)973 void init_chebyquad (int ndim, SerialDenseVector<int,double>& x)
974 {
975 
976   for (int i=0; i<ndim; i++) x(i)   = double(i+1)/double(ndim+1);
977 
978 }void chebyquad(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
979 {
980 //  Chebyquad
981 //     parameters
982 //
983 //     input
984 //        ndim          dimension of the problem to be solved
985 //        x             point at which the function is to be evaluated
986 //     output
987 //        f             function value at x
988 
989 
990   int i, i2, j;
991   int m2   = ndim/2;
992   int mdim = ndim;
993 
994   SerialDenseVector<int,double> fi(ndim);
995   double coeff = 1.0/double(ndim);
996   double t, t0, t1, tj;
997 
998   fx = 0.0;
999   fi = 0.0;
1000 
1001   for (i=0; i<m2; i++) {
1002     i2 = 2*(i+1);
1003     fi(i2-1) = 1.0 / double((i2*i2) - 1.0);
1004   }
1005 
1006   for (j=0; j<ndim; j++) {
1007     t0 = 1.0;
1008     t1 = 2.0 * x(j) - 1.0;
1009     tj = t1;
1010     fi(0) = fi(0) + coeff*t1;
1011 
1012     for (i=1; i<mdim; i++) {
1013       t = 2.0 * tj * t1 - t0;
1014       fi(i) = fi(i) + coeff*t;
1015       t0 = t1;
1016       t1 = t;
1017     }
1018   }
1019   for (i=0; i<ndim; i++) {
1020     fx += fi(i)*fi(i);
1021   }
1022 
1023   result = NLPFunction;
1024 }
init_nelder(int ndim,SerialDenseVector<int,double> & x)1025 void init_nelder (int ndim, SerialDenseVector<int,double>& x)
1026 {
1027 
1028   for (int i=0; i<ndim; i++) x(i)   = double(i+1)/double(ndim+1);
1029 
1030 }
nelder(int ndim,const SerialDenseVector<int,double> & x,double & fx,int & result)1031 void nelder(int ndim, const SerialDenseVector<int,double>& x, double& fx, int& result)
1032 {
1033 //  Nelder
1034 //     parameters
1035 //
1036 //     input
1037 //        ndim          dimension of the problem to be solved
1038 //        x             point at which the function is to be evaluated
1039 //     output
1040 //        f             function value at x
1041 
1042   double f1, f2, f3;
1043 
1044   if (ndim != 2)
1045   {
1046     exit (1);
1047   }
1048 
1049   f1 = x(3) - x(0)*x(1)*x(2);
1050   f2 = x(2) - x(0)*x(1);
1051   f3 = x(1) - x(0);
1052 
1053   fx = f1*f1 + f2*f2 + f3*f3 + x(0)*x(0);
1054 
1055   result = NLPFunction;
1056 }
1057