MISolver::solve_multiple_rounding() in O(n*log(n)) (#5)

* Implement the MISolve::mutliple_rounding() method with O(n*log(n)) complexity.
* The CF compute is now ~94% faster for a 1M triangles mesh.
* Clean up the MISolver declaration

Solver/ConstrainedSolver.hh

@@ -258,7 +258,7 @@ public:
   void set_noisy( int _noisy) { noisy_ = _noisy;}
 
   // Get/Set whether the constraint reordering is used (default true)
-  bool& use_constraint_reordering() { return miso_.use_constraint_reordering(); }
+  bool use_constraint_reordering = true;
 
   /// Access the MISolver (e.g. to change settings)
   COMISO::MISolver& misolver() { return miso_;}

Solver/ConstrainedSolverT.cc

@@ -167,7 +167,7 @@ ConstrainedSolver::solve(
 
   DEB_out_if(_show_timings, 1, "Initital dimension: " << nrows << " x " << ncols
     << ", number of constraints: " << ncons << ", number of integer variables: " << _idx_to_round.size()
-    << ", use reordering: " << (use_constraint_reordering() ? "yes" : "no") << "\n")
+    << ", use reordering: " << (use_constraint_reordering ? "yes" : "no") << "\n")
 
     // StopWatch for Timings
     Base::StopWatch sw, sw2; sw.start(); sw2.start();
@@ -177,7 +177,7 @@ ConstrainedSolver::solve(
   std::vector<int> c_elim(ncons);
 
   // apply sparse gauss elimination to make subsequent _conditions independent
-  if (use_constraint_reordering())
+  if (use_constraint_reordering)
     make_constraints_independent_reordering(_constraints, _idx_to_round, c_elim);
   else
     make_constraints_independent(_constraints, _idx_to_round, c_elim);

Solver/IterativeSolverT.cc

@@ -9,6 +9,7 @@
 //== INCLUDES =================================================================
 
 #include "IterativeSolverT.hh"
+#include <Base/Debug/DebOut.hh>
 
 //== NAMESPACES ===============================================================
 
@@ -19,8 +20,8 @@ namespace COMISO
 
 template <class RealT>
 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)
+    const Vector& _rhs, const IndexVector& _idxs, const int _max_iter,
+    const Real& _tolerance)
 {
   if (_max_iter == 0)
     return false;
@@ -28,21 +29,20 @@ bool IterativeSolverT<RealT>::gauss_seidel_local(const Matrix& _A, Vector& _x,
   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();
+  updt_vrbl_indcs_.clear();
+  indx_queue_.clear(); 
 
   for (unsigned int i = 0; i < _idxs.size(); ++i)
-    q.push_back(_idxs[i]);
+    indx_queue_.push_back(_idxs[i]);
 
   int it_count = 0;
 
-  while (!q.empty() && it_count < _max_iter)
+  while (!indx_queue_.empty() && it_count < _max_iter)
   {
     ++it_count;
-    const auto i = q.front();
-    q.pop_front();
-    i_temp.clear();
+    const auto i = indx_queue_.front();
+    indx_queue_.pop_front();
+    indx_temp_.clear();
 
     double res_i = -_rhs[i];
     double x_i_new = _rhs[i];
@@ -56,7 +56,7 @@ bool IterativeSolverT<RealT>::gauss_seidel_local(const Matrix& _A, Vector& _x,
       res_i += (*it) * _x[j];
       x_i_new -= (*it) * _x[j];
       if (j != i)
-        i_temp.push_back(j);
+        indx_temp_.push_back(j);
       else
         diag = *it;
     }
@@ -68,82 +68,13 @@ bool IterativeSolverT<RealT>::gauss_seidel_local(const Matrix& _A, Vector& _x,
     if (fabs(res_i * diag) > _tolerance)
     {
       _x[i] += x_i_new * diag;
-
-      for (unsigned int j = 0; j < i_temp.size(); ++j)
-        q.push_back(i_temp[j]);
+      updt_vrbl_indcs_.push_back(i);
+      for (unsigned int j = 0; j < indx_temp_.size(); ++j)
+        indx_queue_.push_back(indx_temp_[j]);
     }
   }
 
-  return q.empty(); // converged?
-}
-
-//-----------------------------------------------------------------------------
-
-template <class RealT>
-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<Matrix>::const_sub_col_type ColT;
-  typedef typename gmm::linalg_traits<ColT>::const_iterator CIter;
-
-  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]);
-
-  int it_count = 0;
-
-  bool finished = false;
-
-  while (!finished && it_count < _max_iter)
-  {
-    finished = true;
-    std::set<int>::iterator s_it = s.begin();
-    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);
-
-      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)
-      {
-        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)
-      {
-        _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]);
-        }
-      }
-    }
-  }
-
-  // converged?
-  return finished;
+  return indx_queue_.empty(); // converged?
 }
 
 //-----------------------------------------------------------------------------
@@ -152,6 +83,7 @@ template <class RealT>
 bool IterativeSolverT<RealT>::conjugate_gradient(const Matrix& _A, Vector& _x,
     const Vector& _rhs, int& _max_iter, Real& _tolerance)
 {
+  DEB_enter_func;
   Real rho, rho_1(0), a;
 
   // initialize vectors
@@ -179,7 +111,7 @@ bool IterativeSolverT<RealT>::conjugate_gradient(const Matrix& _A, Vector& _x,
   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;
+    DEB_line(11, "iter " << cur_iter << "  res " << res_norm);
 
     if (cur_iter != 0)
     {
@@ -200,7 +132,7 @@ bool IterativeSolverT<RealT>::conjugate_gradient(const Matrix& _A, Vector& _x,
   _max_iter = cur_iter;
   _tolerance = res_norm;
 
-  return (!not_converged);
+  return !not_converged;
 }
 
 //-----------------------------------------------------------------------------
@@ -210,19 +142,11 @@ 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)
-  {
     res = std::max(fabs(_v[i] * _diag[i]), res);
-
-    //     Real cur = fabs(_v[i]*_diag[i]);
-    //     if(cur > res)
-    //       res = cur;
-  }
   return res;
 }
 
-//-----------------------------------------------------------------------------
 
 //=============================================================================
 } // namespace COMISO

Solver/IterativeSolverT.hh

@@ -11,7 +11,6 @@
 
 #include <CoMISo/Utils/gmm.hh>
 #include <deque>
-#include <set>
 
 //== FORWARDDECLARATIONS ======================================================
 
@@ -32,19 +31,21 @@ namespace COMISO
 template <class RealT> class IterativeSolverT
 {
 public:
+  typedef unsigned int uint;
   typedef RealT Real;
   typedef std::vector<Real> Vector;
+  typedef std::vector<uint> IndexVector;
   typedef gmm::csc_matrix<Real> Matrix;
 
-  // local gauss_seidel
+  // local Gauss-Seidel
   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);
+      const IndexVector& _idxs, const int _max_iter, const Real& _tolerance);
 
-  // local gauss_seidel
-  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);
+  // get the indices of any variables updated during the last local Gauss-Seidel
+  const IndexVector& updated_variable_indices() const
+  {
+    return updt_vrbl_indcs_;
+  }
 
   // conjugate gradient
   bool conjugate_gradient(const Matrix& _A, Vector& _x, const Vector& _rhs,
@@ -55,16 +56,16 @@ private:
   Real vect_norm_rel(const Vector& _v, const Vector& _diag) const;
 
 private:
-  // helper for conjugate gradient
+  // context  for Conjugate Gradient
   Vector p_;
   Vector q_;
   Vector r_;
   Vector d_;
 
-  //  helper for local gauss seidel
-  std::vector<unsigned int> i_temp;
-  std::deque<unsigned int> q;
-  std::set<int> s;
+  // context for local Gauss-Seidel
+  IndexVector indx_temp_;
+  std::deque<uint> indx_queue_;
+  IndexVector updt_vrbl_indcs_; // updated variable indices
 };
 
 //=============================================================================

Solver/MISolver.cc

@@ -46,6 +46,16 @@ ILOSTLBEGIN
 #include "EigenLDLTSolver.hh"
 #endif
 
+#if COMISO_SUITESPARSE_AVAILABLE
+#include "CholmodSolver.hh"
+#elif COMISO_EIGEN3_AVAILABLE
+#include "EigenLDLTSolver.hh"
+#else
+#error "MISolver requires Suitesparse or Eigen3 support"
+#endif
+
+#include "IterativeSolverT.hh"
+
 #include <CoMISo/Utils/gmm.hh>
 #include <CoMISo/Utils/Tools.hh>
 
@@ -59,62 +69,65 @@ namespace COMISO
 {
 namespace
 {
+typedef unsigned int uint;
 
-class RoundingQueue : protected std::priority_queue<std::pair<double, int>>
-{
+// gmm Column and ColumnIterator of CSC matrix
+typedef gmm::linalg_traits<MISolver::CSCMatrix>::const_sub_col_type Col;
+typedef gmm::linalg_traits<Col>::const_iterator ColIter;
+
+using ToRoundSet = std::set<std::pair<double, uint>>;
+using ToRoundSetIter = ToRoundSet::iterator;
+
+// extend ToRoundSetIter with a flag to tell if it has been set or not
+class ToRoundSetIterExt
+{ 
 public:
-  RoundingQueue(const double _threshold) : threshold_(_threshold), cur_sum_(0.0)
+  void set(ToRoundSetIter _iter)
   {
-    c.reserve(128);
+    iter_ = _iter;
+    null_ = false;
   }
 
-  void clear()
-  {
-    c.clear();
-    cur_sum_ = 0.0;
-  }
+  void clear() { null_ = true; }
 
-  void add(const int _id, const double _rd_val)
-  {
-    if (empty() || cur_sum_ + _rd_val <= threshold_)
-    { // empty? -> always add one element
-      emplace(_rd_val, _id);
-      cur_sum_ += _rd_val;
-      return;
-    }
-    if (top().first > _rd_val) // replace the last element if worse
-    {
-      cur_sum_ -= top().first;
-      pop();
-      emplace(_rd_val, _id);
-      cur_sum_ += _rd_val;
-    }
-  }
+  bool is_null() const { return null_; }
 
-  void get_ids(std::vector<int>& _ids)
+  const ToRoundSetIter& get() const
   {
-    _ids.clear();
-    _ids.reserve(size());
-    std::sort_heap(c.begin(), c.end());
-    for (auto s_it = c.begin(); s_it != c.end(); ++s_it)
-      _ids.push_back(s_it->second);
+    DEB_error_if(is_null(), "accessing an null iterator");
+    return iter_;
   }
 
+  auto operator->() const { return get().operator->(); }
+  auto operator*() const { return get().operator*(); }
+  operator ToRoundSetIter() { return get(); }
+
 private:
-  const double threshold_;
-  double cur_sum_;
+  ToRoundSetIter iter_;
+  bool null_ = true;
 };
 
-// gmm Column and ColumnIterator of CSC matrix
-typedef gmm::linalg_traits<MISolver::CSCMatrix>::const_sub_col_type Col;
-typedef gmm::linalg_traits<Col>::const_iterator ColIter;
-
 } // namespace
 
+// base class selected based on the available packages
+class MISolver::DirectSolver : public
+#if COMISO_SUITESPARSE_AVAILABLE
+                               CholmodSolver
+#elif COMISO_EIGEN3_AVAILABLE
+                               EigenLDLTSolver
+#else
+#error "MISolver requires Suitesparse or Eigen3 support"
+#endif
+{
+};
+
+class MISolver::IterativeSolver : public IterativeSolverT<double>
+{
+};
+
 // Constructor
 MISolver::MISolver()
-{
-  // default parameters
+{// default parameters
   initial_full_solution_ = true;
   iter_full_solution_ = true;
   final_full_solution_ = true;
@@ -134,22 +147,23 @@ MISolver::MISolver()
 
   gurobi_max_time_ = 60;
 
-  noisy_ = 0;
-  stats_ = true;
-
-  use_constraint_reordering_ = true;
+  direct_solver_ = new DirectSolver;
+  iter_solver_ = new IterativeSolver;
 }
 
+MISolver::~MISolver()
+{
+  delete direct_solver_;
+  delete iter_solver_;
+}
 //-----------------------------------------------------------------------------
 
 void MISolver::solve(
     CSCMatrix& _A, Vecd& _x, Vecd& _rhs, Veci& _to_round, bool _fixed_order)
 {
   DEB_enter_func;
-
-  DEB_out(6, "# integer    variables: "
-                 << _to_round.size() << "\n# continuous variables: "
-                 << _x.size() - _to_round.size() << "\n");
+  DEB_line(2, "integer variables #: " << _to_round.size());
+  DEB_line(2, "continuous variables #: " << _x.size() - _to_round.size());
 
   // nothing to solve?
   if (gmm::mat_ncols(_A) == 0 || gmm::mat_nrows(_A) == 0)
@@ -272,13 +286,16 @@ void MISolver::solve_cplex(CSCMatrix& _A, Vecd& _x, Vecd& _rhs, Veci& _to_round)
 void MISolver::solve_no_rounding(CSCMatrix& _A, Vecd& _x, Vecd& _rhs)
 {
   COMISO_THROW_if(
-      !direct_solver_.calc_system_gmm(_A), UNSPECIFIED_EIGEN_FAILURE);
-  COMISO_THROW_if(!direct_solver_.solve(_x, _rhs), UNSPECIFIED_EIGEN_FAILURE);
+      !direct_solver_->calc_system_gmm(_A), UNSPECIFIED_EIGEN_FAILURE);
+  COMISO_THROW_if(!direct_solver_->solve(_x, _rhs), UNSPECIFIED_EIGEN_FAILURE);
 }
 
 //-----------------------------------------------------------------------------
 
-void MISolver::resolve(Vecd& _x, Vecd& _rhs) { direct_solver_.solve(_x, _rhs); }
+void MISolver::resolve(Vecd& _x, Vecd& _rhs)
+{
+  direct_solver_->solve(_x, _rhs);
+}
 
 //-----------------------------------------------------------------------------
 
@@ -298,8 +315,8 @@ void MISolver::solve_direct_rounding(
   Veci old_idx(_rhs.size());
   for (unsigned int i = 0; i < old_idx.size(); ++i)
     old_idx[i] = i;
-  direct_solver_.calc_system_gmm(_A);
-  direct_solver_.solve(_x, _rhs);
+  direct_solver_->calc_system_gmm(_A);
+  direct_solver_->solve(_x, _rhs);
 
 #ifdef COMISO_MISOLVER_PERFORMANCE_TEST
   // check solver performance (only for testing!!!)
@@ -403,9 +420,9 @@ void MISolver::solve_direct_rounding(
   // final full solution
   if (gmm::mat_ncols(_A) > 0)
   {
-    //    direct_solver_.update_system_gmm(_A);
-    direct_solver_.calc_system_gmm(_A);
-    direct_solver_.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
@@ -450,9 +467,9 @@ void MISolver::solve_iterative(
 
   if (initial_full_solution_)
   {
-    DEB_out_if(noisy_ > 2, 2, "initial full solution\n")
-        direct_solver_.calc_system_gmm(_A);
-    direct_solver_.solve(_x, _rhs);
+    DEB_line(2, "initial full solution");
+    direct_solver_->calc_system_gmm(_A);
+    direct_solver_->solve(_x, _rhs);
 
     cholmod_step_done_ = true;
 
@@ -468,10 +485,8 @@ void MISolver::solve_iterative(
   // loop until solution computed
   for (unsigned int i = 0; i < to_round.size(); ++i)
   {
-    DEB_out_if(noisy_ > 0, 2,
-        "Integer DOF's left: " << to_round.size() - (i + 1)
-                               << " ") DEB_out_if(noisy_ > 1, 1,
-        "residuum_norm: " << COMISO_GMM::residuum_norm(_A, xr, _rhs) << "\n");
+    DEB_line(11, "Integer DOF's left: " << to_round.size() - (i + 1) << " ");
+    DEB_line(11, "residuum_norm: " << COMISO_GMM::residuum_norm(_A, xr, _rhs));
 
     // index to eliminate
     unsigned int i_best = 0;
@@ -525,22 +540,22 @@ void MISolver::solve_iterative(
 
     // 3-stage update of solution w.r.t. roundings
     // local GS / CG / SparseCholesky
-    update_solution(_A, xr, _rhs, neigh_i);
+    update_solution_is_local(_A, xr, _rhs, neigh_i);
   }
 
   // final full solution?
   if (final_full_solution_)
   {
-    DEB_line_if(noisy_ > 2, 2, "final full solution");
+    DEB_line(2, "final full solution");
 
     if (gmm::mat_ncols(_A) > 0)
     {
       if (cholmod_step_done_)
-        direct_solver_.update_system_gmm(_A);
+        direct_solver_->update_system_gmm(_A);
       else
-        direct_solver_.calc_system_gmm(_A);
+        direct_solver_->calc_system_gmm(_A);
 
-      direct_solver_.solve(xr, _rhs);
+      direct_solver_->solve(xr, _rhs);
       ++n_full_;
     }
   }
@@ -550,19 +565,18 @@ void MISolver::solve_iterative(
     _x[old_idx[i]] = xr[i];
 
   // output statistics
-  DEB_out_if(stats_, 2,
-      " *** Statistics of MiSo Solver ***"
-          << "\n Number of CG    iterations  = " << n_cg_
-          << "\n Number of LOCAL iterations  = " << n_local_
-          << "\n Number of FULL  iterations  = " << n_full_
-          << "\n Number of ROUNDING          = " << _to_round.size()
-          << "\n time searching next integer = "
-          << time_search_next_integer / 1000.0 << "s\n\n");
+  DEB_line(2, " *** Statistics of MiSo Solver ***");
+  DEB_line(2, "Number of CG    iterations  = " << n_cg_);
+  DEB_line(2, "Number of LOCAL iterations  = " << n_local_);
+  DEB_line(2, "Number of FULL  iterations  = " << n_full_);
+  DEB_line(2, "Number of ROUNDING          = " << _to_round.size());
+  DEB_line(2, "time searching next integer = "
+                  << time_search_next_integer / 1000.0 << "s");
 }
 
 //-----------------------------------------------------------------------------
 
-void MISolver::update_solution(
+bool MISolver::update_solution_is_local(
     const CSCMatrix& _A, Vecd& _x, const Vecd& _rhs, const Vecui& _neigh_i)
 {
   DEB_enter_func;
@@ -570,29 +584,35 @@ void MISolver::update_solution(
 
   if (max_local_iters_ > 0) // compute new solution
   {
-    DEB_out(6, "use local iteration ");
-
-    int n_its = max_local_iters_;
-    double tolerance = max_local_error_;
-    converged =
-        siter_.gauss_seidel_local(_A, _x, _rhs, _neigh_i, n_its, tolerance);
+    DEB_out(11, "use local iteration ");
+    converged = iter_solver_->gauss_seidel_local(
+        _A, _x, _rhs, _neigh_i, max_local_iters_, max_local_error_);
 
     ++n_local_;
   }
 
-  // conjugate gradient
-  if (!converged && max_cg_iters_ > 0)
+  if (converged)
   {
+    DEB_line(11, "Local iteration converged");
+    return true;
+  }
+
+  bool only_local_updates = true; // set to false if we do any global updates
+
+  DEB_line(6, "Local iterations failed to converge, needs global update");
+  if (max_cg_iters_ > 0) 
+  {// conjugate gradient iterations
     DEB_out(6, ", cg ");
 
     int max_cg_iters = max_cg_iters_;
     double tolerance = max_cg_error_;
     converged =
-        siter_.conjugate_gradient(_A, _x, _rhs, max_cg_iters, tolerance);
+        iter_solver_->conjugate_gradient(_A, _x, _rhs, max_cg_iters, tolerance);
 
     DEB_out(6, "( converged " << converged << " "
                               << " iters " << max_cg_iters << " "
                               << " res_norm " << tolerance << "\n");
+    only_local_updates = false;
     ++n_cg_;
   }
 
@@ -600,19 +620,21 @@ void MISolver::update_solution(
   {
     DEB_out(6, ", full ");
     if (cholmod_step_done_)
-      direct_solver_.update_system_gmm(_A);
+      direct_solver_->update_system_gmm(_A);
     else
     {
-      direct_solver_.calc_system_gmm(_A);
+      direct_solver_->calc_system_gmm(_A);
       cholmod_step_done_ = true;
     }
     // const_cast<> to work around broken const correctness in Eigen
-    direct_solver_.solve(_x, const_cast<Vecd&>(_rhs)); 
+    direct_solver_->solve(_x, const_cast<Vecd&>(_rhs)); 
 
+    only_local_updates = false;
     ++n_full_;
   }
 
   DEB_line(6, "");
+  return only_local_updates;
 }
 
 //-----------------------------------------------------------------------------
@@ -629,117 +651,112 @@ void MISolver::solve_multiple_rounding(
   // reset cholmod step flag
   cholmod_step_done_ = false;
 
-  Veci to_round(_to_round);
-  { // copy to round vector and make it unique
-    std::sort(to_round.begin(), to_round.end());
-    const auto last_unique = std::unique(to_round.begin(), to_round.end());
-    to_round.resize(last_unique - to_round.begin());
-  }
-
   if (initial_full_solution_)
   {
-    DEB_out_if(noisy_ > 2, 2, "initial full solution\n");
+    DEB_line(2, "initial full solution");
     // TODO: we can throw more specific outcomes in the body of the
     // functions below
     COMISO_THROW_if(
-        !direct_solver_.calc_system_gmm(_A), UNSPECIFIED_EIGEN_FAILURE);
-    COMISO_THROW_if(!direct_solver_.solve(_x, _rhs), UNSPECIFIED_EIGEN_FAILURE);
+        !direct_solver_->calc_system_gmm(_A), UNSPECIFIED_EIGEN_FAILURE);
+    COMISO_THROW_if(!direct_solver_->solve(_x, _rhs), UNSPECIFIED_EIGEN_FAILURE);
 
     cholmod_step_done_ = true;
 
     ++n_full_;