Commit 25d7dace authored by David Bommes's avatar David Bommes

-added arpack support

git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@123 1355f012-dd97-4b2f-ae87-10fa9f823a57
parent 6d2a6077
......@@ -162,6 +162,17 @@ else ()
set (COMISO_GUROBI_CONFIG_FILE_SETTINGS "#define COMISO_GUROBI_AVAILABLE 0" )
endif ()
find_package (ARPACK)
if (ARPACK_FOUND )
set (COMISO_ARPACK_CONFIG_FILE_SETTINGS "#define COMISO_ARPACK_AVAILABLE 1" )
list( APPEND COMISO_INCLUDE_DIRECTORIES ${ARPACK_INCLUDE_DIR} )
# list( APPEND COMISO_LINK_DIRECTORIES ${ARPACK_LIBRARY_DIR} )
list( APPEND COMISO_LINK_LIBRARIES ${ARPACK_LIBRARY} )
else ()
message (STATUS "ARPACK not found!")
set (COMISO_ARPACK_CONFIG_FILE_SETTINGS "#define COMISO_ARPACK_AVAILABLE 0" )
endif ()
include_directories (
..
${CMAKE_CURRENT_SOURCE_DIR}
......@@ -184,6 +195,7 @@ set (directories
.
Solver
NSolver
EigenSolver
Config
Utils
QtWidgets
......@@ -274,5 +286,9 @@ endif()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/super_sparse_matrix/CMakeLists.txt" )
add_subdirectory (Examples/super_sparse_matrix)
endif()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/eigen_solver/CMakeLists.txt" )
add_subdirectory (Examples/eigen_solver)
endif()
......@@ -11,5 +11,6 @@
@COMISO_MUMPS_CONFIG_FILE_SETTINGS@
@COMISO_TAUCS_CONFIG_FILE_SETTINGS@
@COMISO_GUROBI_CONFIG_FILE_SETTINGS@
@COMISO_ARPACK_CONFIG_FILE_SETTINGS@
//=============================================================================
//
// CLASS ArpackSolver - IMPLEMENTATION
//
//=============================================================================
#define COMISO_ARPACKSOLVER_C
//== INCLUDES =================================================================
#include "ArpackSolver.hh"
//== NAMESPACES ===============================================================
namespace COMISO {
//== IMPLEMENTATION ==========================================================
template<class MatrixT,class MatrixT2>
void
ArpackSolver::
solve(const MatrixT& _A,
std::vector<double>& _eigenvalues,
MatrixT2& _eigenvectors,
const int _n_eigvalues,
const char* _which_eigs )
{
#if COMISO_ARPACK_AVAILABLE
Matrix A(_A);
ARSymStdEig<double, Matrix> eig_prob(A.matrix().cols(), _n_eigvalues, &A, &Matrix::mult_Mv, (char*)_which_eigs,
0, 0.0, 2000);
int n_converged = eig_prob.FindEigenvectors();
// store result
_eigenvalues.resize(n_converged);
_eigenvectors.resize(A.matrix().rows(),n_converged);
for( int i=0; i<n_converged; ++i)
{
_eigenvalues[i] = eig_prob.Eigenvalue(i);
for(int j = 0; j<A.matrix().rows(); ++j)
_eigenvectors.coeffRef(j,i) = eig_prob.RawEigenvector(i)[j];
}
#endif
}
//-----------------------------------------------------------------------------
template<class MatrixT,class MatrixT2>
void
ArpackSolver::
solve_inverse(const MatrixT& _A,
std::vector<double>& _eigenvalues,
MatrixT2& _eigenvectors,
const int _n_eigvalues,
const char* _which_eigs)
{
#if COMISO_ARPACK_AVAILABLE
Matrix A(_A,true);
ARSymStdEig<double, Matrix> eig_prob(A.matrix().cols(), _n_eigvalues, &A, &Matrix::mult_M_inv_v, (char*)_which_eigs,
0, 0.0, 2000);
// ARSymStdEig(int np, int nevp, ARFOP* objOPp,
// void (ARFOP::* MultOPxp)(ARFLOAT[], ARFLOAT[]),
// char* whichp = "LM", int ncvp = 0, ARFLOAT tolp = 0.0,
// int maxitp = 0, ARFLOAT* residp = NULL, bool ishiftp = true);
int n_converged = eig_prob.FindEigenvectors();
// store result
_eigenvalues.resize(n_converged);
_eigenvectors.resize(A.matrix().rows(),n_converged);
for( int i=0; i<n_converged; ++i)
{
_eigenvalues[i] = eig_prob.Eigenvalue(i);
for(int j = 0; j<A.matrix().rows(); ++j)
_eigenvectors.coeffRef(j,i) = eig_prob.RawEigenvector(i)[j];
}
#endif
}
//-----------------------------------------------------------------------------
template<class MatrixT,class MatrixT2>
void
ArpackSolver::
check_result(const MatrixT& _A, std::vector<double>& _eigenvalues, MatrixT2& _eigenvectors)
{
int n=_eigenvectors.rows();
if(n<20)
std::cerr << _A << std::endl;
for(unsigned int i=0; i<_eigenvalues.size(); ++i)
{
std::cerr << "eigenvalue " << i << ": " << _eigenvalues[i] << ", ";
if(n < 20)
{
std::cerr << "eigenvector: ";
for(int j=0; j<n; ++j)
{
std::cerr << _eigenvectors.coeffRef(j,i) << ", ";
}
}
// compute residuum
Eigen::Matrix<double, Eigen::Dynamic, 1> v = _eigenvectors.block(0,i,n,1);
std::cerr << "residuum norm: " << (_A*v - _eigenvalues[i]*v).norm() << std::endl;
}
}
//-----------------------------------------------------------------------------
//=============================================================================
} // namespace COMISO
//=============================================================================
//=============================================================================
//
// CLASS ArpackSolver
//
//=============================================================================
#ifndef COMISO_ARPACKSOLVER_HH
#define COMISO_ARPACKSOLVER_HH
//== COMPILE-TIME PACKAGE REQUIREMENTS ========================================
#include <CoMISo/Config/config.hh>
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include "EigenArpackMatrixT.hh"
#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
#include <Eigen/Sparse>
#if COMISO_ARPACK_AVAILABLE
#include <arpack++/arssym.h>
#endif
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
/** \class ArpackSolver ArpackSolver.hh <COMISO/.../ArpackSolver.hh>
Brief Description.
A more elaborate description follows.
*/
class COMISODLLEXPORT ArpackSolver
{
public:
// sparse matrix type
typedef EigenArpackMatrixT<double,Eigen::SparseMatrix<double,Eigen::ColMajor> > Matrix;
/// Constructor
ArpackSolver() {}
/// Destructor
~ArpackSolver() {}
// solve eigenproblem
// number of desired eigenvalues -> _n_eigenvalues
// which eigenvalues -> one of {LA (largest algebraic), SA (smalles algebraic), LM (largest magnitude), SM(smallest magnitued), BE(both ends)}
template<class MatrixT,class MatrixT2>
void solve(const MatrixT& _A,
std::vector<double>& _eigenvalues,
MatrixT2& _eigenvectors,
const int _n_eigvalues = 1,
const char* _which_eigs = "SM");
// solve eigenproblem
// number of desired eigenvalues -> _n_eigenvalues
// which eigenvalues -> one of {LA (largest algebraic), SA (smalles algebraic), LM (largest magnitude), SM(smallest magnitued), BE(both ends)}
template<class MatrixT,class MatrixT2>
void solve_inverse(const MatrixT& _A,
std::vector<double>& _eigenvalues,
MatrixT2& _eigenvectors,
const int _n_eigvalues = 1,
const char* _which_eigs = "LM");
// check resulting eigenvalues/eigenvectors
template<class MatrixT,class MatrixT2>
void check_result(const MatrixT& _A, std::vector<double>& _eigenvalues, MatrixT2& _eigenvectors);
private:
};
//=============================================================================
} // namespace ACG
//=============================================================================
#if defined(INCLUDE_TEMPLATES) && !defined(COMISO_ARPACKSOLVER_C)
#define COMISO_ARPACKSOLVER_TEMPLATES
#include "ArpackSolver.cc"
#endif
//=============================================================================
#endif // ACG_ARPACKSOLVER_HH defined
//=============================================================================
//=============================================================================
//
// CLASS EigenArpackMatrixT - IMPLEMENTATION
//
//=============================================================================
#define COMISO_EIGENARPACKMATRIXT_C
//== INCLUDES =================================================================
#include "EigenArpackMatrixT.hh"
//== NAMESPACES ===============================================================
namespace COMISO {
//== IMPLEMENTATION ==========================================================
//-----------------------------------------------------------------------------
//=============================================================================
} // namespace ACG
//=============================================================================
//=============================================================================
//
// CLASS EigenArpackMatrixT
//
//=============================================================================
#ifndef COMISO_EIGENARPACKMATRIXT_HH
#define COMISO_EIGENARPACKMATRIXT_HH
//== INCLUDES =================================================================
#include <iostream>
#include <Eigen/Eigen>
#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
#include <Eigen/Sparse>
#include <unsupported/Eigen/CholmodSupport>
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
/** \class EigenArpackMatrixT EigenArpackMatrixT.hh <COMISO/.../EigenArpackMatrixT.hh>
Brief Description.
A more elaborate description follows.
*/
template <class RealT,class MatrixT>
class EigenArpackMatrixT
{
public:
typedef MatrixT Matrix;
typedef RealT Real;
/// Constructor
template<class MatrixT2>
EigenArpackMatrixT(const MatrixT2& _m, bool _use_inverse = false)
{
mat_ = _m;
if(_use_inverse)
{
sllt_.compute(mat_);
if ( !sllt_.succeeded() )
std::cout << "[ERROR] EigenArpackMatrix(): Could not compute llt factorization." << std::endl;
}
}
/// Destructor
~EigenArpackMatrixT() {}
// get reference on matrix
Matrix& matrix() { return mat_; }
// matrix-vector multiplication _w = mat_*_v
void mult_Mv(Real* _v, Real* _w)
{
Eigen::Map<Eigen::Matrix<Real, Eigen::Dynamic, 1> > v(_v,mat_.rows()); // uses v as a ArrayXf object
Eigen::Map<Eigen::Matrix<Real, Eigen::Dynamic, 1> > w(_w,mat_.cols()); // uses w as a ArrayXf object
w = mat_*v;
}
// matrix-vector multiplication _w = mat_*_v
void mult_M_inv_v(Real* _v, Real* _w)
{
Eigen::Map<Eigen::Matrix<Real, Eigen::Dynamic, 1> > v(_v,mat_.rows()); // uses v as a ArrayXf object
Eigen::Map<Eigen::Matrix<Real, Eigen::Dynamic, 1> > w(_w,mat_.cols()); // uses w as a ArrayXf object
w = sllt_.solve(v);
// std::cerr << "input:" << std::endl;
// std::cerr << v << std::endl;
// std::cerr << "output:" << std::endl;
// std::cerr << w << std::endl;
}
private:
Matrix mat_;
Eigen::SparseLLT<Eigen::SparseMatrix<Real>, Eigen::Cholmod> sllt_;
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#if defined(INCLUDE_TEMPLATES) && !defined(COMISO_EIGENARPACKMATRIXT_C)
#define COMISO_EIGENARPACKMATRIXT_TEMPLATES
#include "EigenArpackMatrixT.cc"
#endif
//=============================================================================
#endif // COMISO_EIGENARPACKMATRIXT_HH defined
//=============================================================================
......@@ -434,7 +434,7 @@ bool NProblemIPOPT::eval_h(Index n, const Number* x, bool new_x,
for(unsigned int j=0; j<constraints_.size(); ++j)
{
SMatrixNC H;
constraints_[j]->eval_hessian(x, H);
constraints_[j]->eval_hessian(&(x_rnd[0]), H);
SMatrixNC::iterator m_it = H.begin();
SMatrixNC::iterator m_end = H.end();
......
......@@ -116,11 +116,10 @@ public:
else ++ n_ok;
}
}
if (n_errors > 0) {
std::cerr << "############## Gradient Checker #############\n";
std::cerr << "#ok : " << n_ok << std::endl;
std::cerr << "#error: " << n_errors << std::endl;
}
std::cerr << "############## Gradient Checker #############\n";
std::cerr << "#ok : " << n_ok << std::endl;
std::cerr << "#error: " << n_errors << std::endl;
return (n_errors == 0);
}
......@@ -174,11 +173,9 @@ public:
}
}
if (n_errors > 0) {
std::cerr << "############## Hessian Checker #############\n";
std::cerr << "#ok : " << n_ok << std::endl;
std::cerr << "#error: " << n_errors << std::endl;
}
std::cerr << "############## Hessian Checker #############\n";
std::cerr << "#ok : " << n_ok << std::endl;
std::cerr << "#error: " << n_errors << std::endl;
return (n_errors == 0);
}
......
......@@ -27,11 +27,26 @@ if( GUROBI_FOUND)
list( APPEND COMISO_LINK_LIBRARIES ${GUROBI_LIBRARY} )
endif()
find_package(ARPACK)
if( ARPACK_FOUND)
list( APPEND COMISO_INCLUDE_DIRECTORIES ${ARPACK_INCLUDE_DIR} )
list( APPEND COMISO_LINK_LIBRARIES ${ARPACK_LIBRARIES} )
endif()
FIND_PACKAGE( Boost 1.42.0 COMPONENTS system filesystem regex)
if( Boost_FOUND)
list( APPEND COMISO_INCLUDE_DIRECTORIES ${Boost_INCLUDE_DIR} )
list( APPEND COMISO_LINK_DIRECTORIES ${Boost_LIBRARY_DIRS} )
list( APPEND COMISO_LINK_LIBRARIES ${Boost_LIBRARIES} )
endif()
list( APPEND COMISO_LINK_LIBRARIES ${SUITESPARSE_LIBRARIES} )
list( APPEND COMISO_LINK_LIBRARIES ${MPI_CXX_LIBRARIES} )
list( APPEND COMISO_LINK_LIBRARIES ${PETSC_LIBRARY} )
list( APPEND COMISO_LINK_LIBRARIES ${TAO_LIBRARY} )
#MESSAGE( ${COMISO_LINK_LIBRARIES})
include_directories (
..
${CMAKE_SOURCE_DIR}
......
# - Try to find ARPACK
# Once done this will define
#
# ARPACK_FOUND - system has ARPACK
# ARPACK_INCLUDE_DIR - the ARPACK include directories
# ARPACK_LIBRARY - Link this to use ARPACK
# ARPACKpp_LIBRARY - Link this to use ARPACK++
# ARPACK_LIBRARIES - Link this to use ARPACK together with ARPACK++
find_path (ARPACK_INCLUDE_DIR NAMES arpack++/argeig.h
HINTS ENV ARPACK_INCLUDE_DIR
PATHS /usr/include/arpack++ "C:\\libs\\arpack++\\include"
DOC "ARPACK Include Directory")
IF ( WIN32 )
find_library( ARPACK_LIBRARY arpack.lib
PATHS "C:\\libs\\arpack++\\lib"
)
ELSE( WIN32 )
find_library( ARPACK_LIBRARY arpack arpack++
PATHS /usr/lib
)
find_library( ARPACKpp_LIBRARY arpack++
PATHS /usr/lib
)
list( APPEND ARPACK_LIBRARIES ${ARPACK_LIBRARY} )
list( APPEND ARPACK_LIBRARIES ${ARPACKpp_LIBRARY} )
ENDIF( WIN32 )
IF (ARPACK_INCLUDE_DIR AND ARPACK_LIBRARY)
SET(ARPACK_FOUND TRUE)
ELSE ()
SET(ARPACK_FOUND FALSE)
ENDIF ()
......@@ -121,6 +121,20 @@ if ( COMISO_INCLUDE_DIR )
endif()
STRING(REGEX MATCH "\#define COMISO_ARPACK_AVAILABLE 1" COMISO_ARPACK_BUILD_TIME_AVAILABLE ${CURRENT_COMISO_CONFIG} )
if ( COMISO_ARPACK_BUILD_TIME_AVAILABLE )
find_package(ARPACK)
if ( NOT ARPACK_FOUND )
message(ERROR "COMISO configured with ARPACK but ARPACK not available")
endif()
list (APPEND COMISO_OPT_DEPS "ARPACK")
endif()
add_definitions (-DCOMISODLL -DUSECOMISO )
endif(COMISO_INCLUDE_DIR)
......
......@@ -39,21 +39,33 @@ ELSE( WIN32 )
# MESSAGE(STATUS "$ENV{GUROBI_HOME}/include")
IF(GUROBI_INCLUDE_DIR)
SET(GUROBI_FOUND TRUE)
SET(GUROBI_INCLUDE_DIR ${GUROBI_INCLUDE_DIR})
SET(GUROBI_LIBRARY_DIR "$ENV{GUROBI_HOME}/lib/" CACHE PATH "Path to GUROBI Library")
SET(GUROBI_LIBRARY "gurobi45;gurobi_c++;pthread" CACHE STRING "GUROBI Libraries")
MESSAGE(STATUS "${GUROBI_LIBRARY_DIR}")
MESSAGE(STATUS "${GUROBI_LIBRARY}")
ELSE(GUROBI_INCLUDE_DIR)
SET(GUROBI_FOUND FALSE)
SET(GUROBI_INCLUDE_DIR ${GUROBI_INCLUDE_DIR})
ENDIF(GUROBI_INCLUDE_DIR)
SET(GUROBI_FOUND TRUE)
SET(GUROBI_INCLUDE_DIR ${GUROBI_INCLUDE_DIR})
ELSE(GUROBI_INCLUDE_DIR)
SET(GUROBI_FOUND FALSE)
SET(GUROBI_INCLUDE_DIR ${GUROBI_INCLUDE_DIR})
ENDIF(GUROBI_INCLUDE_DIR)
#find_library( GUROBI_LIBRARY
# gurobi
# PATHS "${GUROBI_HOME}/lib" )
ENDIF()
endif(GUROBI_INCLUDE_DIR)
endif(GUROBI_INCLUDE_DIR)
\ No newline at end of file
IF(UNIX)
IF(GUROBI_INCLUDE_DIR AND NOT GUROBI_LIBRARY)
IF (NOT DEFINED ENV{GUROBI_HOME})
message(FATAL_ERROR "Environment variable GUROBI_HOME not set. Set it so I can find the gurobi libs in \${GUROBI_HOME}/lib/.")
ENDIF()
SET(GUROBI_LIBRARY_DIR "$ENV{GUROBI_HOME}/lib/" CACHE PATH "Path to GUROBI Library")
IF (EXISTS "${GUROBI_LIBRARY_DIR}/libgurobi45.so")
SET(GUROBI_LIBRARY "gurobi45;gurobi_c++;pthread" CACHE STRING "GUROBI Libraries")
ELSEIF(EXISTS "${GUROBI_LIBRARY_DIR}/libgurobi46.so")
SET(GUROBI_LIBRARY "gurobi46;gurobi_c++;pthread" CACHE STRING "GUROBI Libraries")
ELSE()
message(FATAL_ERROR "Couldn't find a gurobi lib in ${GUROBI_LIBRARY_DIR}. Maybe it's a version I don't know about, yet.")
ENDIF()
ENDIF()
ENDIF(UNIX)
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment