Merge commit '2fe50f60' into marinom/merge-from-ReForm

Solver/GMM_ToolsT.cc

@@ -816,56 +816,52 @@ 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;
+  _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)
-  {
-    if(_A.ir[iv] == _i)
+  // 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)
     {
       _A.pr[iv] = 1.0;
-      _rhs[_i]  = _xi;
+      _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;
+        _A.pr[j] = 0.0;
+        // move to next
+        break;
       }
+    }
   }
 }
 

Solver/IterativeSolverT.cc

@@ -12,144 +12,132 @@
 
 //== 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[i];
+    double x_i_new = _rhs[i];
+    double diag = 1.0;
 
-    double res_i   = -_rhs[cur_i];
-    double x_i_new = _rhs[cur_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;
+        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 );
-      
-      double res_i   = -_rhs[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 diag    = 1.0;
-      for ( ; it!=ite; ++it )
+      double diag = 1.0;
+      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()) );
-	else
-	  diag = *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()));
+        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 );
-	
-	for ( unsigned int j=0; j<i_temp.size(); ++j )
-	{
-	  finished = false;
-	  s.insert( i_temp[j] );
-	}
+        _x[cur_i] += x_i_new / _A(cur_i, cur_i);
+
+        for (unsigned int j = 0; j < i_temp.size(); ++j)
+        {
+          finished = false;
+          s.insert(i_temp[j]);
+        }
       }
     }
   }
@@ -157,17 +145,12 @@ gauss_seidel_local2( typename gmm::csc_matrix<Real>&  _A,
   // converged?
   return finished;
 }
-    
+
 //-----------------------------------------------------------------------------
 
 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,78 +159,71 @@ 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)) && 
-	 cur_iter < _max_iter)
+  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_);
+    if (cur_iter != 0)
+    {
+      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;
 
     ++cur_iter;
   }
-  
-  _max_iter  = cur_iter;
+
+  _max_iter = cur_iter;
   _tolerance = res_norm;
 
   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,85 +4,69 @@
 //
 //=============================================================================
 
-
 #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.
-  
+
     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:
-
   // helper for conjugate gradient
-  std::vector<Real> p_;
-  std::vector<Real> q_;
-  std::vector<Real> r_;
-  std::vector<Real> d_;
+  Vector p_;
+  Vector q_;
+  Vector r_;
+  Vector d_;
 
   //  helper for local gauss seidel
-  std::vector<unsigned int>  i_temp;
+  std::vector<unsigned int> i_temp;
   std::deque<unsigned int> q;
   std::set<int> s;
 };
 
-
 //=============================================================================
 } // namespace COMISO
 //=============================================================================
@@ -93,4 +77,3 @@ private:
 //=============================================================================
 #endif // COMISO_ITERATIVESOLVERT_HH defined
 //=============================================================================
-

Solver/MISolver.hh

@@ -88,11 +88,6 @@ public:
   typedef std::vector< int >              Veci;
   typedef std::vector< unsigned int >     Vecui;
 
-  // gmm Column and ColumnIterator of CSC matrix
-  typedef gmm::linalg_traits< CSCMatrix >::const_sub_col_type Col;
-  typedef gmm::linalg_traits< Col >::const_iterator           ColIter;
-
-
   /// default Constructor
   MISolver();
 
@@ -222,66 +217,6 @@ public:
   /// Set/Get use_constraint_reordering for constraint solver (default = true)
   bool& use_constraint_reordering() { return use_constraint_reordering_;}
 
- private:
-
-  // find set of variables for simultaneous rounding
-  class RoundingSet
-  {
-  public:
-    typedef std::pair<double,int> PairDI;
-    
-    RoundingSet() : threshold_(0.5), cur_sum_(0.0) {}
-
-    void clear() { rset_.clear(); cur_sum_ = 0.0;}
-
-    bool add( int _id, double _rd_val)
-    {
-      // empty set? -> always add one element
-      if( rset_.empty() || cur_sum_+_rd_val <= threshold_)
-      {
-	rset_.insert( PairDI(_rd_val,_id) );
-	cur_sum_ += _rd_val;
-	return true;
-      }
-      else
-      {
-	// move to last element
-	std::set<PairDI>::iterator s_it = rset_.end();
-	--s_it;
-
-	if( s_it->first > _rd_val)
-	{
-	  cur_sum_ -= s_it->first;
-	  rset_.erase(s_it);
-	  rset_.insert( PairDI(_rd_val,_id) );
-	  cur_sum_ += _rd_val;
-	  return true;
-	}
-      }
-      return false;
-    }
-    
-    void set_threshold( double _thres) { threshold_ = _thres; }
-
-    void get_ids( std::vector<int>& _ids )
-    {
-      _ids.clear();
-      _ids.reserve( rset_.size());
-      std::set<PairDI>::iterator s_it = rset_.begin();
-      for(; s_it != rset_.end(); ++s_it)
-	_ids.push_back( s_it->second);
-    }
-
-  private:
-
-    double threshold_;
-    double cur_sum_;
-
-    std::set<PairDI> rset_;
-
-    std::set<PairDI> test_;
-  };
-
 private:
 
   void solve_no_rounding( 
@@ -295,11 +230,8 @@ private:
     Vecd&      _rhs, 
     Veci&      _to_round);
 
-  void solve_multiple_rounding( 
-    CSCMatrix& _A, 
-    Vecd&      _x, 
-    Vecd&      _rhs, 
-    Veci&      _to_round );
+  void solve_multiple_rounding(
+      CSCMatrix& _A, Vecd& _x, Vecd& _rhs, const Veci& _to_round);
 
   void solve_iterative(
     CSCMatrix& _A, 
@@ -320,12 +252,8 @@ private:
     Vecd&      _rhs,
     Veci&      _to_round );
 
-
-void update_solution( 
-    CSCMatrix& _A, 
-    Vecd&      _x,