Commit d946c11a authored by David Bommes's avatar David Bommes
Browse files

cholmod no longer required in MISolver, EigenLDLT sufficient instead

git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@171 1355f012-dd97-4b2f-ae87-10fa9f823a57
parent 6011647e
/*===========================================================================*\
* *
* CoMISo *
* Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
* www.rwth-graphics.de *
* *
*---------------------------------------------------------------------------*
* This file is part of CoMISo. *
* *
* CoMISo is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* CoMISo is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
#include "EigenLDLTSolver.hh"
namespace COMISO {
EigenLDLTSolver::EigenLDLTSolver() : n_(0)
{
}
//-----------------------------------------------------------------------------
EigenLDLTSolver::~EigenLDLTSolver()
{
}
//-----------------------------------------------------------------------------
bool EigenLDLTSolver::calc_system( const std::vector<int>& _colptr,
const std::vector<int>& _rowind,
const std::vector<double>& _values)
{
std::cerr << "Info: EigenLDLTSolver::calc_system( const std::vector<int> ...) not implemented yet..." << std::endl;
return false;
}
//-----------------------------------------------------------------------------
bool EigenLDLTSolver::update_system( const std::vector<int>& _colptr,
const std::vector<int>& _rowind,
const std::vector<double>& _values )
{
std::cerr << "Info: EigenLDLTSolver::update_system( const std::vector<int> ...) not implemented yet..." << std::endl;
return false;
}
//-----------------------------------------------------------------------------
bool EigenLDLTSolver::solve( double * _x, double * _b)
{
// map arrays to Eigen-Vectors
Eigen::Map<Eigen::VectorXd> x(_x,n_);
Eigen::Map<Eigen::VectorXd> b(_b,n_);
// solve for another right hand side:
x = ldlt_.solve(b);
return (ldlt_.info()==Eigen::Success);
}
//-----------------------------------------------------------------------------
int EigenLDLTSolver::dimension()
{
return n_;
}
//-----------------------------------------------------------------------------
bool EigenLDLTSolver::
solve ( std::vector<double>& _x0, std::vector<double>& _b)
{
return solve( &(_x0[0]), &(_b[0]));
}
//-----------------------------------------------------------------------------
bool& EigenLDLTSolver::
show_timings()
{
return show_timings_;
}
}
/*===========================================================================*\
* *
* CoMISo *
* Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
* www.rwth-graphics.de *
* *
*---------------------------------------------------------------------------*
* This file is part of CoMISo. *
* *
* CoMISo is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* CoMISo is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
//=============================================================================
//
// CLASS CholmodSolver
//
//=============================================================================
#ifndef COMISO_EIGEN_LDLT_SOLVER_HH
#define COMISO_EIGEN_LDLT_SOLVER_HH
//== INCLUDES =================================================================
#include <CoMISo/Config/CoMISoDefines.hh>
#include <CoMISo/Utils/StopWatch.hh>
#include <iostream>
#include <vector>
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
//== NAMESPACES ===============================================================
namespace COMISO {
//== CLASS DEFINITION =========================================================
class COMISODLLEXPORT EigenLDLTSolver
{
public:
// _size is maximal size this instance can handle (smaller problems are possible!!!)
EigenLDLTSolver();
~EigenLDLTSolver();
bool calc_system( const std::vector<int>& _colptr,
const std::vector<int>& _rowind,
const std::vector<double>& _values );
template< class GMM_MatrixT>
bool calc_system_gmm( const GMM_MatrixT& _mat);
template< class Eigen_MatrixT>
bool calc_system_eigen( const Eigen_MatrixT& _mat);
bool update_system( const std::vector<int>& _colptr,
const std::vector<int>& _rowind,
const std::vector<double>& _values );
template< class GMM_MatrixT>
bool update_system_gmm( const GMM_MatrixT& _mat);
template< class Eigen_MatrixT>
bool update_system_eigen( const Eigen_MatrixT& _mat);
bool solve ( double * _x0, double * _b);
bool solve ( std::vector<double>& _x0, std::vector<double>& _b);
bool& show_timings();
int dimension();
private:
// dimension n_
unsigned int n_;
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > ldlt_;
bool show_timings_;
StopWatch sw_;
};
//=============================================================================
} // namespace COMISO
//=============================================================================
#if defined(INCLUDE_TEMPLATES) && !defined(COMISO_EIGEN_LDLT_SOLVER_TEMPLATES_C)
#define COMISO_EIGEN_LDLT_SOLVER_TEMPLATES
#include "EigenLDLTSolverT.cc"
#endif
//=============================================================================
#endif // COMISO_EIGEN_LDLT_SOLVER_HH defined
//=============================================================================
/*===========================================================================*\
* *
* CoMISo *
* Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
* www.rwth-graphics.de *
* *
*---------------------------------------------------------------------------*
* This file is part of CoMISo. *
* *
* CoMISo is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* CoMISo is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
#define COMISO_EIGEN_LDLT_SOLVER_TEMPLATES_C
#include <CoMISo/Solver/GMM_Tools.hh>
#include <CoMISo/Solver/Eigen_Tools.hh>
#include <CoMISo/Solver/EigenLDLTSolver.hh>
namespace COMISO {
template< class GMM_MatrixT>
bool EigenLDLTSolver::calc_system_gmm( const GMM_MatrixT& _mat)
{
if(show_timings_) sw_.start();
Eigen::SparseMatrix<double> E;
COMISO_EIGEN::gmm_to_eigen(_mat, E);
if(show_timings_)
{
std::cerr << "EigenLDLT Timing GMM convert: " << sw_.stop()/1000.0 << "s\n";
std::cerr << "#nnz: " << E.nonZeros() << std::endl;
}
return calc_system_eigen( E);
}
//-----------------------------------------------------------------------------
template< class GMM_MatrixT>
bool EigenLDLTSolver::update_system_gmm( const GMM_MatrixT& _mat)
{
if(show_timings_) sw_.start();
Eigen::SparseMatrix<double> E;
COMISO_EIGEN::gmm_to_eigen(_mat, E);
if(show_timings_)
{
std::cerr << "EigenLDLT Timing GMM convert: " << sw_.stop()/1000.0 << "s\n";
std::cerr << "#nnz: " << E.nonZeros() << std::endl;
}
return update_system_eigen( E);
}
//-----------------------------------------------------------------------------
template< class Eigen_MatrixT>
bool EigenLDLTSolver::calc_system_eigen( const Eigen_MatrixT& _mat)
{
n_ = _mat.rows();
if(show_timings_) sw_.start();
ldlt_.compute(_mat);
if(show_timings_)
{
std::cerr << "EigenLDLT Timing EIGEN compute: " << sw_.stop()/1000.0 << "s\n";
}
return (ldlt_.info()==Eigen::Success);
}
//-----------------------------------------------------------------------------
template< class Eigen_MatrixT>
bool EigenLDLTSolver::update_system_eigen( const Eigen_MatrixT& _mat)
{
if(show_timings_) sw_.start();
ldlt_.factorize(_mat);
if(show_timings_)
{
std::cerr << "EigenLDLT Timing EIGEN factorize: " << sw_.stop()/1000.0 << "s\n";
}
return (ldlt_.info()==Eigen::Success );
}
}
......@@ -637,7 +637,36 @@ void eigen_to_cholmod_dense( const MatrixT& _A, cholmod_dense* &_AC, cholmod_com
}*/
// convert a gmm col-sparse matrix into an eigen sparse matrix
template<class GMM_MatrixT, class EIGEN_MatrixT>
void gmm_to_eigen( const GMM_MatrixT& _G, EIGEN_MatrixT& _E)
{
#ifdef COMISO_Eigen3_AVAILABLE
typedef typename EIGEN_MatrixT::Scalar Scalar;
typedef typename gmm::linalg_traits<GMM_MatrixT>::const_sub_col_type ColT;
typedef typename gmm::linalg_traits<ColT>::const_iterator CIter;
// build matrix triplets
typedef Eigen::Triplet< Scalar > Triplet;
std::vector< Triplet > triplets;
triplets.reserve(gmm::nnz(_G));
for(unsigned int i=0; i<gmm::mat_ncols(_G); ++i)
{
ColT col = mat_const_col( _G, i );
CIter it = gmm::vect_const_begin( col );
CIter ite = gmm::vect_const_end( col );
for ( ; it!=ite; ++it )
triplets.push_back( Triplet( it.index(), i, *it));
}
// generate eigen matrix
_E = EIGEN_MatrixT( gmm::mat_nrows(_G), gmm::mat_ncols(_G));
_E.setFromTriplets( triplets.begin(), triplets.end());
#endif
}
#endif
......
......@@ -37,12 +37,9 @@
#ifdef COMISO_Eigen3_AVAILABLE
#include <Eigen/Dense>
#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
#include <Eigen/Sparse>
#endif
//#include <Eigen/Eigen>
//#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
//#include <Eigen/Sparse>
//#endif
#ifndef COMISO_NCHOLMOD
#include <cholmod.h>
......@@ -108,6 +105,10 @@ void eigen_to_cholmod( const MatrixT& _A,
bool _long_int = false);
#endif
// convert a gmm column-sparse matrix into an eigen sparse matrix
template<class GMM_MatrixT, class EIGEN_MatrixT>
void gmm_to_eigen( const GMM_MatrixT& _G, EIGEN_MatrixT& _E);
//=============================================================================
} // namespace COMISO_Eigen
......
......@@ -42,6 +42,7 @@
// hack for testing only
#include "SparseQRSolver.hh"
#include "UMFPACKSolver.hh"
#include "EigenLDLTSolver.hh"
#define ROUND(x) ((x)<0?int((x)-0.5):int((x)+0.5))
......@@ -88,8 +89,8 @@ MISolver::solve_no_rounding(
Vecd& _x,
Vecd& _rhs )
{
chol_.calc_system_gmm(_A);
chol_.solve(_x, _rhs);
direct_solver_.calc_system_gmm(_A);
direct_solver_.solve(_x, _rhs);
}
......@@ -101,7 +102,7 @@ MISolver::resolve(
Vecd& _x,
Vecd& _rhs )
{
chol_.solve(_x, _rhs);
direct_solver_.solve(_x, _rhs);
}
......@@ -127,13 +128,14 @@ MISolver::solve_direct_rounding(
Veci old_idx(_rhs.size());
for(unsigned int i=0; i<old_idx.size(); ++i)
old_idx[i] = i;
chol_.calc_system_gmm(_A);
chol_.solve(_x, _rhs);
direct_solver_.calc_system_gmm(_A);
direct_solver_.solve(_x, _rhs);
// check solver performance (only for testing!!!)
{
StopWatch sw;
// hack
const bool enable_performance_test = false;
// performance comparison code
......@@ -155,6 +157,7 @@ MISolver::solve_direct_rounding(
#endif
// performance comparison code
#if(COMISO_SUITESPARSE_AVAILABLE)
if(enable_performance_test)
{
sw.start();
......@@ -185,6 +188,25 @@ MISolver::solve_direct_rounding(
gmm::add(_x,gmm::scaled(x4,-1.0),res);
std::cerr << "DIFFERENCE IN RESULT: " << gmm::vect_norm2(res) << std::endl;
}
#endif
#if(COMISO_Eigen3_AVAILABLE)
// performance comparison code
if(enable_performance_test)
{
sw.start();
COMISO::EigenLDLTSolver ldlt;
ldlt.calc_system_gmm(_A);
std::cerr << "Eigen LDLT factor took: " << sw.stop()/1000.0 << "s\n";
Vecd x5(_x);
sw.start();
ldlt.solve(x5,_rhs);
std::cerr << "Eigen LDLT solve took: " << sw.stop()/1000.0 << "s\n";
Vecd res(_x);
gmm::add(_x,gmm::scaled(x5,-1.0),res);
std::cerr << "DIFFERENCE IN RESULT: " << gmm::vect_norm2(res) << std::endl;
}
#endif
}
// round and eliminate variables
......@@ -213,9 +235,9 @@ MISolver::solve_direct_rounding(
// final full solution
if( gmm::mat_ncols( _A) > 0)
{
// chol_.update_system_gmm(_A);
chol_.calc_system_gmm(_A);
chol_.solve( xr, _rhs);
// direct_solver_.update_system_gmm(_A);
direct_solver_.calc_system_gmm(_A);
direct_solver_.solve( xr, _rhs);
}
// store solution values to result vector
......@@ -269,8 +291,8 @@ MISolver::solve_iterative(
if( initial_full_solution_)
{
if( noisy_ > 2) std::cerr << "initial full solution" << std::endl;
chol_.calc_system_gmm(_A);
chol_.solve(_x, _rhs);
direct_solver_.calc_system_gmm(_A);
direct_solver_.solve(_x, _rhs);
cholmod_step_done_ = true;
......@@ -355,11 +377,11 @@ MISolver::solve_iterative(
if( gmm::mat_ncols( _A) > 0)
{
if(cholmod_step_done_)
chol_.update_system_gmm(_A);
direct_solver_.update_system_gmm(_A);
else
chol_.calc_system_gmm(_A);
direct_solver_.calc_system_gmm(_A);
chol_.solve( xr, _rhs);
direct_solver_.solve( xr, _rhs);
++n_full_;
}
}
......@@ -434,13 +456,13 @@ MISolver::update_solution(
if( gmm::mat_ncols( _A) > 0)
{
if(cholmod_step_done_)
chol_.update_system_gmm(_A);
direct_solver_.update_system_gmm(_A);
else
{
chol_.calc_system_gmm(_A);
direct_solver_.calc_system_gmm(_A);
cholmod_step_done_ = true;
}
chol_.solve(_x,_rhs);
direct_solver_.solve(_x,_rhs);
++n_full_;
}
......@@ -487,8 +509,8 @@ MISolver::solve_multiple_rounding(
if( initial_full_solution_)
{
if( noisy_ > 2) std::cerr << "initial full solution" << std::endl;
chol_.calc_system_gmm(_A);
chol_.solve(_x, _rhs);
direct_solver_.calc_system_gmm(_A);
direct_solver_.solve(_x, _rhs);
cholmod_step_done_ = true;
......@@ -579,11 +601,11 @@ MISolver::solve_multiple_rounding(
if( gmm::mat_ncols( _A) > 0)
{
if(cholmod_step_done_)
chol_.update_system_gmm(_A);
direct_solver_.update_system_gmm(_A);
else
chol_.calc_system_gmm(_A);
direct_solver_.calc_system_gmm(_A);
chol_.solve( xr, _rhs);
direct_solver_.solve( xr, _rhs);
++n_full_;
}
}
......
......@@ -38,8 +38,15 @@
#include <CoMISo/Config/CoMISoDefines.hh>
#include <CoMISo/Config/config.hh>
#if COMISO_SUITESPARSE_AVAILABLE
#include "CholmodSolver.hh"
#elif COMISO_Eigen3_AVAILABLE
#include "EigenLDLTSolver.hh"
#else
#print "Warning: MISolver requires Suitesparse or Eigen3 support"
#endif
#include "GMM_Tools.hh"
#include "CholmodSolver.hh"
#include "IterativeSolverT.hh"
#include <vector>
......@@ -360,7 +367,14 @@ private:
// flag
bool cholmod_step_done_;
COMISO::CholmodSolver chol_;
// declar direct solver depending on availability
#if COMISO_SUITESPARSE_AVAILABLE
COMISO::CholmodSolver direct_solver_;
#elif COMISO_Eigen3_AVAILABLE
COMISO::EigenLDLTSolver direct_solver_;
#else
#print "Warning: MISolver requires Suitesparse or Eigen3 support"
#endif
IterativeSolverT<double> siter_;
......
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