Commit 282232b5 authored by Martin Marinov's avatar Martin Marinov Committed by GitHub Enterprise
Browse files

Optimize MISolver::solve_multiple_rounding (#3)

* Avoid processing twice indices of values to round
* Replace the rounding set with a priority queue
* Remove unnecessary "old index" map
* Reformat MISolver.cc to conform to the style guidelines and add missing ; after DEB_* statements
* Remove MISolverT.hh
* Move private symbols to the cc file
* Improve code quality, formatting and const-correctness of IterativeSolverT
* Improve code quality and formatting for MISolver::update_solution
parent 74c18843
...@@ -20,7 +20,6 @@ SET(my_t_impls ...@@ -20,7 +20,6 @@ SET(my_t_impls
${CMAKE_CURRENT_SOURCE_DIR}/GMM_ToolsT.cc ${CMAKE_CURRENT_SOURCE_DIR}/GMM_ToolsT.cc
${CMAKE_CURRENT_SOURCE_DIR}/IterativeSolverT.cc ${CMAKE_CURRENT_SOURCE_DIR}/IterativeSolverT.cc
${CMAKE_CURRENT_SOURCE_DIR}/SparseQRSolverT.cc ${CMAKE_CURRENT_SOURCE_DIR}/SparseQRSolverT.cc
${CMAKE_CURRENT_SOURCE_DIR}/MISolverT.cc
${CMAKE_CURRENT_SOURCE_DIR}/TaucsSolverT.cc ${CMAKE_CURRENT_SOURCE_DIR}/TaucsSolverT.cc
${CMAKE_CURRENT_SOURCE_DIR}/UMFPACKSolverT.cc ${CMAKE_CURRENT_SOURCE_DIR}/UMFPACKSolverT.cc
PARENT_SCOPE PARENT_SCOPE
......
...@@ -816,56 +816,52 @@ void eliminate_var_idx( const IntegerT _evar, ...@@ -816,56 +816,52 @@ void eliminate_var_idx( const IntegerT _evar,
//----------------------------------------------------------------------------- //-----------------------------------------------------------------------------
template<class ScalarT, class VectorT, class RealT> template <class ScalarT, class VectorT, class RealT>
void fix_var_csc_symmetric( const unsigned int _i, void fix_var_csc_symmetric(const unsigned int _i, const ScalarT _xi,
const ScalarT _xi, typename gmm::csc_matrix<RealT>& _A, VectorT& _x, VectorT& _rhs)
typename gmm::csc_matrix<RealT>& _A,
VectorT& _x,
VectorT& _rhs )
{ {
// GMM CSC FORMAT // GMM CSC FORMAT
// T *pr; // values. // T *pr; // values.
// IND_TYPE *ir; // row indices. // IND_TYPE *ir; // row indices.
// IND_TYPE *jc; // column repartition on pr and ir. // IND_TYPE *jc; // column repartition on pr and ir.
gmm::size_type n = _A.nc; gmm::size_type n = _A.nc;
// update x // update x
_x[_i] = _xi; _x[_i] = _xi;
std::vector<unsigned int> idx; idx.reserve(16); std::vector<unsigned int> idx;
idx.reserve(16);
// clear i-th column and collect nonzeros // clear i-th column and collect non-zeros
for( unsigned int iv=_A.jc[_i]; iv<_A.jc[_i+1]; ++iv) for (unsigned int iv = _A.jc[_i]; iv < _A.jc[_i + 1]; ++iv)
{ {
if(_A.ir[iv] == _i) if (_A.ir[iv] == _i)
{ {
_A.pr[iv] = 1.0; _A.pr[iv] = 1.0;
_rhs[_i] = _xi; _rhs[_i] = _xi;
} }
else else
{ {
// update rhs _rhs[_A.ir[iv]] -= _A.pr[iv] * _xi; // update rhs
_rhs[_A.ir[iv]] -= _A.pr[iv]*_xi; _A.pr[iv] = 0; // clear entry
// clear entry idx.push_back(_A.ir[iv]); // store index
_A.pr[iv] = 0;
// store index
idx.push_back( _A.ir[iv]);
} }
} }
for(std::size_t i=0; i<idx.size(); ++i) for (std::size_t i = 0; i < idx.size(); ++i)
{ {
unsigned int col = idx[i]; unsigned int col = idx[i];
for( unsigned int j=_A.jc[col]; j<_A.jc[col+1]; ++j) for (unsigned int j = _A.jc[col]; j < _A.jc[col + 1]; ++j)
if(_A.ir[j] == _i) {
if (_A.ir[j] == _i)
{ {
_A.pr[j] = 0.0; _A.pr[j] = 0.0;
// move to next // move to next
break; break;
} }
}
} }
} }
......
...@@ -12,144 +12,132 @@ ...@@ -12,144 +12,132 @@
//== NAMESPACES =============================================================== //== NAMESPACES ===============================================================
namespace COMISO{ namespace COMISO
{
//== IMPLEMENTATION ========================================================== //== IMPLEMENTATION ==========================================================
template <class RealT> template <class RealT>
bool bool IterativeSolverT<RealT>::gauss_seidel_local(const Matrix& _A, Vector& _x,
IterativeSolverT<RealT>:: const Vector& _rhs, const std::vector<unsigned int>& _idxs,
gauss_seidel_local( typename gmm::csc_matrix<Real>& _A, const int _max_iter, const Real& _tolerance)
std::vector<Real>& _x,
std::vector<Real>& _rhs,
std::vector<unsigned int>& _idxs,
int& _max_iter,
Real& _tolerance )
{ {
if( _max_iter == 0) return false; if (_max_iter == 0)
return false;
typedef typename gmm::linalg_traits< gmm::csc_matrix<Real> >::const_sub_col_type ColT; typedef typename gmm::linalg_traits<Matrix>::const_sub_col_type ColT;
typedef typename gmm::linalg_traits<ColT>::const_iterator CIter; typedef typename gmm::linalg_traits<ColT>::const_iterator CIter;
// clear old data // clear old data
i_temp.clear(); i_temp.clear();
q.clear(); q.clear();
for ( unsigned int i=0; i<_idxs.size(); ++i ) for (unsigned int i = 0; i < _idxs.size(); ++i)
q.push_back( _idxs[i] ); q.push_back(_idxs[i]);
int it_count = 0; int it_count = 0;
while ( !q.empty() && it_count < _max_iter ) while (!q.empty() && it_count < _max_iter)
{ {
++it_count; ++it_count;
unsigned int cur_i = q.front(); const auto i = q.front();
q.pop_front(); q.pop_front();
i_temp.clear(); i_temp.clear();
ColT col = mat_const_col( _A, cur_i ); double res_i = -_rhs[i];
double x_i_new = _rhs[i];
CIter it = gmm::vect_const_begin( col ); double diag = 1.0;
CIter ite = gmm::vect_const_end( col );
double res_i = -_rhs[cur_i]; const ColT col = mat_const_col(_A, i);
double x_i_new = _rhs[cur_i]; for (auto it = gmm::vect_const_begin(col), ite = gmm::vect_const_end(col);
double diag = 1.0; it != ite; ++it)
for ( ; it!=ite; ++it )
{ {
res_i += ( *it ) * _x[it.index()]; const auto j = static_cast<unsigned>(it.index());
x_i_new -= ( *it ) * _x[it.index()]; res_i += (*it) * _x[j];
if( it.index() != cur_i) x_i_new -= (*it) * _x[j];
i_temp.push_back( static_cast<int>(it.index()) ); if (j != i)
i_temp.push_back(j);
else else
diag = *it; diag = *it;
} }
// take inverse of diag // take inverse of diag
diag = 1.0/diag; diag = 1.0 / diag;
// compare relative residuum normalized by diagonal entry // compare relative residuum normalized by diagonal entry
if ( fabs(res_i*diag) > _tolerance ) if (fabs(res_i * diag) > _tolerance)
{ {
_x[cur_i] += x_i_new*diag; _x[i] += x_i_new * diag;
for ( unsigned int j=0; j<i_temp.size(); ++j ) for (unsigned int j = 0; j < i_temp.size(); ++j)
q.push_back( i_temp[j] ); q.push_back(i_temp[j]);
} }
} }
// converged? return q.empty(); // converged?
return q.empty();
} }
//----------------------------------------------------------------------------- //-----------------------------------------------------------------------------
template <class RealT> template <class RealT>
bool bool IterativeSolverT<RealT>::gauss_seidel_local2(const Matrix& _A, Vector& _x,
IterativeSolverT<RealT>:: const Vector& _rhs, const std::vector<unsigned int>& _idxs,
gauss_seidel_local2( typename gmm::csc_matrix<Real>& _A, const int _max_iter, const Real& _tolerance)
std::vector<Real>& _x,
std::vector<Real>& _rhs,
std::vector<unsigned int>& _idxs,
int& _max_iter,
Real& _tolerance )
{ {
typedef typename gmm::linalg_traits< gmm::csc_matrix<Real> >::const_sub_col_type ColT; typedef typename gmm::linalg_traits<Matrix>::const_sub_col_type ColT;
typedef typename gmm::linalg_traits<ColT>::const_iterator CIter; typedef typename gmm::linalg_traits<ColT>::const_iterator CIter;
double t2 = _tolerance*_tolerance; double t2 = _tolerance * _tolerance;
// clear old data // clear old data
i_temp.clear(); i_temp.clear();
s.clear(); s.clear();
for ( unsigned int i=0; i<_idxs.size(); ++i ) for (unsigned int i = 0; i < _idxs.size(); ++i)
s.insert( _idxs[i] ); s.insert(_idxs[i]);
int it_count = 0; int it_count = 0;
bool finished = false; bool finished = false;
while ( !finished && it_count < _max_iter ) while (!finished && it_count < _max_iter)
{ {
finished = true; finished = true;
std::set<int>::iterator s_it = s.begin(); std::set<int>::iterator s_it = s.begin();
for(; s_it != s.end(); ++s_it) for (; s_it != s.end(); ++s_it)
{ {
++it_count; ++it_count;
unsigned int cur_i = *s_it; unsigned int cur_i = *s_it;
i_temp.clear(); i_temp.clear();
ColT col = mat_const_col( _A, cur_i ); ColT col = mat_const_col(_A, cur_i);
CIter it = gmm::vect_const_begin( col ); CIter it = gmm::vect_const_begin(col);
CIter ite = gmm::vect_const_end( col ); CIter ite = gmm::vect_const_end(col);
double res_i = -_rhs[cur_i]; double res_i = -_rhs[cur_i];
double x_i_new = _rhs[cur_i]; double x_i_new = _rhs[cur_i];
double diag = 1.0; double diag = 1.0;
for ( ; it!=ite; ++it ) for (; it != ite; ++it)
{ {
res_i += ( *it ) * _x[it.index()]; res_i += (*it) * _x[it.index()];
x_i_new -= ( *it ) * _x[it.index()]; x_i_new -= (*it) * _x[it.index()];
if( it.index() != cur_i) if (it.index() != cur_i)
i_temp.push_back( static_cast<int>(it.index()) ); i_temp.push_back(static_cast<int>(it.index()));
else else
diag = *it; diag = *it;
} }
// compare relative residuum normalized by diagonal entry // compare relative residuum normalized by diagonal entry
if ( res_i*res_i/diag > t2 ) if (res_i * res_i / diag > t2)
{ {
_x[cur_i] += x_i_new/_A( cur_i, cur_i ); _x[cur_i] += x_i_new / _A(cur_i, cur_i);
for ( unsigned int j=0; j<i_temp.size(); ++j ) for (unsigned int j = 0; j < i_temp.size(); ++j)
{ {
finished = false; finished = false;
s.insert( i_temp[j] ); s.insert(i_temp[j]);
} }
} }
} }
} }
...@@ -157,17 +145,12 @@ gauss_seidel_local2( typename gmm::csc_matrix<Real>& _A, ...@@ -157,17 +145,12 @@ gauss_seidel_local2( typename gmm::csc_matrix<Real>& _A,
// converged? // converged?
return finished; return finished;
} }
//----------------------------------------------------------------------------- //-----------------------------------------------------------------------------
template <class RealT> template <class RealT>
bool bool IterativeSolverT<RealT>::conjugate_gradient(const Matrix& _A, Vector& _x,
IterativeSolverT<RealT>:: const Vector& _rhs, int& _max_iter, Real& _tolerance)
conjugate_gradient( typename gmm::csc_matrix<Real>& _A,
std::vector<Real>& _x,
std::vector<Real>& _rhs,
int& _max_iter,
Real& _tolerance )
{ {
Real rho, rho_1(0), a; Real rho, rho_1(0), a;
...@@ -176,78 +159,71 @@ conjugate_gradient( typename gmm::csc_matrix<Real>& _A, ...@@ -176,78 +159,71 @@ conjugate_gradient( typename gmm::csc_matrix<Real>& _A,
q_.resize(_x.size()); q_.resize(_x.size());
r_.resize(_x.size()); r_.resize(_x.size());
d_.resize(_x.size()); d_.resize(_x.size());
gmm::copy( _x, p_); gmm::copy(_x, p_);
// initialize diagonal (for relative norm) // initialize diagonal (for relative norm)
for(unsigned int i=0; i<_x.size(); ++i) for (unsigned int i = 0; i < _x.size(); ++i)
d_[i] = 1.0/_A(i,i); d_[i] = 1.0 / _A(i, i);
// start with iteration 0 // start with iteration 0
int cur_iter(0); int cur_iter(0);
gmm::mult(_A, gmm::scaled(_x, Real(-1)), _rhs, r_); gmm::mult(_A, gmm::scaled(_x, Real(-1)), _rhs, r_);
rho = gmm::vect_sp( r_, r_); rho = gmm::vect_sp(r_, r_);
gmm::copy(r_, p_); gmm::copy(r_, p_);
bool not_converged = true; bool not_converged = true;
Real res_norm(0); Real res_norm(0);
// while not converged // while not converged
while( (not_converged = ( (res_norm=vect_norm_rel(r_, d_)) > _tolerance)) && while ((not_converged = ((res_norm = vect_norm_rel(r_, d_)) > _tolerance)) &&
cur_iter < _max_iter) cur_iter < _max_iter)
{ {
// std::cerr << "iter " << cur_iter << " res " << res_norm << std::endl; // std::cerr << "iter " << cur_iter << " res " << res_norm << std::endl;
if (cur_iter != 0) if (cur_iter != 0)
{ {
rho = gmm::vect_sp( r_, r_); rho = gmm::vect_sp(r_, r_);
gmm::add(r_, gmm::scaled(p_, rho / rho_1), p_); gmm::add(r_, gmm::scaled(p_, rho / rho_1), p_);
} }
gmm::mult(_A, p_, q_); gmm::mult(_A, p_, q_);
a = rho / gmm::vect_sp( q_, p_); a = rho / gmm::vect_sp(q_, p_);
gmm::add(gmm::scaled(p_, a), _x); gmm::add(gmm::scaled(p_, a), _x);
gmm::add(gmm::scaled(q_, -a), r_); gmm::add(gmm::scaled(q_, -a), r_);
rho_1 = rho; rho_1 = rho;
++cur_iter; ++cur_iter;
} }
_max_iter = cur_iter; _max_iter = cur_iter;
_tolerance = res_norm; _tolerance = res_norm;
return (!not_converged); return (!not_converged);
} }
//----------------------------------------------------------------------------- //-----------------------------------------------------------------------------
template <class RealT> template <class RealT>
typename IterativeSolverT<RealT>::Real typename IterativeSolverT<RealT>::Real IterativeSolverT<RealT>::vect_norm_rel(
IterativeSolverT<RealT>:: const Vector& _v, const Vector& _diag) const
vect_norm_rel(const std::vector<Real>& _v, const std::vector<Real>& _diag) const
{ {
Real res = 0.0; Real res = 0.0;
for(unsigned int i=0; i<_v.size(); ++i) for (unsigned int i = 0; i < _v.size(); ++i)
{ {
res = std::max(fabs(_v[i]*_diag[i]), res); res = std::max(fabs(_v[i] * _diag[i]), res);
// Real cur = fabs(_v[i]*_diag[i]); // Real cur = fabs(_v[i]*_diag[i]);
// if(cur > res) // if(cur > res)
// res = cur; // res = cur;
} }
return res; return res;
} }
//----------------------------------------------------------------------------- //-----------------------------------------------------------------------------
//============================================================================= //=============================================================================
} // namespace COMISO } // namespace COMISO
//============================================================================= //=============================================================================
...@@ -4,85 +4,69 @@ ...@@ -4,85 +4,69 @@
// //
//============================================================================= //=============================================================================
#ifndef COMISO_ITERATIVESOLVERT_HH #ifndef COMISO_ITERATIVESOLVERT_HH
#define COMISO_ITERATIVESOLVERT_HH #define COMISO_ITERATIVESOLVERT_HH
//== INCLUDES ================================================================= //== INCLUDES =================================================================
#include <CoMISo/Utils/gmm.hh> #include <CoMISo/Utils/gmm.hh>
#include <deque> #include <deque>
#include <queue>
#include <set> #include <set>
//== FORWARDDECLARATIONS ====================================================== //== FORWARDDECLARATIONS ======================================================
//== NAMESPACES =============================================================== //== NAMESPACES ===============================================================
namespace COMISO { namespace COMISO
{
//== CLASS DEFINITION ========================================================= //== CLASS DEFINITION =========================================================
/** \class IterativeSolverT IterativeSolverT.hh <COMISO/.../IterativeSolverT.hh> /** \class IterativeSolverT IterativeSolverT.hh <COMISO/.../IterativeSolverT.hh>
Brief Description. Brief Description.
A more elaborate description follows. A more elaborate description follows.
*/ */
template <class RealT> template <class RealT> class IterativeSolverT
class IterativeSolverT
{ {
public: public:
typedef RealT Real; typedef RealT Real;
typedef std::vector<Real> Vector;
typedef gmm::csc_matrix<Real> Matrix;
// local gauss_seidel // local gauss_seidel
bool gauss_seidel_local( typename gmm::csc_matrix<Real>& _A,