Commit 22522835 authored by David Bommes's avatar David Bommes

- fixed cmake issues

- added ipopt and mumps support

git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@65 1355f012-dd97-4b2f-ae87-10fa9f823a57
parent 2fd223c9
......@@ -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 ()
......
......@@ -7,6 +7,6 @@
@COMISO_PETSC_CONFIG_FILE_SETTINGS@
@COMISO_TAO_CONFIG_FILE_SETTINGS@
@COMISO_IPOPT_CONFIG_FILE_SETTINGS@
@COMISO_MUMPS_CONFIG_FILE_SETTINGS@
......@@ -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)
......
......@@ -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)
......
......@@ -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 ()
......@@ -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 ()
This diff is collapsed.
//=============================================================================
//
// CLASS IPOPTSolver
//
//=============================================================================
#ifndef COMISO_IPOPTSOLVER_HH
#define COMISO_IPOPTSOLVER_HH
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
#if COMISO_IPOPT_AVAILABLE
//== INCLUDES =================================================================
#include <vector>
#include <cstddef>
#include <gmm/gmm.h>
#include "NProblemGmmInterface.hh"
#include "NConstraintGmmInterface.hh"
#include <IpTNLP.hpp>
#include <IpIpoptApplication.hpp>
#include <IpSolveStatistics.hpp>
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
/** \class NewtonSolver NewtonSolver.hh <ACG/.../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<NConstraintGmmInterface*>& _constraints);
protected:
double* P(std::vector<double>& _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<SVectorNP>::const_iterator SVectorNP_citer;
typedef typename gmm::linalg_traits<SVectorNP>::iterator SVectorNP_iter;
/** default constructor */
NProblemIPOPT(NProblemGmmInterface* _problem, std::vector<NConstraintGmmInterface*>& _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<NConstraintGmmInterface*>& constraints_;
int nnz_jac_g_;
int nnz_h_lag_;
// constant structure of jacobian_of_constraints and hessian_of_lagrangian
std::vector<Index> jac_g_iRow_;
std::vector<Index> jac_g_jCol_;
std::vector<Index> h_lag_iRow_;
std::vector<Index> h_lag_jCol_;
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#endif // COMISO_IPOPT_AVAILABLE
//=============================================================================
#endif // ACG_IPOPTSOLVER_HH defined
//=============================================================================
//=============================================================================
//
// CLASS NConstraintGmmInterface
//
//=============================================================================
#ifndef COMISO_LINEARCONSTRAINT_HH
#define COMISO_LINEARCONSTRAINT_HH
//== INCLUDES =================================================================
#include <gmm/gmm.h>
#include "NConstraintGmmInterface.hh"
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
/** \class NProblemGmmInterface NProblemGmmInterface.hh <ACG/.../NPRoblemGmmInterface.hh>
Brief Description.
A more elaborate description follows.
*/
class LinearConstraint : public NConstraintGmmInterface
{
public:
// ToDo: appropriate Vector/MatrixType ???
typedef gmm::wsvector<double> SVectorNP;
typedef gmm::row_matrix< SVectorNP > SMatrixNP;
// use c-arrays as vectors for gmm
typedef gmm::array1D_reference<double*> 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
//=============================================================================
......@@ -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<double>& _x, std::vector<double>& _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<RMatrix,std::vector<double> >(A_orig_,Vtemp_, b_orig_) << std::endl;
}
}
......
......@@ -85,6 +85,10 @@ public:
void project_x( const std::vector<double>& _x, std::vector<double>& _xp);
void project_x( double* _x, double* _xp);
// project dx to closest one that A(x0+dx) = b
void project_dx( const std::vector<double>& _x, std::vector<double>& _xp);
void project_dx( double* _x, double* _xp);
// transform gradient
void transform_gradient( const std::vector<double>& _g, std::vector<double>& _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<double> Vtemp_;
// hack -> store initial linear system
RMatrix A_orig_;
std::vector<double> b_orig_;
};
//=============================================================================
......
......@@ -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_));
......
//=============================================================================
//
// CLASS NConstraintGmmInterface
//
//=============================================================================
#ifndef COMISO_NCONSTRAINTGMMINTERFACE_HH
#define COMISO_NCONSTRAINTGMMINTERFACE_HH
//== INCLUDES =================================================================
#include <gmm/gmm.h>
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
/** \class NProblemGmmInterface NProblemGmmInterface.hh <ACG/.../NPRoblemGmmInterface.hh>
Brief Description.
A more elaborate description follows.
*/
class NConstraintGmmInterface
{
public:
// ToDo: appropriate Vector/MatrixType ???
typedef gmm::wsvector<double> 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) {}