diff --git a/src/math/lp/CMakeLists.txt b/src/math/lp/CMakeLists.txt index 5ccc023ab..ef363509f 100644 --- a/src/math/lp/CMakeLists.txt +++ b/src/math/lp/CMakeLists.txt @@ -41,7 +41,6 @@ z3_add_component(lp permutation_matrix.cpp random_updater.cpp row_eta_matrix.cpp - square_sparse_matrix.cpp static_matrix.cpp COMPONENT_DEPENDENCIES util diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp index f890b27fa..9bc850848 100644 --- a/src/math/lp/lar_solver.cpp +++ b/src/math/lp/lar_solver.cpp @@ -783,13 +783,6 @@ namespace lp { return r; } - - template - void lar_solver::add_last_rows_to_lu(lp_primal_core_solver& s) { - lp_assert(false); - - } - bool lar_solver::x_is_correct() const { if (m_mpq_lar_core_solver.m_r_x.size() != A_r().column_count()) { return false; diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h index be255dd65..b925b5ef0 100644 --- a/src/math/lp/lar_solver.h +++ b/src/math/lp/lar_solver.h @@ -222,8 +222,6 @@ class lar_solver : public column_namer { void update_x_and_inf_costs_for_columns_with_changed_bounds_tableau(); void solve_with_core_solver(); numeric_pair get_basic_var_value_from_row(unsigned i); - template - void add_last_rows_to_lu(lp_primal_core_solver & s); bool x_is_correct() const; void fill_last_row_of_A_r(static_matrix> & A, const lar_term * ls); template diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h index 2c7360712..10de624f5 100644 --- a/src/math/lp/lp_primal_core_solver.h +++ b/src/math/lp/lp_primal_core_solver.h @@ -926,7 +926,7 @@ public: m_lower_bounds_dummy.resize(A.column_count(), zero_of_type()); m_enter_price_eps = numeric_traits::precise() ? numeric_traits::zero() : T(1e-5); #ifdef Z3DEBUG - // check_correctness(); + lp_assert(false); #endif } diff --git a/src/math/lp/square_sparse_matrix.cpp b/src/math/lp/square_sparse_matrix.cpp deleted file mode 100644 index 3ec88f47d..000000000 --- a/src/math/lp/square_sparse_matrix.cpp +++ /dev/null @@ -1,118 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#include -#include "util/vector.h" -#include "math/lp/lp_settings.h" -#include "math/lp/square_sparse_matrix_def.h" -#include "math/lp/dense_matrix.h" -namespace lp { -template double square_sparse_matrix::dot_product_with_row(unsigned int, vector const&) const; -template void square_sparse_matrix::add_new_element(unsigned int, unsigned int, const double&); -template void square_sparse_matrix::divide_row_by_constant(unsigned int, const double&, lp_settings&); -template bool square_sparse_matrix::fill_eta_matrix(unsigned int, eta_matrix**); -template const double & square_sparse_matrix::get(unsigned int, unsigned int) const; -template unsigned square_sparse_matrix::get_number_of_nonzeroes() const; -template bool square_sparse_matrix::get_pivot_for_column(unsigned int&, unsigned int&, int, unsigned int); -template unsigned square_sparse_matrix::lowest_row_in_column(unsigned int); -template bool square_sparse_matrix::pivot_row_to_row(unsigned int, const double&, unsigned int, lp_settings&); -template bool square_sparse_matrix::pivot_with_eta(unsigned int, eta_matrix*, lp_settings&); -template void square_sparse_matrix::prepare_for_factorization(); -template void square_sparse_matrix::remove_element(vector >&, indexed_value&); -template void square_sparse_matrix::replace_column(unsigned int, indexed_vector&, lp_settings&); -template void square_sparse_matrix::set(unsigned int, unsigned int, double); -template void square_sparse_matrix::set_max_in_row(vector >&); -template bool square_sparse_matrix::set_row_from_work_vector_and_clean_work_vector_not_adjusted(unsigned int, indexed_vector&, lp_settings&); -template bool square_sparse_matrix::shorten_active_matrix(unsigned int, eta_matrix*); -template void square_sparse_matrix::solve_y_U(vector&) const; -template square_sparse_matrix::square_sparse_matrix(unsigned int, unsigned); -template void square_sparse_matrix::add_new_element(unsigned int, unsigned int, const mpq&); -template void square_sparse_matrix::divide_row_by_constant(unsigned int, const mpq&, lp_settings&); -template bool square_sparse_matrix::fill_eta_matrix(unsigned int, eta_matrix**); -template mpq const & square_sparse_matrix::get(unsigned int, unsigned int) const; -template unsigned square_sparse_matrix::get_number_of_nonzeroes() const; -template bool square_sparse_matrix::get_pivot_for_column(unsigned int&, unsigned int&, int, unsigned int); -template unsigned square_sparse_matrix::lowest_row_in_column(unsigned int); -template bool square_sparse_matrix::pivot_with_eta(unsigned int, eta_matrix*, lp_settings&); -template void square_sparse_matrix::prepare_for_factorization(); -template void square_sparse_matrix::remove_element(vector> &, indexed_value&); -template void square_sparse_matrix::replace_column(unsigned int, indexed_vector&, lp_settings&); -template void square_sparse_matrix::set_max_in_row(vector>&); -template bool square_sparse_matrix::set_row_from_work_vector_and_clean_work_vector_not_adjusted(unsigned int, indexed_vector&, lp_settings&); -template bool square_sparse_matrix::shorten_active_matrix(unsigned int, eta_matrix*); -template void square_sparse_matrix::solve_y_U(vector&) const; -template void square_sparse_matrix>::add_new_element(unsigned int, unsigned int, const mpq&); -template void square_sparse_matrix>::divide_row_by_constant(unsigned int, const mpq&, lp_settings&); -template bool square_sparse_matrix>::fill_eta_matrix(unsigned int, eta_matrix >**); -template const mpq & square_sparse_matrix>::get(unsigned int, unsigned int) const; -template unsigned square_sparse_matrix>::get_number_of_nonzeroes() const; -template bool square_sparse_matrix>::get_pivot_for_column(unsigned int&, unsigned int&, int, unsigned int); -template unsigned square_sparse_matrix>::lowest_row_in_column(unsigned int); -template bool square_sparse_matrix>::pivot_with_eta(unsigned int, eta_matrix >*, lp_settings&); -template void square_sparse_matrix>::prepare_for_factorization(); -template void square_sparse_matrix>::remove_element(vector>&, indexed_value&); -template void square_sparse_matrix>::replace_column(unsigned int, indexed_vector&, lp_settings&); -template void square_sparse_matrix>::set_max_in_row(vector>&); -template bool square_sparse_matrix>::set_row_from_work_vector_and_clean_work_vector_not_adjusted(unsigned int, indexed_vector&, lp_settings&); -template bool square_sparse_matrix>::shorten_active_matrix(unsigned int, eta_matrix >*); -template void square_sparse_matrix>::solve_y_U(vector&) const; -template void square_sparse_matrix::double_solve_U_y(indexed_vector&, const lp_settings &); -template void square_sparse_matrix::double_solve_U_y(indexed_vector&, const lp_settings&); -template void square_sparse_matrix>::double_solve_U_y(indexed_vector&, const lp_settings&); -template void square_sparse_matrix >::double_solve_U_y >(indexed_vector>&, const lp_settings&); -template void square_sparse_matrix::solve_U_y_indexed_only(indexed_vector&, const lp_settings&, vector &); -template void square_sparse_matrix::solve_U_y_indexed_only(indexed_vector&, const lp_settings &, vector &); -#ifdef Z3DEBUG -template bool square_sparse_matrix::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const; -template bool square_sparse_matrix::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const; -template bool square_sparse_matrix >::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const; -#endif - -template void square_sparse_matrix >::solve_U_y_indexed_only(indexed_vector&, const lp_settings &, vector &); -template void square_sparse_matrix::solve_U_y(vector&); -template void square_sparse_matrix::double_solve_U_y(vector&); -template void square_sparse_matrix::solve_U_y(vector&); -template void square_sparse_matrix::double_solve_U_y(vector&); -template void square_sparse_matrix >::solve_U_y >(vector >&); -template void square_sparse_matrix >::double_solve_U_y >(vector >&); -template void square_sparse_matrix::find_error_in_solution_U_y_indexed(indexed_vector&, indexed_vector&, const vector &); -template double square_sparse_matrix::dot_product_with_row(unsigned int, indexed_vector const&) const; -template void square_sparse_matrix::find_error_in_solution_U_y_indexed(indexed_vector&, indexed_vector&, const vector &); -template mpq square_sparse_matrix::dot_product_with_row(unsigned int, indexed_vector const&) const; -template void square_sparse_matrix >::find_error_in_solution_U_y_indexed(indexed_vector&, indexed_vector&, const vector &); -template mpq square_sparse_matrix >::dot_product_with_row(unsigned int, indexed_vector const&) const; -template void square_sparse_matrix >::find_error_in_solution_U_y_indexed >(indexed_vector >&, indexed_vector >&, const vector &); -template numeric_pair square_sparse_matrix >::dot_product_with_row >(unsigned int, indexed_vector > const&) const; -template void square_sparse_matrix::extend_and_sort_active_rows(vector const&, vector&); - -template void square_sparse_matrix >::extend_and_sort_active_rows(vector const&, vector&); - -template void square_sparse_matrix >::solve_U_y(vector&); -template void square_sparse_matrix >::double_solve_U_y(vector&); -template void square_sparse_matrix< mpq,numeric_pair< mpq> >::set(unsigned int,unsigned int, mpq); -template void square_sparse_matrix::solve_y_U_indexed(indexed_vector&, const lp_settings & ); -template void square_sparse_matrix::solve_y_U_indexed(indexed_vector&, const lp_settings &); -template void square_sparse_matrix >::solve_y_U_indexed(indexed_vector&, const lp_settings &); - -template square_sparse_matrix::square_sparse_matrix(static_matrix const&, vector&); -template square_sparse_matrix::square_sparse_matrix (static_matrix const&, vector&); -template square_sparse_matrix >::square_sparse_matrix(static_matrix > const&, vector&); -} -template void lp::square_sparse_matrix::copy_from_input_on_basis >(lp::static_matrix const&, vector&); -template void lp::square_sparse_matrix::copy_from_input_on_basis >(lp::static_matrix const&, vector&); diff --git a/src/math/lp/square_sparse_matrix.h b/src/math/lp/square_sparse_matrix.h deleted file mode 100644 index 5637af99f..000000000 --- a/src/math/lp/square_sparse_matrix.h +++ /dev/null @@ -1,433 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ - -#pragma once -#include "util/vector.h" -#include "math/lp/permutation_matrix.h" -#include "math/lp/static_matrix.h" -#include -#include -#include -#include -#include -#include "math/lp/indexed_value.h" -#include "math/lp/indexed_vector.h" -#include -#include "math/lp/lp_settings.h" -#include "math/lp/eta_matrix.h" -#include "math/lp/binary_heap_upair_queue.h" -#include "math/lp/numeric_pair.h" -#include "math/lp/u_set.h" -namespace lp { -// it is a square matrix -template -class square_sparse_matrix - : public matrix -{ - struct col_header { - unsigned m_shortened_markovitz; - vector> m_values; // the actual column values - - col_header(): m_shortened_markovitz(0) {} - - void shorten_markovich_by_one() { - m_shortened_markovitz++; - } - - void zero_shortened_markovitz() { - m_shortened_markovitz = 0; - } - }; - - unsigned m_n_of_active_elems; - binary_heap_upair_queue m_pivot_queue; -public: - vector>> m_rows; - vector m_columns; - permutation_matrix m_row_permutation; - permutation_matrix m_column_permutation; - // m_work_pivot_vector[j] = offset of elementh of j-th column in the row we are pivoting to - // if the column is not present then m_work_pivot_vector[j] is -1 - vector m_work_pivot_vector; - vector m_processed; - unsigned get_n_of_active_elems() const { return m_n_of_active_elems; } - -#ifdef Z3DEBUG - // dense_matrix m_dense; -#endif - /* - the rule is: row i is mapped to m_row_permutation[i] and - column j is mapped to m_column_permutation.apply_reverse(j) - */ - - unsigned adjust_row(unsigned row) const{ - return m_row_permutation[row]; - } - - unsigned adjust_column(unsigned col) const{ - return m_column_permutation.apply_reverse(col); - } - - unsigned adjust_row_inverse(unsigned row) const{ - return m_row_permutation.apply_reverse(row); - } - - unsigned adjust_column_inverse(unsigned col) const{ - return m_column_permutation[col]; - } - - template - void copy_column_from_input(unsigned input_column, const M& A, unsigned j); - template - void copy_column_from_input_with_possible_zeros(const M& A, unsigned j); - - template - void copy_from_input(const M& A); - template - void copy_from_input_on_basis(const M& A, vector & basis); - -public: - - // constructors - template - square_sparse_matrix(const M &A, vector& basis); - - template - square_sparse_matrix(const M &A); - - square_sparse_matrix(unsigned dim, unsigned); // the second parameter is needed to distinguish this - // constructor from the one above - - - - class ref_matrix_element { - square_sparse_matrix & m_matrix; - unsigned m_row; - unsigned m_col; - public: - ref_matrix_element(square_sparse_matrix & m, unsigned row, unsigned col):m_matrix(m), m_row(row), m_col(col) {} - ref_matrix_element & operator=(T const & v) { m_matrix.set( m_row, m_col, v); return *this; } - ref_matrix_element & operator=(ref_matrix_element const & v) { m_matrix.set(m_row, m_col, v.m_matrix.get(v.m_row, v.m_col)); return *this; } - operator T () const { return m_matrix.get(m_row, m_col); } - }; - - class ref_row { - square_sparse_matrix & m_matrix; - unsigned m_row; - public: - ref_row(square_sparse_matrix & m, unsigned row) : m_matrix(m), m_row(row) {} - ref_matrix_element operator[](unsigned col) const { return ref_matrix_element(m_matrix, m_row, col); } - }; - - void set_with_no_adjusting_for_row(unsigned row, unsigned col, T val); - void set_with_no_adjusting_for_col(unsigned row, unsigned col, T val); - - void set_with_no_adjusting(unsigned row, unsigned col, T val); - - void set(unsigned row, unsigned col, T val); - - T const & get_not_adjusted(unsigned row, unsigned col) const; - T const & get(unsigned row, unsigned col) const; - - ref_row operator[](unsigned row) { return ref_row(*this, row); } - - ref_matrix_element operator()(unsigned row, unsigned col) { return ref_matrix_element(*this, row, col); } - - T operator() (unsigned row, unsigned col) const { return get(row, col); } - - vector> & get_row_values(unsigned row) { - return m_rows[row]; - } - - vector> const & get_row_values(unsigned row) const { - return m_rows[row]; - } - - vector> & get_column_values(unsigned col) { - return m_columns[col].m_values; - } - - vector> const & get_column_values(unsigned col) const { - return m_columns[col].m_values; - } - - unsigned dimension() const {return static_cast(m_row_permutation.size());} - - unsigned row_count() const override {return dimension();} - unsigned column_count() const override {return dimension();} - - void init_row_headers(); - - void init_column_headers(); - - unsigned lowest_row_in_column(unsigned j); - - indexed_value & column_iv_other(indexed_value & iv) { - return m_rows[iv.m_index][iv.m_other]; - } - - indexed_value & row_iv_other(indexed_value & iv) { - return m_columns[iv.m_index].m_values[iv.m_other]; - } - - void remove_element(vector> & row_vals, unsigned row_offset, vector> & column_vals, unsigned column_offset); - - void remove_element(vector> & row_chunk, indexed_value & row_el_iv); - - void put_max_index_to_0(vector> & row_vals, unsigned max_index); - - void set_max_in_row(unsigned row) { - set_max_in_row(m_rows[row]); - } - - - void set_max_in_row(vector> & row_vals); - - bool pivot_with_eta(unsigned i, eta_matrix *eta_matrix, lp_settings & settings); - - void scan_row_to_work_vector_and_remove_pivot_column(unsigned row, unsigned pivot_column); - - // This method pivots row i to row i0 by muliplying row i by - // alpha and adding it to row i0. - // After pivoting the row i0 has a max abs value set correctly at the beginning of m_start, - // Returns false if the resulting row is all zeroes, and true otherwise - bool pivot_row_to_row(unsigned i, const T& alpha, unsigned i0, lp_settings & settings ); - - // set the max val as well - // returns false if the resulting row is all zeroes, and true otherwise - bool set_row_from_work_vector_and_clean_work_vector_not_adjusted(unsigned i0, indexed_vector & work_vec, - lp_settings & settings); - - - // set the max val as well - // returns false if the resulting row is all zeroes, and true otherwise - bool set_row_from_work_vector_and_clean_work_vector(unsigned i0); - - void remove_zero_elements_and_set_data_on_existing_elements(unsigned row); - - // work_vec here has not adjusted column indices - void remove_zero_elements_and_set_data_on_existing_elements_not_adjusted(unsigned row, indexed_vector & work_vec, lp_settings & settings); - - void multiply_from_right(permutation_matrix& p) { - // m_dense = m_dense * p; - m_column_permutation.multiply_by_permutation_from_right(p); - // lp_assert(*this == m_dense); - } - - void multiply_from_left(permutation_matrix& p) { - // m_dense = p * m_dense; - m_row_permutation.multiply_by_permutation_from_left(p); - // lp_assert(*this == m_dense); - } - - void multiply_from_left_with_reverse(permutation_matrix& p) { - // m_dense = p * m_dense; - m_row_permutation.multiply_by_permutation_reverse_from_left(p); - // lp_assert(*this == m_dense); - } - - // adding delta columns at the end of the matrix - void add_columns_at_the_end(unsigned delta); - - void delete_column(int i); - - void swap_columns(unsigned a, unsigned b) { - m_column_permutation.transpose_from_left(a, b); - } - - void swap_rows(unsigned a, unsigned b) { - m_row_permutation.transpose_from_right(a, b); - // m_dense.swap_rows(a, b); - // lp_assert(*this == m_dense); - } - - void divide_row_by_constant(unsigned i, const T & t, lp_settings & settings); - - bool close(T a, T b) { - return // (numeric_traits::precise() && numeric_traits::is_zero(a - b)) - // || - fabs(numeric_traits::get_double(a - b)) < 0.0000001; - } - - // solving x * this = y, and putting the answer into y - // the matrix here has to be upper triangular - void solve_y_U(vector & y) const; - - // solving x * this = y, and putting the answer into y - // the matrix here has to be upper triangular - void solve_y_U_indexed(indexed_vector & y, const lp_settings &); - - // fills the indices for such that y[i] can be not a zero - // sort them so the smaller indices come first - void fill_reachable_indices(std::set & rset, T *y); - - template - void find_error_in_solution_U_y(vector& y_orig, vector & y); - - template - void find_error_in_solution_U_y_indexed(indexed_vector& y_orig, indexed_vector & y, const vector& sorted_active_rows); - - template - void add_delta_to_solution(const vector& del, vector & y); - - template - void add_delta_to_solution(const indexed_vector& del, indexed_vector & y); - - template - void double_solve_U_y(indexed_vector& y, const lp_settings & settings); - - template - void double_solve_U_y(vector& y); - // solving this * x = y, and putting the answer into y - // the matrix here has to be upper triangular - template - void solve_U_y(vector & y); - // solving this * x = y, and putting the answer into y - // the matrix here has to be upper triangular - template - void solve_U_y_indexed_only(indexed_vector & y, const lp_settings&, vector & sorted_active_rows ); - - T get_elem(unsigned i, unsigned j) const override { return get(i, j); } - unsigned get_number_of_rows() const { return dimension(); } - unsigned get_number_of_columns() const { return dimension(); } - void set_number_of_rows(unsigned /*m*/) override { } - void set_number_of_columns(unsigned /*n*/) override { } - template - L dot_product_with_row (unsigned row, const vector & y) const; - - template - L dot_product_with_row (unsigned row, const indexed_vector & y) const; - - unsigned get_number_of_nonzeroes() const; - - bool get_non_zero_column_in_row(unsigned i, unsigned *j) const; - - void remove_element_that_is_not_in_w(vector> & column_vals, indexed_value & col_el_iv); - - - // w contains the new column - // the old column inside of the matrix has not been changed yet - void remove_elements_that_are_not_in_w_and_update_common_elements(unsigned column_to_replace, indexed_vector & w); - - void add_new_element(unsigned row, unsigned col, const T& val); - - // w contains the "rest" of the new column; all common elements of w and the old column has been zeroed - // the old column inside of the matrix has not been changed yet - void add_new_elements_of_w_and_clear_w(unsigned column_to_replace, indexed_vector & w, lp_settings & settings); - - void replace_column(unsigned column_to_replace, indexed_vector & w, lp_settings &settings); - - unsigned pivot_score(unsigned i, unsigned j); - - void enqueue_domain_into_pivot_queue(); - - void set_max_in_rows(); - - void zero_shortened_markovitz_numbers(); - - void prepare_for_factorization(); - - void recover_pivot_queue(vector & rejected_pivots); - - int elem_is_too_small(unsigned i, unsigned j, int c_partial_pivoting); - - bool remove_row_from_active_pivots_and_shorten_columns(unsigned row); - - void remove_pivot_column(unsigned row); - - void update_active_pivots(unsigned row); - - bool shorten_active_matrix(unsigned row, eta_matrix *eta_matrix); - - unsigned pivot_score_without_shortened_counters(unsigned i, unsigned j, unsigned k); -#ifdef Z3DEBUG - bool can_improve_score_for_row(unsigned row, unsigned score, T const & c_partial_pivoting, unsigned k); - bool really_best_pivot(unsigned i, unsigned j, T const & c_partial_pivoting, unsigned k); - void print_active_matrix(unsigned k, std::ostream & out); -#endif - bool pivot_queue_is_correct_for_row(unsigned i, unsigned k); - - bool pivot_queue_is_correct_after_pivoting(int k); - - bool get_pivot_for_column(unsigned &i, unsigned &j, int c_partial_pivoting, unsigned k); - - bool elem_is_too_small(vector> & row_chunk, indexed_value & iv, int c_partial_pivoting); - - unsigned number_of_non_zeroes_in_row(unsigned row) const { - return static_cast(m_rows[row].size()); - } - - unsigned number_of_non_zeroes_in_column(unsigned col) const { - return m_columns[col].m_values.size(); - } - - bool shorten_columns_by_pivot_row(unsigned i, unsigned pivot_column); - - bool col_is_active(unsigned j, unsigned pivot) { - return adjust_column_inverse(j) > pivot; - } - - bool row_is_active(unsigned i, unsigned pivot) { - return adjust_row_inverse(i) > pivot; - } - - bool fill_eta_matrix(unsigned j, eta_matrix ** eta); -#ifdef Z3DEBUG - bool is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings & settings) const; - - bool is_upper_triangular_until(unsigned k) const; - void check_column_vs_rows(unsigned col); - - void check_row_vs_columns(unsigned row); - - void check_rows_vs_columns(); - - void check_columns_vs_rows(); - - void check_matrix(); -#endif - void create_graph_G(const vector & active_rows, vector & sorted_active_rows); - void process_column_recursively(unsigned i, vector & sorted_rows); - void extend_and_sort_active_rows(const vector & active_rows, vector & sorted_active_rows); - void process_index_recursively_for_y_U(unsigned j, vector & sorted_rows); - void resize(unsigned new_dim) { - unsigned old_dim = dimension(); - lp_assert(new_dim >= old_dim); - for (unsigned j = old_dim; j < new_dim; j++) { - m_rows.push_back(vector>()); - m_columns.push_back(col_header()); - } - m_pivot_queue.resize(new_dim); - m_row_permutation.resize(new_dim); - m_column_permutation.resize(new_dim); - m_work_pivot_vector.resize(new_dim); - m_processed.resize(new_dim); - for (unsigned j = old_dim; j < new_dim; j++) { - add_new_element(j, j, numeric_traits::one()); - } - } -#ifdef Z3DEBUG -vector get_full_row(unsigned i) const; -#endif - unsigned pivot_queue_size() const { return m_pivot_queue.size(); } -}; -}; - - diff --git a/src/math/lp/square_sparse_matrix_def.h b/src/math/lp/square_sparse_matrix_def.h deleted file mode 100644 index 192634627..000000000 --- a/src/math/lp/square_sparse_matrix_def.h +++ /dev/null @@ -1,1303 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once - -#include "util/vector.h" -#include "math/lp/square_sparse_matrix.h" -#include "math/lp/dense_matrix.h" - -#include -#include -namespace lp { -template -template -void square_sparse_matrix::copy_column_from_input(unsigned input_column, const M& A, unsigned j) { - vector> & new_column_vector = m_columns[j].m_values; - for (auto c : A.column(input_column)) { - unsigned col_offset = static_cast(new_column_vector.size()); - vector> & row_vector = m_rows[c.var()]; - unsigned row_offset = static_cast(row_vector.size()); - new_column_vector.push_back(indexed_value(c.coeff(), c.var(), row_offset)); - row_vector.push_back(indexed_value(c.coeff(), j, col_offset)); - m_n_of_active_elems++; - } -} - -template -template -void square_sparse_matrix::copy_column_from_input_with_possible_zeros(const M& A, unsigned j) { - vector> & new_column_vector = m_columns[j].m_values; - for (auto c : A.column(j)) { - if (is_zero(c.coeff())) - continue; - unsigned col_offset = static_cast(new_column_vector.size()); - vector> & row_vector = m_rows[c.var()]; - unsigned row_offset = static_cast(row_vector.size()); - new_column_vector.push_back(indexed_value(c.coeff(), c.var(), row_offset)); - row_vector.push_back(indexed_value(c.coeff(), j, col_offset)); - m_n_of_active_elems++; - } -} - -template -template -void square_sparse_matrix::copy_from_input_on_basis(const M & A, vector & basis) { - unsigned m = A.row_count(); - for (unsigned j = m; j-- > 0;) { - copy_column_from_input(basis[j], A, j); - } -} -template -template -void square_sparse_matrix::copy_from_input(const M & A) { - unsigned m = A.row_count(); - for (unsigned j = m; j-- > 0;) { - copy_column_from_input_with_possible_zeros(A, j); - } -} - -// constructor that copies columns of the basis from A -template -template -square_sparse_matrix::square_sparse_matrix(const M &A, vector & basis) : - m_n_of_active_elems(0), - m_pivot_queue(A.row_count()), - m_row_permutation(A.row_count()), - m_column_permutation(A.row_count()), - m_work_pivot_vector(A.row_count(), -1), - m_processed(A.row_count()) { - init_row_headers(); - init_column_headers(); - copy_from_input_on_basis(A, basis); -} - -template -template -square_sparse_matrix::square_sparse_matrix(const M &A) : - m_n_of_active_elems(0), - m_pivot_queue(A.row_count()), - m_row_permutation(A.row_count()), - m_column_permutation(A.row_count()), - m_work_pivot_vector(A.row_count(), -1), - m_processed(A.row_count()) { - init_row_headers(); - init_column_headers(); - copy_from_input(A); -} - - -template -void square_sparse_matrix::set_with_no_adjusting_for_row(unsigned row, unsigned col, T val) { // should not be used in efficient code - vector> & row_vec = m_rows[row]; - for (auto & iv : row_vec) { - if (iv.m_index == col) { - iv.set_value(val); - return; - } - } - // have not found the column between the indices - row_vec.push_back(indexed_value(val, col, -1)); -} - -template -void square_sparse_matrix::set_with_no_adjusting_for_col(unsigned row, unsigned col, T val) { // should not be used in efficient code - vector> & col_vec = m_columns[col].m_values; - for (auto & iv : col_vec) { - if (iv.m_index == row) { - iv.set_value(val); - return; - } - } - // have not found the column between the indices - col_vec.push_back(indexed_value(val, row, -1)); -} - - -template -void square_sparse_matrix::set_with_no_adjusting(unsigned row, unsigned col, T val) { // should not be used in efficient code - set_with_no_adjusting_for_row(row, col, val); - set_with_no_adjusting_for_col(row, col, val); -} - -template -void square_sparse_matrix::set(unsigned row, unsigned col, T val) { // should not be used in efficient code - lp_assert(row < dimension() && col < dimension()); - // m_dense.set_elem(row, col, val); - row = adjust_row(row); - col = adjust_column(col); - set_with_no_adjusting(row, col, val); - // lp_assert(*this == m_dense); -} - -template -T const & square_sparse_matrix::get_not_adjusted(unsigned row, unsigned col) const { - for (indexed_value const & iv : m_rows[row]) { - if (iv.m_index == col) { - return iv.m_value; - } - } - return numeric_traits::zero(); -} - -template -T const & square_sparse_matrix::get(unsigned row, unsigned col) const { // should not be used in efficient code - row = adjust_row(row); - auto & row_chunk = m_rows[row]; - col = adjust_column(col); - for (indexed_value const & iv : row_chunk) { - if (iv.m_index == col) { - return iv.m_value; - } - } - return numeric_traits::zero(); -} - -// constructor creating a zero matrix of dim*dim -template -square_sparse_matrix::square_sparse_matrix(unsigned dim, unsigned ) : - m_pivot_queue(dim), // dim will be the initial size of the queue - m_row_permutation(dim), - m_column_permutation(dim), - m_work_pivot_vector(dim, -1), - m_processed(dim) { - init_row_headers(); - init_column_headers(); - } - -template -void square_sparse_matrix::init_row_headers() { - for (unsigned l = 0; l < m_row_permutation.size(); l++) { - m_rows.push_back(vector>()); - } -} - -template -void square_sparse_matrix::init_column_headers() { // we always have only square square_sparse_matrix - for (unsigned l = 0; l < m_row_permutation.size(); l++) { - m_columns.push_back(col_header()); - } -} - -template -unsigned square_sparse_matrix::lowest_row_in_column(unsigned j) { - auto & mc = get_column_values(adjust_column(j)); - unsigned ret = 0; - for (auto & iv : mc) { - unsigned row = adjust_row_inverse(iv.m_index); - if (row > ret) { - ret = row; - } - } - return ret; -} - -template -void square_sparse_matrix::remove_element(vector> & row_vals, unsigned row_offset, vector> & column_vals, unsigned column_offset) { - if (column_offset != column_vals.size() - 1) { - auto & column_iv = column_vals[column_offset] = column_vals.back(); // copy from the tail - column_iv_other(column_iv).m_other = column_offset; - if (row_offset != row_vals.size() - 1) { - auto & row_iv = row_vals[row_offset] = row_vals.back(); // copy from the tail - row_iv_other(row_iv).m_other = row_offset; - } - } else if (row_offset != row_vals.size() - 1) { - auto & row_iv = row_vals[row_offset] = row_vals.back(); // copy from the tail - row_iv_other(row_iv).m_other = row_offset; - } - // do nothing - just decrease the sizes - column_vals.pop_back(); - row_vals.pop_back(); - m_n_of_active_elems--; // the value is correct only when refactoring -} - -template -void square_sparse_matrix::remove_element(vector> & row_chunk, indexed_value & row_el_iv) { - auto & column_chunk = get_column_values(row_el_iv.m_index); - indexed_value & col_el_iv = column_chunk[row_el_iv.m_other]; - remove_element(row_chunk, col_el_iv.m_other, column_chunk, row_el_iv.m_other); -} - -template -void square_sparse_matrix::put_max_index_to_0(vector> & row_vals, unsigned max_index) { - if (max_index == 0) return; - indexed_value * max_iv = & row_vals[max_index]; - indexed_value * start_iv = & row_vals[0]; - // update the "other" columns elements which are bound to the start_iv and max_iv - m_columns[max_iv->m_index].m_values[max_iv->m_other].m_other = 0; - m_columns[start_iv->m_index].m_values[start_iv->m_other].m_other = max_index; - - // swap the elements - indexed_value t = * max_iv; - * max_iv = * start_iv; - * start_iv = t; -} - -template -void square_sparse_matrix::set_max_in_row(vector> & row_vals) { - if (row_vals.empty()) - return; - T max_val = abs(row_vals[0].m_value); - unsigned max_index = 0; - for (unsigned i = 1; i < row_vals.size(); i++) { - T iabs = abs(row_vals[i].m_value); - if (iabs > max_val) { - max_val = iabs; - max_index = i; - } - } - put_max_index_to_0(row_vals, max_index); -} - -template -bool square_sparse_matrix::pivot_with_eta(unsigned i, eta_matrix *eta_matrix, lp_settings & settings) { - const T& pivot = eta_matrix->get_diagonal_element(); - for (auto & it : eta_matrix->m_column_vector.m_data) { - if (!pivot_row_to_row(i, it.second, it.first, settings)) { - return false; - } - } - - divide_row_by_constant(i, pivot, settings); - if (!shorten_active_matrix(i, eta_matrix)) { - return false; - } - - return true; -} - -// returns the offset of the pivot column in the row -template -void square_sparse_matrix::scan_row_to_work_vector_and_remove_pivot_column(unsigned row, unsigned pivot_column) { - auto & rvals = m_rows[row]; - unsigned size = rvals.size(); - for (unsigned j = 0; j < size; j++) { - auto & iv = rvals[j]; - if (iv.m_index != pivot_column) { - m_work_pivot_vector[iv.m_index] = j; - } else { - remove_element(rvals, iv); - j--; - size--; - } - } -} - -#ifdef Z3DEBUG -template -vector square_sparse_matrix::get_full_row(unsigned i) const { - vector r; - for (unsigned j = 0; j < column_count(); j++) - r.push_back(get(i, j)); - return r; -} -#endif - - - -// This method pivots row i to row i0 by muliplying row i by -// alpha and adding it to row i0. -// After pivoting the row i0 has a max abs value set correctly at the beginning of m_start, -// Returns false if the resulting row is all zeroes, and true otherwise -template -bool square_sparse_matrix::pivot_row_to_row(unsigned i, const T& alpha, unsigned i0, lp_settings & settings ) { - lp_assert(i < dimension() && i0 < dimension()); - lp_assert(i != i0); - unsigned pivot_col = adjust_column(i); - i = adjust_row(i); - i0 = adjust_row(i0); - vector became_zeros; - // the offset of element of the pivot column in row i0 - scan_row_to_work_vector_and_remove_pivot_column(i0, pivot_col); - auto & i0_row_vals = m_rows[i0]; - // run over the pivot row and update row i0 - unsigned prev_size_i0 = i0_row_vals.size(); - for (const auto & iv : m_rows[i]) { - unsigned j = iv.m_index; - if (j == pivot_col) continue; - T alv = alpha * iv.m_value; - int j_offs = m_work_pivot_vector[j]; - if (j_offs == -1) { // it is a new element - if (!settings.abs_val_is_smaller_than_drop_tolerance(alv)) { - add_new_element(i0, j, alv); - } - } - else { - auto & row_el_iv = i0_row_vals[j_offs]; - row_el_iv.m_value += alv; - if (settings.abs_val_is_smaller_than_drop_tolerance(row_el_iv.m_value)) { - became_zeros.push_back(j_offs); - ensure_increasing(became_zeros); - } - else { - m_columns[j].m_values[row_el_iv.m_other].set_value(row_el_iv.m_value); - } - } - } - - - // clean the work vector - for (unsigned k = 0; k < prev_size_i0; k++) { - m_work_pivot_vector[i0_row_vals[k].m_index] = -1; - } - - for (unsigned k = became_zeros.size(); k-- > 0; ) { - unsigned j = became_zeros[k]; - remove_element(i0_row_vals, i0_row_vals[j]); - if (i0_row_vals.empty()) - return false; - } - - if (numeric_traits::precise() == false) - set_max_in_row(i0_row_vals); - - return !i0_row_vals.empty(); -} - - - -// set the max val as well -// returns false if the resulting row is all zeroes, and true otherwise -template -bool square_sparse_matrix::set_row_from_work_vector_and_clean_work_vector_not_adjusted(unsigned i0, indexed_vector & work_vec, - lp_settings & settings) { - remove_zero_elements_and_set_data_on_existing_elements_not_adjusted(i0, work_vec, settings); - // all non-zero elements in m_work_pivot_vector are new - for (unsigned j : work_vec.m_index) { - if (numeric_traits::is_zero(work_vec[j])) { - continue; - } - lp_assert(!settings.abs_val_is_smaller_than_drop_tolerance(work_vec[j])); - add_new_element(i0, adjust_column(j), work_vec[j]); - work_vec[j] = numeric_traits::zero(); - } - work_vec.m_index.clear(); - auto & row_vals = m_rows[i0]; - if (row_vals.empty()) { - return false; - } - set_max_in_row(row_vals); // it helps to find larger pivots - return true; -} - - - -template -void square_sparse_matrix::remove_zero_elements_and_set_data_on_existing_elements(unsigned row) { - auto & row_vals = m_rows[row]; - for (unsigned k = static_cast(row_vals.size()); k-- > 0;) { // we cannot simply run the iterator since we are removing - // elements from row_vals - auto & row_el_iv = row_vals[k]; - unsigned j = row_el_iv.m_index; - T & wj = m_work_pivot_vector[j]; - if (is_zero(wj)) { - remove_element(row_vals, row_el_iv); - } else { - m_columns[j].m_values[row_el_iv.m_other].set_value(wj); - row_el_iv.set_value(wj); - wj = zero_of_type(); - } - } -} - -// work_vec here has not adjusted column indices -template -void square_sparse_matrix::remove_zero_elements_and_set_data_on_existing_elements_not_adjusted(unsigned row, indexed_vector & work_vec, lp_settings & settings) { - auto & row_vals = m_rows[row]; - for (unsigned k = static_cast(row_vals.size()); k-- > 0;) { // we cannot simply run the iterator since we are removing - // elements from row_vals - auto & row_el_iv = row_vals[k]; - unsigned j = row_el_iv.m_index; - unsigned rj = adjust_column_inverse(j); - T val = work_vec[rj]; - if (settings.abs_val_is_smaller_than_drop_tolerance(val)) { - remove_element(row_vals, row_el_iv); - lp_assert(numeric_traits::is_zero(val)); - } else { - m_columns[j].m_values[row_el_iv.m_other].set_value(row_el_iv.m_value = val); - work_vec[rj] = numeric_traits::zero(); - } - } -} - - -// adding delta columns at the end of the matrix -template -void square_sparse_matrix::add_columns_at_the_end(unsigned delta) { - for (unsigned i = 0; i < delta; i++) { - col_header col_head; - m_columns.push_back(col_head); - } - m_column_permutation.enlarge(delta); -} - -template -void square_sparse_matrix::delete_column(int i) { - lp_assert(i < dimension()); - for (auto cell = m_columns[i].m_head; cell != nullptr;) { - auto next_cell = cell->m_down; - kill_cell(cell); - cell = next_cell; - } -} - -template -void square_sparse_matrix::divide_row_by_constant(unsigned i, const T & t, lp_settings & settings) { - lp_assert(!settings.abs_val_is_smaller_than_zero_tolerance(t)); - i = adjust_row(i); - for (auto & iv : m_rows[i]) { - T &v = iv.m_value; - v /= t; - if (settings.abs_val_is_smaller_than_drop_tolerance(v)){ - v = numeric_traits::zero(); - } - m_columns[iv.m_index].m_values[iv.m_other].set_value(v); - } -} - - -// solving x * this = y, and putting the answer into y -// the matrix here has to be upper triangular -template -void square_sparse_matrix::solve_y_U(vector & y) const { // works by rows -#ifdef Z3DEBUG - // T * rs = clone_vector(y, dimension()); -#endif - unsigned end = dimension(); - for (unsigned i = 0; i + 1 < end; i++) { - // all y[i] has correct values already - const T & yv = y[i]; - if (numeric_traits::is_zero(yv)) continue; - auto & mc = get_row_values(adjust_row(i)); - for (auto & c : mc) { - unsigned col = adjust_column_inverse(c.m_index); - if (col != i) { - y[col] -= c.m_value * yv; - } - } - } -#ifdef Z3DEBUG - // dense_matrix deb(*this); - // T * clone_y = clone_vector(y, dimension()); - // deb.apply_from_right(clone_y); - // lp_assert(vectors_are_equal(rs, clone_y, dimension())); - // delete [] clone_y; - // delete [] rs; -#endif -} - -// solving x * this = y, and putting the answer into y -// the matrix here has to be upper triangular -template -void square_sparse_matrix::solve_y_U_indexed(indexed_vector & y, const lp_settings & settings) { -#if 0 && Z3DEBUG - vector ycopy(y.m_data); - if (numeric_traits::precise() == false) - solve_y_U(ycopy); -#endif - vector sorted_active_columns; - extend_and_sort_active_rows(y.m_index, sorted_active_columns); - for (unsigned k = sorted_active_columns.size(); k-- > 0; ) { - unsigned j = sorted_active_columns[k]; - auto & yj = y[j]; - for (auto & c: m_columns[adjust_column(j)].m_values) { - unsigned i = adjust_row_inverse(c.m_index); - if (i == j) continue; - yj -= y[i] * c.m_value; - } - } - y.m_index.clear(); - for (auto j : sorted_active_columns) { - if (!settings.abs_val_is_smaller_than_drop_tolerance(y[j])) - y.m_index.push_back(j); - else if (!numeric_traits::precise()) - y.m_data[j] = zero_of_type(); - } - - lp_assert(y.is_OK()); -#if 0 && Z3DEBUG - if (numeric_traits::precise() == false) - lp_assert(vectors_are_equal(ycopy, y.m_data)); -#endif -} - - -// fills the indices for such that y[i] can be not a zero -// sort them so the smaller indices come first -// void fill_reachable_indices(std::set & rset, T *y) { -// std::queue q; -// int m = dimension(); -// for (int i = m - 1; i >= 0; i--) { -// if (!numeric_traits::is_zero(y[i])){ -// for (cell * c = m_columns[adjust_column(i)].m_head; c != nullptr; c = c->m_down) { -// unsigned row = adjust_row_inverse(c->m_i); -// q.push(row); -// } -// } -// } -// while (!q.empty()) { -// unsigned i = q.front(); -// q.pop(); -// for (cell * c = m_columns[adjust_column(i)].m_head; c != nullptr; c = c->m_down) { -// unsigned row = adjust_row_inverse(c->m_i); -// if (rset.find(row) == rset.end()){ -// rset.insert(row); -// q.push(row); -// } -// } -// } -// } - -template -template -void square_sparse_matrix::find_error_in_solution_U_y(vector& y_orig, vector & y) { - unsigned i = dimension(); - while (i--) { - y_orig[i] -= dot_product_with_row(i, y); - } -} - -template -template -void square_sparse_matrix::find_error_in_solution_U_y_indexed(indexed_vector& y_orig, indexed_vector & y, const vector& sorted_active_rows) { - for (unsigned i: sorted_active_rows) - y_orig.add_value_at_index(i, -dot_product_with_row(i, y)); // cannot round up here!!! - // y_orig can contain very small values -} - - -template -template -void square_sparse_matrix::add_delta_to_solution(const vector& del, vector & y) { - unsigned i = dimension(); - while (i--) { - y[i] += del[i]; - } -} -template -template -void square_sparse_matrix::add_delta_to_solution(const indexed_vector& del, indexed_vector & y) { -// lp_assert(del.is_OK()); - // lp_assert(y.is_OK()); - for (auto i : del.m_index) { - y.add_value_at_index(i, del[i]); - } -} -template -template -void square_sparse_matrix::double_solve_U_y(indexed_vector& y, const lp_settings & settings){ - lp_assert(y.is_OK()); - indexed_vector y_orig(y); // copy y aside - vector active_rows; - solve_U_y_indexed_only(y, settings, active_rows); - lp_assert(y.is_OK()); - find_error_in_solution_U_y_indexed(y_orig, y, active_rows); - // y_orig contains the error now - if (y_orig.m_index.size() * ratio_of_index_size_to_all_size() < 32 * dimension()) { - active_rows.clear(); - solve_U_y_indexed_only(y_orig, settings, active_rows); - add_delta_to_solution(y_orig, y); - y.clean_up(); - } else { // the dense version - solve_U_y(y_orig.m_data); - add_delta_to_solution(y_orig.m_data, y.m_data); - y.restore_index_and_clean_from_data(); - } - lp_assert(y.is_OK()); -} -template -template -void square_sparse_matrix::double_solve_U_y(vector& y){ - vector y_orig(y); // copy y aside - solve_U_y(y); - find_error_in_solution_U_y(y_orig, y); - // y_orig contains the error now - solve_U_y(y_orig); - add_delta_to_solution(y_orig, y); -} - -// solving this * x = y, and putting the answer into y -// the matrix here has to be upper triangular -template -template -void square_sparse_matrix::solve_U_y(vector & y) { // it is a column wise version -#ifdef Z3DEBUG - // T * rs = clone_vector(y, dimension()); -#endif - - for (unsigned j = dimension(); j--; ) { - const L & yj = y[j]; - if (is_zero(yj)) continue; - for (const auto & iv : m_columns[adjust_column(j)].m_values) { - unsigned i = adjust_row_inverse(iv.m_index); - if (i != j) { - y[i] -= iv.m_value * yj; - } - } - } -#ifdef Z3DEBUG - // dense_matrix deb(*this); - // T * clone_y = clone_vector(y, dimension()); - // deb.apply_from_left(clone_y); - // lp_assert(vectors_are_equal(rs, clone_y, dimension())); -#endif -} -template -void square_sparse_matrix::process_index_recursively_for_y_U(unsigned j, vector & sorted_active_rows) { - lp_assert(m_processed[j] == false); - m_processed[j]=true; - auto & row = m_rows[adjust_row(j)]; - for (auto & c : row) { - unsigned i = adjust_column_inverse(c.m_index); - if (i == j) continue; - if (!m_processed[i]) { - process_index_recursively_for_y_U(i, sorted_active_rows); - } - } - sorted_active_rows.push_back(j); -} - -template -void square_sparse_matrix::process_column_recursively(unsigned j, vector & sorted_active_rows) { - lp_assert(m_processed[j] == false); - auto & mc = m_columns[adjust_column(j)].m_values; - for (auto & iv : mc) { - unsigned i = adjust_row_inverse(iv.m_index); - if (i == j) continue; - if (!m_processed[i]) { - process_column_recursively(i, sorted_active_rows); - } - } - m_processed[j]=true; - sorted_active_rows.push_back(j); -} - - -template -void square_sparse_matrix::create_graph_G(const vector & index_or_right_side, vector & sorted_active_rows) { - for (auto i : index_or_right_side) { - if (m_processed[i]) continue; - process_column_recursively(i, sorted_active_rows); - } - - for (auto i : sorted_active_rows) { - m_processed[i] = false; - } -} - - -template -void square_sparse_matrix::extend_and_sort_active_rows(const vector & index_or_right_side, vector & sorted_active_rows) { - for (auto i : index_or_right_side) { - if (m_processed[i]) continue; - process_index_recursively_for_y_U(i, sorted_active_rows); - } - - for (auto i : sorted_active_rows) { - m_processed[i] = false; - } -} - - -template -template -void square_sparse_matrix::solve_U_y_indexed_only(indexed_vector & y, const lp_settings & settings, vector & sorted_active_rows) { // it is a column wise version - create_graph_G(y.m_index, sorted_active_rows); - - for (auto k = sorted_active_rows.size(); k-- > 0;) { - unsigned j = sorted_active_rows[k]; - const L & yj = y[j]; - if (is_zero(yj)) continue; - auto & mc = m_columns[adjust_column(j)].m_values; - for (auto & iv : mc) { - unsigned i = adjust_row_inverse(iv.m_index); - if (i != j) { - y[i] -= iv.m_value * yj; - } - } - } - y.m_index.clear(); - for (auto j : sorted_active_rows) { - if (!settings.abs_val_is_smaller_than_drop_tolerance(y[j])) - y.m_index.push_back(j); - else if (!numeric_traits::precise()) - y[j] = zero_of_type(); - } - - lp_assert(y.is_OK()); -#ifdef Z3DEBUG - // dense_matrix deb(this); - // vector clone_y(y.m_data); - // deb.apply_from_left(clone_y); - // lp_assert(vectors_are_equal(rs, clone_y)); -#endif -} - -template -template -L square_sparse_matrix::dot_product_with_row (unsigned row, const vector & y) const { - L ret = zero_of_type(); - auto & mc = get_row_values(adjust_row(row)); - for (auto & c : mc) { - unsigned col = m_column_permutation[c.m_index]; - ret += c.m_value * y[col]; - } - return ret; -} - -template -template -L square_sparse_matrix::dot_product_with_row (unsigned row, const indexed_vector & y) const { - L ret = zero_of_type(); - auto & mc = get_row_values(adjust_row(row)); - for (auto & c : mc) { - unsigned col = m_column_permutation[c.m_index]; - ret += c.m_value * y[col]; - } - return ret; -} - - -template -unsigned square_sparse_matrix::get_number_of_nonzeroes() const { - unsigned ret = 0; - for (unsigned i = dimension(); i--; ) { - ret += number_of_non_zeroes_in_row(i); - } - return ret; -} - -template -bool square_sparse_matrix::get_non_zero_column_in_row(unsigned i, unsigned *j) const { - // go over the i-th row - auto & mc = get_row_values(adjust_row(i)); - if (mc.size() > 0) { - *j = m_column_permutation[mc[0].m_index]; - return true; - } - return false; -} - -template -void square_sparse_matrix::remove_element_that_is_not_in_w(vector> & column_vals, indexed_value & col_el_iv) { - auto & row_chunk = m_rows[col_el_iv.m_index]; - indexed_value & row_el_iv = row_chunk[col_el_iv.m_other]; - unsigned index_in_row = col_el_iv.m_other; - remove_element(row_chunk, col_el_iv.m_other, column_vals, row_el_iv.m_other); - if (index_in_row == 0) - set_max_in_row(row_chunk); -} - - -// w contains the new column -// the old column inside of the matrix has not been changed yet -template -void square_sparse_matrix::remove_elements_that_are_not_in_w_and_update_common_elements(unsigned column_to_replace, indexed_vector & w) { - // -------------------------------- - // column_vals represents the old column - auto & column_vals = m_columns[column_to_replace].m_values; - for (unsigned k = static_cast(column_vals.size()); k-- > 0;) { - indexed_value & col_el_iv = column_vals[k]; - unsigned i = col_el_iv.m_index; - T &w_data_at_i = w[adjust_row_inverse(i)]; - if (numeric_traits::is_zero(w_data_at_i)) { - remove_element_that_is_not_in_w(column_vals, col_el_iv); - } else { - auto& row_chunk = m_rows[i]; - unsigned index_in_row = col_el_iv.m_other; - if (index_in_row == 0) { - bool look_for_max = abs(w_data_at_i) < abs(row_chunk[0].m_value); - row_chunk[0].set_value(col_el_iv.m_value = w_data_at_i); - if (look_for_max) - set_max_in_row(i); - } else { - row_chunk[index_in_row].set_value(col_el_iv.m_value = w_data_at_i); - if (abs(w_data_at_i) > abs(row_chunk[0].m_value)) - put_max_index_to_0(row_chunk, index_in_row); - } - w_data_at_i = numeric_traits::zero(); - } - } -} - -template -void square_sparse_matrix::add_new_element(unsigned row, unsigned col, const T& val) { - auto & row_vals = m_rows[row]; - auto & col_vals = m_columns[col].m_values; - unsigned row_el_offs = static_cast(row_vals.size()); - unsigned col_el_offs = static_cast(col_vals.size()); - row_vals.push_back(indexed_value(val, col, col_el_offs)); - col_vals.push_back(indexed_value(val, row, row_el_offs)); - m_n_of_active_elems++; -} - -// w contains the "rest" of the new column; all common elements of w and the old column has been zeroed -// the old column inside of the matrix has not been changed yet -template -void square_sparse_matrix::add_new_elements_of_w_and_clear_w(unsigned column_to_replace, indexed_vector & w, lp_settings & settings) { - for (unsigned i : w.m_index) { - T w_at_i = w[i]; - if (numeric_traits::is_zero(w_at_i)) continue; // was dealt with already - if (!settings.abs_val_is_smaller_than_drop_tolerance(w_at_i)) { - unsigned ai = adjust_row(i); - add_new_element(ai, column_to_replace, w_at_i); - auto & row_chunk = m_rows[ai]; - lp_assert(row_chunk.size() > 0); - if (abs(w_at_i) > abs(row_chunk[0].m_value)) - put_max_index_to_0(row_chunk, static_cast(row_chunk.size()) - 1); - } - w[i] = numeric_traits::zero(); - } - w.m_index.clear(); -} - -template -void square_sparse_matrix::replace_column(unsigned column_to_replace, indexed_vector & w, lp_settings &settings) { - column_to_replace = adjust_column(column_to_replace); - remove_elements_that_are_not_in_w_and_update_common_elements(column_to_replace, w); - add_new_elements_of_w_and_clear_w(column_to_replace, w, settings); -} - -template -unsigned square_sparse_matrix::pivot_score(unsigned i, unsigned j) { - // It goes like this (rnz-1)(cnz-1) is the Markovitz number, that is the max number of - // new non zeroes we can obtain after the pivoting. - // In addition we will get another cnz - 1 elements in the eta matrix created for this pivot, - // which gives rnz(cnz-1). For example, is 0 for a column singleton, but not for - // a row singleton ( which is not a column singleton). - - auto col_header = m_columns[j]; - - return static_cast(get_row_values(i).size() * (col_header.m_values.size() - col_header.m_shortened_markovitz - 1)); -} - -template -void square_sparse_matrix::enqueue_domain_into_pivot_queue() { - lp_assert(m_pivot_queue.size() == 0); - for (unsigned i = 0; i < dimension(); i++) { - auto & rh = m_rows[i]; - unsigned rnz = static_cast(rh.size()); - for (auto iv : rh) { - unsigned j = iv.m_index; - m_pivot_queue.enqueue(i, j, rnz * (static_cast(m_columns[j].m_values.size()) - 1)); - } - } -} - -template -void square_sparse_matrix::set_max_in_rows() { - unsigned i = dimension(); - while (i--) - set_max_in_row(i); -} - - -template -void square_sparse_matrix::zero_shortened_markovitz_numbers() { - for (auto & ch : m_columns) - ch.zero_shortened_markovitz(); -} - -template -void square_sparse_matrix::prepare_for_factorization() { - zero_shortened_markovitz_numbers(); - set_max_in_rows(); - enqueue_domain_into_pivot_queue(); -} - -template -void square_sparse_matrix::recover_pivot_queue(vector & rejected_pivots) { - for (auto p : rejected_pivots) { - m_pivot_queue.enqueue(p.first, p.second, pivot_score(p.first, p.second)); - } -} - -template -int square_sparse_matrix::elem_is_too_small(unsigned i, unsigned j, int c_partial_pivoting) { - auto & row_chunk = m_rows[i]; - - if (j == row_chunk[0].m_index) { - return 0; // the max element is at the head - } - T max = abs(row_chunk[0].m_value); - for (unsigned k = 1; k < row_chunk.size(); k++) { - auto &iv = row_chunk[k]; - if (iv.m_index == j) - return abs(iv.m_value) * c_partial_pivoting < max ? 1: 0; - } - return 2; // the element became zero but it still sits in the active pivots? -} - -template -bool square_sparse_matrix::remove_row_from_active_pivots_and_shorten_columns(unsigned row) { - unsigned arow = adjust_row(row); - for (auto & iv : m_rows[arow]) { - m_pivot_queue.remove(arow, iv.m_index); - m_n_of_active_elems--; // the value is correct only when refactoring - if (adjust_column_inverse(iv.m_index) <= row) - continue; // this column will be removed anyway - auto & col = m_columns[iv.m_index]; - - col.shorten_markovich_by_one(); - if (col.m_values.size() <= col.m_shortened_markovitz) - return false; // got a zero column - } - return true; -} - -template -void square_sparse_matrix::remove_pivot_column(unsigned row) { - unsigned acol = adjust_column(row); - for (const auto & iv : m_columns[acol].m_values) - if (adjust_row_inverse(iv.m_index) >= row) - m_pivot_queue.remove(iv.m_index, acol); -} - -template -void square_sparse_matrix::update_active_pivots(unsigned row) { - unsigned arow = adjust_row(row); - for (const auto & iv : m_rows[arow]) { - col_header & ch = m_columns[iv.m_index]; - int cols = static_cast(ch.m_values.size()) - ch.m_shortened_markovitz - 1; - lp_assert(cols >= 0); - for (const auto &ivc : ch.m_values) { - unsigned i = ivc.m_index; - if (adjust_row_inverse(i) <= row) continue; // the i is not an active row - m_pivot_queue.enqueue(i, iv.m_index, static_cast(m_rows[i].size())*cols); - } - } -} - -template -bool square_sparse_matrix::shorten_active_matrix(unsigned row, eta_matrix *eta_matrix) { - if (!remove_row_from_active_pivots_and_shorten_columns(row)) - return false; - remove_pivot_column(row); - // need to know the max priority of the queue here - update_active_pivots(row); - if (eta_matrix == nullptr) return true; - // it looks like double work, but the pivot scores have changed for all rows - // touched by eta_matrix - for (auto & it : eta_matrix->m_column_vector.m_data) { - unsigned row = adjust_row(it.first); - const auto & row_values = m_rows[row]; - unsigned rnz = static_cast(row_values.size()); - for (auto & iv : row_values) { - const col_header& ch = m_columns[iv.m_index]; - int cnz = static_cast(ch.m_values.size()) - ch.m_shortened_markovitz - 1; - lp_assert(cnz >= 0); - m_pivot_queue.enqueue(row, iv.m_index, rnz * cnz); - } - } - - return true; -} - -template -unsigned square_sparse_matrix::pivot_score_without_shortened_counters(unsigned i, unsigned j, unsigned k) { - auto &cols = m_columns[j].m_values; - unsigned cnz = cols.size(); - for (auto & iv : cols) { - if (adjust_row_inverse(iv.m_index) < k) - cnz--; - } - lp_assert(cnz > 0); - return m_rows[i].m_values.size() * (cnz - 1); -} -#ifdef Z3DEBUG -template -bool square_sparse_matrix::can_improve_score_for_row(unsigned row, unsigned score, T const & c_partial_pivoting, unsigned k) { - unsigned arow = adjust_row(row); - auto & row_vals = m_rows[arow].m_values; - auto & begin_iv = row_vals[0]; - T row_max = abs(begin_iv.m_value); - lp_assert(adjust_column_inverse(begin_iv.m_index) >= k); - if (pivot_score_without_shortened_counters(arow, begin_iv.m_index, k) < score) { - print_active_matrix(k); - return true; - } - for (unsigned jj = 1; jj < row_vals.size(); jj++) { - auto & iv = row_vals[jj]; - lp_assert(adjust_column_inverse(iv.m_index) >= k); - lp_assert(abs(iv.m_value) <= row_max); - if (c_partial_pivoting * abs(iv.m_value) < row_max) continue; - if (pivot_score_without_shortened_counters(arow, iv.m_index, k) < score) { - print_active_matrix(k); - return true; - } - } - return false; -} - -template -bool square_sparse_matrix::really_best_pivot(unsigned i, unsigned j, T const & c_partial_pivoting, unsigned k) { - unsigned queue_pivot_score = pivot_score_without_shortened_counters(i, j, k); - for (unsigned ii = k; ii < dimension(); ii++) { - lp_assert(!can_improve_score_for_row(ii, queue_pivot_score, c_partial_pivoting, k)); - } - return true; -} -template -void square_sparse_matrix::print_active_matrix(unsigned k, std::ostream & out) { - out << "active matrix for k = " << k << std::endl; - if (k >= dimension()) { - out << "empty" << std::endl; - return; - } - unsigned dim = dimension() - k; - dense_matrix b(dim, dim); - for (unsigned i = 0; i < dim; i++) - for (unsigned j = 0; j < dim; j++ ) - b.set_elem(i, j, zero_of_type()); - for (int i = k; i < dimension(); i++) { - unsigned col = adjust_column(i); - for (auto &iv : get_column_values(col)) { - unsigned row = iv.m_index; - unsigned row_ex = this->adjust_row_inverse(row); - if (row_ex < k) continue; - auto v = this->get_not_adjusted(row, col); - b.set_elem(row_ex - k, i -k, v); - } - } - print_matrix(b, out); -} - -template -bool square_sparse_matrix::pivot_queue_is_correct_for_row(unsigned i, unsigned k) { - unsigned arow = adjust_row(i); - for (auto & iv : m_rows[arow].m_values) { - lp_assert(pivot_score_without_shortened_counters(arow, iv.m_index, k + 1) == - m_pivot_queue.get_priority(arow, iv.m_index)); - } - return true; -} - -template -bool square_sparse_matrix::pivot_queue_is_correct_after_pivoting(int k) { - for (unsigned i = k + 1; i < dimension(); i++ ) - lp_assert(pivot_queue_is_correct_for_row(i, k)); - lp_assert(m_pivot_queue.is_correct()); - return true; -} -#endif - -template -bool square_sparse_matrix::get_pivot_for_column(unsigned &i, unsigned &j, int c_partial_pivoting, unsigned k) { - vector pivots_candidates_that_are_too_small; - while (!m_pivot_queue.is_empty()) { - m_pivot_queue.dequeue(i, j); - unsigned i_inv = adjust_row_inverse(i); - if (i_inv < k) continue; - unsigned j_inv = adjust_column_inverse(j); - if (j_inv < k) continue; - int _small = elem_is_too_small(i, j, c_partial_pivoting); - if (!_small) { -#ifdef Z3DEBUG - // if (!really_best_pivot(i, j, c_partial_pivoting, k)) { - // print_active_matrix(k); - // lp_assert(false); - // } -#endif - recover_pivot_queue(pivots_candidates_that_are_too_small); - i = i_inv; - j = j_inv; - return true; - } - if (_small != 2) { // 2 means that the pair is not in the matrix - pivots_candidates_that_are_too_small.push_back(std::make_pair(i, j)); - } - } - recover_pivot_queue(pivots_candidates_that_are_too_small); - return false; -} - -template -bool square_sparse_matrix::elem_is_too_small(vector> & row_chunk, indexed_value & iv, int c_partial_pivoting) { - if (&iv == &row_chunk[0]) { - return false; // the max element is at the head - } - T val = abs(iv.m_value); - T max = abs(row_chunk[0].m_value); - return val * c_partial_pivoting < max; -} - -template -bool square_sparse_matrix::shorten_columns_by_pivot_row(unsigned i, unsigned pivot_column) { - vector> & row_chunk = get_row_values(i); - - for (indexed_value & iv : row_chunk) { - unsigned j = iv.m_index; - if (j == pivot_column) { - lp_assert(!col_is_active(j)); - continue; - } - m_columns[j].shorten_markovich_by_one(); - - if (m_columns[j].m_shortened_markovitz >= get_column_values(j).size()) { // got the zero column under the row! - return false; - } - } - return true; -} - -template -bool square_sparse_matrix::fill_eta_matrix(unsigned j, eta_matrix ** eta) { - const vector> & col_chunk = get_column_values(adjust_column(j)); - bool is_unit = true; - for (const auto & iv : col_chunk) { - unsigned i = adjust_row_inverse(iv.m_index); - if (i > j) { - is_unit = false; - break; - } - if (i == j && iv.m_value != 1) { - is_unit = false; - break; - } - } - - if (is_unit) { - *eta = nullptr; - return true; - } - -#ifdef Z3DEBUG - *eta = new eta_matrix(j, dimension()); -#else - *eta = new eta_matrix(j); -#endif - for (const auto & iv : col_chunk) { - unsigned i = adjust_row_inverse(iv.m_index); - if (i < j) { - continue; - } - if (i > j) { - (*eta)->push_back(i, - iv.m_value); - } else { // i == j - if ( !(*eta)->set_diagonal_element(iv.m_value)) { - delete *eta; - *eta = nullptr; - return false; - } - - } - } - - (*eta)->divide_by_diagonal_element(); - return true; -} -#ifdef Z3DEBUG -template -bool square_sparse_matrix::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings & settings) const { - for (unsigned i = 0; i < dimension(); i++) { - vector> const & row_chunk = get_row_values(i); - lp_assert(row_chunk.size()); - T const & max = abs(row_chunk[0].m_value); - unsigned ai = adjust_row_inverse(i); - for (auto & iv : row_chunk) { - lp_assert(abs(iv.m_value) <= max); - unsigned aj = adjust_column_inverse(iv.m_index); - if (!(ai <= aj || numeric_traits::is_zero(iv.m_value))) - return false; - if (aj == ai) { - if (iv.m_value != 1) { - return false; - } - } - if (settings.abs_val_is_smaller_than_drop_tolerance(iv.m_value) && (!is_zero(iv.m_value))) - return false; - } - } - return true; -} - -template -bool square_sparse_matrix::is_upper_triangular_until(unsigned k) const { - for (unsigned j = 0; j < dimension() && j < k; j++) { - unsigned aj = adjust_column(j); - auto & col = get_column_values(aj); - for (auto & iv : col) { - unsigned row = adjust_row_inverse(iv.m_index); - if (row > j) - return false; - } - } - return true; -} - -template -void square_sparse_matrix::check_column_vs_rows(unsigned col) { - auto & mc = get_column_values(col); - for (indexed_value & column_iv : mc) { - indexed_value & row_iv = column_iv_other(column_iv); - if (row_iv.m_index != col) { - lp_assert(false); - } - - if (& row_iv_other(row_iv) != &column_iv) { - lp_assert(false); - } - - if (row_iv.m_value != column_iv.m_value) { - lp_assert(false); - } - } -} - -template -void square_sparse_matrix::check_row_vs_columns(unsigned row) { - auto & mc = get_row_values(row); - for (indexed_value & row_iv : mc) { - indexed_value & column_iv = row_iv_other(row_iv); - - if (column_iv.m_index != row) { - lp_assert(false); - } - - if (& row_iv != & column_iv_other(column_iv)) { - lp_assert(false); - } - - if (row_iv.m_value != column_iv.m_value) { - lp_assert(false); - } - } -} - -template -void square_sparse_matrix::check_rows_vs_columns() { - for (unsigned i = 0; i < dimension(); i++) { - check_row_vs_columns(i); - } -} - -template -void square_sparse_matrix::check_columns_vs_rows() { - for (unsigned i = 0; i < dimension(); i++) { - check_column_vs_rows(i); - } -} -template -void square_sparse_matrix::check_matrix() { - check_rows_vs_columns(); - check_columns_vs_rows(); -} -#endif -} - diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 2e8faeaed..b997b38f7 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -576,7 +576,7 @@ void test_lp_primal_core_solver() { - +#ifdef Z3DEBUG void fill_uniformly(dense_matrix & m, unsigned dim) { int v = 0; @@ -619,7 +619,7 @@ void test_dense_matrix() { auto c2 = d * p2; } - +#endif vector> vector_of_permutations() {