Commit 547ce967 authored by Martin Marinov's avatar Martin Marinov

Merge commit 'f5b324f7' into marinom/merge-from-ReForm

parents db2564f9 f5b324f7
Pipeline #15420 failed with stages
in 1 minute and 15 seconds
......@@ -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_;}
......
......@@ -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);
......
......@@ -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
......
......@@ -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
};
//=============================================================================
......
......@@ -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;