source: trunk/ADOL-C/examples/additional_examples/ipopt/LuksanVlcek1_sparse/ADOL-C_sparseNLP.cpp @ 344

Last change on this file since 344 was 344, checked in by awalther, 7 years ago

set derivative values array pointer to NULL in LuksanVlek1_sparse

File size: 13.3 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     ADOL-C_sparseNLP.cpp
4 Revision: $$
5 Contents: class myADOLC_sparseNPL for interfacing with Ipopt
6 
7 Copyright (c) Andrea Walther
8   
9 This file is part of ADOL-C. This software is provided as open source.
10 Any use, reproduction, or distribution of the software constitutes
11 recipient's acceptance of the terms of the accompanying license file.
12 
13 This code is based on the file  MyNLP.cpp contained in the Ipopt package
14 with the authors:  Carl Laird, Andreas Waechter   
15----------------------------------------------------------------------------*/
16
17/** C++ Example NLP for interfacing a problem with IPOPT and ADOL-C.
18 *  MyADOLC_sparseNLP implements a C++ example showing how to interface
19 *  with IPOPT and ADOL-C through the TNLP interface. This class
20 *  implements the Example 5.1 from "Sparse and Parially Separable
21 *  Test Problems for Unconstrained and Equality Constrained
22 *  Optimization" by L. Luksan and J. Vlcek taking sparsity
23 *  into account.
24 *
25 *  exploitation of sparsity !!
26 *
27 */
28
29#include <cassert>
30
31#include "ADOL-C_sparseNLP.hpp"
32
33using namespace Ipopt;
34
35/* Constructor. */
36MyADOLC_sparseNLP::MyADOLC_sparseNLP()
37{}
38
39MyADOLC_sparseNLP::~MyADOLC_sparseNLP()
40{}
41
42bool MyADOLC_sparseNLP::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g,
43                         Index& nnz_h_lag, IndexStyleEnum& index_style)
44{
45  n = 4000;
46
47  m = n-2;
48
49  generate_tapes(n, m, nnz_jac_g, nnz_h_lag);
50
51  // use the C style indexing (0-based)
52  index_style = C_STYLE;
53
54  return true;
55}
56
57bool MyADOLC_sparseNLP::get_bounds_info(Index n, Number* x_l, Number* x_u,
58                            Index m, Number* g_l, Number* g_u)
59{
60  // none of the variables have bounds
61  for (Index i=0; i<n; i++) {
62    x_l[i] = -1e20;
63    x_u[i] =  1e20;
64  }
65
66  // Set the bounds for the constraints
67  for (Index i=0; i<m; i++) {
68    g_l[i] = 0;
69    g_u[i] = 0;
70  }
71
72  return true;
73}
74
75bool MyADOLC_sparseNLP::get_starting_point(Index n, bool init_x, Number* x,
76                               bool init_z, Number* z_L, Number* z_U,
77                               Index m, bool init_lambda,
78                               Number* lambda)
79{
80  // Here, we assume we only have starting values for x, if you code
81  // your own NLP, you can provide starting values for the others if
82  // you wish.
83  assert(init_x == true);
84  assert(init_z == false);
85  assert(init_lambda == false);
86
87  // set the starting point
88  for (Index i=0; i<n/2; i++) {
89    x[2*i] = -1.2;
90    x[2*i+1] = 1.;
91  }
92  if (n != 2*(n/2)) {
93    x[n-1] = -1.2;
94  }
95
96  return true;
97}
98
99template<class T> bool  MyADOLC_sparseNLP::eval_obj(Index n, const T *x, T& obj_value)
100{
101  T a1, a2;
102  obj_value = 0.;
103  for (Index i=0; i<n-1; i++) {
104    a1 = x[i]*x[i]-x[i+1];
105    a2 = x[i] - 1.;
106    obj_value += 100.*a1*a1 + a2*a2;
107  }
108
109  return true;
110}
111
112template<class T> bool  MyADOLC_sparseNLP::eval_constraints(Index n, const T *x, Index m, T* g)
113{
114  for (Index i=0; i<m; i++) {
115    g[i] = 3.*pow(x[i+1],3.) + 2.*x[i+2] - 5.
116           + sin(x[i+1]-x[i+2])*sin(x[i+1]+x[i+2]) + 4.*x[i+1]
117           - x[i]*exp(x[i]-x[i+1]) - 3.;
118  }
119
120  return true;
121}
122
123//*************************************************************************
124//
125//
126//         Nothing has to be changed below this point !!
127//
128//
129//*************************************************************************
130
131
132bool MyADOLC_sparseNLP::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
133{
134  eval_obj(n,x,obj_value);
135
136  return true;
137}
138
139bool MyADOLC_sparseNLP::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
140{
141
142  gradient(tag_f,n,x,grad_f);
143
144  return true;
145}
146
147bool MyADOLC_sparseNLP::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
148{
149
150  eval_constraints(n,x,m,g);
151
152  return true;
153}
154
155bool MyADOLC_sparseNLP::eval_jac_g(Index n, const Number* x, bool new_x,
156                       Index m, Index nele_jac, Index* iRow, Index *jCol,
157                       Number* values)
158{
159
160  if (values == NULL) {
161    // return the structure of the jacobian
162
163    for(Index idx=0; idx<nnz_jac; idx++)
164      {
165        iRow[idx] = rind_g[idx];
166        jCol[idx] = cind_g[idx];
167      }
168  }
169  else {
170    // return the values of the jacobian of the constraints
171
172    sparse_jac(tag_g, m, n, 1, x, &nnz_jac, &rind_g, &cind_g, &jacval, options_g); 
173
174    for(Index idx=0; idx<nnz_jac; idx++)
175      {
176        values[idx] = jacval[idx];
177
178      }
179  }
180  return true;
181}
182
183bool MyADOLC_sparseNLP::eval_h(Index n, const Number* x, bool new_x,
184                   Number obj_factor, Index m, const Number* lambda,
185                   bool new_lambda, Index nele_hess, Index* iRow,
186                   Index* jCol, Number* values)
187{
188
189  if (values == NULL) {
190    // return the structure. This is a symmetric matrix, fill the lower left
191    // triangle only.
192
193    for(Index idx=0; idx<nnz_L; idx++)
194      {
195        iRow[idx] = rind_L[idx];
196        jCol[idx] = cind_L[idx];
197      }
198  }
199  else {
200    // return the values. This is a symmetric matrix, fill the lower left
201    // triangle only
202
203    for(Index idx = 0; idx<n ; idx++)
204      x_lam[idx] = x[idx];
205    for(Index idx = 0; idx<m ; idx++)
206      x_lam[n+idx] = lambda[idx];
207    x_lam[n+m] = obj_factor;
208
209    sparse_hess(tag_L, n+m+1, 1, x_lam, &nnz_L_total, &rind_L_total, &cind_L_total, &hessval, options_L);
210     
211    Index idx = 0;
212    for(Index idx_total = 0; idx_total <nnz_L_total ; idx_total++)
213      {
214        if((rind_L_total[idx_total] < (unsigned int) n) && (cind_L_total[idx_total] < (unsigned int) n))
215          {
216            values[idx] = hessval[idx_total];
217            idx++;
218          }
219      }
220  }
221
222  return true;
223}
224
225void MyADOLC_sparseNLP::finalize_solution(SolverReturn status,
226                              Index n, const Number* x, const Number* z_L, const Number* z_U,
227                              Index m, const Number* g, const Number* lambda,
228                              Number obj_value,
229                              const IpoptData* ip_data,
230                              IpoptCalculatedQuantities* ip_cq)
231{
232
233  printf("\n\nObjective value\n");
234  printf("f(x*) = %e\n", obj_value);
235
236// memory deallocation of ADOL-C variables
237
238  delete x_lam;
239  delete rind_g;
240  delete cind_g;
241  delete rind_L;
242  delete cind_L;
243  delete rind_L_total;
244  delete cind_L_total;
245  delete jacval;
246  delete hessval;
247
248  for (int i=0;i<n+m+1;i++) {
249     free(HP_t[i]);
250   }
251  free(HP_t);
252}
253
254
255//***************    ADOL-C part ***********************************
256
257void MyADOLC_sparseNLP::generate_tapes(Index n, Index m, Index& nnz_jac_g, Index& nnz_h_lag)
258{
259  Number *xp    = new double[n];
260  Number *lamp  = new double[m];
261  Number *zl    = new double[m];
262  Number *zu    = new double[m];
263
264  adouble *xa   = new adouble[n];
265  adouble *g    = new adouble[m];
266  adouble *lam  = new adouble[m];
267  adouble sig;
268  adouble obj_value;
269 
270  double dummy;
271  double *jacval;
272
273  int i,j,k,l,ii;
274
275  x_lam   = new double[n+m+1];
276
277  get_starting_point(n, 1, xp, 0, zl, zu, m, 0, lamp);
278
279  trace_on(tag_f);
280   
281    for(Index idx=0;idx<n;idx++)
282      xa[idx] <<= xp[idx];
283
284    eval_obj(n,xa,obj_value);
285
286    obj_value >>= dummy;
287
288  trace_off();
289 
290  trace_on(tag_g);
291   
292    for(Index idx=0;idx<n;idx++)
293      xa[idx] <<= xp[idx];
294
295    eval_constraints(n,xa,m,g);
296
297
298    for(Index idx=0;idx<m;idx++)
299      g[idx] >>= dummy;
300
301  trace_off();
302
303  trace_on(tag_L);
304   
305    for(Index idx=0;idx<n;idx++)
306      xa[idx] <<= xp[idx];
307    for(Index idx=0;idx<m;idx++)
308      lam[idx] <<= 1.0;
309    sig <<= 1.0;
310
311    eval_obj(n,xa,obj_value);
312
313    obj_value *= sig;
314    eval_constraints(n,xa,m,g);
315 
316    for(Index idx=0;idx<m;idx++)
317      obj_value += g[idx]*lam[idx];
318
319    obj_value >>= dummy;
320
321  trace_off();
322
323  rind_g = NULL; 
324  cind_g = NULL;
325
326  options_g[0] = 0;          /* sparsity pattern by index domains (default) */ 
327  options_g[1] = 0;          /*                         safe mode (default) */ 
328  options_g[2] = -1;         /*                     &jacval is not computed */ 
329  options_g[3] = 0;          /*                column compression (default) */ 
330 
331  this->jacval=NULL;
332  this->hessval=NULL;
333  sparse_jac(tag_g, m, n, 0, xp, &nnz_jac, &rind_g, &cind_g, &jacval, options_g); 
334
335  options_g[2] = 0;
336  nnz_jac_g = nnz_jac;
337
338  unsigned int  **JP_f=NULL;                /* compressed block row storage */
339  unsigned int  **JP_g=NULL;                /* compressed block row storage */
340  unsigned int  **HP_f=NULL;                /* compressed block row storage */
341  unsigned int  **HP_g=NULL;                /* compressed block row storage */
342  unsigned int  *HP_length=NULL;            /* length of arrays */
343  unsigned int  *temp=NULL;                 /* help array */
344
345  int ctrl_H;
346
347  JP_f = (unsigned int **) malloc(sizeof(unsigned int*));
348  JP_g = (unsigned int **) malloc(m*sizeof(unsigned int*));
349  HP_f = (unsigned int **) malloc(n*sizeof(unsigned int*)); 
350  HP_g = (unsigned int **) malloc(n*sizeof(unsigned int*));
351  HP_t = (unsigned int **) malloc((n+m+1)*sizeof(unsigned int*));
352  HP_length = (unsigned int *) malloc((n)*sizeof(unsigned int));
353  ctrl_H = 0;
354
355  hess_pat(tag_f, n, xp, HP_f, ctrl_H);
356
357  indopro_forward_safe(tag_f, 1, n, xp, JP_f);
358  indopro_forward_safe(tag_g, m, n, xp, JP_g);
359  nonl_ind_forward_safe(tag_g, m, n, xp, HP_g);
360
361
362  for (i=0;i<n;i++) 
363    {
364      if (HP_f[i][0]+HP_g[i][0]!=0)
365        {
366          if (HP_f[i][0]==0)
367            {
368              HP_t[i] = (unsigned int *) malloc((HP_g[i][0]+HPOFF)*sizeof(unsigned int));
369              for(j=0;j<=(int) HP_g[i][0];j++)
370                {
371                  HP_t[i][j] = HP_g[i][j];
372                }
373              HP_length[i] = HP_g[i][0]+HPOFF;
374            }
375          else
376            {
377              if (HP_g[i][0]==0)
378                {
379                  HP_t[i] = (unsigned int *) malloc((HP_f[i][0]+HPOFF)*sizeof(unsigned int));
380                  for(j=0;j<=(int) HP_f[i][0];j++)
381                    {
382                      HP_t[i][j] = HP_f[i][j];
383                    }
384                  HP_length[i] = HP_f[i][0]+HPOFF;
385                }
386              else
387                {
388                  HP_t[i] = (unsigned int *) malloc((HP_f[i][0]+HP_g[i][0]+HPOFF)*sizeof(unsigned int));
389                  k = l = j = 1;
390                  while ((k<=(int) HP_f[i][0]) && (l <= (int) HP_g[i][0]))
391                    {
392                      if (HP_f[i][k] < HP_g[i][l])
393                        {
394                          HP_t[i][j]=HP_f[i][k];
395                          j++; k++;
396                        }
397                      else
398                        {
399                          if (HP_f[i][k] == HP_g[i][l])
400                            {
401                              HP_t[i][j]=HP_f[i][k];
402                              l++;j++;k++;
403                            }
404                          else
405                            {
406                              HP_t[i][j]=HP_g[i][l];
407                              j++;l++;               
408                            }
409                        }
410                    } // end while
411                  for(ii=k;ii<=(int) HP_f[i][0];ii++)
412                    {
413                      HP_t[i][j] = HP_f[i][ii];
414                      j++;
415                    }
416                  for(ii=l;ii<=(int) HP_g[i][0];ii++)
417                    {
418                      HP_t[i][j] = HP_g[i][ii];
419                      j++;
420                    }
421                 
422                }
423            }
424          HP_t[i][0]=j-1;
425          HP_length[i] = HP_f[i][0]+HP_g[i][0]+HPOFF;
426        }
427      else
428        {
429          HP_t[i] = (unsigned int *) malloc((HPOFF+1)*sizeof(unsigned int));
430          HP_t[i][0]=0;
431          HP_length[i]=HPOFF;
432        }
433    }   
434
435  for (i=0;i<m;i++) 
436    {
437      HP_t[n+i] = (unsigned int *) malloc((JP_g[i][0]+1)*sizeof(unsigned int));
438      HP_t[n+i][0]=JP_g[i][0];
439      for(j=1;j<= (int) JP_g[i][0];j++)
440        {
441          HP_t[n+i][j]=JP_g[i][j];
442          if (HP_length[JP_g[i][j]] < HP_t[JP_g[i][j]][0]+1)
443            {
444              temp = (unsigned int *) malloc((HP_t[JP_g[i][j]][0])*sizeof(unsigned int));
445              for(l=0;l<=(int)HP_t[JP_g[i][j]][0];l++)
446                temp[l] = HP_t[JP_g[i][j]][l];
447              free(HP_t[JP_g[i][j]]);
448              HP_t[JP_g[i][j]] = (unsigned int *) malloc(2*HP_length[JP_g[i][j]]*sizeof(unsigned int));
449              HP_length[JP_g[i][j]] = 2*HP_length[JP_g[i][j]];
450              for(l=0;l<=(int)temp[0];l++)
451                HP_t[JP_g[i][j]][l] =temp[l];
452              free(temp);
453            }
454          HP_t[JP_g[i][j]][0] = HP_t[JP_g[i][j]][0]+1;
455          HP_t[JP_g[i][j]][HP_t[JP_g[i][j]][0]] = i+n;
456        }
457    }
458
459
460  for(j=1;j<= (int) JP_f[0][0];j++)
461    {
462      if (HP_length[JP_f[0][j]] < HP_t[JP_f[0][j]][0]+1)
463        {
464          temp = (unsigned int *) malloc((HP_t[JP_f[0][j]][0])*sizeof(unsigned int));
465          for(l=0;l<=(int)HP_t[JP_f[0][j]][0];l++)
466            temp[l] = HP_t[JP_f[0][j]][l];
467          free(HP_t[JP_f[0][j]]);
468          HP_t[JP_f[0][j]] = (unsigned int *) malloc(2*HP_length[JP_f[0][j]]*sizeof(unsigned int));
469          HP_length[JP_f[0][j]] = 2*HP_length[JP_f[0][j]];
470          for(l=0;l<=(int)temp[0];l++)
471            HP_t[JP_f[0][j]][l] =temp[l];
472          free(temp);
473        }
474      HP_t[JP_f[0][j]][0] = HP_t[JP_f[0][j]][0]+1;
475      HP_t[JP_f[0][j]][HP_t[JP_f[0][j]][0]] = n+m;
476    }
477
478
479  HP_t[n+m] = (unsigned int *) malloc((JP_f[0][0]+2)*sizeof(unsigned int));
480  HP_t[n+m][0]=JP_f[0][0]+1;
481  for(j=1;j<= (int) JP_f[0][0];j++)
482    HP_t[n+m][j]=JP_f[0][j];
483  HP_t[n+m][JP_f[0][0]+1]=n+m;
484
485
486
487  set_HP(tag_L,n+m+1,HP_t);
488
489  nnz_h_lag = 0;
490   for (i=0;i<n;i++) {
491    for (j=1;j<=(int) HP_t[i][0];j++)
492      if ((int) HP_t[i][j] <= i)
493        nnz_h_lag++;
494     free(HP_f[i]);
495     free(HP_g[i]);
496   }
497  nnz_L = nnz_h_lag;
498
499  options_L[0] = 0;         
500  options_L[1] = 1;       
501
502  sparse_hess(tag_L, n+m+1, -1, xp, &nnz_L_total, &rind_L_total, &cind_L_total, &hessval, options_L);
503
504  rind_L = new unsigned int[nnz_L];
505  cind_L = new unsigned int[nnz_L];
506  rind_L_total = new unsigned int[nnz_L_total];
507  cind_L_total = new unsigned int[nnz_L_total];
508
509  unsigned int ind = 0;
510
511  for (int i=0;i<n;i++) 
512    for (unsigned int j=1;j<=HP_t[i][0];j++)
513      {
514        if (((int) HP_t[i][j]>=i) &&((int) HP_t[i][j]<n)) 
515          {
516            rind_L[ind] = i;
517            cind_L[ind++] = HP_t[i][j];
518          }
519      }
520
521   ind = 0;
522   for (int i=0;i<n+m+1;i++) 
523     for (unsigned int j=1;j<=HP_t[i][0];j++)
524       {
525        if ((int) HP_t[i][j]>=i) 
526          {
527            rind_L_total[ind] = i;
528            cind_L_total[ind++] = HP_t[i][j];
529          }
530       }
531
532  for (i=0;i<m;i++) {
533     free(JP_g[i]);
534   }
535
536  free(JP_f[0]);
537  free(JP_f);
538  free(JP_g);
539  free(HP_f);
540  free(HP_g);
541  free(HP_length);
542
543  delete[] lam;
544  delete[] g;
545  delete[] xa;
546  delete[] zu;
547  delete[] zl;
548  delete[] lamp;
549  delete[] xp;
550}
Note: See TracBrowser for help on using the repository browser.