From 22522835d54f97d54730b41afacc7ff9dc1d5ce8 Mon Sep 17 00:00:00 2001 From: David Bommes Date: Sun, 4 Sep 2011 13:04:27 +0000 Subject: [PATCH] - fixed cmake issues - added ipopt and mumps support git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@65 1355f012-dd97-4b2f-ae87-10fa9f823a57 --- CoMISo/CMakeLists.txt | 17 +- CoMISo/Config/config.hh.in | 2 +- .../Examples/factored_solver/CMakeLists.txt | 6 + .../Examples/quadratic_solver/CMakeLists.txt | 7 + .../small_factored_example/CMakeLists.txt | 9 +- .../small_quadratic_example/CMakeLists.txt | 9 +- CoMISo/NSolver/IPOPTSolver.cc | 401 ++++++++++++++++++ CoMISo/NSolver/IPOPTSolver.hh | 229 ++++++++++ CoMISo/NSolver/LinearConstraint.hh | 117 +++++ .../LinearConstraintHandlerElimination.cc | 51 ++- .../LinearConstraintHandlerElimination.hh | 8 + .../LinearConstraintHandlerEliminationT.cc | 16 + CoMISo/NSolver/NConstraintGmmInterface.hh | 97 +++++ CoMISo/Solver/GMM_Tools.cc | 10 +- CoMISo/cmake/FindCoMISo.cmake | 17 +- CoMISo/cmake/FindMUMPS.cmake | 26 ++ 16 files changed, 1005 insertions(+), 17 deletions(-) create mode 100644 CoMISo/NSolver/IPOPTSolver.cc create mode 100644 CoMISo/NSolver/IPOPTSolver.hh create mode 100644 CoMISo/NSolver/LinearConstraint.hh create mode 100644 CoMISo/NSolver/NConstraintGmmInterface.hh create mode 100644 CoMISo/cmake/FindMUMPS.cmake diff --git a/CoMISo/CMakeLists.txt b/CoMISo/CMakeLists.txt index 15858a0..b17d7b6 100644 --- a/CoMISo/CMakeLists.txt +++ b/CoMISo/CMakeLists.txt @@ -65,7 +65,7 @@ if (PETSC_FOUND AND MPI_FOUND) set (COMISO_PETSC_CONFIG_FILE_SETTINGS "#define COMISO_PETSC_AVAILABLE 1" ) list( APPEND COMISO_INCLUDE_DIRECTORIES ${PETSC_INCLUDE_DIRS} ) else () - message (STATUS "PETSC not found!") + message (STATUS "PETSC or dependency not found!") set (COMISO_PETSC_CONFIG_FILE_SETTINGS "#define COMISO_PETSC_AVAILABLE 0" ) endif () @@ -75,16 +75,25 @@ if (TAO_FOUND AND PETSC_FOUND AND MPI_FOUND) set (COMISO_TAO_CONFIG_FILE_SETTINGS "#define COMISO_TAO_AVAILABLE 1" ) list( APPEND COMISO_INCLUDE_DIRECTORIES ${TAO_INCLUDE_DIRS} ) else () - message (STATUS "TAO not found!") + message (STATUS "TAO or dependency not found!") set (COMISO_TAO_CONFIG_FILE_SETTINGS "#define COMISO_TAO_AVAILABLE 0" ) endif () +find_package (MUMPS) +if (MUMPS_FOUND ) + set (COMISO_MUMPS_CONFIG_FILE_SETTINGS "#define COMISO_MUMPS_AVAILABLE 1" ) + list( APPEND COMISO_INCLUDE_DIRECTORIES ${MUMPS_INCLUDE_DIR} ) +else () + message (STATUS "MUMPS not found!") + set (COMISO_MUMPS_CONFIG_FILE_SETTINGS "#define COMISO_MUMPS_AVAILABLE 0" ) +endif () + find_package (IPOPT) -if (IPOPT_FOUND ) +if (IPOPT_FOUND AND MUMPS_FOUND) set (COMISO_IPOPT_CONFIG_FILE_SETTINGS "#define COMISO_IPOPT_AVAILABLE 1" ) list( APPEND COMISO_INCLUDE_DIRECTORIES ${IPOPT_INCLUDE_DIR} ) else () - message (STATUS "IPOPT not found!") + message (STATUS "IPOPT or dependency not found!") set (COMISO_IPOPT_CONFIG_FILE_SETTINGS "#define COMISO_IPOPT_AVAILABLE 0" ) endif () diff --git a/CoMISo/Config/config.hh.in b/CoMISo/Config/config.hh.in index 0495229..fecc0a2 100644 --- a/CoMISo/Config/config.hh.in +++ b/CoMISo/Config/config.hh.in @@ -7,6 +7,6 @@ @COMISO_PETSC_CONFIG_FILE_SETTINGS@ @COMISO_TAO_CONFIG_FILE_SETTINGS@ @COMISO_IPOPT_CONFIG_FILE_SETTINGS@ - +@COMISO_MUMPS_CONFIG_FILE_SETTINGS@ diff --git a/CoMISo/Examples/factored_solver/CMakeLists.txt b/CoMISo/Examples/factored_solver/CMakeLists.txt index 613a9e2..45222b2 100644 --- a/CoMISo/Examples/factored_solver/CMakeLists.txt +++ b/CoMISo/Examples/factored_solver/CMakeLists.txt @@ -2,6 +2,7 @@ include (ACGCommon) find_package(CoMISo) find_package(SUITESPARSE) find_package(MPI) +find_package(IPOPT) include_directories ( @@ -11,12 +12,15 @@ include_directories ( ${CMAKE_CURRENT_BINARY_DIR} ${COMISO_INCLUDE_DIR} ${SUITESPARSE_INCLUDE_DIRS} + ${MUMPS_INCLUDE_DIR} + ${IPOPT_INCLUDE_DIR} ) link_directories ( ${SUITESPARSE_LIBRARY_DIRS} ${TAO_LIBRARY_DIR} ${PETSC_LIBRARY_DIR} + ${IPOPT_HSL_LIBRARY_DIR} ) # source code directories @@ -49,6 +53,8 @@ target_link_libraries (factored_solver ${MPI_LIBRARY} ${PETSC_LIBRARY} ${TAO_LIBRARY} + ${MUMPS_LIBRARY} + ${IPOPT_LIBRARY} ) if (APPLE) diff --git a/CoMISo/Examples/quadratic_solver/CMakeLists.txt b/CoMISo/Examples/quadratic_solver/CMakeLists.txt index e7a87a7..71e3229 100644 --- a/CoMISo/Examples/quadratic_solver/CMakeLists.txt +++ b/CoMISo/Examples/quadratic_solver/CMakeLists.txt @@ -2,6 +2,8 @@ include (ACGCommon) find_package(CoMISo) find_package(SUITESPARSE) find_package(MPI) +find_package(IPOPT) + include_directories ( .. @@ -10,12 +12,15 @@ include_directories ( ${CMAKE_CURRENT_BINARY_DIR} ${COMISO_INCLUDE_DIR} ${SUITESPARSE_INCLUDE_DIRS} + ${MUMPS_INCLUDE_DIR} + ${IPOPT_INCLUDE_DIR} ) link_directories ( ${SUITESPARSE_LIBRARY_DIRS} ${TAO_LIBRARY_DIR} ${PETSC_LIBRARY_DIR} + ${IPOPT_HSL_LIBRARY_DIR} ) # source code directories @@ -48,6 +53,8 @@ target_link_libraries (quadratic_solver ${MPI_LIBRARY} ${PETSC_LIBRARY} ${TAO_LIBRARY} + ${MUMPS_LIBRARY} + ${IPOPT_LIBRARY} ) if (APPLE) diff --git a/CoMISo/Examples/small_factored_example/CMakeLists.txt b/CoMISo/Examples/small_factored_example/CMakeLists.txt index 56b8405..121ffed 100644 --- a/CoMISo/Examples/small_factored_example/CMakeLists.txt +++ b/CoMISo/Examples/small_factored_example/CMakeLists.txt @@ -2,6 +2,8 @@ include (ACGCommon) find_package(CoMISo) find_package(SUITESPARSE) find_package(MPI) +find_package(IPOPT) + include_directories ( .. @@ -10,12 +12,15 @@ include_directories ( ${CMAKE_CURRENT_BINARY_DIR} ${COMISO_INCLUDE_DIR} ${SUITESPARSE_INCLUDE_DIRS} + ${MUMPS_INCLUDE_DIR} + ${IPOPT_INCLUDE_DIR} ) link_directories ( ${SUITESPARSE_LIBRARY_DIRS} ${TAO_LIBRARY_DIR} ${PETSC_LIBRARY_DIR} + ${IPOPT_HSL_LIBRARY_DIR} ) # source code directories @@ -48,6 +53,8 @@ target_link_libraries (small_factored_solver ${MPI_LIBRARY} ${PETSC_LIBRARY} ${TAO_LIBRARY} + ${MUMPS_LIBRARY} + ${IPOPT_LIBRARY} ) if (APPLE) @@ -60,5 +67,3 @@ if (APPLE) MACOSX_BUNDLE_INFO_STRING "CoMISo small_factored_solver" ) endif () - - diff --git a/CoMISo/Examples/small_quadratic_example/CMakeLists.txt b/CoMISo/Examples/small_quadratic_example/CMakeLists.txt index 7ed6aec..4fd4f80 100644 --- a/CoMISo/Examples/small_quadratic_example/CMakeLists.txt +++ b/CoMISo/Examples/small_quadratic_example/CMakeLists.txt @@ -2,6 +2,8 @@ include (ACGCommon) find_package(CoMISo) find_package(SUITESPARSE) find_package(MPI) +find_package(IPOPT) + include_directories ( .. @@ -10,12 +12,15 @@ include_directories ( ${CMAKE_CURRENT_BINARY_DIR} ${COMISO_INCLUDE_DIR} ${SUITESPARSE_INCLUDE_DIRS} + ${MUMPS_INCLUDE_DIR} + ${IPOPT_INCLUDE_DIR} ) link_directories ( ${SUITESPARSE_LIBRARY_DIRS} ${TAO_LIBRARY_DIR} ${PETSC_LIBRARY_DIR} + ${IPOPT_HSL_LIBRARY_DIR} ) # source code directories @@ -48,6 +53,8 @@ target_link_libraries (small_quadratic_solver ${MPI_LIBRARY} ${PETSC_LIBRARY} ${TAO_LIBRARY} + ${MUMPS_LIBRARY} + ${IPOPT_LIBRARY} ) if (APPLE) @@ -60,5 +67,3 @@ if (APPLE) MACOSX_BUNDLE_INFO_STRING "CoMISo small_quadratic_solver" ) endif () - - diff --git a/CoMISo/NSolver/IPOPTSolver.cc b/CoMISo/NSolver/IPOPTSolver.cc new file mode 100644 index 0000000..79cdf35 --- /dev/null +++ b/CoMISo/NSolver/IPOPTSolver.cc @@ -0,0 +1,401 @@ +//============================================================================= +// +// CLASS IPOPTSolver - IMPLEMENTATION +// +//============================================================================= + +//== INCLUDES ================================================================= + +//== COMPILE-TIME PACKAGE REQUIREMENTS ======================================== +#include +#if COMISO_IPOPT_AVAILABLE +//============================================================================= + + +#include "IPOPTSolver.hh" + +//== NAMESPACES =============================================================== + +namespace COMISO { + +//== IMPLEMENTATION ========================================================== + + + +int +IPOPTSolver:: +solve(NProblemGmmInterface* _problem, std::vector& _constraints) +{ + //---------------------------------------------------------------------------- + // 1. Create an instance of IPOPT NLP + //---------------------------------------------------------------------------- + Ipopt::SmartPtr np = new NProblemIPOPT(_problem, _constraints); + + //---------------------------------------------------------------------------- + // 2. solve problem + //---------------------------------------------------------------------------- + // Create an instance of the IpoptApplication + Ipopt::SmartPtr app = IpoptApplicationFactory(); + + app->Options()->SetStringValue("linear_solver", "ma57"); + // app->Options()->SetStringValue("derivative_test", "second-order"); + // app->Options()->SetIntegerValue("print_level", 0); + + // Initialize the IpoptApplication and process the options + Ipopt::ApplicationReturnStatus status; + status = app->Initialize(); + if (status != Ipopt::Solve_Succeeded) + { + printf("\n\n*** Error IPOPT during initialization!\n"); + } + + //---------------------------------------------------------------------------- + // 3. solve problem + //---------------------------------------------------------------------------- + status = app->OptimizeTNLP(np); + + //---------------------------------------------------------------------------- + // 4. output statistics + //---------------------------------------------------------------------------- + if (status == Ipopt::Solve_Succeeded || status == Ipopt::Solved_To_Acceptable_Level) + { + // Retrieve some statistics about the solve + Ipopt::Index iter_count = app->Statistics()->IterationCount(); + printf("\n\n*** IPOPT: The problem solved in %d iterations!\n", iter_count); + + Ipopt::Number final_obj = app->Statistics()->FinalObjective(); + printf("\n\n*** IPOPT: The final value of the objective function is %e.\n", final_obj); + } + + return status; +} + + +//== IMPLEMENTATION PROBLEM INSTANCE========================================================== + + +bool NProblemIPOPT::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, + Index& nnz_h_lag, IndexStyleEnum& index_style) +{ + // number of variables + n = problem_->n_unknowns(); + + // number of constraints + m = constraints_.size(); + + // get nonzero structure + std::vector x(n); + problem_->initial_x(&(x[0])); + // ToDo: perturb x + + // nonzeros in the jacobian of C_ and the hessian of the lagrangian + SVectorNP g; + SMatrixNP H; + problem_->eval_hessian(&(x[0]), H); + nnz_jac_g = 0; + nnz_h_lag = gmm::nnz(H); + + // clear old data + jac_g_iRow_.clear(); + jac_g_jCol_.clear(); + h_lag_iRow_.clear(); + h_lag_jCol_.clear(); + + // get non-zero structure of initial hessian + // iterate over rows + for( int i=0; i= (int)v_it.index()) + { + h_lag_iRow_.push_back(i); + h_lag_jCol_.push_back(v_it.index()); + } + } + } + + + // get nonzero structure of constraints + for( int i=0; ieval_gradient(&(x[0]),g); + constraints_[i]->eval_hessian (&(x[0]),H); + nnz_jac_g += gmm::nnz(g); + nnz_h_lag += gmm::nnz(H); + + // build table + SVectorNP_citer v_it = gmm::vect_const_begin(g); + SVectorNP_citer v_end = gmm::vect_const_end (g); + + for(; v_it != v_end; ++v_it) + { + jac_g_iRow_.push_back(i); + jac_g_jCol_.push_back(v_it.index()); + } + + for( int i=0; i= (int)v_it.index()) + { + h_lag_iRow_.push_back(i); + h_lag_jCol_.push_back(v_it.index()); + } + } + } + } + + // store for error checking... + nnz_jac_g_ = nnz_jac_g; + nnz_h_lag_ = nnz_h_lag; + + // We use the standard fortran index style for row/col entries + index_style = C_STYLE; + + return true; +} + + +//----------------------------------------------------------------------------- + + +bool NProblemIPOPT::get_bounds_info(Index n, Number* x_l, Number* x_u, + Index m, Number* g_l, Number* g_u) +{ + // first clear all variable bounds + for( int i=0; iconstraint_type()) + { + case NConstraintGmmInterface::NC_EQUAL : g_u[i] = 0.0 ; g_l[i] = 0.0 ; break; + case NConstraintGmmInterface::NC_LESS_EQUAL : g_u[i] = 0.0 ; g_l[i] = -1.0e19; break; + case NConstraintGmmInterface::NC_GREATER_EQUAL : g_u[i] = 1.0e19; g_l[i] = 0.0 ; break; + default : g_u[i] = 1.0e19; g_l[i] = -1.0e19; break; + } + } + + return true; +} + + +//----------------------------------------------------------------------------- + + +bool NProblemIPOPT::get_starting_point(Index n, bool init_x, Number* x, + bool init_z, Number* z_L, Number* z_U, + Index m, bool init_lambda, + Number* lambda) +{ + // get initial value of problem instance + problem_->initial_x(x); + + return true; +} + + +//----------------------------------------------------------------------------- + + +bool NProblemIPOPT::eval_f(Index n, const Number* x, bool new_x, Number& obj_value) +{ + // return the value of the objective function + obj_value = problem_->eval_f(x); + return true; +} + + +//----------------------------------------------------------------------------- + + +bool NProblemIPOPT::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f) +{ + problem_->eval_gradient(x, grad_f); + + return true; +} + + +//----------------------------------------------------------------------------- + + +bool NProblemIPOPT::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g) +{ + // evaluate all constraint functions + for( int i=0; ieval_constraint(x); + + return true; +} + + +//----------------------------------------------------------------------------- + + +bool NProblemIPOPT::eval_jac_g(Index n, const Number* x, bool new_x, + Index m, Index nele_jac, Index* iRow, Index *jCol, + Number* values) +{ + if (values == NULL) + { + // return the (cached) structure of the jacobian of the constraints + gmm::copy(jac_g_iRow_, VectorPTi(iRow, jac_g_iRow_.size())); + gmm::copy(jac_g_jCol_, VectorPTi(jCol, jac_g_jCol_.size())); + } + else + { + // return the values of the jacobian of the constraints + + // return the structure of the jacobian of the constraints + // global index + int gi = 0; + SVectorNP g; + + for( int i=0; ieval_gradient(x, g); + + // iterate over non-zero values + SVectorNP_citer it = gmm::vect_const_begin(g); + SVectorNP_citer ite = gmm::vect_const_end(g); + + for (; it != ite; ++it) + { + if(gi < nele_jac) + values[gi] = *it; + ++gi; + } + } + + if( gi != nele_jac) + std::cerr << "Warning: number of non-zeros in Jacobian of C is incorrect: " + << gi << " vs " << nele_jac << std::endl; + } + + return true; +} + + +//----------------------------------------------------------------------------- + + +bool NProblemIPOPT::eval_h(Index n, const Number* x, bool new_x, + Number obj_factor, Index m, const Number* lambda, + bool new_lambda, Index nele_hess, Index* iRow, + Index* jCol, Number* values) +{ + if (values == NULL) + { + // return the (cached) structure of the hessian + gmm::copy(h_lag_iRow_, VectorPTi(iRow, h_lag_iRow_.size())); + gmm::copy(h_lag_jCol_, VectorPTi(jCol, h_lag_jCol_.size())); + } + else + { + // return values. + + // global index + int gi = 0; + + // get hessian of problem + SMatrixNP H; + problem_->eval_hessian(x, H); + + for( int i=0; i= (int)v_it.index()) + { + if( gi < nele_hess) + values[gi] = obj_factor*(*v_it); + ++gi; + } + } + } + + // Hessians of Constraints + for(unsigned int j=0; jeval_hessian(x, H); + + for( int i=0; i= (int)v_it.index()) + { + if( gi < nele_hess) + values[gi] = lambda[j]*(*v_it); + ++gi; + } + } + } + } + + // error check + if( gi != nele_hess) + std::cerr << "Warning: number of non-zeros in Hessian of Lagrangian is incorrect: " + << gi << " vs " << nele_hess << std::endl; + } + return true; +} + + +//----------------------------------------------------------------------------- + + +void NProblemIPOPT::finalize_solution(SolverReturn status, + Index n, const Number* x, const Number* z_L, const Number* z_U, + Index m, const Number* g, const Number* lambda, + Number obj_value, + const IpoptData* ip_data, + IpoptCalculatedQuantities* ip_cq) +{ + // problem knows what to do + problem_->store_result(x); +} + + + +//============================================================================= +} // namespace COMISO +//============================================================================= +#endif // COMISO_IPOPT_AVAILABLE +//============================================================================= \ No newline at end of file diff --git a/CoMISo/NSolver/IPOPTSolver.hh b/CoMISo/NSolver/IPOPTSolver.hh new file mode 100644 index 0000000..1f15959 --- /dev/null +++ b/CoMISo/NSolver/IPOPTSolver.hh @@ -0,0 +1,229 @@ +//============================================================================= +// +// CLASS IPOPTSolver +// +//============================================================================= + + +#ifndef COMISO_IPOPTSOLVER_HH +#define COMISO_IPOPTSOLVER_HH + + +//== COMPILE-TIME PACKAGE REQUIREMENTS ======================================== +#include +#if COMISO_IPOPT_AVAILABLE + +//== INCLUDES ================================================================= + +#include +#include +#include +#include "NProblemGmmInterface.hh" +#include "NConstraintGmmInterface.hh" +#include +#include +#include + +//== FORWARDDECLARATIONS ====================================================== + +//== NAMESPACES =============================================================== + +namespace COMISO { + +//== CLASS DEFINITION ========================================================= + + + +/** \class NewtonSolver NewtonSolver.hh + + Brief Description. + + A more elaborate description follows. +*/ +class IPOPTSolver +{ +public: + + /// Default constructor + IPOPTSolver() {} + + /// Destructor + ~IPOPTSolver() {} + + // solve -> returns ipopt status code +//------------------------------------------------------ +// enum ApplicationReturnStatus +// { +// Solve_Succeeded=0, +// Solved_To_Acceptable_Level=1, +// Infeasible_Problem_Detected=2, +// Search_Direction_Becomes_Too_Small=3, +// Diverging_Iterates=4, +// User_Requested_Stop=5, +// Feasible_Point_Found=6, +// +// Maximum_Iterations_Exceeded=-1, +// Restoration_Failed=-2, +// Error_In_Step_Computation=-3, +// Maximum_CpuTime_Exceeded=-4, +// Not_Enough_Degrees_Of_Freedom=-10, +// Invalid_Problem_Definition=-11, +// Invalid_Option=-12, +// Invalid_Number_Detected=-13, +// +// Unrecoverable_Exception=-100, +// NonIpopt_Exception_Thrown=-101, +// Insufficient_Memory=-102, +// Internal_Error=-199 +// }; +//------------------------------------------------------ + + int solve(NProblemGmmInterface* _problem, std::vector& _constraints); + +protected: + double* P(std::vector& _v) + { + if( !_v.empty()) + return ((double*)&_v[0]); + else + return 0; + } + +private: + +}; + + +//== CLASS DEFINITION PROBLEM INSTANCE========================================================= + + +class NProblemIPOPT : public Ipopt::TNLP +{ +public: + + // Ipopt Types + typedef Ipopt::Number Number; + typedef Ipopt::Index Index; + typedef Ipopt::SolverReturn SolverReturn; + typedef Ipopt::IpoptData IpoptData; + typedef Ipopt::IpoptCalculatedQuantities IpoptCalculatedQuantities; + + typedef NConstraintGmmInterface::SVectorNP SVectorNP; + typedef NConstraintGmmInterface::SMatrixNP SMatrixNP; + + typedef gmm::array1D_reference< double* > VectorPT; + typedef gmm::array1D_reference< const double* > VectorPTC; + + typedef gmm::array1D_reference< Index* > VectorPTi; + typedef gmm::array1D_reference< const Index* > VectorPTCi; + + typedef typename gmm::linalg_traits::const_iterator SVectorNP_citer; + typedef typename gmm::linalg_traits::iterator SVectorNP_iter; + + /** default constructor */ + NProblemIPOPT(NProblemGmmInterface* _problem, std::vector& _constraints) + : problem_(_problem), constraints_(_constraints) {} + + /** default destructor */ + virtual ~NProblemIPOPT() {}; + + /**@name Overloaded from TNLP */ + //@{ + /** Method to return some info about the nlp */ + virtual bool get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, + Index& nnz_h_lag, IndexStyleEnum& index_style); + + /** Method to return the bounds for my problem */ + virtual bool get_bounds_info(Index n, Number* x_l, Number* x_u, + Index m, Number* g_l, Number* g_u); + + /** Method to return the starting point for the algorithm */ + virtual bool get_starting_point(Index n, bool init_x, Number* x, + bool init_z, Number* z_L, Number* z_U, + Index m, bool init_lambda, + Number* lambda); + + /** Method to return the objective value */ + virtual bool eval_f(Index n, const Number* x, bool new_x, Number& obj_value); + + /** Method to return the gradient of the objective */ + virtual bool eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f); + + /** Method to return the constraint residuals */ + virtual bool eval_g(Index n, const Number* x, bool new_x, Index m, Number* g); + + /** Method to return: + * 1) The structure of the jacobian (if "values" is NULL) + * 2) The values of the jacobian (if "values" is not NULL) + */ + virtual bool eval_jac_g(Index n, const Number* x, bool new_x, + Index m, Index nele_jac, Index* iRow, Index *jCol, + Number* values); + + /** Method to return: + * 1) The structure of the hessian of the lagrangian (if "values" is NULL) + * 2) The values of the hessian of the lagrangian (if "values" is not NULL) + */ + virtual bool eval_h(Index n, const Number* x, bool new_x, + Number obj_factor, Index m, const Number* lambda, + bool new_lambda, Index nele_hess, Index* iRow, + Index* jCol, Number* values); + + //@} + + /** @name Solution Methods */ + //@{ + /** This method is called when the algorithm is complete so the TNLP can store/write the solution */ + virtual void finalize_solution(SolverReturn status, + Index n, const Number* x, const Number* z_L, const Number* z_U, + Index m, const Number* g, const Number* lambda, + Number obj_value, + const IpoptData* ip_data, + IpoptCalculatedQuantities* ip_cq); + //@} + +private: + /**@name Methods to block default compiler methods. + * The compiler automatically generates the following three methods. + * Since the default compiler implementation is generally not what + * you want (for all but the most simple classes), we usually + * put the declarations of these methods in the private section + * and never implement them. This prevents the compiler from + * implementing an incorrect "default" behavior without us + * knowing. (See Scott Meyers book, "Effective C++") + * + */ + //@{ + // MyNLP(); + NProblemIPOPT(const NProblemIPOPT&); + NProblemIPOPT& operator=(const NProblemIPOPT&); + //@} + + +private: + + // pointer to problem instance + NProblemGmmInterface* problem_; + // reference to constraints vector + std::vector& constraints_; + + int nnz_jac_g_; + int nnz_h_lag_; + + // constant structure of jacobian_of_constraints and hessian_of_lagrangian + std::vector jac_g_iRow_; + std::vector jac_g_jCol_; + std::vector h_lag_iRow_; + std::vector h_lag_jCol_; +}; + + +//============================================================================= +} // namespace COMISO + +//============================================================================= +#endif // COMISO_IPOPT_AVAILABLE +//============================================================================= +#endif // ACG_IPOPTSOLVER_HH defined +//============================================================================= + diff --git a/CoMISo/NSolver/LinearConstraint.hh b/CoMISo/NSolver/LinearConstraint.hh new file mode 100644 index 0000000..9cfdc69 --- /dev/null +++ b/CoMISo/NSolver/LinearConstraint.hh @@ -0,0 +1,117 @@ +//============================================================================= +// +// CLASS NConstraintGmmInterface +// +//============================================================================= + + +#ifndef COMISO_LINEARCONSTRAINT_HH +#define COMISO_LINEARCONSTRAINT_HH + + +//== INCLUDES ================================================================= + +#include +#include "NConstraintGmmInterface.hh" + +//== FORWARDDECLARATIONS ====================================================== + +//== NAMESPACES =============================================================== + +namespace COMISO { + +//== CLASS DEFINITION ========================================================= + + + +/** \class NProblemGmmInterface NProblemGmmInterface.hh + + Brief Description. + + A more elaborate description follows. +*/ +class LinearConstraint : public NConstraintGmmInterface +{ +public: + + // ToDo: appropriate Vector/MatrixType ??? + typedef gmm::wsvector SVectorNP; + typedef gmm::row_matrix< SVectorNP > SMatrixNP; + + // use c-arrays as vectors for gmm + typedef gmm::array1D_reference VectorPT; + + // different types of constraints +// enum ConstraintType {NC_EQUAL, NC_LESS_EQUAL, NC_GREATER_EQUAL}; + + /// Default constructor + LinearConstraint(const ConstraintType _type = NC_EQUAL) : NConstraintGmmInterface(_type) + {} + + // linear equation of the form -> coeffs_^T * (x,1) =_type= 0 + LinearConstraint(const SVectorNP& _coeffs_plus_const, const ConstraintType _type = NC_EQUAL) : NConstraintGmmInterface(_type) + { + int n = gmm::vect_size(_coeffs_plus_const)-1; + + gmm::resize(coeffs_, n+1); + gmm::copy(_coeffs_plus_const, coeffs_); + + if( n >= 0) + b_ = coeffs_[n]; + else + b_ = 0.0; + + gmm::resize(coeffs_, n); + } + + // linear equation of the form -> coeffs_^T * (x,1) =_type= 0 + LinearConstraint(const SVectorNP& _coeffs, const double _b, const ConstraintType _type = NC_EQUAL) : NConstraintGmmInterface(_type) + { + int n = gmm::vect_size(_coeffs); + gmm::resize(coeffs_, n); + gmm::copy(_coeffs, coeffs_); + b_ = _b; + } + + /// Destructor + ~LinearConstraint() {} + + virtual int n_unknowns() + { + return gmm::vect_size(coeffs_); + } + + virtual double eval_constraint ( const double* _x ) + { + return (gmm::vect_sp(coeffs_, VectorPT((double*)_x, gmm::vect_size(coeffs_))) + b_); + } + + virtual void eval_gradient( const double* _x, SVectorNP& _g ) + { + gmm::resize(_g, gmm::vect_size(coeffs_)); + gmm::copy (coeffs_, _g); + } + + virtual void eval_hessian ( const double* _x, SMatrixNP& _h ) + { + gmm::resize(_h, gmm::vect_size(coeffs_), gmm::vect_size(coeffs_)); + gmm::clear(_h); + } + + // inherited from base +// virtual ConstraintType constraint_type ( ) { return type_; } + +private: + + // linear equation of the form -> coeffs_^T * x + b_ + SVectorNP coeffs_; + double b_; +}; + + +//============================================================================= +} // namespace COMISO +//============================================================================= +#endif // ACG_LINEARCONSTRAINT_HH defined +//============================================================================= + diff --git a/CoMISo/NSolver/LinearConstraintHandlerElimination.cc b/CoMISo/NSolver/LinearConstraintHandlerElimination.cc index a34d85a..0594103 100644 --- a/CoMISo/NSolver/LinearConstraintHandlerElimination.cc +++ b/CoMISo/NSolver/LinearConstraintHandlerElimination.cc @@ -152,8 +152,55 @@ project_x( double* _x, double* _xp) { Vtemp_.resize( n_red_); - transform_x (_x , &Vtemp_[0]); - inv_transform_x( &Vtemp_[0], _xp ); + transform_x (_x , &(Vtemp_[0])); + inv_transform_x( &(Vtemp_[0]), _xp ); + } +} + + +//----------------------------------------------------------------------------- + + +void +LinearConstraintHandlerElimination:: +project_dx( const std::vector& _x, std::vector& _xp) +{ + _xp.resize(n_); + if( _x.size()) + project_dx((double*)&(_x[0]), &(_xp[0])); +} + + +//----------------------------------------------------------------------------- + + +void +LinearConstraintHandlerElimination:: +project_dx( double* _x, double* _xp) +{ + if( n_ == n_red_) // special case of no constraints + { + // just copy + gmm::copy(VectorPT(_x, n_), VectorPT(_xp, n_)); + } + else + { + // idea: C_ is an orthonormal basis of the kernel + // -> simply project out the non-constraint part + + Vtemp_.resize( n_red_); + + gmm::mult(Ct_, VectorPT(_x, n_), Vtemp_); + gmm::mult(C_, Vtemp_, VectorPT(_xp, n_)); + +// // check result +// std::cerr << "check result..." << std::endl; +// Vtemp_.resize(b_orig_.size()); +// gmm::mult(A_orig_, VectorPT(_xp, n_), Vtemp_); +//// gmm::add(gmm::scaled(b_orig_, -1.0), Vtemp_); +// std::cerr << "check constraint update... " << gmm::vect_norm2(Vtemp_) << std::endl; +// +//// std::cerr << "check constraint update... " << COMISO_GMM::residuum_norm >(A_orig_,Vtemp_, b_orig_) << std::endl; } } diff --git a/CoMISo/NSolver/LinearConstraintHandlerElimination.hh b/CoMISo/NSolver/LinearConstraintHandlerElimination.hh index 0f65f22..f76ced8 100644 --- a/CoMISo/NSolver/LinearConstraintHandlerElimination.hh +++ b/CoMISo/NSolver/LinearConstraintHandlerElimination.hh @@ -85,6 +85,10 @@ public: void project_x( const std::vector& _x, std::vector& _xp); void project_x( double* _x, double* _xp); + // project dx to closest one that A(x0+dx) = b + void project_dx( const std::vector& _x, std::vector& _xp); + void project_dx( double* _x, double* _xp); + // transform gradient void transform_gradient( const std::vector& _g, std::vector& _gC); void transform_gradient( double* _g, double* _gC); @@ -110,6 +114,10 @@ private: // temp matrix to transform hessian and temp vectors RMatrix Mtemp_; std::vector Vtemp_; + + // hack -> store initial linear system + RMatrix A_orig_; + std::vector b_orig_; }; //============================================================================= diff --git a/CoMISo/NSolver/LinearConstraintHandlerEliminationT.cc b/CoMISo/NSolver/LinearConstraintHandlerEliminationT.cc index 55c39cf..335fdf0 100644 --- a/CoMISo/NSolver/LinearConstraintHandlerEliminationT.cc +++ b/CoMISo/NSolver/LinearConstraintHandlerEliminationT.cc @@ -59,6 +59,22 @@ initialize( const MatrixT& _C, const VectorT& _c) n_ = gmm::mat_nrows(C_); n_red_ = gmm::mat_ncols(C_); m_ = n_ - n_red_; + + // hack -> store initial system + if(0) + { + gmm::resize(A_orig_, gmm::mat_nrows(_C), gmm::mat_ncols(_C)); + gmm::resize(b_orig_, _c.size()); + gmm::copy(_C, A_orig_); + gmm::copy(_c, b_orig_); + } + +// RMatrix CtC(n_red_, n_red_); +// gmm::mult(Ct_,C_, CtC); +// std::cerr << "CtC\n"; +// std::cerr << CtC << std::endl; + + /* // set up least squares problem gmm::resize(Mtemp_, gmm::mat_ncols(C_), gmm::mat_ncols(C_)); diff --git a/CoMISo/NSolver/NConstraintGmmInterface.hh b/CoMISo/NSolver/NConstraintGmmInterface.hh new file mode 100644 index 0000000..14ca21c --- /dev/null +++ b/CoMISo/NSolver/NConstraintGmmInterface.hh @@ -0,0 +1,97 @@ +//============================================================================= +// +// CLASS NConstraintGmmInterface +// +//============================================================================= + + +#ifndef COMISO_NCONSTRAINTGMMINTERFACE_HH +#define COMISO_NCONSTRAINTGMMINTERFACE_HH + + +//== INCLUDES ================================================================= + +#include + +//== FORWARDDECLARATIONS ====================================================== + +//== NAMESPACES =============================================================== + +namespace COMISO { + +//== CLASS DEFINITION ========================================================= + + + +/** \class NProblemGmmInterface NProblemGmmInterface.hh + + Brief Description. + + A more elaborate description follows. +*/ +class NConstraintGmmInterface +{ +public: + + // ToDo: appropriate Vector/MatrixType ??? + typedef gmm::wsvector SVectorNP; + typedef gmm::row_matrix< SVectorNP > SMatrixNP; + + // different types of constraints + enum ConstraintType {NC_EQUAL, NC_LESS_EQUAL, NC_GREATER_EQUAL}; + + /// Default constructor + NConstraintGmmInterface(const ConstraintType _type = NC_EQUAL) : type_(_type) {} + + /// Destructor + ~NConstraintGmmInterface() {} + + virtual int n_unknowns ( ) = 0; + virtual double eval_constraint ( const double* _x ) = 0; + virtual void eval_gradient ( const double* _x, SVectorNP& _g ) = 0; + virtual void eval_hessian ( const double* _x, SMatrixNP& _h ) = 0; + + virtual ConstraintType constraint_type ( ) { return type_; } + + virtual bool is_satisfied ( const double* _x, double _eps = 1e-6 ) + { + switch( type_) + { + case NC_EQUAL : return (fabs(eval_constraint(_x)) <= _eps); break; + case NC_LESS_EQUAL : return ( eval_constraint(_x) <= _eps); break; + case NC_GREATER_EQUAL: return ( eval_constraint(_x) >= -_eps); break; + } + return false; + } + + virtual double gradient_update_factor( const double* _x, double _eps = 1e-6 ) + { + double val = eval_constraint(_x); + bool upper_bound_ok = ( val <= _eps); + bool lower_bound_ok = ( val >= -_eps); + + if(upper_bound_ok) + { + if(lower_bound_ok || type_ == NC_LESS_EQUAL) return 0.0; + else return 1.0; + } + else + { + if(lower_bound_ok && type_ == NC_GREATER_EQUAL) return 0.0; + else return -1.0; + } + } + + +private: + // constraint type + ConstraintType type_; +}; + + +//============================================================================= +} // namespace COMISO +//============================================================================= +#endif // ACG_NCONSTRAINTGMMINTERFACE_HH defined +//============================================================================= + diff --git a/CoMISo/Solver/GMM_Tools.cc b/CoMISo/Solver/GMM_Tools.cc index ba2df15..99ee375 100644 --- a/CoMISo/Solver/GMM_Tools.cc +++ b/CoMISo/Solver/GMM_Tools.cc @@ -963,13 +963,13 @@ template double residuum_norm( MatrixT& _A, VectorT& _x, VectorT& _rhs ) { if ( gmm::mat_ncols( _A ) != _x.size() ) - std::cerr << "DIM ERROR (residuum_norm): " << gmm::mat_ncols( _A ) << " " << _x.size() << std::endl; - if ( _rhs.size() != _x.size() ) - std::cerr << "DIM ERROR 2 (residuum_norm): " << _rhs.size() << " " << _x.size() << std::endl; + std::cerr << "DIM ERROR (residuum_norm): " << gmm::mat_ncols( _A ) << " vs " << _x.size() << std::endl; + if ( gmm::mat_nrows( _A) !=_rhs.size() ) + std::cerr << "DIM ERROR 2 (residuum_norm): " << gmm::mat_nrows( _A) << " vs " << _rhs.size() << std::endl; // temp vectors - VectorT Ax( _x.size() ); - VectorT res( _x.size() ); + VectorT Ax( _rhs.size()); + VectorT res( _rhs.size() ); gmm::mult( _A,_x, Ax ); gmm::add( Ax, gmm::scaled( _rhs, -1.0 ), res ); diff --git a/CoMISo/cmake/FindCoMISo.cmake b/CoMISo/cmake/FindCoMISo.cmake index 187c567..ce60015 100644 --- a/CoMISo/cmake/FindCoMISo.cmake +++ b/CoMISo/cmake/FindCoMISo.cmake @@ -14,7 +14,7 @@ ENDIF (COMISO_INCLUDE_DIR) # Find CoMISo config file FIND_PATH( COMISO_INCLUDE_DIR CoMISo/Config/config.hh - PATHS "${CMAKE_SOURCE_DIR}/libs/" ) + PATHS "${CMAKE_SOURCE_DIR}/../" "${CMAKE_SOURCE_DIR}/libs/" ) FILE(READ ${COMISO_INCLUDE_DIR}/CoMISo/Config/config.hh CURRENT_COMISO_CONFIG) @@ -63,6 +63,21 @@ if ( COMISO_IPOPT_BUILD_TIME_AVAILABLE ) endif() +STRING(REGEX MATCH "\#define COMISO_MUMPS_AVAILABLE 1" COMISO_MUMPS_BUILD_TIME_AVAILABLE ${CURRENT_COMISO_CONFIG} ) + +if ( COMISO_MUMPS_BUILD_TIME_AVAILABLE ) + + find_package(MUMPS) + + if ( NOT MUMPS_FOUND ) + message(ERROR "COMISO configured with mumps but mumps not available") + endif() + + list (APPEND COMISO_OPT_DEPS "MUMPS") + +endif() + + STRING(REGEX MATCH "\#define COMISO_TAO_AVAILABLE 1" COMISO_TAO_BUILD_TIME_AVAILABLE ${CURRENT_COMISO_CONFIG} ) if ( COMISO_TAO_BUILD_TIME_AVAILABLE ) diff --git a/CoMISo/cmake/FindMUMPS.cmake b/CoMISo/cmake/FindMUMPS.cmake new file mode 100644 index 0000000..58db0c4 --- /dev/null +++ b/CoMISo/cmake/FindMUMPS.cmake @@ -0,0 +1,26 @@ +if (MUMPS_INCLUDE_DIR) + # in cache already + SET(MUMPS_FIND_QUIETLY TRUE) +endif (MUMPS_INCLUDE_DIR) + +if (WIN32) +#TODO +ELSEIF(APPLE) +#TODO +ELSE( WIN32 ) + find_path(MUMPS_INCLUDE_DIR NAMES dmumps_c.h + PATHS "/usr/include/" + ) + + IF(MUMPS_INCLUDE_DIR) + SET(MUMPS_FOUND TRUE) + SET(MUMPS_INCLUDE_DIR ${MUMPS_INCLUDE_DIR}) + ELSE(MUMPS_INCLUDE_DIR) + SET(MUMPS_FOUND FALSE) + SET(MUMPS_INCLUDE_DIR ${IPOPT_INCLUDE_DIR}) + ENDIF(MUMPS_INCLUDE_DIR) + + find_library( MUMPS_LIBRARY + dmumps + PATHS "/usr/lib" ) +ENDIF() -- GitLab