source: branches/devel/Bonmin/src/Interfaces/BonCurvatureEstimator.cpp @ 142

Last change on this file since 142 was 142, checked in by andreasw, 13 years ago

committed changes in curvature branching heuristic

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 19.5 KB
Line 
1// Copyright (C) 2006 International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: BonCurvatureEstimator.cpp 142 2006-11-29 20:54:55Z andreasw $
6//
7// Author:   Andreas Waechter                 IBM    2006-10-11
8
9#include "BonCurvatureEstimator.hpp"
10#include "IpGenTMatrix.hpp"
11#include "IpIdentityMatrix.hpp"
12#include "IpZeroMatrix.hpp"
13#include "IpDenseVector.hpp"
14#include "IpBlas.hpp"
15#include <cmath>
16#ifdef HAVE_MA27
17# include "IpMa27TSolverInterface.hpp"
18#endif
19#ifdef HAVE_MA57
20# include "IpMa57TSolverInterface.hpp"
21#endif
22#ifdef HAVE_MC19
23# include "IpMc19TSymScalingMethod.hpp"
24#endif
25#ifdef HAVE_PARDISO
26# include "IpPardisoSolverInterface.hpp"
27#endif
28#ifdef HAVE_TAUCS
29# include "IpTAUCSSolverInterface.hpp"
30#endif
31#ifdef HAVE_WSMP
32# include "IpWsmpSolverInterface.hpp"
33#endif
34#ifdef HAVE_MUMPS
35# include "IpMumpsSolverInterface.hpp"
36#endif
37
38namespace Bonmin
39{
40  using namespace Ipopt;
41
42  // ToDo: Consider NLP scaling?
43
44  CurvatureEstimator::CurvatureEstimator
45    (SmartPtr<Journalist> jnlst,
46     SmartPtr<OptionsList> options,
47     SmartPtr<TNLP> tnlp)
48      :
49      jnlst_(jnlst),
50      options_(options),
51      prefix_(""),
52      tnlp_(tnlp),
53      grad_f_(NULL),
54      irows_jac_(NULL),
55      jcols_jac_(NULL),
56      jac_vals_(NULL),
57      irows_hess_(NULL),
58      jcols_hess_(NULL),
59      hess_vals_(NULL),
60      x_l_(NULL),
61      x_u_(NULL),
62      g_l_(NULL),
63      g_u_(NULL),
64      x_free_map_(NULL),
65      g_fixed_map_(NULL),
66      lambda_(NULL),
67      initialized_(false)
68  {
69    DBG_ASSERT(IsValid(jnlst));
70    DBG_ASSERT(IsValid(options));
71    DBG_ASSERT(IsValid(tnlp));
72
73    ////////////////////////////////////////////////////////////////////
74    // Create a strategy object for solving the linear system for the //
75    // projection matrix                                              //
76    ////////////////////////////////////////////////////////////////////
77
78    // The following linear are taken from AlgBuilder in Ipopt
79    SmartPtr<SparseSymLinearSolverInterface> SolverInterface;
80    std::string linear_solver;
81    options->GetStringValue("linear_solver", linear_solver, prefix_);
82    if (linear_solver=="ma27") {
83#ifdef HAVE_MA27
84      SolverInterface = new Ma27TSolverInterface();
85#else
86
87      THROW_EXCEPTION(OPTION_INVALID,
88                      "Selected linear solver MA27 not available.");
89#endif
90
91    }
92    else if (linear_solver=="ma57") {
93#ifdef HAVE_MA57
94      SolverInterface = new Ma57TSolverInterface();
95#else
96
97      THROW_EXCEPTION(OPTION_INVALID,
98                      "Selected linear solver MA57 not available.");
99#endif
100
101    }
102    else if (linear_solver=="pardiso") {
103#ifdef HAVE_PARDISO
104      SolverInterface = new PardisoSolverInterface();
105#else
106
107      THROW_EXCEPTION(OPTION_INVALID,
108                      "Selected linear solver Pardiso not available.");
109#endif
110
111    }
112    else if (linear_solver=="taucs") {
113#ifdef HAVE_TAUCS
114      SolverInterface = new TAUCSSolverInterface();
115#else
116
117      THROW_EXCEPTION(OPTION_INVALID,
118                      "Selected linear solver TAUCS not available.");
119#endif
120
121    }
122    else if (linear_solver=="wsmp") {
123#ifdef HAVE_WSMP
124      SolverInterface = new WsmpSolverInterface();
125#else
126
127      THROW_EXCEPTION(OPTION_INVALID,
128                      "Selected linear solver WSMP not available.");
129#endif
130
131    }
132    else if (linear_solver=="mumps") {
133#ifdef HAVE_MUMPS
134      SolverInterface = new MumpsSolverInterface();
135#else
136
137      THROW_EXCEPTION(OPTION_INVALID,
138                      "Selected linear solver MUMPS not available.");
139#endif
140
141    }
142
143    SmartPtr<TSymScalingMethod> ScalingMethod;
144    std::string linear_system_scaling;
145    if (!options->GetStringValue("linear_system_scaling",
146                                 linear_system_scaling, prefix_)) {
147      // By default, don't use mc19 for non-HSL solvers
148      if (linear_solver!="ma27" && linear_solver!="ma57") {
149        linear_system_scaling="none";
150      }
151    }
152    if (linear_system_scaling=="mc19") {
153#ifdef HAVE_MC19
154      ScalingMethod = new Mc19TSymScalingMethod();
155#else
156
157      THROW_EXCEPTION(OPTION_INVALID,
158                      "Selected linear system scaling method MC19 not available.");
159#endif
160
161    }
162
163    tsymlinearsolver_ = new TSymLinearSolver(SolverInterface, ScalingMethod);
164    // End of lines from AlgBuilder
165  }
166
167  CurvatureEstimator::~CurvatureEstimator()
168  {
169    if (initialized_) {
170      delete [] grad_f_;
171      delete [] irows_jac_;
172      delete [] jcols_jac_;
173      delete [] jac_vals_;
174      delete [] irows_hess_;
175      delete [] jcols_hess_;
176      delete [] hess_vals_;
177      delete [] x_l_;
178      delete [] x_u_;
179      delete [] g_l_;
180      delete [] g_u_;
181      delete [] x_free_map_;
182      delete [] g_fixed_map_;
183      delete [] lambda_;
184    }
185  }
186
187  bool CurvatureEstimator::Initialize()
188  {
189    DBG_ASSERT(!initialized_);
190
191    //////////////////////////////////////
192    // Prepare internal data structures //
193    //////////////////////////////////////
194
195    // Get sizes
196    TNLP::IndexStyleEnum index_style;
197    if (!tnlp_->get_nlp_info(n_, m_, nnz_jac_, nnz_hess_, index_style)) {
198      return false;
199    }
200
201    // Space for gradient
202    grad_f_ = new Number[n_];
203
204    // Get nonzero entries in the matrices
205    irows_jac_ = new Index[nnz_jac_];
206    jcols_jac_ = new Index[nnz_jac_];
207    if (!tnlp_->eval_jac_g(n_, NULL, false, m_, nnz_jac_,
208                           irows_jac_, jcols_jac_, NULL)) {
209      return false;
210    }
211    if (index_style == TNLP::FORTRAN_STYLE) {
212      for (Index i=0; i<nnz_jac_; i++) {
213        irows_jac_[i]--;
214        jcols_jac_[i]--;
215      }
216    }
217    jac_vals_ = new Number[nnz_jac_];
218
219    irows_hess_ = new Index[nnz_hess_];
220    jcols_hess_ = new Index[nnz_hess_];
221    if (!tnlp_->eval_h(n_, NULL, false, 1., m_, NULL, false, nnz_hess_,
222                       irows_hess_, jcols_hess_, NULL)) {
223      return false;
224    }
225    if (index_style == TNLP::FORTRAN_STYLE) {
226      for (Index i=0; i<nnz_hess_; i++) {
227        irows_hess_[i]--;
228        jcols_hess_[i]--;
229      }
230    }
231    hess_vals_ = NULL; // We set it to NULL, so that we know later
232    // that we still need to compute the values
233
234    // Space for bounds
235    x_l_ = new Number[n_];
236    x_u_ = new Number[n_];
237    g_l_ = new Number[m_];
238    g_u_ = new Number[m_];
239
240    // Get space for the activities maps
241    x_free_map_ = new Index[n_];
242    g_fixed_map_ = new Index[m_];
243
244    // Get space for the multipliers
245    lambda_ = new Number[m_];
246
247    initialized_ = true;
248    return true;
249  }
250
251  bool
252  CurvatureEstimator::ComputeNullSpaceCurvature(
253    bool new_bounds,
254    int n,
255    const Number* x,
256    bool new_x,
257    const Number* z_L,
258    const Number* z_U,
259    int m,
260    const Number* lam,
261    bool new_mults,
262    const Number* orig_d,
263    Number* projected_d,
264    Number& gradLagTd,
265    Number& dTHLagd)
266  {
267    if (!initialized_) {
268      Initialize();
269      new_bounds = true;
270      new_mults = true;
271    }
272    //DELETEME
273    new_bounds = true;
274
275    DBG_ASSERT(n == n_);
276    DBG_ASSERT(m == m_);
277
278    // If necessary, get new Jacobian values (for the original matrix)
279    if (new_x) {
280      if (!tnlp_->eval_jac_g(n_, x, new_x, m_, nnz_jac_,
281                             NULL, NULL, jac_vals_)) {
282        return false;
283      }
284    }
285
286    if (new_bounds) {
287      // Get bounds
288      if (!tnlp_->get_bounds_info(n_, x_l_, x_u_, m_, g_l_, g_u_)) {
289        return false;
290      }
291    }
292
293    // If necessary, determine new activities
294    bool new_activities = false;
295    if (new_bounds || new_mults) {
296      new_activities = true;
297      active_x_.clear();
298      const Number zTol = 1e300; //1e-2; // ToDo: make this more elaborate
299      jnlst_->Printf(J_MOREDETAILED, J_NLP,
300                     "List of variables considered fixed (with z_L and z_U)\n");
301      for (Index i=0; i<n; i++) {
302        if(z_L[i] > zTol || z_U[i] > zTol) {
303          if (z_L[i] > z_U[i]) {
304            active_x_.push_back(-(i+1));
305            DBG_ASSERT(x_l_[i] > -1e19);
306          }
307          else {
308            active_x_.push_back(i+1);
309            DBG_ASSERT(x_u_[i] < 1e19);
310          }
311          jnlst_->Printf(J_MOREDETAILED, J_NLP,
312                         "x[%5d] (%e,%e)\n", i, z_L[i], z_U[i]);
313          DBG_ASSERT(orig_d[i] == 0.);
314        }
315      }
316
317      active_g_.clear();
318      // Compute the product of the direction with the constraint Jacobian
319      // This could be done more efficient if we have d in sparse format
320      Number* jacTd = new Number[m];
321      const Number zero = 0.;
322      IpBlasDcopy(m, &zero, 0, jacTd, 1);
323      for (Index i=0; i<nnz_jac_; i++) {
324        const Index& irow = irows_jac_[i];
325        const Index& jcol = jcols_jac_[i];
326        jacTd[irow] += jac_vals_[i]*orig_d[jcol];
327      }
328
329      const Number lamTol = 1e-4;
330      jnlst_->Printf(J_MOREDETAILED, J_NLP,
331                     "List of constraints considered fixed (with lam and jacTd)\n");
332      for (Index j=0; j<m; j++) {
333        if (g_l_[j] < g_u_[j] && fabs(lam[j]) > lamTol) {
334          if (lam[j]*jacTd[j] > 0) {
335            if (lam[j] < 0.) {
336              active_g_.push_back(-(j+1));
337              DBG_ASSERT(g_l_[j] > -1e19);
338            }
339            else {
340              active_g_.push_back(j+1);
341              DBG_ASSERT(g_u_[j] < 1e19);
342            }
343            //      active_g_.push_back(j+1);
344            jnlst_->Printf(J_MOREDETAILED, J_NLP,
345                           "g[%5d] (%e,%e)\n", j, lam[j], jacTd[j]);
346          }
347        }
348      }
349      delete [] jacTd;
350    }
351
352    // Check if the structure of the matrix has changed
353    if (new_activities) {
354      if (!PrepareNewMatrixStructure(new_activities)) {
355        return false;
356      }
357    }
358
359    bool new_lambda = false;
360    if (new_x || new_activities) {
361      if (!PrepareNewMatrixValues(new_activities, x, new_x)) {
362        return false;
363      }
364
365      // Compute least square multipliers for the given activities
366      if (!tnlp_->eval_grad_f(n_, x, new_x, grad_f_)) {
367        return false;
368      }
369#ifdef lambdas
370      if (!SolveSystem(grad_f_, NULL, NULL, lambda_)) {
371        return false;
372      }
373      IpBlasDscal(m_, -1., lambda_, 1);
374      if (jnlst_->ProduceOutput(J_MOREVECTOR, J_NLP)) {
375        jnlst_->Printf(J_MOREVECTOR, J_NLP,
376                       "Curvature Estimator: least square multiplier:\n");
377        for (Index i=0; i<m_; i++) {
378          jnlst_->Printf(J_MOREVECTOR, J_NLP, "lambda[%5d] = %23.16e\n",
379                         i, lambda_[i]);
380        }
381      }
382      new_lambda = true;
383#endif
384    }
385
386    // Compute the projection of the direction
387    if (!SolveSystem(orig_d, NULL, projected_d, NULL)) {
388      return false;
389    }
390
391    // Sanity check to see if the thing is indeed in the null space
392    // (if the constraint gradients are rank-deficient, the solver
393    // might not have done a good job)
394    Number* trash = new Number[m_];
395    for (Index j=0; j<m_; j++) {
396      trash[j] = 0.;
397    }
398    for (Index i=0; i<nnz_jac_; i++) {
399      const Index &irow = irows_jac_[i];
400      const Index &jcol = jcols_jac_[i];
401      if (x_free_map_[jcol] >= 0 && g_fixed_map_[irow] >= 0) {
402        trash[irow] += jac_vals_[i]*projected_d[jcol];
403      }
404    }
405    if (jnlst_->ProduceOutput(J_MOREVECTOR, J_NLP)) {   
406      for (Index j=0; j<m_; j++) {
407        jnlst_->Printf(J_MOREVECTOR, J_NLP,
408                       "nullspacevector[%5d] = %e\n", j, trash[j]);
409      }
410    }
411    Index imax = IpBlasIdamax(m_, trash, 1);
412    Number max_trash = trash[imax];
413    delete [] trash;
414    const Number max_trash_tol = 1e-6;
415    if (max_trash > max_trash_tol) {
416      return false;
417    }
418
419    if (jnlst_->ProduceOutput(J_MOREVECTOR, J_NLP)) {
420      jnlst_->Printf(J_MOREVECTOR, J_NLP,
421                     "Curvature Estimator: original and projected directions are:\n");
422      for (Index i=0; i<n_; i++) {
423        jnlst_->Printf(J_MOREVECTOR, J_NLP,
424                       "orig_d[%5d] = %23.16e proj_d[%5d] = %23.16e\n",
425                       i, orig_d[i], i, projected_d[i]);
426      }
427    }
428
429    // Compute the product with the gradient of the Lagrangian
430    if (true) {
431      gradLagTd = 0.;
432    }
433    else {
434      gradLagTd = IpBlasDdot(n, projected_d, 1, grad_f_, 1);
435      for (Index i=0; i<nnz_jac_; i++) {
436        const Index &irow = irows_jac_[i];
437        const Index &jcol = jcols_jac_[i];
438        gradLagTd += lambda_[irow]*jac_vals_[i]*projected_d[jcol];
439      }
440    }
441
442    // Compute the product with the Hessian of the Lagrangian
443    //    if (!Compute_dTHLagd(projected_d, x, new_x, lambda_, new_lambda, dTHLagd)) {
444    if (!Compute_dTHLagd(projected_d, x, new_x, lam, new_lambda, dTHLagd)) {
445      return false;
446    }
447
448    return true;
449  }
450
451  bool CurvatureEstimator::PrepareNewMatrixStructure(
452    bool new_activities)
453  {
454    if (new_activities) {
455      // Deterimine which variables are free
456      for (Index i=0; i<n_; i++) {
457        x_free_map_[i] = 0;
458      }
459      // fixed by activities
460      for (std::vector<int>::iterator i=active_x_.begin();
461           i != active_x_.end(); i++) {
462        DBG_ASSERT(*i != 0 && *i<=n_ && *i>=-n_);
463        if (*i<0) {
464          x_free_map_[(-*i)-1] = -1;
465          DBG_ASSERT(x_l_[(-*i)-1] > -1e19);
466        }
467        else {
468          x_free_map_[(*i)-1] = -1;
469          DBG_ASSERT(x_u_[(*i)-1] < 1e19);
470        }
471      }
472      // fixed by bounds
473      for (Index i=0; i<n_; i++) {
474        if (x_l_[i] == x_u_[i]) {
475          x_free_map_[i] = -1;
476        }
477      }
478      // Correct the numbering in the x map and determine number of
479      // free variables
480      Index nx_free_ = 0;
481      for (Index i=0; i<n_; i++) {
482        if (x_free_map_[i] == 0) {
483          x_free_map_[i] = nx_free_;
484          nx_free_++;
485        }
486      }
487
488      // Determine which constraints are fixed
489      for (Index j=0; j<m_; j++) {
490        g_fixed_map_[j] = -1;
491      }
492      // fixed by activities
493      for (std::vector<int>::iterator i=active_g_.begin();
494           i != active_g_.end(); i++) {
495        DBG_ASSERT(*i != 0 && *i<=m_ && *i>=-m_);
496        if (*i<0) {
497          g_fixed_map_[(-*i)-1] = 0;
498          DBG_ASSERT(g_l_[(-*i)-1] > -1e19); //ToDo look at option?
499        }
500        else {
501          g_fixed_map_[(*i)-1] = 0;
502          DBG_ASSERT(g_u_[(*i)-1] < 1e19);
503        }
504      }
505      // fixed by bounds
506      for (Index j=0; j<m_; j++) {
507        if (g_l_[j] == g_u_[j]) {
508          g_fixed_map_[j] = 0;
509        }
510      }
511      // Correct the numbering in the g map and determine number of
512      // fixed constraints
513      Index ng_fixed_ = 0;
514      for (Index j=0; j<m_; j++) {
515        if (g_fixed_map_[j] == 0) {
516          g_fixed_map_[j] = ng_fixed_;
517          ng_fixed_++;
518        }
519      }
520
521      // Determine the row and column indices for the Jacobian of the fixed
522      // constraints in the space of the free variables
523      Index* iRows = new Index[nnz_jac_];
524      Index* jCols = new Index[nnz_jac_];
525      Index nnz_proj_jac = 0;
526      for (Index i=0; i<nnz_jac_; i++) {
527        const Index &irow = irows_jac_[i];
528        const Index &jcol = jcols_jac_[i];
529        if (x_free_map_[jcol] >= 0 && g_fixed_map_[irow] >= 0) {
530          iRows[nnz_proj_jac] = g_fixed_map_[irow] + 1;
531          jCols[nnz_proj_jac] = x_free_map_[jcol] + 1;
532          nnz_proj_jac++;
533        }
534      }
535
536      // Create the matrix space for the Jacobian matrices
537      SmartPtr<GenTMatrixSpace> proj_jac_space =
538        new GenTMatrixSpace(ng_fixed_, nx_free_, nnz_proj_jac, iRows, jCols);
539      delete [] iRows;
540      delete [] jCols;
541
542      // Create Matrix space for the projection matrix
543      comp_proj_matrix_space_ =
544        new CompoundSymMatrixSpace(2, nx_free_+ng_fixed_);
545      comp_proj_matrix_space_->SetBlockDim(0, nx_free_);
546      comp_proj_matrix_space_->SetBlockDim(1, ng_fixed_);
547      SmartPtr<SymMatrixSpace> identity_space =
548        new IdentityMatrixSpace(nx_free_);
549      comp_proj_matrix_space_->SetCompSpace(0, 0, *identity_space, true);
550      comp_proj_matrix_space_->SetCompSpace(1, 0, *proj_jac_space, true);
551      //      SmartPtr<MatrixSpace> zero_space =
552      //        new ZeroMatrixSpace(ng_fixed_, ng_fixed_);
553      //      comp_proj_matrix_space_->SetCompSpace(1, 1, *zero_space, true);
554
555      // Create a Vector space for the rhs and sol
556      comp_vec_space_ = new CompoundVectorSpace(2, nx_free_+ng_fixed_);
557      SmartPtr<DenseVectorSpace> x_space = new DenseVectorSpace(nx_free_);
558      comp_vec_space_->SetCompSpace(0, *x_space);
559      SmartPtr<DenseVectorSpace> g_space = new DenseVectorSpace(ng_fixed_);
560      comp_vec_space_->SetCompSpace(1, *g_space);
561    } 
562
563    return true;
564  }
565
566  bool CurvatureEstimator::PrepareNewMatrixValues(
567    bool new_activities,
568    const Number* x,
569    bool new_x)
570  {
571    if (new_x || new_activities) {
572      comp_proj_matrix_ = comp_proj_matrix_space_->MakeNewCompoundSymMatrix();
573      SmartPtr<Matrix> jac = comp_proj_matrix_->GetCompNonConst(1, 0);
574      SmartPtr<GenTMatrix> tjac = dynamic_cast<GenTMatrix*> (GetRawPtr(jac));
575      Number* vals = tjac->Values();
576      Index inz=0;
577      for (Index i=0; i<nnz_jac_; i++) {
578        Index irow = irows_jac_[i];
579        Index jcol = jcols_jac_[i];
580        if (x_free_map_[jcol] >= 0 && g_fixed_map_[irow] >= 0) {
581          vals[inz++] = jac_vals_[i];
582        }
583      }
584      //DELETEME
585      // comp_proj_matrix_->Print(*jnlst_,J_WARNING,J_MAIN,"compmat");
586      DBG_ASSERT(inz == tjac->Nonzeros());
587
588      // We need to reset the linear solver object, so that it knows
589      // that the structure of the linear system has changed
590      tsymlinearsolver_->ReducedInitialize(*jnlst_, *options_, prefix_);
591    }
592
593    return true;
594  }
595
596  bool CurvatureEstimator::SolveSystem(const Number* rhs_x,
597                                       const Number* rhs_g,
598                                       Number* sol_x, Number* sol_g)
599  {
600    // Create a vector for the RHS
601    SmartPtr<CompoundVector> rhs = comp_vec_space_->MakeNewCompoundVector();
602    SmartPtr<Vector> vrhs_x = rhs->GetCompNonConst(0);
603    SmartPtr<Vector> vrhs_g = rhs->GetCompNonConst(1);
604    // Now fill this vector with the values, extracting the relevant entries
605    if (rhs_x) {
606      SmartPtr<DenseVector> drhs_x =
607        dynamic_cast<DenseVector*> (GetRawPtr(vrhs_x));
608      Number* xvals = drhs_x->Values();
609      for (Index i=0; i<n_; i++) {
610        if (x_free_map_[i]>=0) {
611          xvals[x_free_map_[i]] = rhs_x[i];
612        }
613      }
614    }
615    else {
616      vrhs_x->Set(0.);
617    }
618    if (rhs_g) {
619      SmartPtr<DenseVector> drhs_g =
620        dynamic_cast<DenseVector*> (GetRawPtr(vrhs_g));
621      Number* gvals = drhs_g->Values();
622      for (Index j=0; j<m_; j++) {
623        if (g_fixed_map_[j]>=0) {
624          gvals[g_fixed_map_[j]] = rhs_g[j];
625        }
626      }
627    }
628    else {
629      vrhs_g->Set(0.);
630    }
631
632    // Solve the linear system
633    SmartPtr<CompoundVector> sol = comp_vec_space_->MakeNewCompoundVector();
634    ESymSolverStatus solver_retval =
635      tsymlinearsolver_->Solve(*comp_proj_matrix_, *rhs, *sol, false, 0);
636    // For now we try, if the system is reported to be singular, to
637    // still get a solution; ToDo
638    if (solver_retval == SYMSOLVER_SINGULAR) {
639      jnlst_->Printf(J_DETAILED, J_NLP,
640                     "Projection matrix reported to be singular, but we try to obtain a solution anyway.\n");
641      solver_retval =
642        tsymlinearsolver_->Solve(*comp_proj_matrix_, *rhs, *sol, false, 0);
643    }
644
645    if (solver_retval != SYMSOLVER_SUCCESS) {
646      // DELETEME
647      printf("Return code from Solver is %d\n", solver_retval);
648      return false;
649    }
650
651    // and get the solution out of it
652    if (sol_x) {
653      SmartPtr<Vector> vsol_x = sol->GetCompNonConst(0);
654      SmartPtr<const DenseVector> dsol_x =
655        dynamic_cast<const DenseVector*> (GetRawPtr(vsol_x));
656      const Number* xvals = dsol_x->Values();
657      for (Index i=0; i<n_; i++) {
658        if (x_free_map_[i]>=0) {
659          sol_x[i] = xvals[x_free_map_[i]];
660        }
661        else {
662          sol_x[i] = 0.;
663        }
664      }
665    }
666    if (sol_g) {
667      SmartPtr<Vector> vsol_g = sol->GetCompNonConst(1);
668      SmartPtr<const DenseVector> dsol_g =
669        dynamic_cast<const DenseVector*> (GetRawPtr(vsol_g));
670      const Number* gvals = dsol_g->Values();
671      for (Index j=0; j<m_; j++) {
672        if (g_fixed_map_[j]>=0) {
673          sol_g[j] = gvals[g_fixed_map_[j]];
674        }
675        else {
676          sol_g[j] = 0.;
677        }
678      }
679    }
680
681    return true;
682  }
683
684  bool CurvatureEstimator::Compute_dTHLagd(
685    const Number* d, const Number* x, bool new_x, const Number* lambda,
686    bool new_lambda,  Number& dTHLagd)
687  {
688    if (new_x || new_lambda || !hess_vals_) {
689      if (!hess_vals_) {
690        hess_vals_ = new Number[nnz_hess_];
691      }
692      if (!tnlp_->eval_h(n_, x, new_x, 1., m_, lambda, new_lambda, nnz_hess_,
693                         NULL, NULL, hess_vals_)) {
694        return false;
695      }
696    }
697    dTHLagd = 0.;
698    for (Index i=0; i<nnz_hess_; i++) {
699      Index irow = irows_hess_[i];
700      Index jcol = jcols_hess_[i];
701      if (irow == jcol) {
702        dTHLagd += d[irow]*d[irow]*hess_vals_[i];
703      }
704      else {
705        dTHLagd += 2.*d[irow]*d[jcol]*hess_vals_[i];
706      }
707    }
708    return true;
709  }
710
711
712} // namespace Bonmin
Note: See TracBrowser for help on using the repository browser.