From 377ceba6d58e436a05092b2f57ebe5f068cc78eb Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Mon, 6 Mar 2023 09:33:53 -0800 Subject: [PATCH] rm lu --- src/math/lp/CMakeLists.txt | 1 - src/math/lp/lp_core_solver_base.h | 14 +- src/math/lp/lp_primal_core_solver.h | 4 - src/math/lp/lp_primal_core_solver_def.h | 13 +- .../lp/lp_primal_core_solver_tableau_def.h | 6 +- src/math/lp/lp_settings.h | 1 - src/math/lp/lp_settings_def.h | 2 - src/math/lp/lu.cpp | 81 -- src/math/lp/lu.h | 325 -------- src/math/lp/lu_def.h | 720 ------------------ src/math/lp/row_eta_matrix.cpp | 1 - src/math/lp/square_sparse_matrix.cpp | 1 - src/math/lp/square_sparse_matrix_def.h | 2 + src/smt/theory_lra.cpp | 2 +- 14 files changed, 22 insertions(+), 1151 deletions(-) delete mode 100644 src/math/lp/lu.cpp delete mode 100644 src/math/lp/lu.h delete mode 100644 src/math/lp/lu_def.h diff --git a/src/math/lp/CMakeLists.txt b/src/math/lp/CMakeLists.txt index cc58c9610fc..97073736717 100644 --- a/src/math/lp/CMakeLists.txt +++ b/src/math/lp/CMakeLists.txt @@ -21,7 +21,6 @@ z3_add_component(lp lp_core_solver_base.cpp lp_primal_core_solver.cpp lp_settings.cpp - lu.cpp lp_utils.cpp matrix.cpp mon_eq.cpp diff --git a/src/math/lp/lp_core_solver_base.h b/src/math/lp/lp_core_solver_base.h index c9b42f28be0..0e8a2abdf3a 100644 --- a/src/math/lp/lp_core_solver_base.h +++ b/src/math/lp/lp_core_solver_base.h @@ -25,11 +25,21 @@ Revision History: #include "math/lp/core_solver_pretty_printer.h" #include "math/lp/numeric_pair.h" #include "math/lp/static_matrix.h" -#include "math/lp/lu.h" #include "math/lp/permutation_matrix.h" #include "math/lp/column_namer.h" +#include "math/lp/u_set.h" + namespace lp { +template +X dot_product(const vector & a, const vector & b) { + lp_assert(a.size() == b.size()); + auto r = zero_of_type(); + for (unsigned i = 0; i < a.size(); i++) { + r += a[i] * b[i]; + } + return r; +} template // X represents the type of the x variable and the bounds class lp_core_solver_base { @@ -156,6 +166,8 @@ class lp_core_solver_base { void pretty_print(std::ostream & out); + + X get_cost() const { return dot_product(m_costs, m_x); } diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h index abee0cc6efd..38a2aa5f5e0 100644 --- a/src/math/lp/lp_primal_core_solver.h +++ b/src/math/lp/lp_primal_core_solver.h @@ -29,7 +29,6 @@ Revision History: #include #include #include -#include "math/lp/lu.h" #include "math/lp/static_matrix.h" #include "math/lp/core_solver_pretty_printer.h" #include "math/lp/lp_core_solver_base.h" @@ -418,9 +417,6 @@ class lp_primal_core_solver:public lp_core_solver_base { // returns the number of iterations unsigned solve(); - lu> * factorization() {return nullptr;} - - void delete_factorization(); // according to Swietanowski, " A new steepest edge approximation for the simplex method for linear programming" void init_column_norms(); diff --git a/src/math/lp/lp_primal_core_solver_def.h b/src/math/lp/lp_primal_core_solver_def.h index 7764ce6d2d9..3967b6e6698 100644 --- a/src/math/lp/lp_primal_core_solver_def.h +++ b/src/math/lp/lp_primal_core_solver_def.h @@ -26,6 +26,7 @@ Revision History: #include #include #include "math/lp/lp_primal_core_solver.h" +#include "math/lp/dense_matrix.h" namespace lp { // This core solver solves (Ax=b, lower_bound_values \leq x \leq upper_bound_values, maximize costs*x ) // The right side b is given implicitly by x and the basis @@ -733,8 +734,7 @@ template unsigned lp_primal_core_solver::solve() default: break; // do nothing } - } while (this->get_status() != lp_status::FLOATING_POINT_ERROR - && + } while ( this->get_status() != lp_status::UNBOUNDED && this->get_status() != lp_status::OPTIMAL @@ -745,18 +745,13 @@ template unsigned lp_primal_core_solver::solve() && !(this->current_x_is_feasible() && this->m_look_for_feasible_solution_only)); - lp_assert(this->get_status() == lp_status::FLOATING_POINT_ERROR - || + lp_assert( this->current_x_is_feasible() == false || this->calc_current_x_is_feasible_include_non_basis()); return this->total_iterations(); } -template void lp_primal_core_solver::delete_factorization() { - lp_assert(false); -} - // according to Swietanowski, " A new steepest edge approximation for the simplex method for linear programming" template void lp_primal_core_solver::init_column_norms() { lp_assert(numeric_traits::precise() == false); @@ -860,7 +855,7 @@ template void lp_primal_core_solver::fill_breakpo template bool lp_primal_core_solver::done() { - if (this->get_status() == lp_status::OPTIMAL || this->get_status() == lp_status::FLOATING_POINT_ERROR) return true; + if (this->get_status() == lp_status::OPTIMAL) return true; if (this->get_status() == lp_status::INFEASIBLE) { return true; } diff --git a/src/math/lp/lp_primal_core_solver_tableau_def.h b/src/math/lp/lp_primal_core_solver_tableau_def.h index 736b79db5ad..e742b0cd01a 100644 --- a/src/math/lp/lp_primal_core_solver_tableau_def.h +++ b/src/math/lp/lp_primal_core_solver_tableau_def.h @@ -157,8 +157,7 @@ unsigned lp_primal_core_solver::solve_with_tableau() { this->set_status(lp_status::CANCELLED); break; // from the loop } - } while (this->get_status() != lp_status::FLOATING_POINT_ERROR - && + } while ( this->get_status() != lp_status::UNBOUNDED && this->get_status() != lp_status::OPTIMAL @@ -168,8 +167,7 @@ unsigned lp_primal_core_solver::solve_with_tableau() { !(this->current_x_is_feasible() && this->m_look_for_feasible_solution_only) ); - lp_assert(this->get_status() == lp_status::FLOATING_POINT_ERROR - || + lp_assert( this->get_status() == lp_status::CANCELLED || this->current_x_is_feasible() == false diff --git a/src/math/lp/lp_settings.h b/src/math/lp/lp_settings.h index 360ef99bf1d..cb60eebc9c5 100644 --- a/src/math/lp/lp_settings.h +++ b/src/math/lp/lp_settings.h @@ -69,7 +69,6 @@ enum class lp_status { DUAL_UNBOUNDED, OPTIMAL, FEASIBLE, - FLOATING_POINT_ERROR, TIME_EXHAUSTED, EMPTY, UNSTABLE, diff --git a/src/math/lp/lp_settings_def.h b/src/math/lp/lp_settings_def.h index 58b37a19dcf..bb87109f6fb 100644 --- a/src/math/lp/lp_settings_def.h +++ b/src/math/lp/lp_settings_def.h @@ -45,7 +45,6 @@ const char* lp_status_to_string(lp_status status) { case lp_status::DUAL_UNBOUNDED: return "DUAL_UNBOUNDED"; case lp_status::OPTIMAL: return "OPTIMAL"; case lp_status::FEASIBLE: return "FEASIBLE"; - case lp_status::FLOATING_POINT_ERROR: return "FLOATING_POINT_ERROR"; case lp_status::TIME_EXHAUSTED: return "TIME_EXHAUSTED"; case lp_status::EMPTY: return "EMPTY"; case lp_status::UNSTABLE: return "UNSTABLE"; @@ -62,7 +61,6 @@ lp_status lp_status_from_string(std::string status) { if (status == "UNBOUNDED") return lp_status::UNBOUNDED; if (status == "OPTIMAL") return lp_status::OPTIMAL; if (status == "FEASIBLE") return lp_status::FEASIBLE; - if (status == "FLOATING_POINT_ERROR") return lp_status::FLOATING_POINT_ERROR; if (status == "TIME_EXHAUSTED") return lp_status::TIME_EXHAUSTED; if (status == "EMPTY") return lp_status::EMPTY; lp_unreachable(); diff --git a/src/math/lp/lu.cpp b/src/math/lp/lu.cpp deleted file mode 100644 index 313fa990151..00000000000 --- a/src/math/lp/lu.cpp +++ /dev/null @@ -1,81 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#include -#include -#include -#include "util/vector.h" -#include "util/debug.h" -#include "math/lp/lu_def.h" -namespace lp { -template double dot_product(vector const&, vector const&); -template lu>::lu(static_matrix const&, vector&, lp_settings&); -template void lu>::push_matrix_to_tail(tail_matrix*); -template void lu>::replace_column(double, indexed_vector&, unsigned); -template lu>::~lu(); -template void lu>::push_matrix_to_tail(tail_matrix*); -template lu>::~lu(); -template void lu>::push_matrix_to_tail(tail_matrix*); -template lu>::~lu(); -template mpq dot_product(vector const&, vector const&); -template void init_factorization> - (lu>*&, static_matrix&, vector&, lp_settings&); -template void init_factorization> - (lu>*&, static_matrix&, vector&, lp_settings&); -template void init_factorization>(lu >*&, static_matrix&, vector&, lp_settings&); -template void print_matrix>(square_sparse_matrix&, std::ostream & out); -template void print_matrix>(static_matrix&, std::ostream&); -template void print_matrix >(static_matrix&, std::ostream&); -template void print_matrix>(static_matrix&, std::ostream & out); -#ifdef Z3DEBUG -template bool lu>::is_correct(const vector& basis); -template bool lu>::is_correct( vector const &); -template dense_matrix get_B>(lu>&, const vector& basis); -template dense_matrix get_B>(lu>&, vector const&); - -#endif - -template bool lu>::pivot_the_row(int); // NOLINT -template void lu>::init_vector_w(unsigned int, indexed_vector&); -template void lu>::solve_By(vector&); -template void lu>::solve_By_when_y_is_ready_for_X(vector&); -template void lu>::solve_yB_with_error_check(vector&, const vector& basis); -template void lu>::solve_yB_with_error_check_indexed(indexed_vector&, vector const&, const vector & basis, const lp_settings&); -template void lu>::replace_column(mpq, indexed_vector&, unsigned); -template void lu>::solve_By(vector&); -template void lu>::solve_By_when_y_is_ready_for_X(vector&); -template void lu>::solve_yB_with_error_check(vector&, const vector& basis); -template void lu>::solve_yB_with_error_check_indexed(indexed_vector&, vector< int > const&, const vector & basis, const lp_settings&); -template void lu >::solve_yB_with_error_check_indexed(indexed_vector&, vector< int > const&, const vector & basis, const lp_settings&); -template void lu >::init_vector_w(unsigned int, indexed_vector&); -template void lu >::replace_column(mpq, indexed_vector&, unsigned); -template void lu >::solve_Bd_faster(unsigned int, indexed_vector&); -template void lu >::solve_By(vector&); -template void lu >::solve_By_when_y_is_ready_for_X(vector&); -template void lu >::solve_yB_with_error_check(vector&, const vector& basis); -template void lu>::solve_By(indexed_vector&); -template void lu>::solve_By(indexed_vector&); -template void lu>::solve_yB_indexed(indexed_vector&); -template void lu >::solve_yB_indexed(indexed_vector&); -template void lu>::solve_By_for_T_indexed_only(indexed_vector&, lp_settings const&); -template void lu>::solve_By_for_T_indexed_only(indexed_vector&, lp_settings const&); -#ifdef Z3DEBUG -template void print_matrix>(tail_matrix&, std::ostream&); -#endif -} diff --git a/src/math/lp/lu.h b/src/math/lp/lu.h deleted file mode 100644 index 18251e3ada5..00000000000 --- a/src/math/lp/lu.h +++ /dev/null @@ -1,325 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - for matrix B we have - t0*...*tn-1*B = Q*U*R - here ti are matrices corresponding to pivot operations, - including columns and rows swaps, - or a multiplication matrix row by a number - Q, R - permutations and U is an upper triangular matrix -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ - -#pragma once - -#include "util/vector.h" -#include "util/debug.h" -#include -#include -#include "math/lp/square_sparse_matrix.h" -#include "math/lp/static_matrix.h" -#include -#include "math/lp/numeric_pair.h" -#include -#include -#include "math/lp/row_eta_matrix.h" -#include "math/lp/square_dense_submatrix.h" -#include "math/lp/dense_matrix.h" -namespace lp { -template // print the nr x nc submatrix at the top left corner -void print_submatrix(square_sparse_matrix & m, unsigned mr, unsigned nc); - -template -void print_matrix(M &m, std::ostream & out); - -template -X dot_product(const vector & a, const vector & b) { - lp_assert(a.size() == b.size()); - auto r = zero_of_type(); - for (unsigned i = 0; i < a.size(); i++) { - r += a[i] * b[i]; - } - return r; -} - - -template -class one_elem_on_diag: public tail_matrix { - unsigned m_i; - T m_val; -public: - one_elem_on_diag(unsigned i, T val) : m_i(i), m_val(val) { -#ifdef Z3DEBUG - m_one_over_val = numeric_traits::one() / m_val; -#endif - } - - bool is_dense() const override { return false; } - - one_elem_on_diag(const one_elem_on_diag & o); - -#ifdef Z3DEBUG - unsigned m_m; - unsigned m_n; - void set_number_of_rows(unsigned m) override { m_m = m; m_n = m; } - void set_number_of_columns(unsigned n) override { m_m = n; m_n = n; } - T m_one_over_val; - - T get_elem (unsigned i, unsigned j) const override; - - unsigned row_count() const override { return m_m; } // not defined } - unsigned column_count() const override { return m_m; } // not defined } -#endif - void apply_from_left(vector & w, lp_settings &) override { - w[m_i] /= m_val; - } - - void apply_from_right(vector & w) override { - w[m_i] /= m_val; - } - - void apply_from_right(indexed_vector & w) override { - if (is_zero(w.m_data[m_i])) - return; - auto & v = w.m_data[m_i] /= m_val; - if (lp_settings::is_eps_small_general(v, 1e-14)) { - w.erase_from_index(m_i); - v = zero_of_type(); - } - } - - - void apply_from_left_to_T(indexed_vector & w, lp_settings & settings) override; - - void conjugate_by_permutation(permutation_matrix & p) { - // this = p * this * p(-1) -#ifdef Z3DEBUG - // auto rev = p.get_reverse(); - // auto deb = ((*this) * rev); - // deb = p * deb; -#endif - m_i = p.apply_reverse(m_i); - -#ifdef Z3DEBUG - // lp_assert(*this == deb); -#endif - } -}; // end of one_elem_on_diag - -enum class LU_status { OK, Degenerated}; - -// This class supports updates of the columns of B, and solves systems Bx=b,and yB=c -// Using Suhl-Suhl method described in the dissertation of Achim Koberstein, Chapter 5 -template -class lu { - LU_status m_status; -public: - typedef typename M::coefftype T; - typedef typename M::argtype X; - - // the fields - unsigned m_dim; - const M & m_A; - permutation_matrix m_Q; - permutation_matrix m_R; - permutation_matrix m_r_wave; - square_sparse_matrix m_U; - square_dense_submatrix* m_dense_LU; - - vector *> m_tail; - lp_settings & m_settings; - bool m_failure; - indexed_vector m_row_eta_work_vector; - indexed_vector m_w_for_extension; - indexed_vector m_y_copy; - indexed_vector m_ii; //to optimize the work with the m_index fields - unsigned m_refactor_counter; - // constructor - // if A is an m by n matrix then basis has length m and values in [0,n); the values are all different - // they represent the set of m columns - lu(const M & A, - vector& basis, - lp_settings & settings); - lu(const M & A, lp_settings&); - void debug_test_of_basis(const M & A, vector & basis); - void solve_Bd_when_w_is_ready(vector & d, indexed_vector& w ); - void solve_By(indexed_vector & y); - - void solve_By(vector & y); - - void solve_By_for_T_indexed_only(indexed_vector& y, const lp_settings &); - - template - void solve_By_when_y_is_ready(indexed_vector & y); - void solve_By_when_y_is_ready_for_X(vector & y); - void solve_By_when_y_is_ready_for_T(vector & y, vector & index); - void print_indexed_vector(indexed_vector & w, std::ofstream & f); - - void print_matrix_compact(std::ostream & f); - - void print(indexed_vector & w, const vector& basis); - void solve_Bd_faster(unsigned a_column, indexed_vector & d); // d is the right side on the input and the solution at the exit - - - void solve_yB_indexed(indexed_vector& y); - - void add_delta_to_solution_indexed(indexed_vector& y); - - void add_delta_to_solution(const vector& yc, vector& y); - - - void find_error_of_yB(vector& yc, const vector& y, - const vector& basis); - - void find_error_of_yB_indexed(const indexed_vector& y, - const vector& heading, const lp_settings& settings); - - - void solve_yB_with_error_check(vector & y, const vector& basis); - - void solve_yB_with_error_check_indexed(indexed_vector & y, const vector& heading, const vector & basis, const lp_settings &); - - void apply_Q_R_to_U(permutation_matrix & r_wave); - - - LU_status get_status() { return m_status; } - - void set_status(LU_status status) { - m_status = status; - } - - ~lu(); - - void init_vector_y(vector & y); - - void perform_transformations_on_w(indexed_vector& w); - - void init_vector_w(unsigned entering, indexed_vector & w); - void apply_lp_list_to_w(indexed_vector & w); - void apply_lp_list_to_y(vector& y); - - void swap_rows(int j, int k); - - void swap_columns(int j, int pivot_column); - - void push_matrix_to_tail(tail_matrix* tm) { - m_tail.push_back(tm); - } - - bool pivot_the_row(int row); - - eta_matrix * get_eta_matrix_for_pivot(unsigned j); - // we're processing the column j now - eta_matrix * get_eta_matrix_for_pivot(unsigned j, square_sparse_matrix& copy_of_U); - - // see page 407 of Chvatal - unsigned transform_U_to_V_by_replacing_column(indexed_vector & w, unsigned leaving_column_of_U); - -#ifdef Z3DEBUG - void check_vector_w(unsigned entering); - - void check_apply_matrix_to_vector(matrix *lp, T *w); - - void check_apply_lp_lists_to_w(T * w); - - // provide some access operators for testing - permutation_matrix & Q() { return m_Q; } - permutation_matrix & R() { return m_R; } - matrix & U() { return m_U; } - unsigned tail_size() { return m_tail.size(); } - - tail_matrix * get_lp_matrix(unsigned i) { - return m_tail[i]; - } - - T B_(unsigned i, unsigned j, const vector& basis) { - return m_A[i][basis[j]]; - } - - unsigned dimension() { return m_dim; } - -#endif - - - unsigned get_number_of_nonzeroes() { - return m_U.get_number_of_nonzeroes(); - } - - - void process_column(int j); - - bool is_correct(const vector& basis); - bool is_correct(); - - - - // needed for debugging purposes - void copy_w(T *buffer, indexed_vector & w); - - // needed for debugging purposes - void restore_w(T *buffer, indexed_vector & w); - bool all_columns_and_rows_are_active(); - - bool too_dense(unsigned j) const; - - void pivot_in_dense_mode(unsigned i); - - void create_initial_factorization(); - - void calculate_r_wave_and_update_U(unsigned bump_start, unsigned bump_end, permutation_matrix & r_wave); - - void scan_last_row_to_work_vector(unsigned lowest_row_of_the_bump); - - bool diagonal_element_is_off(T /* diag_element */) { return false; } - - void pivot_and_solve_the_system(unsigned replaced_column, unsigned lowest_row_of_the_bump); - // see Achim Koberstein's thesis page 58, but here we solve the system and pivot to the last - // row at the same time - - void replace_column(T pivot_elem, indexed_vector & w, unsigned leaving_column_of_U); - - void calculate_Lwave_Pwave_for_bump(unsigned replaced_column, unsigned lowest_row_of_the_bump); - - void calculate_Lwave_Pwave_for_last_row(unsigned lowest_row_of_the_bump, T diagonal_element); - - void prepare_entering(unsigned entering, indexed_vector & w) { - lp_assert(false); - } - bool need_to_refactor() { lp_assert(false); - return m_refactor_counter >= 200; } - - void adjust_dimension_with_matrix_A() { - lp_assert(false); - } - - - - - - -}; // end of lu - -template -void init_factorization(lu* & factorization, M & m_A, vector & m_basis, lp_settings &m_settings); - -#ifdef Z3DEBUG -template -dense_matrix get_B(lu& f, const vector& basis); - -template -dense_matrix get_B(lu& f); -#endif -} diff --git a/src/math/lp/lu_def.h b/src/math/lp/lu_def.h deleted file mode 100644 index ca34481cca8..00000000000 --- a/src/math/lp/lu_def.h +++ /dev/null @@ -1,720 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once - -#include -#include -#include -#include "util/vector.h" -#include -#include "util/debug.h" -#include "math/lp/lu.h" -namespace lp { -template // print the nr x nc submatrix at the top left corner -void print_submatrix(square_sparse_matrix & m, unsigned mr, unsigned nc, std::ostream & out) { - vector> A; - vector widths; - for (unsigned i = 0; i < m.row_count() && i < mr ; i++) { - A.push_back(vector()); - for (unsigned j = 0; j < m.column_count() && j < nc; j++) { - A[i].push_back(T_to_string(static_cast(m(i, j)))); - } - } - - for (unsigned j = 0; j < m.column_count() && j < nc; j++) { - widths.push_back(get_width_of_column(j, A)); - } - - print_matrix_with_widths(A, widths, out); -} - -template -void print_matrix(M &m, std::ostream & out) { - vector> A; - vector widths; - for (unsigned i = 0; i < m.row_count(); i++) { - A.push_back(vector()); - for (unsigned j = 0; j < m.column_count(); j++) { - A[i].push_back(T_to_string(m[i][j])); - } - } - - for (unsigned j = 0; j < m.column_count(); j++) { - widths.push_back(get_width_of_column(j, A)); - } - - print_matrix_with_widths(A, widths, out); -} - -template -one_elem_on_diag::one_elem_on_diag(const one_elem_on_diag & o) { - m_i = o.m_i; - m_val = o.m_val; -#ifdef Z3DEBUG - m_m = m_n = o.m_m; - m_one_over_val = numeric_traits::one() / o.m_val; -#endif -} - -#ifdef Z3DEBUG -template -T one_elem_on_diag::get_elem(unsigned i, unsigned j) const { - if (i == j){ - if (j == m_i) { - return m_one_over_val; - } - return numeric_traits::one(); - } - - return numeric_traits::zero(); -} -#endif -template -void one_elem_on_diag::apply_from_left_to_T(indexed_vector & w, lp_settings & settings) { - T & t = w[m_i]; - if (numeric_traits::is_zero(t)) { - return; - } - t /= m_val; - if (numeric_traits::precise()) return; - if (settings.abs_val_is_smaller_than_drop_tolerance(t)) { - w.erase_from_index(m_i); - t = numeric_traits::zero(); - } -} - -// This class supports updates of the columns of B, and solves systems Bx=b,and yB=c -// Using Suhl-Suhl method described in the dissertation of Achim Koberstein, Chapter 5 -template -lu::lu(const M& A, - vector& basis, - lp_settings & settings): - m_status(LU_status::OK), - m_dim(A.row_count()), - m_A(A), - m_Q(m_dim), - m_R(m_dim), - m_r_wave(m_dim), - m_U(A, basis), // create the square matrix that eventually will be factorized - m_settings(settings), - m_failure(false), - m_row_eta_work_vector(A.row_count()), - m_refactor_counter(0) { - lp_assert(!(numeric_traits::precise() )); -#ifdef Z3DEBUG - debug_test_of_basis(A, basis); -#endif - ++m_settings.stats().m_num_factorizations; - create_initial_factorization(); -#ifdef Z3DEBUG - // lp_assert(check_correctness()); -#endif -} -template -lu::lu(const M& A, - lp_settings & settings): - m_status(LU_status::OK), - m_dim(A.row_count()), - m_A(A), - m_Q(m_dim), - m_R(m_dim), - m_r_wave(m_dim), - m_U(A), // create the square matrix that eventually will be factorized - m_settings(settings), - m_failure(false), - m_row_eta_work_vector(A.row_count()), - m_refactor_counter(0) { - lp_assert(A.row_count() == A.column_count()); - create_initial_factorization(); -#ifdef Z3DEBUG - lp_assert(is_correct()); -#endif -} -template -void lu::debug_test_of_basis( M const & A, vector & basis) { - std::set set; - for (unsigned i = 0; i < A.row_count(); i++) { - lp_assert(basis[i]< A.column_count()); - set.insert(basis[i]); - } - lp_assert(set.size() == A.row_count()); -} - -template -void lu::solve_By(indexed_vector & y) { - lp_assert(false); // not implemented - // init_vector_y(y); - // solve_By_when_y_is_ready(y); - } - - -template -void lu::solve_By(vector & y) { - init_vector_y(y); - solve_By_when_y_is_ready_for_X(y); -} - -template -void lu::solve_By_when_y_is_ready_for_X(vector & y) { - if (numeric_traits::precise()) { - m_U.solve_U_y(y); - m_R.apply_reverse_from_left_to_X(y); // see 24.3 from Chvatal - return; - } - m_U.double_solve_U_y(y); - m_R.apply_reverse_from_left_to_X(y); // see 24.3 from Chvatal - unsigned i = m_dim; - while (i--) { - if (is_zero(y[i])) continue; - if (m_settings.abs_val_is_smaller_than_drop_tolerance(y[i])){ - y[i] = zero_of_type(); - } - } -} - -template -void lu::solve_By_when_y_is_ready_for_T(vector & y, vector & index) { - if (numeric_traits::precise()) { - m_U.solve_U_y(y); - m_R.apply_reverse_from_left_to_T(y); // see 24.3 from Chvatal - unsigned j = m_dim; - while (j--) { - if (!is_zero(y[j])) - index.push_back(j); - } - return; - } - m_U.double_solve_U_y(y); - m_R.apply_reverse_from_left_to_T(y); // see 24.3 from Chvatal - unsigned i = m_dim; - while (i--) { - if (is_zero(y[i])) continue; - if (m_settings.abs_val_is_smaller_than_drop_tolerance(y[i])){ - y[i] = zero_of_type(); - } else { - index.push_back(i); - } - } -} - -template -void lu::solve_By_for_T_indexed_only(indexed_vector & y, const lp_settings & settings) { - if (numeric_traits::precise()) { - vector active_rows; - m_U.solve_U_y_indexed_only(y, settings, active_rows); - m_R.apply_reverse_from_left(y); // see 24.3 from Chvatal - return; - } - m_U.double_solve_U_y(y, m_settings); - m_R.apply_reverse_from_left(y); // see 24.3 from Chvatal -} - -template -void lu::print_matrix_compact(std::ostream & f) { - f << "matrix_start" << std::endl; - f << "nrows " << m_A.row_count() << std::endl; - f << "ncolumns " << m_A.column_count() << std::endl; - for (unsigned i = 0; i < m_A.row_count(); i++) { - auto & row = m_A.m_rows[i]; - f << "row " << i << std::endl; - for (auto & t : row) { - f << "column " << t.m_j << " value " << t.m_value << std::endl; - } - f << "row_end" << std::endl; - } - f << "matrix_end" << std::endl; -} -template -void lu< M>::print(indexed_vector & w, const vector& basis) { - std::string dump_file_name("/tmp/lu"); - remove(dump_file_name.c_str()); - std::ofstream f(dump_file_name); - if (!f.is_open()) { - LP_OUT(m_settings, "cannot open file " << dump_file_name << std::endl); - return; - } - LP_OUT(m_settings, "writing lu dump to " << dump_file_name << std::endl); - print_matrix_compact(f); - print_vector(basis, f); - print_indexed_vector(w, f); - f.close(); -} - -template -void lu< M>::solve_Bd_faster(unsigned a_column, indexed_vector & d) { // puts the a_column into d - init_vector_w(a_column, d); - solve_By_for_T_indexed_only(d, m_settings); -} - -template -void lu< M>::solve_yB_indexed(indexed_vector& y) { - lp_assert(y.is_OK()); - // first solve yU = cb*R(-1) - m_R.apply_reverse_from_right_to_T(y); // got y = cb*R(-1) - lp_assert(y.is_OK()); - m_U.solve_y_U_indexed(y, m_settings); // got y*U=cb*R(-1) - lp_assert(y.is_OK()); - m_Q.apply_reverse_from_right_to_T(y); - lp_assert(y.is_OK()); - for (auto e = m_tail.rbegin(); e != m_tail.rend(); ++e) { -#ifdef Z3DEBUG - (*e)->set_number_of_columns(m_dim); -#endif - (*e)->apply_from_right(y); - lp_assert(y.is_OK()); - } -} - -template -void lu< M>::add_delta_to_solution(const vector& yc, vector& y){ - unsigned i = static_cast(y.size()); - while (i--) - y[i]+=yc[i]; -} - -template -void lu< M>::add_delta_to_solution_indexed(indexed_vector& y) { - // the delta sits in m_y_copy, put result into y - lp_assert(y.is_OK()); - lp_assert(m_y_copy.is_OK()); - m_ii.clear(); - m_ii.resize(y.data_size()); - for (unsigned i : y.m_index) - m_ii.set_value(1, i); - for (unsigned i : m_y_copy.m_index) { - y.m_data[i] += m_y_copy[i]; - if (m_ii[i] == 0) - m_ii.set_value(1, i); - } - lp_assert(m_ii.is_OK()); - y.m_index.clear(); - - for (unsigned i : m_ii.m_index) { - T & v = y.m_data[i]; - if (!lp_settings::is_eps_small_general(v, 1e-14)) - y.m_index.push_back(i); - else if (!numeric_traits::is_zero(v)) - v = zero_of_type(); - } - - lp_assert(y.is_OK()); -} - -template -void lu< M>::find_error_of_yB(vector& yc, const vector& y, const vector& m_basis) { - unsigned i = m_dim; - while (i--) { - yc[i] -= m_A.dot_product_with_column(y, m_basis[i]); - } -} - -template -void lu< M>::find_error_of_yB_indexed(const indexed_vector& y, const vector& heading, const lp_settings& settings) { -#if 0 == 1 - // it is a non efficient version - indexed_vector yc = m_y_copy; - yc.m_index.clear(); - lp_assert(!numeric_traits::precise()); - { - - vector d_basis(y.m_data.size()); - for (unsigned j = 0; j < heading.size(); j++) { - if (heading[j] >= 0) { - d_basis[heading[j]] = j; - } - } - - - unsigned i = m_dim; - while (i--) { - T & v = yc.m_data[i] -= m_A.dot_product_with_column(y.m_data, d_basis[i]); - if (settings.abs_val_is_smaller_than_drop_tolerance(v)) - v = zero_of_type(); - else - yc.m_index.push_back(i); - } - } -#endif - lp_assert(m_ii.is_OK()); - m_ii.clear(); - m_ii.resize(y.data_size()); - lp_assert(m_y_copy.is_OK()); - // put the error into m_y_copy - for (auto k : y.m_index) { - auto & row = m_A.m_rows[k]; - const T & y_k = y.m_data[k]; - for (auto & c : row) { - unsigned j = c.var(); - int hj = heading[j]; - if (hj < 0) continue; - if (m_ii.m_data[hj] == 0) - m_ii.set_value(1, hj); - m_y_copy.m_data[hj] -= c.coeff() * y_k; - } - } - // add the index of m_y_copy to m_ii - for (unsigned i : m_y_copy.m_index) { - if (m_ii.m_data[i] == 0) - m_ii.set_value(1, i); - } - - // there is no guarantee that m_y_copy is OK here, but its index - // is contained in m_ii index - m_y_copy.m_index.clear(); - // setup the index of m_y_copy - for (auto k : m_ii.m_index) { - T& v = m_y_copy.m_data[k]; - if (settings.abs_val_is_smaller_than_drop_tolerance(v)) - v = zero_of_type(); - else { - m_y_copy.set_value(v, k); - } - } - lp_assert(m_y_copy.is_OK()); - -} - - - - -// solves y*B = y -// y is the input -template -void lu< M>::solve_yB_with_error_check_indexed(indexed_vector & y, const vector& heading, const vector & basis, const lp_settings & settings) { - lp_assert(false); - } - - -// solves y*B = y -// y is the input -template -void lu< M>::solve_yB_with_error_check(vector & y, const vector& basis) { - lp_assert(false); -} -template -void lu< M>::apply_Q_R_to_U(permutation_matrix & r_wave) { - m_U.multiply_from_right(r_wave); - m_U.multiply_from_left_with_reverse(r_wave); -} - - -// Solving yB = cb to find the entering variable, -// where cb is the cost vector projected to B. -// The result is stored in cb. - -// solving Bd = a ( to find the column d of B^{-1} A_N corresponding to the entering -// variable -template -lu< M>::~lu(){ - for (auto t : m_tail) { - delete t; - } -} -template -void lu< M>::init_vector_y(vector & y) { - apply_lp_list_to_y(y); - m_Q.apply_reverse_from_left_to_X(y); -} - -template -void lu< M>::perform_transformations_on_w(indexed_vector& w) { - apply_lp_list_to_w(w); - m_Q.apply_reverse_from_left(w); - // TBD does not compile: lp_assert(numeric_traits::precise() || check_vector_for_small_values(w, m_settings)); -} - -// see Chvatal 24.3 -template -void lu< M>::init_vector_w(unsigned entering, indexed_vector & w) { - w.clear(); - m_A.copy_column_to_indexed_vector(entering, w); // w = a, the column - perform_transformations_on_w(w); -} -template -void lu< M>::apply_lp_list_to_w(indexed_vector & w) { - for (unsigned i = 0; i < m_tail.size(); i++) { - m_tail[i]->apply_from_left_to_T(w, m_settings); - // TBD does not compile: lp_assert(check_vector_for_small_values(w, m_settings)); - } -} -template -void lu< M>::apply_lp_list_to_y(vector& y) { - for (unsigned i = 0; i < m_tail.size(); i++) { - m_tail[i]->apply_from_left(y, m_settings); - } -} -template -void lu< M>::swap_rows(int j, int k) { - if (j != k) { - m_Q.transpose_from_left(j, k); - m_U.swap_rows(j, k); - } -} - -template -void lu< M>::swap_columns(int j, int pivot_column) { - if (j == pivot_column) - return; - m_R.transpose_from_right(j, pivot_column); - m_U.swap_columns(j, pivot_column); -} -template -bool lu< M>::pivot_the_row(int row) { - eta_matrix * eta_matrix = get_eta_matrix_for_pivot(row); - if (get_status() != LU_status::OK) { - return false; - } - - if (eta_matrix == nullptr) { - m_U.shorten_active_matrix(row, nullptr); - return true; - } - if (!m_U.pivot_with_eta(row, eta_matrix, m_settings)) - return false; - eta_matrix->conjugate_by_permutation(m_Q); - push_matrix_to_tail(eta_matrix); - return true; -} -// we're processing the column j now -template -eta_matrix * lu< M>::get_eta_matrix_for_pivot(unsigned j) { - eta_matrix *ret; - if(!m_U.fill_eta_matrix(j, &ret)) { - set_status(LU_status::Degenerated); - } - return ret; -} -// we're processing the column j now -template -eta_matrix * lu::get_eta_matrix_for_pivot(unsigned j, square_sparse_matrix& copy_of_U) { - eta_matrix *ret; - copy_of_U.fill_eta_matrix(j, &ret); - return ret; -} - -// see page 407 of Chvatal -template -unsigned lu::transform_U_to_V_by_replacing_column(indexed_vector & w, - unsigned leaving_column) { - unsigned column_to_replace = m_R.apply_reverse(leaving_column); - m_U.replace_column(column_to_replace, w, m_settings); - return column_to_replace; -} - -#ifdef Z3DEBUG -template -void lu::check_vector_w(unsigned entering) { - T * w = new T[m_dim]; - m_A.copy_column_to_vector(entering, w); - check_apply_lp_lists_to_w(w); - delete [] w; -} -template -void lu::check_apply_matrix_to_vector(matrix *lp, T *w) { - if (lp != nullptr) { - lp -> set_number_of_rows(m_dim); - lp -> set_number_of_columns(m_dim); - apply_to_vector(*lp, w); - } -} - -template -void lu::check_apply_lp_lists_to_w(T * w) { - for (unsigned i = 0; i < m_tail.size(); i++) { - check_apply_matrix_to_vector(m_tail[i], w); - } - permutation_matrix qr = m_Q.get_reverse(); - apply_to_vector(qr, w); - for (int i = m_dim - 1; i >= 0; i--) { - lp_assert(abs(w[i] - w[i]) < 0.0000001); - } -} - -#endif -template -void lu::process_column(int j) { - unsigned pi, pj; - bool success = m_U.get_pivot_for_column(pi, pj, m_settings.c_partial_pivoting, j); - if (!success) { - // LP_OUT(m_settings, "get_pivot returned false: cannot find the pivot for column " << j << std::endl); - m_failure = true; - return; - } - - if (static_cast(pi) == -1) { - // LP_OUT(m_settings, "cannot find the pivot for column " << j << std::endl); - m_failure = true; - return; - } - swap_columns(j, pj); - swap_rows(j, pi); - if (!pivot_the_row(j)) { - // LP_OUT(m_settings, "pivot_the_row(" << j << ") failed" << std::endl); - m_failure = true; - } -} -template -bool lu::is_correct(const vector& basis) { -return true; -} - -template -bool lu::is_correct() { - return true; -} - - -// needed for debugging purposes -template -void lu::copy_w(T *buffer, indexed_vector & w) { - -} - -// needed for debugging purposes -template -void lu::restore_w(T *buffer, indexed_vector & w) { - -} -template -bool lu::all_columns_and_rows_are_active() { - return true; -} -template -bool lu::too_dense(unsigned j) const { - return false; -} -template -void lu::pivot_in_dense_mode(unsigned i) { - -} -template -void lu::create_initial_factorization(){ - -} - -template -void lu::calculate_r_wave_and_update_U(unsigned bump_start, unsigned bump_end, permutation_matrix & r_wave) { - if (bump_start > bump_end) { - set_status(LU_status::Degenerated); - return; - } - if (bump_start == bump_end) { - return; - } - - r_wave[bump_start] = bump_end; // sending the offensive column to the end of the bump - - for ( unsigned i = bump_start + 1 ; i <= bump_end; i++ ) { - r_wave[i] = i - 1; - } - - m_U.multiply_from_right(r_wave); - m_U.multiply_from_left_with_reverse(r_wave); -} -template -void lu::scan_last_row_to_work_vector(unsigned lowest_row_of_the_bump) { - vector> & last_row_vec = m_U.get_row_values(m_U.adjust_row(lowest_row_of_the_bump)); - for (auto & iv : last_row_vec) { - if (is_zero(iv.m_value)) continue; - lp_assert(!m_settings.abs_val_is_smaller_than_drop_tolerance(iv.m_value)); - unsigned adjusted_col = m_U.adjust_column_inverse(iv.m_index); - if (adjusted_col < lowest_row_of_the_bump) { - m_row_eta_work_vector.set_value(-iv.m_value, adjusted_col); - } else { - m_row_eta_work_vector.set_value(iv.m_value, adjusted_col); // preparing to calculate the real value in the matrix - } - } -} - -template -void lu::pivot_and_solve_the_system(unsigned replaced_column, unsigned lowest_row_of_the_bump) { - // we have the system right side at m_row_eta_work_vector now - // solve the system column wise - for (unsigned j = replaced_column; j < lowest_row_of_the_bump; j++) { - T v = m_row_eta_work_vector[j]; - if (numeric_traits::is_zero(v)) continue; // this column does not contribute to the solution - unsigned aj = m_U.adjust_row(j); - vector> & row = m_U.get_row_values(aj); - for (auto & iv : row) { - unsigned col = m_U.adjust_column_inverse(iv.m_index); - lp_assert(col >= j || numeric_traits::is_zero(iv.m_value)); - if (col == j) continue; - if (numeric_traits::is_zero(iv.m_value)) { - continue; - } - // the -v is for solving the system ( to zero the last row), and +v is for pivoting - T delta = col < lowest_row_of_the_bump? -v * iv.m_value: v * iv.m_value; - lp_assert(numeric_traits::is_zero(delta) == false); - - - - // m_row_eta_work_vector.add_value_at_index_with_drop_tolerance(col, delta); - if (numeric_traits::is_zero(m_row_eta_work_vector[col])) { - if (!m_settings.abs_val_is_smaller_than_drop_tolerance(delta)){ - m_row_eta_work_vector.set_value(delta, col); - } - } else { - T t = (m_row_eta_work_vector[col] += delta); - if (m_settings.abs_val_is_smaller_than_drop_tolerance(t)){ - m_row_eta_work_vector[col] = numeric_traits::zero(); - auto it = std::find(m_row_eta_work_vector.m_index.begin(), m_row_eta_work_vector.m_index.end(), col); - if (it != m_row_eta_work_vector.m_index.end()) - m_row_eta_work_vector.m_index.erase(it); - } - } - } - } -} - -template -void lu::replace_column(T pivot_elem_for_checking, indexed_vector & w, unsigned leaving_column_of_U){ - lp_assert(false); -} -template -void lu::calculate_Lwave_Pwave_for_bump(unsigned replaced_column, unsigned lowest_row_of_the_bump){ - lp_assert(false);// lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings)); -} - -template -void lu::calculate_Lwave_Pwave_for_last_row(unsigned lowest_row_of_the_bump, T diagonal_element) { - lp_assert(false); -} - -template -void init_factorization(lu* & factorization, M & m_A, vector & m_basis, lp_settings &m_settings) { - lp_assert(false); -} - -#ifdef Z3DEBUG -template -dense_matrix get_B(lu& f, const vector& basis) { - lp_assert(false); - - dense_matrix B(0, 0); - return B; -} -template -dense_matrix get_B(lu& f) { - lp_assert(false); - dense_matrix B(0,0); - return B; -} -#endif -} diff --git a/src/math/lp/row_eta_matrix.cpp b/src/math/lp/row_eta_matrix.cpp index 6fafb83ed6a..356f80b3c01 100644 --- a/src/math/lp/row_eta_matrix.cpp +++ b/src/math/lp/row_eta_matrix.cpp @@ -20,7 +20,6 @@ Revision History: #include #include "util/vector.h" #include "math/lp/row_eta_matrix_def.h" -#include "math/lp/lu.h" namespace lp { template void row_eta_matrix::conjugate_by_permutation(permutation_matrix&); template void row_eta_matrix >::conjugate_by_permutation(permutation_matrix >&); diff --git a/src/math/lp/square_sparse_matrix.cpp b/src/math/lp/square_sparse_matrix.cpp index 35d38e52944..3ec88f47d85 100644 --- a/src/math/lp/square_sparse_matrix.cpp +++ b/src/math/lp/square_sparse_matrix.cpp @@ -20,7 +20,6 @@ Revision History: #include #include "util/vector.h" #include "math/lp/lp_settings.h" -#include "math/lp/lu.h" #include "math/lp/square_sparse_matrix_def.h" #include "math/lp/dense_matrix.h" namespace lp { diff --git a/src/math/lp/square_sparse_matrix_def.h b/src/math/lp/square_sparse_matrix_def.h index 3533ba066b5..19263462717 100644 --- a/src/math/lp/square_sparse_matrix_def.h +++ b/src/math/lp/square_sparse_matrix_def.h @@ -21,6 +21,8 @@ Revision History: #include "util/vector.h" #include "math/lp/square_sparse_matrix.h" +#include "math/lp/dense_matrix.h" + #include #include namespace lp { diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index 2aa988282d2..3eda1ddff43 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -3126,7 +3126,7 @@ class theory_lra::imp { return l_false; TRACE("arith", tout << "status treated as inconclusive: " << status << "\n";); // TENTATIVE_UNBOUNDED, UNBOUNDED, TENTATIVE_DUAL_UNBOUNDED, DUAL_UNBOUNDED, - // FLOATING_POINT_ERROR, TIME_EXAUSTED, EMPTY, UNSTABLE + // TIME_EXAUSTED, EMPTY, UNSTABLE return l_undef; }