Merge branch 'marinom/merge-from-ReForm' into 'master'

Merge from ReForm: Optimize MISolver::solve_multiple_rounding()

See merge request !56

Solver/GMM_ToolsT.cc

@@ -816,57 +816,53 @@ void eliminate_var_idx( const IntegerT           _evar,
 //-----------------------------------------------------------------------------
 
 
-template<class ScalarT, class VectorT, class RealT>
-void fix_var_csc_symmetric( const unsigned int                _i,
-			    const ScalarT                     _xi,
-			    typename gmm::csc_matrix<RealT>&  _A,
-			    VectorT&                          _x,
-			    VectorT&                          _rhs )
+template <class ScalarT, class VectorT, class RealT>
+void fix_var_csc_symmetric(const unsigned int _i, const ScalarT _xi,
+    typename gmm::csc_matrix<RealT>& _A, VectorT& _x, VectorT& _rhs)
 {
-// GMM CSC FORMAT
-//     T *pr;        // values.
-//     IND_TYPE *ir; // row indices.
-//     IND_TYPE *jc; // column repartition on pr and ir.
+  // GMM CSC FORMAT
+  //     T *pr;        // values.
+  //     IND_TYPE *ir; // row indices.
+  //     IND_TYPE *jc; // column repartition on pr and ir.
 
   gmm::size_type n = _A.nc;
 
   // update x
   _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
-  for( unsigned int iv=_A.jc[_i]; iv<_A.jc[_i+1]; ++iv)
+  // clear i-th column and collect non-zeros
+  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;
       _rhs[_i] = _xi;
     }
     else
     {
-      // update rhs
-      _rhs[_A.ir[iv]] -= _A.pr[iv]*_xi;
-      // clear entry
-      _A.pr[iv] = 0;
-      // store index
-      idx.push_back( _A.ir[iv]);
+      _rhs[_A.ir[iv]] -= _A.pr[iv] * _xi; // update rhs
+      _A.pr[iv] = 0; // clear entry
+      idx.push_back(_A.ir[iv]); // store index
     }
-
   }
 
-  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];
 
-    for( unsigned int j=_A.jc[col]; j<_A.jc[col+1]; ++j)
-      if(_A.ir[j] == _i)
+    for (unsigned int j = _A.jc[col]; j < _A.jc[col + 1]; ++j)
+    {
+      if (_A.ir[j] == _i)
       {
         _A.pr[j] = 0.0;
         // move to next
         break;
       }
     }
+  }
 }
 
 

Solver/IterativeSolverT.cc

@@ -12,143 +12,131 @@
 
 //== NAMESPACES ===============================================================
 
-namespace COMISO{
+namespace COMISO
+{
 
 //== IMPLEMENTATION ==========================================================
 
 template <class RealT>
-bool
-IterativeSolverT<RealT>::
-gauss_seidel_local( typename gmm::csc_matrix<Real>&  _A,
-		    std::vector<Real>&               _x, 
-		    std::vector<Real>&               _rhs, 
-		    std::vector<unsigned int>&       _idxs, 
-		    int&                             _max_iter, 
-		    Real&                            _tolerance )
+bool IterativeSolverT<RealT>::gauss_seidel_local(const Matrix& _A, Vector& _x,
+    const Vector& _rhs, const std::vector<unsigned int>& _idxs,
+    const int _max_iter, const 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;
 
   // clear old data
   i_temp.clear();
   q.clear();
 
-  for ( unsigned int i=0; i<_idxs.size(); ++i )
-    q.push_back( _idxs[i] );
+  for (unsigned int i = 0; i < _idxs.size(); ++i)
+    q.push_back(_idxs[i]);
 
   int it_count = 0;
 
-  while ( !q.empty() && it_count < _max_iter )
+  while (!q.empty() && it_count < _max_iter)
   {
     ++it_count;
-    unsigned int cur_i = q.front();
+    const auto i = q.front();
     q.pop_front();
     i_temp.clear();
 
-    ColT col = mat_const_col( _A, cur_i );
-
-    CIter it  = gmm::vect_const_begin( col );
-    CIter ite = gmm::vect_const_end( col );
-
-    double res_i   = -_rhs[cur_i];
-    double x_i_new = _rhs[cur_i];
+    double res_i = -_rhs[i];
+    double x_i_new = _rhs[i];
     double diag = 1.0;
-    for ( ; it!=ite; ++it )
+
+    const ColT col = mat_const_col(_A, i);
+    for (auto it = gmm::vect_const_begin(col), ite = gmm::vect_const_end(col);
+         it != ite; ++it)
     {
-      res_i   += ( *it ) * _x[it.index()];
-      x_i_new -= ( *it ) * _x[it.index()];
-      if( it.index() != cur_i)
-  	i_temp.push_back( static_cast<int>(it.index()) );
+      const auto j = static_cast<unsigned>(it.index());
+      res_i += (*it) * _x[j];
+      x_i_new -= (*it) * _x[j];
+      if (j != i)
+        i_temp.push_back(j);
       else
         diag = *it;
     }
 
     // take inverse of diag
-    diag = 1.0/diag;
+    diag = 1.0 / diag;
 
     // 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 )
-	q.push_back( i_temp[j] );
+      for (unsigned int j = 0; j < i_temp.size(); ++j)
+        q.push_back(i_temp[j]);
     }
   }
 
-  // converged?
-  return q.empty();
+  return q.empty(); // converged?
 }
 
-
 //-----------------------------------------------------------------------------
 
-
 template <class RealT>
-bool
-IterativeSolverT<RealT>::
-gauss_seidel_local2( typename gmm::csc_matrix<Real>&  _A,
-		     std::vector<Real>&               _x, 
-		     std::vector<Real>&               _rhs, 
-		     std::vector<unsigned int>&       _idxs, 
-		     int&                             _max_iter, 
-		     Real&                            _tolerance )
+bool IterativeSolverT<RealT>::gauss_seidel_local2(const Matrix& _A, Vector& _x,
+    const Vector& _rhs, const std::vector<unsigned int>& _idxs,
+    const int _max_iter, const 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;
 
-  double t2 = _tolerance*_tolerance;
+  double t2 = _tolerance * _tolerance;
 
   // clear old data
   i_temp.clear();
   s.clear();
 
-  for ( unsigned int i=0; i<_idxs.size(); ++i )
-    s.insert( _idxs[i] );
+  for (unsigned int i = 0; i < _idxs.size(); ++i)
+    s.insert(_idxs[i]);
 
   int it_count = 0;
 
   bool finished = false;
 
-  while ( !finished && it_count < _max_iter )
+  while (!finished && it_count < _max_iter)
   {
     finished = true;
     std::set<int>::iterator s_it = s.begin();
-    for(; s_it != s.end(); ++s_it)
+    for (; s_it != s.end(); ++s_it)
     {
       ++it_count;
       unsigned int cur_i = *s_it;
       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 ite = gmm::vect_const_end( col );
+      CIter it = gmm::vect_const_begin(col);
+      CIter ite = gmm::vect_const_end(col);
 
       double res_i = -_rhs[cur_i];
       double x_i_new = _rhs[cur_i];
       double diag = 1.0;
-      for ( ; it!=ite; ++it )
+      for (; it != ite; ++it)
       {
-	res_i   += ( *it ) * _x[it.index()];
-	x_i_new -= ( *it ) * _x[it.index()];
-	if( it.index() != cur_i)
-	  i_temp.push_back( static_cast<int>(it.index()) );
+        res_i += (*it) * _x[it.index()];
+        x_i_new -= (*it) * _x[it.index()];
+        if (it.index() != cur_i)
+          i_temp.push_back(static_cast<int>(it.index()));
         else
           diag = *it;
       }
 
       // 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;
-	  s.insert( i_temp[j] );
+          s.insert(i_temp[j]);
         }
       }
     }
@@ -161,13 +149,8 @@ gauss_seidel_local2( typename gmm::csc_matrix<Real>&  _A,
 //-----------------------------------------------------------------------------
 
 template <class RealT>
-bool
-IterativeSolverT<RealT>::
-conjugate_gradient( typename gmm::csc_matrix<Real>&  _A,
-		    std::vector<Real>&               _x, 
-		    std::vector<Real>&               _rhs, 
-		    int&                             _max_iter, 
-		    Real&                            _tolerance )
+bool IterativeSolverT<RealT>::conjugate_gradient(const Matrix& _A, Vector& _x,
+    const Vector& _rhs, int& _max_iter, Real& _tolerance)
 {
   Real rho, rho_1(0), a;
 
@@ -176,37 +159,37 @@ conjugate_gradient( typename gmm::csc_matrix<Real>&  _A,
   q_.resize(_x.size());
   r_.resize(_x.size());
   d_.resize(_x.size());
-  gmm::copy( _x, p_);
+  gmm::copy(_x, p_);
 
   // initialize diagonal (for relative norm)
-  for(unsigned int i=0; i<_x.size(); ++i)
-    d_[i] = 1.0/_A(i,i);
+  for (unsigned int i = 0; i < _x.size(); ++i)
+    d_[i] = 1.0 / _A(i, i);
 
   // start with iteration 0
   int cur_iter(0);
 
   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_);
 
   bool not_converged = true;
   Real res_norm(0);
 
   // 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)
   {
     //    std::cerr << "iter " << cur_iter << "  res " << res_norm << std::endl;
 
     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::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(q_, -a), r_);
     rho_1 = rho;
@@ -220,34 +203,27 @@ conjugate_gradient( typename gmm::csc_matrix<Real>&  _A,
   return (!not_converged);
 }
 
-
-
 //-----------------------------------------------------------------------------
 
-
 template <class RealT>
-typename IterativeSolverT<RealT>::Real
-IterativeSolverT<RealT>::
-vect_norm_rel(const std::vector<Real>& _v, const std::vector<Real>& _diag) const
+typename IterativeSolverT<RealT>::Real IterativeSolverT<RealT>::vect_norm_rel(
+    const Vector& _v, const Vector& _diag) const
 {
   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]);
-//     if(cur > res)
-//       res = cur;
+    //     Real cur = fabs(_v[i]*_diag[i]);
+    //     if(cur > res)
+    //       res = cur;
   }
   return res;
 }
 
-
 //-----------------------------------------------------------------------------
 
-
-
 //=============================================================================
 } // namespace COMISO
 //=============================================================================

Solver/IterativeSolverT.hh

@@ -4,29 +4,24 @@
 //
 //=============================================================================
 
-
 #ifndef COMISO_ITERATIVESOLVERT_HH
 #define COMISO_ITERATIVESOLVERT_HH
 
-
 //== INCLUDES =================================================================
 
 #include <CoMISo/Utils/gmm.hh>
 #include <deque>
-#include <queue>
 #include <set>
 
 //== FORWARDDECLARATIONS ======================================================
 
 //== NAMESPACES ===============================================================
 
-namespace COMISO {
+namespace COMISO
+{
 
 //== CLASS DEFINITION =========================================================
 
-
-
-	      
 /** \class IterativeSolverT IterativeSolverT.hh <COMISO/.../IterativeSolverT.hh>
 
     Brief Description.
@@ -34,47 +29,37 @@ namespace COMISO {
     A more elaborate description follows.
 */
 
-template <class RealT>
-class IterativeSolverT
+template <class RealT> class IterativeSolverT
 {
 public:
   typedef RealT Real;
+  typedef std::vector<Real> Vector;
+  typedef gmm::csc_matrix<Real> Matrix;
 
   // local gauss_seidel
-  bool gauss_seidel_local( typename gmm::csc_matrix<Real>&  _A,
-   			   std::vector<Real>&               _x, 
-			   std::vector<Real>&               _rhs, 
-  			   std::vector<unsigned int>&       _idxs, 
-  			   int&                             _max_iter, 
-  			   Real&                            _tolerance );
+  bool gauss_seidel_local(const Matrix& _A, Vector& _x, const Vector& _rhs,
+      const std::vector<unsigned int>& _idxs, const int _max_iter,
+      const Real& _tolerance);
 
   // local gauss_seidel
-  bool gauss_seidel_local2( typename gmm::csc_matrix<Real>&  _A,
-			    std::vector<Real>&               _x, 
-			    std::vector<Real>&               _rhs, 
-			    std::vector<unsigned int>&       _idxs, 
-			    int&                             _max_iter, 
-			    Real&                            _tolerance );
+  bool gauss_seidel_local2(const Matrix& _A, Vector& _x, const Vector& _rhs,
+      const std::vector<unsigned int>& _idxs, const int _max_iter,
+      const Real& _tolerance);
 
   // conjugate gradient
-  bool conjugate_gradient( typename gmm::csc_matrix<Real>&  _A,
-			   std::vector<Real>&               _x, 
-			   std::vector<Real>&               _rhs, 
-			   int&                             _max_iter, 
-			   Real&                            _tolerance );
+  bool conjugate_gradient(const Matrix& _A, Vector& _x, const Vector& _rhs,
+      int& _max_iter, Real& _tolerance);
 
 private:
   // compute relative norm
-  Real vect_norm_rel(const std::vector<Real>& _v, const std::vector<Real>& _diag) const;
-
+  Real vect_norm_rel(const Vector& _v, const Vector& _diag) const;
 
 private: