/home/coin/SVN-release/OS-2.4.0/Bonmin/src/Interfaces/BonCurvatureEstimator.cpp

Go to the documentation of this file.
00001 // Copyright (C) 2006, 2007 International Business Machines and others.
00002 // All Rights Reserved.
00003 // This code is published under the Common Public License.
00004 //
00005 // $Id: BonCurvatureEstimator.cpp 1757 2010-12-09 07:21:18Z pbonami $
00006 //
00007 // Author:   Andreas Waechter                 IBM    2006-10-11
00008 
00009 #include "BonCurvatureEstimator.hpp"
00010 #include "IpTSymLinearSolver.hpp"
00011 #include "IpGenTMatrix.hpp"
00012 #include "IpIdentityMatrix.hpp"
00013 #include "IpZeroMatrix.hpp"
00014 #include "IpDenseVector.hpp"
00015 #include "IpBlas.hpp"
00016 #include <cmath>
00017 #ifdef HAVE_MA27
00018 # include "IpMa27TSolverInterface.hpp"
00019 #endif
00020 #ifdef HAVE_MA57
00021 # include "IpMa57TSolverInterface.hpp"
00022 #endif
00023 #ifdef HAVE_MC19
00024 # include "IpMc19TSymScalingMethod.hpp"
00025 #endif
00026 #ifdef HAVE_PARDISO
00027 # include "IpPardisoSolverInterface.hpp"
00028 #endif
00029 #ifdef HAVE_TAUCS
00030 # include "IpTAUCSSolverInterface.hpp"
00031 #endif
00032 #ifdef HAVE_WSMP
00033 # include "IpWsmpSolverInterface.hpp"
00034 #endif
00035 #ifdef COIN_HAS_MUMPS
00036 # include "IpMumpsSolverInterface.hpp"
00037 #endif
00038 
00039 namespace Bonmin
00040 {
00041   using namespace Ipopt;
00042 
00043   // ToDo: Consider NLP scaling?
00044 
00045   CurvatureEstimator::CurvatureEstimator
00046     (SmartPtr<Journalist> jnlst,
00047      SmartPtr<OptionsList> options,
00048      SmartPtr<TNLP> tnlp)
00049       :
00050       jnlst_(jnlst),
00051       options_(options),
00052       prefix_(""),
00053       tnlp_(tnlp),
00054       grad_f_(NULL),
00055       irows_jac_(NULL),
00056       jcols_jac_(NULL),
00057       jac_vals_(NULL),
00058       irows_hess_(NULL),
00059       jcols_hess_(NULL),
00060       hess_vals_(NULL),
00061       eq_x_free_map_(NULL),
00062       eq_g_fixed_map_(NULL),
00063       all_x_free_map_(NULL),
00064       all_g_fixed_map_(NULL),
00065       lambda_(NULL),
00066       eq_projected_d_(NULL),
00067       initialized_(false)
00068   {
00069     DBG_ASSERT(IsValid(jnlst));
00070     DBG_ASSERT(IsValid(options));
00071     DBG_ASSERT(IsValid(tnlp));
00072 
00074     // Create a strategy object for solving the linear system for the //
00075     // projection matrix                                              //
00077 
00078     // The following linear are taken from AlgBuilder in Ipopt
00079     SmartPtr<SparseSymLinearSolverInterface> SolverInterface1;
00080     SmartPtr<SparseSymLinearSolverInterface> SolverInterface2;
00081     std::string linear_solver;
00082     options->GetStringValue("linear_solver", linear_solver, prefix_);
00083     if (linear_solver=="ma27") {
00084 #ifdef HAVE_MA27
00085       SolverInterface1 = new Ma27TSolverInterface();
00086       SolverInterface2 = new Ma27TSolverInterface();
00087 #else
00088 
00089       THROW_EXCEPTION(OPTION_INVALID,
00090                       "Selected linear solver MA27 not available.");
00091 #endif
00092 
00093     }
00094     else if (linear_solver=="ma57") {
00095 #ifdef HAVE_MA57
00096       SolverInterface1 = new Ma57TSolverInterface();
00097       SolverInterface2 = new Ma57TSolverInterface();
00098 #else
00099 
00100       THROW_EXCEPTION(OPTION_INVALID,
00101                       "Selected linear solver MA57 not available.");
00102 #endif
00103 
00104     }
00105     else if (linear_solver=="pardiso") {
00106 #ifdef HAVE_PARDISO
00107       SolverInterface1 = new PardisoSolverInterface();
00108       SolverInterface2 = new PardisoSolverInterface();
00109 #else
00110 
00111       THROW_EXCEPTION(OPTION_INVALID,
00112                       "Selected linear solver Pardiso not available.");
00113 #endif
00114 
00115     }
00116     else if (linear_solver=="taucs") {
00117 #ifdef HAVE_TAUCS
00118       SolverInterface1 = new TAUCSSolverInterface();
00119       SolverInterface2 = new TAUCSSolverInterface();
00120 #else
00121 
00122       THROW_EXCEPTION(OPTION_INVALID,
00123                       "Selected linear solver TAUCS not available.");
00124 #endif
00125 
00126     }
00127     else if (linear_solver=="wsmp") {
00128 #ifdef HAVE_WSMP
00129       SolverInterface1 = new WsmpSolverInterface();
00130       SolverInterface2 = new WsmpSolverInterface();
00131 #else
00132 
00133       THROW_EXCEPTION(OPTION_INVALID,
00134                       "Selected linear solver WSMP not available.");
00135 #endif
00136 
00137     }
00138     else if (linear_solver=="mumps") {
00139 #ifdef COIN_HAS_MUMPS
00140       SolverInterface1 = new MumpsSolverInterface();
00141       SolverInterface2 = new MumpsSolverInterface();
00142 #else
00143 
00144       THROW_EXCEPTION(OPTION_INVALID,
00145                       "Selected linear solver MUMPS not available.");
00146 #endif
00147 
00148     }
00149 
00150     SmartPtr<TSymScalingMethod> ScalingMethod1;
00151     SmartPtr<TSymScalingMethod> ScalingMethod2;
00152     std::string linear_system_scaling;
00153     if (!options->GetStringValue("linear_system_scaling",
00154                                  linear_system_scaling, prefix_)) {
00155       // By default, don't use mc19 for non-HSL solvers
00156       if (linear_solver!="ma27" && linear_solver!="ma57") {
00157         linear_system_scaling="none";
00158       }
00159     }
00160     if (linear_system_scaling=="mc19") {
00161 #ifdef HAVE_MC19
00162       ScalingMethod1 = new Mc19TSymScalingMethod();
00163       ScalingMethod2 = new Mc19TSymScalingMethod();
00164 #else
00165 
00166       THROW_EXCEPTION(OPTION_INVALID,
00167                       "Selected linear system scaling method MC19 not available.");
00168 #endif
00169 
00170     }
00171 
00172     eq_tsymlinearsolver_ =
00173       new TSymLinearSolver(SolverInterface1, ScalingMethod1);
00174     all_tsymlinearsolver_ =
00175       new TSymLinearSolver(SolverInterface2, ScalingMethod2);
00176     // End of lines from AlgBuilder
00177   }
00178 
00179   CurvatureEstimator::~CurvatureEstimator()
00180   {
00181     if (initialized_) {
00182       delete [] grad_f_;
00183       delete [] irows_jac_;
00184       delete [] jcols_jac_;
00185       delete [] jac_vals_;
00186       delete [] irows_hess_;
00187       delete [] jcols_hess_;
00188       delete [] hess_vals_;
00189       delete [] eq_x_free_map_;
00190       delete [] eq_g_fixed_map_;
00191       delete [] all_x_free_map_;
00192       delete [] all_g_fixed_map_;
00193       delete [] eq_projected_d_;
00194       delete [] lambda_;
00195     }
00196   }
00197 
00198   bool CurvatureEstimator::Initialize()
00199   {
00200     DBG_ASSERT(!initialized_);
00201 
00203     // Prepare internal data structures //
00205 
00206     // Get sizes
00207     TNLP::IndexStyleEnum index_style;
00208     if (!tnlp_->get_nlp_info(n_, m_, nnz_jac_, nnz_hess_, index_style)) {
00209       return false;
00210     }
00211 
00212     // Space for gradient
00213     delete [] grad_f_;
00214     grad_f_ = NULL;
00215     grad_f_ = new Number[n_];
00216 
00217     // Get nonzero entries in the matrices
00218     delete [] irows_jac_;
00219     delete [] jcols_jac_;
00220     irows_jac_ = NULL;
00221     jcols_jac_ = NULL;
00222     irows_jac_ = new Index[nnz_jac_];
00223     jcols_jac_ = new Index[nnz_jac_];
00224     if (!tnlp_->eval_jac_g(n_, NULL, false, m_, nnz_jac_,
00225                            irows_jac_, jcols_jac_, NULL)) {
00226       return false;
00227     }
00228     if (index_style == TNLP::FORTRAN_STYLE) {
00229       for (Index i=0; i<nnz_jac_; i++) {
00230         irows_jac_[i]--;
00231         jcols_jac_[i]--;
00232       }
00233     }
00234     delete [] jac_vals_;
00235     jac_vals_ = NULL;
00236     jac_vals_ = new Number[nnz_jac_];
00237 
00238     delete [] irows_hess_;
00239     irows_hess_ = NULL;
00240     irows_hess_ = new Index[nnz_hess_];
00241     delete [] jcols_hess_;
00242     jcols_hess_ = NULL;
00243     jcols_hess_ = new Index[nnz_hess_];
00244     if (!tnlp_->eval_h(n_, NULL, false, 1., m_, NULL, false, nnz_hess_,
00245                        irows_hess_, jcols_hess_, NULL)) {
00246       return false;
00247     }
00248     if (index_style == TNLP::FORTRAN_STYLE) {
00249       for (Index i=0; i<nnz_hess_; i++) {
00250         irows_hess_[i]--;
00251         jcols_hess_[i]--;
00252       }
00253     }
00254     delete [] hess_vals_;
00255     hess_vals_ = NULL; // We set it to NULL, so that we know later
00256     // that we still need to compute the values
00257 
00258     // Get space for the activities maps
00259     delete [] eq_x_free_map_;
00260     delete [] eq_g_fixed_map_;
00261     delete [] all_x_free_map_;
00262     delete [] all_g_fixed_map_;
00263     eq_x_free_map_ = NULL;
00264     eq_g_fixed_map_ = NULL;
00265     all_x_free_map_ = NULL;
00266     all_g_fixed_map_ = NULL;
00267     eq_x_free_map_ = new Index[n_];
00268     eq_g_fixed_map_ = new Index[m_];
00269     all_x_free_map_ = new Index[n_];
00270     all_g_fixed_map_ = new Index[m_];
00271 
00272     // Get space for the multipliers
00273     delete [] lambda_;
00274     lambda_ = NULL;
00275     lambda_ = new Number[m_];
00276 
00277     // Get space for projected d
00278     delete [] eq_projected_d_;
00279     eq_projected_d_ = NULL;
00280     eq_projected_d_ = new Number[n_];
00281 
00282     initialized_ = true;
00283     return true;
00284   }
00285 
00286   bool
00287   CurvatureEstimator::ComputeNullSpaceCurvature(
00288     int n,
00289     const Number* x,
00290     bool new_x,
00291     const Number* x_l,
00292     const Number* x_u,
00293     const Number* g_l,
00294     const Number* g_u,
00295     bool new_bounds,
00296     const Number* z_L,
00297     const Number* z_U,
00298     int m,
00299     const Number* lam,
00300     bool new_mults,
00301     const Number* orig_d,
00302     Number* projected_d,
00303     Number& gradLagTd,
00304     Number& dTHLagd)
00305   {
00306 #if 0
00307     printf("new_bounds = %d new_x = %d new_mults = %d\n", new_bounds, new_x, new_mults);
00308   for (int i=0;  i<n; i++) {
00309     printf("x[%3d] = %15.8e orig_d[%3d] = %15.8e z_L[%3d] = %15.8e z_U[%3d] = %15.8e\n",i,x[i],i,orig_d[i],i,z_L[i],i,z_U[i]);
00310   }
00311   for (int i=0; i<m; i++) {
00312     printf("lam[%3d] = %15.8e\n", i, lam[i]);
00313   }
00314 #endif
00315     if (!initialized_) {
00316       Initialize();
00317       new_bounds = true;
00318       new_mults = true;
00319     }
00320     //DELETEME
00321     new_bounds = true;
00322 
00323     DBG_ASSERT(n == n_);
00324     DBG_ASSERT(m == m_);
00325 
00326     // If necessary, get new Jacobian values (for the original matrix)
00327     if (new_x) {
00328       if (!tnlp_->eval_jac_g(n_, x, new_x, m_, nnz_jac_,
00329                              NULL, NULL, jac_vals_)) {
00330         return false;
00331       }
00332     }
00333 
00334     // First we compute the direction projected into the space of only
00335     // the equality constraints
00336     if (new_x) {
00337       std::vector<int> dummy_active_x;
00338       std::vector<int> dummy_active_g;
00339       if (!PrepareNewMatrixStructure(x_l, x_u, g_l, g_u,
00340                                      dummy_active_x, dummy_active_g,
00341                                      eq_nx_free_, eq_x_free_map_,
00342                                      eq_ng_fixed_, eq_g_fixed_map_,
00343                                      eq_comp_proj_matrix_space_,
00344                                      eq_comp_vec_space_)) {
00345         return false;
00346       }
00347 
00348       if (!PrepareNewMatrixValues(eq_x_free_map_, eq_g_fixed_map_,
00349                                   eq_comp_proj_matrix_space_,
00350                                   eq_comp_proj_matrix_,
00351                                   eq_tsymlinearsolver_)) {
00352         return false;
00353       }
00354     }
00355 
00356     if (eq_ng_fixed_>0) {
00357       // Compute the projection of the direction
00358       if (!SolveSystem(orig_d, NULL, eq_projected_d_, NULL,
00359                        eq_x_free_map_, eq_g_fixed_map_,
00360                        eq_comp_vec_space_, eq_comp_proj_matrix_,
00361                        eq_tsymlinearsolver_)) {
00362         return false;
00363       }
00364       orig_d = eq_projected_d_; // This way we don't need to rememeber
00365                                 // how the direction was projected;
00366     }
00367 
00368     // If necessary, determine new activities
00369     bool new_activities = false;
00370     if (new_bounds || new_mults) {
00371       new_activities = true;
00372       active_x_.clear();
00373       if (eq_ng_fixed_>0) {
00374         const Number zTol = 1e-4;
00375         jnlst_->Printf(J_MOREDETAILED, J_NLP,
00376                        "List of variables considered fixed (with orig_d and z)\n");
00377         for (Index i=0; i<n; i++) {
00378           if (x_l[i] < x_u[i]) {
00379             if (orig_d[i]>0. && z_U[i]*orig_d[i]>zTol) {
00380               active_x_.push_back(i+1);
00381               jnlst_->Printf(J_MOREDETAILED, J_NLP,
00382                              "x[%5d] (%e,%e)\n", i, orig_d[i], z_U[i]);
00383               DBG_ASSERT(x_u[i] < 1e19);
00384             }
00385             else if (orig_d[i]<0. && -z_L[i]*orig_d[i]>zTol) {
00386               active_x_.push_back(-(i+1));
00387               jnlst_->Printf(J_MOREDETAILED, J_NLP,
00388                              "x[%5d] (%e,%e)\n", i, orig_d[i], z_L[i]);
00389               DBG_ASSERT(x_l[i] > -1e19);
00390             }
00391           }
00392         }
00393       }
00394 
00395       active_g_.clear();
00396       // Compute the product of the direction with the constraint Jacobian
00397       // This could be done more efficient if we have d in sparse format
00398       Number* jacTd = new Number[m];
00399       const Number zero = 0.;
00400       IpBlasDcopy(m, &zero, 0, jacTd, 1);
00401       for (Index i=0; i<nnz_jac_; i++) {
00402         const Index& irow = irows_jac_[i];
00403         const Index& jcol = jcols_jac_[i];
00404         jacTd[irow] += jac_vals_[i]*orig_d[jcol];
00405       }
00406 
00407       const Number lamTol = 1e-4;
00408       jnlst_->Printf(J_MOREDETAILED, J_NLP,
00409                      "List of constraints considered fixed (with lam and jacTd)\n");
00410       for (Index j=0; j<m; j++) {
00411         if (g_l[j] < g_u[j] && fabs(lam[j]) > lamTol) {
00412           if (lam[j]*jacTd[j] > 0) {
00413             if (lam[j] < 0.) {
00414               active_g_.push_back(-(j+1));
00415               DBG_ASSERT(g_l[j] > -1e19);
00416             }
00417             else {
00418               active_g_.push_back(j+1);
00419               DBG_ASSERT(g_u[j] < 1e19);
00420             }
00421             //      active_g_.push_back(j+1);
00422             jnlst_->Printf(J_MOREDETAILED, J_NLP,
00423                            "g[%5d] (%e,%e)\n", j, lam[j], jacTd[j]);
00424           }
00425         }
00426       }
00427       delete [] jacTd;
00428     }
00429 
00430     // Check if the structure of the matrix has changed
00431     if (new_activities) {
00432       if (!PrepareNewMatrixStructure(x_l, x_u, g_l, g_u,
00433                                      active_x_, active_g_,
00434                                      all_nx_free_, all_x_free_map_,
00435                                      all_ng_fixed_, all_g_fixed_map_,
00436                                      all_comp_proj_matrix_space_,
00437                                      all_comp_vec_space_)) {
00438         return false;
00439       }
00440     }
00441 
00442     bool new_lambda = false;
00443     if (new_x || new_activities) {
00444       if (!PrepareNewMatrixValues(all_x_free_map_, all_g_fixed_map_,
00445                                   all_comp_proj_matrix_space_,
00446                                   all_comp_proj_matrix_,
00447                                   all_tsymlinearsolver_)) {
00448         return false;
00449       }
00450 
00451 #ifdef lambdas
00452       // Compute least square multipliers for the given activities
00453       if (!tnlp_->eval_grad_f(n_, x, new_x, grad_f_)) {
00454         return false;
00455       }
00456       if (!SolveSystem(grad_f_, NULL, NULL, lambda_)) {
00457         return false;
00458       }
00459       IpBlasDscal(m_, -1., lambda_, 1);
00460       if (jnlst_->ProduceOutput(J_MOREVECTOR, J_NLP)) {
00461         jnlst_->Printf(J_MOREVECTOR, J_NLP,
00462                        "Curvature Estimator: least square multiplier:\n");
00463         for (Index i=0; i<m_; i++) {
00464           jnlst_->Printf(J_MOREVECTOR, J_NLP, "lambda[%5d] = %23.16e\n",
00465                          i, lambda_[i]);
00466         }
00467       }
00468       new_lambda = true;
00469 #endif
00470     }
00471 
00472     // Compute the projection of the direction
00473     if (!SolveSystem(orig_d, NULL, projected_d, NULL,
00474                      all_x_free_map_, all_g_fixed_map_,
00475                      all_comp_vec_space_, all_comp_proj_matrix_,
00476                      all_tsymlinearsolver_)) {
00477       return false;
00478     }
00479 
00480     // Sanity check to see if the thing is indeed in the null space
00481     // (if the constraint gradients are rank-deficient, the solver
00482     // might not have done a good job)
00483     Number* trash = new Number[m_];
00484     for (Index j=0; j<m_; j++) {
00485       trash[j] = 0.;
00486     }
00487     for (Index i=0; i<nnz_jac_; i++) {
00488       const Index &irow = irows_jac_[i];
00489       const Index &jcol = jcols_jac_[i];
00490       if (all_x_free_map_[jcol] >= 0 && all_g_fixed_map_[irow] >= 0) {
00491         trash[irow] += jac_vals_[i]*projected_d[jcol];
00492       }
00493     }
00494     if (jnlst_->ProduceOutput(J_MOREVECTOR, J_NLP)) {    
00495       for (Index j=0; j<m_; j++) {
00496         jnlst_->Printf(J_MOREVECTOR, J_NLP,
00497                        "nullspacevector[%5d] = %e\n", j, trash[j]);
00498       }
00499     }
00500     Index imax = IpBlasIdamax(m_, trash, 1)-1;
00501     Number max_trash = trash[imax];
00502     delete [] trash;
00503     const Number max_trash_tol = 1e-6;
00504     if (max_trash > max_trash_tol) {
00505       jnlst_->Printf(J_WARNING, J_NLP,
00506                      "Curvature Estimator: Bad solution from linear solver with max_red = %e:\n", max_trash);
00507       return false;
00508     }
00509 
00510     if (jnlst_->ProduceOutput(J_MOREVECTOR, J_NLP)) {
00511       jnlst_->Printf(J_MOREVECTOR, J_NLP,
00512                      "Curvature Estimator: original and projected directions are:\n");
00513       for (Index i=0; i<n_; i++) {
00514         jnlst_->Printf(J_MOREVECTOR, J_NLP,
00515                        "orig_d[%5d] = %23.16e proj_d[%5d] = %23.16e\n",
00516                        i, orig_d[i], i, projected_d[i]);
00517       }
00518     }
00519 
00520     gradLagTd = 0.;
00521 #ifdef lambdas
00522     // Compute the product with the gradient of the Lagrangian
00523     gradLagTd = IpBlasDdot(n, projected_d, 1, grad_f_, 1);
00524     for (Index i=0; i<nnz_jac_; i++) {
00525       const Index &irow = irows_jac_[i];
00526       const Index &jcol = jcols_jac_[i];
00527       gradLagTd += lambda_[irow]*jac_vals_[i]*projected_d[jcol];
00528     }
00529 #endif
00530 
00531     // Compute the product with the Hessian of the Lagrangian
00532     //    if (!Compute_dTHLagd(projected_d, x, new_x, lambda_, new_lambda, dTHLagd)) {
00533     if (!Compute_dTHLagd(projected_d, x, new_x, lam, new_lambda, dTHLagd)) {
00534       return false;
00535     }
00536 
00537 #if 0
00538     printf("gradLagTd = %e dTHLagd = %e\n",gradLagTd,dTHLagd);
00539 #endif
00540     return true;
00541   }
00542 
00543   bool CurvatureEstimator::PrepareNewMatrixStructure(
00544     const Number* x_l,
00545     const Number* x_u,
00546     const Number* g_l,
00547     const Number* g_u,
00548     std::vector<int>& active_x,
00549     std::vector<int>& active_g,
00550     Index& nx_free,
00551     Index* x_free_map,
00552     Index& ng_fixed,
00553     Index* g_fixed_map,
00554     SmartPtr<CompoundSymMatrixSpace>& comp_proj_matrix_space,
00555     SmartPtr<CompoundVectorSpace>& comp_vec_space)
00556   {
00557     // Deterimine which variables are free
00558     for (Index i=0; i<n_; i++) {
00559       x_free_map[i] = 0;
00560     }
00561     // fixed by activities
00562     for (std::vector<int>::iterator i=active_x.begin();
00563          i != active_x.end(); i++) {
00564       DBG_ASSERT(*i != 0 && *i<=n_ && *i>=-n_);
00565       if (*i<0) {
00566         x_free_map[(-*i)-1] = -1;
00567         DBG_ASSERT(x_l[(-*i)-1] > -1e19);
00568       }
00569       else {
00570         x_free_map[(*i)-1] = -1;
00571         DBG_ASSERT(x_u[(*i)-1] < 1e19);
00572       }
00573     }
00574     // fixed by bounds
00575     for (Index i=0; i<n_; i++) {
00576       if (x_l[i] == x_u[i]) {
00577         x_free_map[i] = -1;
00578       }
00579     }
00580     // Correct the numbering in the x map and determine number of
00581     // free variables
00582     nx_free = 0;
00583     for (Index i=0; i<n_; i++) {
00584       if (x_free_map[i] == 0) {
00585         x_free_map[i] = nx_free++;
00586       }
00587     }
00588 
00589     // Determine which constraints are fixed
00590     for (Index j=0; j<m_; j++) {
00591       g_fixed_map[j] = -1;
00592     }
00593     // fixed by activities
00594     for (std::vector<int>::iterator i=active_g.begin();
00595          i != active_g.end(); i++) {
00596       DBG_ASSERT(*i != 0 && *i<=m_ && *i>=-m_);
00597       if (*i<0) {
00598         g_fixed_map[(-*i)-1] = 0;
00599         DBG_ASSERT(g_l[(-*i)-1] > -1e19); //ToDo look at option?
00600       }
00601       else {
00602         g_fixed_map[(*i)-1] = 0;
00603         DBG_ASSERT(g_u[(*i)-1] < 1e19);
00604       }
00605     }
00606     // fixed by bounds
00607     for (Index j=0; j<m_; j++) {
00608       if (g_l[j] == g_u[j]) {
00609         g_fixed_map[j] = 0;
00610       }
00611     }
00612     // Correct the numbering in the g map and determine number of
00613     // fixed constraints
00614     ng_fixed = 0;
00615     for (Index j=0; j<m_; j++) {
00616       if (g_fixed_map[j] == 0) {
00617         g_fixed_map[j] = ng_fixed++;
00618       }
00619     }
00620 
00621     // Determine the row and column indices for the Jacobian of the fixed
00622     // constraints in the space of the free variables
00623     Index* iRows = new Index[nnz_jac_];
00624     Index* jCols = new Index[nnz_jac_];
00625     Index nnz_proj_jac = 0;
00626     for (Index i=0; i<nnz_jac_; i++) {
00627       const Index &irow = irows_jac_[i];
00628       const Index &jcol = jcols_jac_[i];
00629       if (x_free_map[jcol] >= 0 && g_fixed_map[irow] >= 0) {
00630         iRows[nnz_proj_jac] = g_fixed_map[irow] + 1;
00631         jCols[nnz_proj_jac] = x_free_map[jcol] + 1;
00632         nnz_proj_jac++;
00633       }
00634     }
00635 
00636     // Create the matrix space for the Jacobian matrices
00637     SmartPtr<GenTMatrixSpace> proj_jac_space =
00638       new GenTMatrixSpace(ng_fixed, nx_free, nnz_proj_jac, iRows, jCols);
00639     delete [] iRows;
00640     delete [] jCols;
00641 
00642     // Create Matrix space for the projection matrix
00643     comp_proj_matrix_space =
00644       new CompoundSymMatrixSpace(2, nx_free+ng_fixed);
00645     comp_proj_matrix_space->SetBlockDim(0, nx_free);
00646     comp_proj_matrix_space->SetBlockDim(1, ng_fixed);
00647     SmartPtr<SymMatrixSpace> identity_space =
00648       new IdentityMatrixSpace(nx_free);
00649     comp_proj_matrix_space->SetCompSpace(0, 0, *identity_space, true);
00650     comp_proj_matrix_space->SetCompSpace(1, 0, *proj_jac_space, true);
00651 
00652     // Create a Vector space for the rhs and sol
00653     comp_vec_space = new CompoundVectorSpace(2, nx_free+ng_fixed);
00654     SmartPtr<DenseVectorSpace> x_space = new DenseVectorSpace(nx_free);
00655     comp_vec_space->SetCompSpace(0, *x_space);
00656     SmartPtr<DenseVectorSpace> g_space = new DenseVectorSpace(ng_fixed);
00657     comp_vec_space->SetCompSpace(1, *g_space);
00658 
00659     return true;
00660   }
00661 
00662   bool CurvatureEstimator::PrepareNewMatrixValues(
00663     const Index* x_free_map,
00664     const Index* g_fixed_map,
00665     SmartPtr<CompoundSymMatrixSpace>& comp_proj_matrix_space,
00666     SmartPtr<CompoundSymMatrix>& comp_proj_matrix,
00667     SmartPtr<TSymLinearSolver>& tsymlinearsolver)
00668   {
00669     comp_proj_matrix = comp_proj_matrix_space->MakeNewCompoundSymMatrix();
00670     SmartPtr<Matrix> jac = comp_proj_matrix->GetCompNonConst(1, 0);
00671     SmartPtr<GenTMatrix> tjac = static_cast<GenTMatrix*> (GetRawPtr(jac));
00672     Number* vals = tjac->Values();
00673     Index inz=0;
00674     for (Index i=0; i<nnz_jac_; i++) {
00675       Index irow = irows_jac_[i];
00676       Index jcol = jcols_jac_[i];
00677       if (x_free_map[jcol] >= 0 && g_fixed_map[irow] >= 0) {
00678         vals[inz++] = jac_vals_[i];
00679       }
00680     }
00681     DBG_ASSERT(inz == tjac->Nonzeros());
00682 
00683     // We need to reset the linear solver object, so that it knows
00684     // that the structure of the linear system has changed
00685     tsymlinearsolver->ReducedInitialize(*jnlst_, *options_, prefix_);
00686 
00687     return true;
00688   }
00689 
00690   bool CurvatureEstimator::SolveSystem(
00691     const Number* rhs_x,
00692     const Number* rhs_g,
00693     Number* sol_x, Number* sol_g,
00694     const Index* x_free_map,
00695     const Index* g_fixed_map,
00696     SmartPtr<CompoundVectorSpace>& comp_vec_space,
00697     SmartPtr<CompoundSymMatrix>& comp_proj_matrix,
00698     SmartPtr<TSymLinearSolver>& tsymlinearsolver)
00699   {
00700     // Create a vector for the RHS
00701     SmartPtr<CompoundVector> rhs = comp_vec_space->MakeNewCompoundVector();
00702     SmartPtr<Vector> vrhs_x = rhs->GetCompNonConst(0);
00703     SmartPtr<Vector> vrhs_g = rhs->GetCompNonConst(1);
00704     // Now fill this vector with the values, extracting the relevant entries
00705     if (rhs_x) {
00706       SmartPtr<DenseVector> drhs_x =
00707         static_cast<DenseVector*> (GetRawPtr(vrhs_x));
00708       Number* xvals = drhs_x->Values();
00709       for (Index i=0; i<n_; i++) {
00710         const Index& ix = x_free_map[i];
00711         if (ix>=0) {
00712           xvals[ix] = rhs_x[i];
00713         }
00714       }
00715     }
00716     else {
00717       vrhs_x->Set(0.);
00718     }
00719     if (rhs_g) {
00720       SmartPtr<DenseVector> drhs_g =
00721         static_cast<DenseVector*> (GetRawPtr(vrhs_g));
00722       Number* gvals = drhs_g->Values();
00723       for (Index j=0; j<m_; j++) {
00724         const Index& jg = g_fixed_map[j];
00725         if (jg>=0) {
00726           gvals[jg] = rhs_g[j];
00727         }
00728       }
00729     }
00730     else {
00731       vrhs_g->Set(0.);
00732     }
00733 
00734     // Solve the linear system
00735     SmartPtr<CompoundVector> sol = comp_vec_space->MakeNewCompoundVector();
00736     ESymSolverStatus solver_retval =
00737       tsymlinearsolver->Solve(*comp_proj_matrix, *rhs, *sol, false, 0);
00738     // For now we try, if the system is reported to be singular, to
00739     // still get a solution; ToDo
00740     if (solver_retval == SYMSOLVER_SINGULAR) {
00741       jnlst_->Printf(J_DETAILED, J_NLP,
00742                      "Projection matrix reported to be singular, but we try to obtain a solution anyway.\n");
00743       solver_retval =
00744         tsymlinearsolver->Solve(*comp_proj_matrix, *rhs, *sol, false, 0);
00745     }
00746 
00747     if (solver_retval != SYMSOLVER_SUCCESS) {
00748       // DELETEME
00749       //printf("Return code from Solver is %d\n", solver_retval);
00750       return false;
00751     }
00752 
00753     // and get the solution out of it
00754     if (sol_x) {
00755       SmartPtr<Vector> vsol_x = sol->GetCompNonConst(0);
00756       SmartPtr<const DenseVector> dsol_x =
00757         static_cast<const DenseVector*> (GetRawPtr(vsol_x));
00758       const Number* xvals = dsol_x->Values();
00759       for (Index i=0; i<n_; i++) {
00760         const Index& ix = x_free_map[i];
00761         if (ix >=0) {
00762           sol_x[i] = xvals[ix];
00763         }
00764         else {
00765           sol_x[i] = 0.;
00766         }
00767       }
00768     }
00769     if (sol_g) {
00770       SmartPtr<Vector> vsol_g = sol->GetCompNonConst(1);
00771       SmartPtr<const DenseVector> dsol_g =
00772         static_cast<const DenseVector*> (GetRawPtr(vsol_g));
00773       const Number* gvals = dsol_g->Values();
00774       for (Index j=0; j<m_; j++) {
00775         const Index& ig = g_fixed_map[j];
00776         if (ig>=0) {
00777           sol_g[j] = gvals[ig];
00778         }
00779         else {
00780           sol_g[j] = 0.;
00781         }
00782       }
00783     }
00784 
00785     return true;
00786   }
00787 
00788   bool CurvatureEstimator::Compute_dTHLagd(
00789     const Number* d, const Number* x, bool new_x, const Number* lambda,
00790     bool new_lambda,  Number& dTHLagd)
00791   {
00792     if (new_x || new_lambda || !hess_vals_) {
00793       if (!hess_vals_) {
00794         hess_vals_ = new Number[nnz_hess_];
00795       }
00796       if (!tnlp_->eval_h(n_, x, new_x, 1., m_, lambda, new_lambda, nnz_hess_,
00797                          NULL, NULL, hess_vals_)) {
00798         return false;
00799       }
00800     }
00801     dTHLagd = 0.;
00802     for (Index i=0; i<nnz_hess_; i++) {
00803       Index irow = irows_hess_[i];
00804       Index jcol = jcols_hess_[i];
00805       if (irow == jcol) {
00806         dTHLagd += d[irow]*d[irow]*hess_vals_[i];
00807       }
00808       else {
00809         dTHLagd += 2.*d[irow]*d[jcol]*hess_vals_[i];
00810       }
00811     }
00812     return true;
00813   }
00814 
00815 
00816 } // namespace Bonmin

Generated on Thu Sep 22 03:05:54 2011 by  doxygen 1.4.7