diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp index eb1199c30..4953cdb0a 100644 --- a/src/smt/theory_lra.cpp +++ b/src/smt/theory_lra.cpp @@ -308,7 +308,7 @@ class theory_lra::imp { m_solver->set_track_pivoted_rows(lp.bprop_on_pivoted_rows()); m_solver->settings().m_int_gomory_cut_period = ctx().get_fparams().m_arith_branch_cut_ratio; m_solver->settings().m_int_cuts_etc_period = ctx().get_fparams().m_arith_branch_cut_ratio; - m_solver->settings().m_int_cut_solver_period = std::max(8u, ctx().get_fparams().m_arith_branch_cut_ratio); + m_solver->settings().m_int_chase_cut_solver_period = std::max(8u, ctx().get_fparams().m_arith_branch_cut_ratio); m_solver->settings().m_int_run_gcd_test = ctx().get_fparams().m_arith_gcd_test; m_solver->settings().set_random_seed(ctx().get_fparams().m_random_seed); @@ -2812,10 +2812,10 @@ public: st.update("arith-make-feasible", m_solver->settings().st().m_make_feasible); st.update("arith-max-columns", m_solver->settings().st().m_max_cols); st.update("arith-max-rows", m_solver->settings().st().m_max_rows); - st.update("cut-solver-calls", m_solver->settings().st().m_cut_solver_calls); - st.update("cut-solver-true", m_solver->settings().st().m_cut_solver_true); - st.update("cut-solver-false", m_solver->settings().st().m_cut_solver_false); - st.update("cut-solver-undef", m_solver->settings().st().m_cut_solver_undef); + st.update("cut-solver-calls", m_solver->settings().st().m_chase_cut_solver_calls); + st.update("cut-solver-true", m_solver->settings().st().m_chase_cut_solver_true); + st.update("cut-solver-false", m_solver->settings().st().m_chase_cut_solver_false); + st.update("cut-solver-undef", m_solver->settings().st().m_chase_cut_solver_undef); st.update("gcd-calls", m_solver->settings().st().m_gcd_calls); st.update("gcd-conflict", m_solver->settings().st().m_gcd_conflicts); st.update("cube-calls", m_solver->settings().st().m_cube_calls); diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 3e6e93f62..ab4611dab 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -1,22 +1,22 @@ /*++ -Copyright (c) 2017 Microsoft Corporation + Copyright (c) 2017 Microsoft Corporation -Module Name: + Module Name: - + -Abstract: + Abstract: - + -Author: + Author: - Lev Nachmanson (levnach) + Lev Nachmanson (levnach) -Revision History: + Revision History: ---*/ + --*/ #include #if _LINUX_ @@ -51,7 +51,11 @@ Revision History: #include "util/lp/stacked_map.h" #include #include "test/lp/gomory_test.h" - +#include "util/lp/matrix.h" +#include "util/lp/hnf.h" +#include "util/lp/square_sparse_matrix_def.h" +#include "util/lp/lu_def.h" +#include "util/lp/general_matrix.h" namespace lp { unsigned seed = 1; @@ -68,11 +72,11 @@ struct simple_column_namer:public column_namer template -void test_matrix(sparse_matrix & a) { +void test_matrix(square_sparse_matrix & a) { auto m = a.dimension(); // copy a to b in the reversed order - sparse_matrix b(m); + square_sparse_matrix b(m, m); std::cout << "copy b to a"<< std::endl; for (int row = m - 1; row >= 0; row--) for (int col = m - 1; col >= 0; col --) { @@ -110,7 +114,7 @@ void test_matrix(sparse_matrix & a) { void tst1() { std::cout << "testing the minimial matrix with 1 row and 1 column" << std::endl; - sparse_matrix m0(1); + square_sparse_matrix m0(1, 1); m0.set(0, 0, 1); // print_matrix(m0); m0.set(0, 0, 0); @@ -118,7 +122,7 @@ void tst1() { test_matrix(m0); unsigned rows = 2; - sparse_matrix m(rows); + square_sparse_matrix m(rows, rows); std::cout << "setting m(0,1)=" << std::endl; m.set(0, 1, 11); @@ -128,7 +132,7 @@ void tst1() { test_matrix(m); - sparse_matrix m1(2); + square_sparse_matrix m1(2, 2); m1.set(0, 0, 2); m1.set(1, 0, 3); // print_matrix(m1); @@ -141,7 +145,7 @@ void tst1() { std::cout << "printing zero matrix 3 by 1" << std::endl; - sparse_matrix m2(3); + square_sparse_matrix m2(3, 3); // print_matrix(m2); m2.set(0, 0, 1); @@ -151,7 +155,7 @@ void tst1() { test_matrix(m2); - sparse_matrix m10by9(10); + square_sparse_matrix m10by9(10, 10); m10by9.set(0, 1, 1); m10by9(0, 1) = 4; @@ -226,6 +230,7 @@ void change_basis(unsigned entering, unsigned leaving, vector& basis, } + #ifdef Z3DEBUG void test_small_lu(lp_settings & settings) { std::cout << " test_small_lu" << std::endl; @@ -245,7 +250,7 @@ void test_small_lu(lp_settings & settings) { vector heading = allocate_basis_heading(m.column_count()); vector non_basic_columns; init_basis_heading_and_non_basic_columns_vector(basis, heading, non_basic_columns); - lu l(m, basis, settings); + lu> l(m, basis, settings); lp_assert(l.is_correct(basis)); indexed_vector w(m.row_count()); std::cout << "entering 2, leaving 0" << std::endl; @@ -317,7 +322,7 @@ void test_small_lu(lp_settings & settings) { #endif -void fill_long_row(sparse_matrix &m, int i) { +void fill_long_row(square_sparse_matrix &m, int i) { int n = m.dimension(); for (int j = 0; j < n; j ++) { m (i, (j + i) % n) = j * j; @@ -332,7 +337,7 @@ void fill_long_row(static_matrix &m, int i) { } -void fill_long_row_exp(sparse_matrix &m, int i) { +void fill_long_row_exp(square_sparse_matrix &m, int i) { int n = m.dimension(); for (int j = 0; j < n; j ++) { @@ -348,23 +353,23 @@ void fill_long_row_exp(static_matrix &m, int i) { } } -void fill_larger_sparse_matrix_exp(sparse_matrix & m){ +void fill_larger_square_sparse_matrix_exp(square_sparse_matrix & m){ for ( unsigned i = 0; i < m.dimension(); i++ ) fill_long_row_exp(m, i); } -void fill_larger_sparse_matrix_exp(static_matrix & m){ +void fill_larger_square_sparse_matrix_exp(static_matrix & m){ for ( unsigned i = 0; i < m.row_count(); i++ ) fill_long_row_exp(m, i); } -void fill_larger_sparse_matrix(sparse_matrix & m){ +void fill_larger_square_sparse_matrix(square_sparse_matrix & m){ for ( unsigned i = 0; i < m.dimension(); i++ ) fill_long_row(m, i); } -void fill_larger_sparse_matrix(static_matrix & m){ +void fill_larger_square_sparse_matrix(static_matrix & m){ for ( unsigned i = 0; i < m.row_count(); i++ ) fill_long_row(m, i); } @@ -385,12 +390,12 @@ void test_larger_lu_exp(lp_settings & settings) { basis[5] = 6; - fill_larger_sparse_matrix_exp(m); + fill_larger_square_sparse_matrix_exp(m); // print_matrix(m); vector heading = allocate_basis_heading(m.column_count()); vector non_basic_columns; init_basis_heading_and_non_basic_columns_vector(basis, heading, non_basic_columns); - lu l(m, basis, settings); + lu> l(m, basis, settings); dense_matrix left_side = l.get_left_side(basis); dense_matrix right_side = l.get_right_side(); @@ -434,13 +439,13 @@ void test_larger_lu_with_holes(lp_settings & settings) { vector heading = allocate_basis_heading(m.column_count()); vector non_basic_columns; init_basis_heading_and_non_basic_columns_vector(basis, heading, non_basic_columns); - lu l(m, basis, settings); + lu> l(m, basis, settings); std::cout << "printing factorization" << std::endl; for (int i = l.tail_size() - 1; i >=0; i--) { auto lp = l.get_lp_matrix(i); lp->set_number_of_columns(m.row_count()); lp->set_number_of_rows(m.row_count()); - print_matrix( lp, std::cout); + print_matrix( *lp, std::cout); } dense_matrix left_side = l.get_left_side(basis); @@ -469,13 +474,13 @@ void test_larger_lu(lp_settings& settings) { basis[5] = 6; - fill_larger_sparse_matrix(m); + fill_larger_square_sparse_matrix(m); print_matrix(m, std::cout); vector heading = allocate_basis_heading(m.column_count()); vector non_basic_columns; init_basis_heading_and_non_basic_columns_vector(basis, heading, non_basic_columns); - auto l = lu (m, basis, settings); + auto l = lu> (m, basis, settings); // std::cout << "printing factorization" << std::endl; // for (int i = lu.tail_size() - 1; i >=0; i--) { // auto lp = lu.get_lp_matrix(i); @@ -517,7 +522,7 @@ void test_lu(lp_settings & settings) { -void init_b(vector & b, sparse_matrix & m, vector& x) { +void init_b(vector & b, square_sparse_matrix & m, vector& x) { for (unsigned i = 0; i < m.dimension(); i++) { b.push_back(m.dot_product_with_row(i, x)); } @@ -627,7 +632,7 @@ void test_lp_primal_core_solver() { #ifdef Z3DEBUG template -void test_swap_rows_with_permutation(sparse_matrix& m){ +void test_swap_rows_with_permutation(square_sparse_matrix& m){ std::cout << "testing swaps" << std::endl; unsigned dim = m.row_count(); dense_matrix original(&m); @@ -648,10 +653,10 @@ void test_swap_rows_with_permutation(sparse_matrix& m){ } #endif template -void fill_matrix(sparse_matrix& m); // forward definition +void fill_matrix(square_sparse_matrix& m); // forward definition #ifdef Z3DEBUG template -void test_swap_cols_with_permutation(sparse_matrix& m){ +void test_swap_cols_with_permutation(square_sparse_matrix& m){ std::cout << "testing swaps" << std::endl; unsigned dim = m.row_count(); dense_matrix original(&m); @@ -673,9 +678,9 @@ void test_swap_cols_with_permutation(sparse_matrix& m){ template -void test_swap_rows(sparse_matrix& m, unsigned i0, unsigned i1){ +void test_swap_rows(square_sparse_matrix& m, unsigned i0, unsigned i1){ std::cout << "test_swap_rows(" << i0 << "," << i1 << ")" << std::endl; - sparse_matrix mcopy(m.dimension()); + square_sparse_matrix mcopy(m.dimension(), 0); for (unsigned i = 0; i < m.dimension(); i++) for (unsigned j = 0; j < m.dimension(); j++) { mcopy(i, j)= m(i, j); @@ -689,9 +694,9 @@ void test_swap_rows(sparse_matrix& m, unsigned i0, unsigned i1){ } } template -void test_swap_columns(sparse_matrix& m, unsigned i0, unsigned i1){ +void test_swap_columns(square_sparse_matrix& m, unsigned i0, unsigned i1){ std::cout << "test_swap_columns(" << i0 << "," << i1 << ")" << std::endl; - sparse_matrix mcopy(m.dimension()); + square_sparse_matrix mcopy(m.dimension(), 0); // the second argument does not matter for (unsigned i = 0; i < m.dimension(); i++) for (unsigned j = 0; j < m.dimension(); j++) { mcopy(i, j)= m(i, j); @@ -714,7 +719,7 @@ void test_swap_columns(sparse_matrix& m, unsigned i0, unsigned i1){ #endif template -void fill_matrix(sparse_matrix& m){ +void fill_matrix(square_sparse_matrix& m){ int v = 0; for (int i = m.dimension() - 1; i >= 0; i--) { for (int j = m.dimension() - 1; j >=0; j--){ @@ -724,7 +729,7 @@ void fill_matrix(sparse_matrix& m){ } void test_pivot_like_swaps_and_pivot(){ - sparse_matrix m(10); + square_sparse_matrix m(10, 10); fill_matrix(m); // print_matrix(m); // pivot at 2,7 @@ -775,7 +780,7 @@ void test_pivot_like_swaps_and_pivot(){ #ifdef Z3DEBUG void test_swap_rows() { - sparse_matrix m(10); + square_sparse_matrix m(10, 10); fill_matrix(m); // print_matrix(m); test_swap_rows(m, 3, 5); @@ -796,7 +801,7 @@ void test_swap_rows() { test_swap_rows(m, 0, 7); // go over some corner cases - sparse_matrix m0(2); + square_sparse_matrix m0(2, 2); test_swap_rows(m0, 0, 1); m0(0, 0) = 3; test_swap_rows(m0, 0, 1); @@ -804,7 +809,7 @@ void test_swap_rows() { test_swap_rows(m0, 0, 1); - sparse_matrix m1(10); + square_sparse_matrix m1(10, 10); test_swap_rows(m1, 0, 1); m1(0, 0) = 3; test_swap_rows(m1, 0, 1); @@ -816,7 +821,7 @@ void test_swap_rows() { test_swap_rows(m1, 0, 1); - sparse_matrix m2(3); + square_sparse_matrix m2(3, 3); test_swap_rows(m2, 0, 1); m2(0, 0) = 3; test_swap_rows(m2, 0, 1); @@ -824,7 +829,7 @@ void test_swap_rows() { test_swap_rows(m2, 0, 2); } -void fill_uniformly(sparse_matrix & m, unsigned dim) { +void fill_uniformly(square_sparse_matrix & m, unsigned dim) { int v = 0; for (unsigned i = 0; i < dim; i++) { for (unsigned j = 0; j < dim; j++) { @@ -842,9 +847,9 @@ void fill_uniformly(dense_matrix & m, unsigned dim) { } } -void sparse_matrix_with_permutaions_test() { +void square_sparse_matrix_with_permutaions_test() { unsigned dim = 4; - sparse_matrix m(dim); + square_sparse_matrix m(dim, dim); fill_uniformly(m, dim); dense_matrix dm(dim, dim); fill_uniformly(dm, dim); @@ -928,7 +933,7 @@ void sparse_matrix_with_permutaions_test() { } void test_swap_columns() { - sparse_matrix m(10); + square_sparse_matrix m(10, 10); fill_matrix(m); // print_matrix(m); @@ -948,7 +953,7 @@ void test_swap_columns() { test_swap_columns(m, 0, 7); // go over some corner cases - sparse_matrix m0(2); + square_sparse_matrix m0(2, 2); test_swap_columns(m0, 0, 1); m0(0, 0) = 3; test_swap_columns(m0, 0, 1); @@ -956,7 +961,7 @@ void test_swap_columns() { test_swap_columns(m0, 0, 1); - sparse_matrix m1(10); + square_sparse_matrix m1(10, 10); test_swap_columns(m1, 0, 1); m1(0, 0) = 3; test_swap_columns(m1, 0, 1); @@ -968,7 +973,7 @@ void test_swap_columns() { test_swap_columns(m1, 0, 1); - sparse_matrix m2(3); + square_sparse_matrix m2(3, 3); test_swap_columns(m2, 0, 1); m2(0, 0) = 3; test_swap_columns(m2, 0, 1); @@ -1846,7 +1851,7 @@ void test_init_U() { basis[1] = 2; basis[2] = 4; - sparse_matrix u(m, basis); + square_sparse_matrix u(m, basis); for (unsigned i = 0; i < 3; i++) { for (unsigned j = 0; j < 3; j ++) { @@ -1859,7 +1864,7 @@ void test_init_U() { } void test_replace_column() { - sparse_matrix m(10); + square_sparse_matrix m(10, 10); fill_matrix(m); m.swap_columns(0, 7); m.swap_columns(6, 3); @@ -1884,7 +1889,9 @@ void test_replace_column() { } + void setup_args_parser(argument_parser & parser) { + parser.add_option_with_help_string("-hnf", "test hermite normal form"); parser.add_option_with_help_string("-gomory", "gomory"); parser.add_option_with_help_string("-intd", "test integer_domain"); parser.add_option_with_help_string("-xyz_sample", "run a small interactive scenario"); @@ -2660,7 +2667,7 @@ void check_lu_from_file(std::string lufile_name) { vector basis_heading; lp_settings settings; vector non_basic_columns; - lu lsuhl(A, basis, settings); + lu> lsuhl(A, basis, settings); indexed_vector d(A.row_count()); unsigned entering = 26; lsuhl.solve_Bd(entering, d, v); @@ -2677,7 +2684,7 @@ void check_lu_from_file(std::string lufile_name) { void test_square_dense_submatrix() { std::cout << "testing square_dense_submatrix" << std::endl; unsigned parent_dim = 7; - sparse_matrix parent(parent_dim); + square_sparse_matrix parent(parent_dim, 0); fill_matrix(parent); unsigned index_start = 3; square_dense_submatrix d; @@ -3156,6 +3163,8 @@ void test_integer_domain_randomly(integer_domain & d) { } void test_integer_domain() { +#ifdef Z3DEBUG + std::cout << "test_integer_domain\n"; unsigned e0 = 0; unsigned e1 = 1; @@ -3197,6 +3206,7 @@ void test_integer_domain() { d.pop(2); d.print(std::cout); lp_assert(d.has_neg_inf() && d.has_pos_inf()); +#endif // integer_domain d; // std::vector> stack; // for (int i = 0; i < 10000; i++) { @@ -3223,10 +3233,10 @@ void test_integer_domain() { -void test_resolve_with_tight_constraint(cut_solver& cs, - lp::cut_solver::polynomial&i , +void test_resolve_with_tight_constraint(chase_cut_solver& cs, + lp::chase_cut_solver::polynomial&i , unsigned int j, - cut_solver::polynomial& ti) { + chase_cut_solver::polynomial& ti) { // std::cout << "resolve constraint "; // cs.print_polynomial(std::cout, i); @@ -3234,172 +3244,453 @@ void test_resolve_with_tight_constraint(cut_solver& cs, // cs.print_polynomial(std::cout, ti); // std::cout << std::endl; // bool j_coeff_is_one = ti.coeff(j) == 1; - // cut_solver::polynomial result; + // chase_cut_solver::polynomial result; // cs.resolve(i, j, j_coeff_is_one, ti); // std::cout << "resolve result is "; // cs.print_polynomial(std::cout, i); // std::cout << std::endl; } -typedef cut_solver::monomial mono; +typedef chase_cut_solver::monomial mono; -void test_resolve(cut_solver& cs, unsigned constraint_index, unsigned i0) { +void test_resolve(chase_cut_solver& cs, unsigned constraint_index, unsigned i0) { var_index x = 0; var_index y = 1; var_index z = 2; std::cout << "test_resolve\n"; - cut_solver::polynomial i; i += mono(2, x);i += mono(-3,y); + chase_cut_solver::polynomial i; i += mono(2, x);i += mono(-3,y); i+= mono(4, z); i.m_a = 5; - cut_solver::polynomial ti; ti += mono(1, x); ti+= mono(1,y);ti.m_a = 3; + chase_cut_solver::polynomial ti; ti += mono(1, x); ti+= mono(1,y);ti.m_a = 3; test_resolve_with_tight_constraint(cs, i, x, ti); test_resolve_with_tight_constraint(cs, i, y ,ti); - } +} void test_gomory_cut_0() { - gomory_test g( - [](unsigned j) { return "v" + T_to_string(j);} // name_function_p - , - [](unsigned j) { //get_value_p - if (j == 1) - return mpq(2730, 1727); - if (j == 2) - return zero_of_type(); - if (j == 3) return mpq(3); - lp_assert(false); + gomory_test g( + [](unsigned j) { return "v" + T_to_string(j);} // name_function_p + , + [](unsigned j) { //get_value_p + if (j == 1) + return mpq(2730, 1727); + if (j == 2) return zero_of_type(); - }, - [](unsigned j) { // at_low_p - if (j == 1) - return false; - if (j == 2) - return true; - if (j == 3) - return true; - lp_assert(false); + if (j == 3) return mpq(3); + lp_assert(false); + return zero_of_type(); + }, + [](unsigned j) { // at_low_p + if (j == 1) return false; - }, - [](unsigned j) { // at_upper - if (j == 1) - return false; - if (j == 2) - return true; - if (j == 3) - return false; - lp_assert(false); + if (j == 2) + return true; + if (j == 3) + return true; + lp_assert(false); + return false; + }, + [](unsigned j) { // at_upper + if (j == 1) return false; - }, - [](unsigned j) { // lower_bound - if (j == 1) { - lp_assert(false); //unlimited from below - return 0; - } - if (j == 2) - return 0; - if (j == 3) - return 3; - lp_assert(false); + if (j == 2) + return true; + if (j == 3) + return false; + lp_assert(false); + return false; + }, + [](unsigned j) { // lower_bound + if (j == 1) { + lp_assert(false); //unlimited from below return 0; - }, - [](unsigned j) { // upper - if (j == 1) { - lp_assert(false); //unlimited from above - return 0; - } - if (j == 2) - return 0; - if (j == 3) - return 10; - lp_assert(false); + } + if (j == 2) return 0; - }, - [] (unsigned) { return 0; }, - [] (unsigned) { return 0; } - ); - lar_term t; - mpq k; - explanation expl; - unsigned inf_col = 1; - vector> row; - row.push_back(std::make_pair(mpq(1), 1)); - row.push_back(std::make_pair(mpq(2731, 1727), 2)); - row.push_back(std::make_pair(mpq(-910, 1727), 3)); - g.mk_gomory_cut(t, k, expl, inf_col, row); + if (j == 3) + return 3; + lp_assert(false); + return 0; + }, + [](unsigned j) { // upper + if (j == 1) { + lp_assert(false); //unlimited from above + return 0; + } + if (j == 2) + return 0; + if (j == 3) + return 10; + lp_assert(false); + return 0; + }, + [] (unsigned) { return 0; }, + [] (unsigned) { return 0; } + ); + lar_term t; + mpq k; + explanation expl; + unsigned inf_col = 1; + vector> row; + row.push_back(std::make_pair(mpq(1), 1)); + row.push_back(std::make_pair(mpq(2731, 1727), 2)); + row.push_back(std::make_pair(mpq(-910, 1727), 3)); + g.mk_gomory_cut(t, k, expl, inf_col, row); } void test_gomory_cut_1() { - gomory_test g( - [](unsigned j) { return "v" + T_to_string(j);} // name_function_p - , - [](unsigned j) { //get_value_p - if (j == 1) - return mpq(-2); - if (j == 2) - return mpq(4363334, 2730001); - if (j == 3) - return mpq(1); - lp_assert(false); - return zero_of_type(); - }, - [](unsigned j) { // at_low_p - if (j == 1) - return false; - if (j == 2) - return false; - if (j == 3) - return true; - lp_assert(false); + gomory_test g( + [](unsigned j) { return "v" + T_to_string(j);} // name_function_p + , + [](unsigned j) { //get_value_p + if (j == 1) + return mpq(-2); + if (j == 2) + return mpq(4363334, 2730001); + if (j == 3) + return mpq(1); + lp_assert(false); + return zero_of_type(); + }, + [](unsigned j) { // at_low_p + if (j == 1) return false; - }, - [](unsigned j) { // at_upper - if (j == 1) - return true; - if (j == 2) - return false; - if (j == 3) - return true; - lp_assert(false); + if (j == 2) return false; - }, - [](unsigned j) { // lower_bound - if (j == 1) { - lp_assert(false); //unlimited from below - return 0; - } - if (j == 2) - return 1; - if (j == 3) - return 1; - lp_assert(false); + if (j == 3) + return true; + lp_assert(false); + return false; + }, + [](unsigned j) { // at_upper + if (j == 1) + return true; + if (j == 2) + return false; + if (j == 3) + return true; + lp_assert(false); + return false; + }, + [](unsigned j) { // lower_bound + if (j == 1) { + lp_assert(false); //unlimited from below return 0; - }, - [](unsigned j) { // upper - if (j == 1) { - return -2; - } - if (j == 2) - return 3333; - if (j == 3) - return 10000; - lp_assert(false); - return 0; - }, - [] (unsigned) { return 0; }, - [] (unsigned) { return 0; } - ); - lar_term t; - mpq k; - explanation expl; - unsigned inf_col = 2; - vector> row; - row.push_back(std::make_pair(mpq(1726667, 2730001), 1)); - row.push_back(std::make_pair(mpq(-910000, 2730001), 3)); - row.push_back(std::make_pair(mpq(1), 2)); - g.mk_gomory_cut(t, k, expl, inf_col, row); + } + if (j == 2) + return 1; + if (j == 3) + return 1; + lp_assert(false); + return 0; + }, + [](unsigned j) { // upper + if (j == 1) { + return -2; + } + if (j == 2) + return 3333; + if (j == 3) + return 10000; + lp_assert(false); + return 0; + }, + [] (unsigned) { return 0; }, + [] (unsigned) { return 0; } + ); + lar_term t; + mpq k; + explanation expl; + unsigned inf_col = 2; + vector> row; + row.push_back(std::make_pair(mpq(1726667, 2730001), 1)); + row.push_back(std::make_pair(mpq(-910000, 2730001), 3)); + row.push_back(std::make_pair(mpq(1), 2)); + g.mk_gomory_cut(t, k, expl, inf_col, row); } +void test_hnf_m_less_than_n() { +#ifdef Z3DEBUG + general_matrix A; + vector v; + // example 4.3 from Nemhauser, Wolsey + v.push_back(mpq(2)); + v.push_back(mpq(6)); + v.push_back(mpq(1)); + v.push_back(mpq(3)); + A.push_row(v); + v.clear(); + v.push_back(mpq(4)); + v.push_back(mpq(7)); + v.push_back(mpq(7)); + v.push_back(mpq(3)); + A.push_row(v); + v.clear(); + v.push_back(mpq(0)); + v.push_back(mpq(0)); + v.push_back(mpq(1)); + v.push_back(mpq(5)); + A.push_row(v); + unsigned r; + mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r); + hnf h(A, d, r); +#endif +} +void test_hnf_m_greater_than_n() { +#ifdef Z3DEBUG + general_matrix A; + vector v; + v.push_back(mpq(2)); + v.push_back(mpq(6)); + A.push_row(v); + v.clear(); + v.push_back(mpq(4)); + v.push_back(mpq(7)); + A.push_row(v); + v.clear(); + v.push_back(mpq(0)); + v.push_back(mpq(0)); + A.push_row(v); + v.clear(); + v.push_back(mpq(12)); + v.push_back(mpq(55)); + A.push_row(v); + unsigned r; + mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r); + hnf h(A, d, r); +#endif +} + + +void cutting_the_mix_example_1() { + mpq sev(7); + mpq nine(9); + mpq d, u, vv; + hnf_calc::extended_gcd_minimal_uv(sev, nine, d, u, vv); + std::cout << "d = " << d << ", u = " << u << ", vv = " << vv << std::endl; + hnf_calc::extended_gcd_minimal_uv(sev, -nine, d, u, vv); + std::cout << "d = " << d << ", u = " << u << ", vv = " << vv << std::endl; + + hnf_calc::extended_gcd_minimal_uv(-nine, -nine, d, u, vv); + std::cout << "d = " << d << ", u = " << u << ", vv = " << vv << std::endl; + + hnf_calc::extended_gcd_minimal_uv(-sev*2, sev, d, u, vv); + std::cout << "d = " << d << ", u = " << u << ", vv = " << vv << std::endl; + + hnf_calc::extended_gcd_minimal_uv(mpq(24), mpq(-7), d, u, vv); + std::cout << "d = " << d << ", u = " << u << ", vv = " << vv << std::endl; + hnf_calc::extended_gcd_minimal_uv(-mpq(24), mpq(7), d, u, vv); + std::cout << "d = " << d << ", u = " << u << ", vv = " << vv << std::endl; + hnf_calc::extended_gcd_minimal_uv(mpq(24), mpq(7), d, u, vv); + std::cout << "d = " << d << ", u = " << u << ", vv = " << vv << std::endl; + hnf_calc::extended_gcd_minimal_uv(-mpq(21), mpq(7), d, u, vv); + std::cout << "d = " << d << ", u = " << u << ", vv = " << vv << std::endl; + + hnf_calc::extended_gcd_minimal_uv(mpq(21), -mpq(7), d, u, vv); + std::cout << "d = " << d << ", u = " << u << ", vv = " << vv << std::endl; +} + +void test_determinant() { +#ifdef Z3DEBUG + { + auto M = general_matrix(4); + M[0][0] = 1; M[0][1] = -1; M[0][2] = 1; M[0][3] = 1; + M[1][0] = 1; M[1][1] = 0; M[1][2] = 0; M[1][3] = 0; + M[2][0] = 0; M[2][1] = 1; M[2][2] = 4; M[2][3] = 0; + M[3][0] = 0; M[3][1] = 0; M[3][2] = 0; M[3][3] = 4; + std::cout << "M = "; M.print(std::cout, 4); endl(std::cout); + mpq d = hnf_calc::determinant(M); + std::cout << "det M = " << d << std::endl; + } + { + auto M = general_matrix(3); + M[0][0] = 3; M[0][1] = -1; M[0][2] = 1; + M[1][0] = 1; M[1][1] = 0; M[1][2] = 0; + M[2][0] = 0; M[2][1] = 1; M[2][2] = 4; + unsigned r; + std::cout << "M = "; M.print(std::cout, 4); endl(std::cout); + mpq d = hnf_calc::determinant_of_rectangular_matrix(M, r); + std::cout << "det M = " << d << std::endl; + std::cout << "rank = " << r << std::endl; + } + { + auto M = general_matrix(4, 6); + M[0][0] = 3; M[0][1] = -1; M[0][2] = 1; M[0][3] = 1; M[0][4] = 3; M[0][5] = -1; + M[1][0] = 1; M[1][1] = 0; M[1][2] = 0; M[1][3] = 0; M[1][4] = 2; M[1][5] = 7; + M[2][0] = 0; M[2][1] = 1; M[2][2] = 4; M[2][3] = 0; M[2][4] = 2; M[2][5] = 8; + M[3][0] = 6; M[3][1] = -2; M[3][2] = 2; M[3][3] = 2; M[3][4] = 6; M[3][5] = -2; + unsigned r; + std::cout << "M = "; M.print(std::cout, 4); endl(std::cout); + mpq d = hnf_calc::determinant_of_rectangular_matrix(M, r); + std::cout << "det M = " << d << std::endl; + std::cout << "rank = " << r << std::endl; + } +#endif +} + +#ifdef Z3DEBUG + +void fill_general_matrix(general_matrix & M) { + unsigned m = M.row_count(); + unsigned n = M.column_count(); + for (unsigned i = 0; i < m; i++) + for (unsigned j = 0; j < n; j++) + M[i][j] = mpq(static_cast(my_random() % 13) - 6); +} + +void call_hnf(general_matrix A) { + unsigned r; + mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r); + if (r == A.row_count()) + hnf h(A, d, r); +} + + +void test_hnf_for_dim(int m) { + general_matrix M(m, m + my_random() % m); + fill_general_matrix(M); + call_hnf(M); +} +void test_hnf_1_2() { + std::cout << "test_hnf_1_2" << std::endl; + general_matrix A; + vector v; + v.push_back(mpq(5)); + v.push_back(mpq(26)); + A.push_row(v); + unsigned r; + mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r); + hnf h(A, d, r); + std::cout << "test_hnf_1_2 passed" << std::endl; +} +void test_hnf_2_2() { + std::cout << "test_hnf_2_2" << std::endl; + general_matrix A; + vector v; + v.push_back(mpq(5)); + v.push_back(mpq(26)); + A.push_row(v); + v.clear(); + v.push_back(mpq(2)); + v.push_back(mpq(11)); + A.push_row(v); + unsigned r; + mpq d = hnf_calc::determinant_of_rectangular_matrix(A, r); + hnf h(A, d, r); + + std::cout << "test_hnf_2_2 passed" << std::endl; +} + +void test_hnf_3_3() { + std::cout << "test_hnf_3_3" << std::endl; + general_matrix A; + vector v; + v.push_back(mpq(-3)); + v.push_back(mpq(0)); + v.push_back(mpq(-1)); + A.push_row(v); + v.clear(); + v.push_back(mpq(-1)); + v.push_back(mpq(0)); + v.push_back(mpq(-6)); + A.push_row(v); + v.clear(); + v.push_back(mpq(-2)); + v.push_back(mpq(-4)); + v.push_back(mpq(-3)); + A.push_row(v); + + call_hnf(A); + std::cout << "test_hnf_3_3 passed" << std::endl; +} +void test_hnf_4_4() { + std::cout << "test_hnf_4_4" << std::endl; + general_matrix A; + vector v; + v.push_back(mpq(4)); + v.push_back(mpq(3)); + v.push_back(mpq(-5)); + v.push_back(mpq(6)); + A.push_row(v); + v.clear(); + v.push_back(mpq(1)); + v.push_back(mpq(-3)); + v.push_back(mpq(1)); + v.push_back(mpq(-4)); + A.push_row(v); + v.clear(); + v.push_back(mpq(4)); + v.push_back(mpq(4)); + v.push_back(mpq(4)); + v.push_back(mpq(4)); + A.push_row(v); + v.clear(); + v.push_back(mpq(2)); + v.push_back(mpq(-2)); + v.push_back(mpq(-5)); + v.push_back(mpq(6)); + A.push_row(v); + call_hnf(A); + std::cout << "test_hnf_4_4 passed" << std::endl; +} +void test_hnf_5_5() { + std::cout << "test_hnf_5_5" << std::endl; + general_matrix A; + vector v; + v.push_back(mpq(-4)); + v.push_back(mpq(5)); + v.push_back(mpq(-5)); + v.push_back(mpq(1)); + v.push_back(mpq(-3)); + A.push_row(v); + v.clear(); + v.push_back(mpq(3)); + v.push_back(mpq(-1)); + v.push_back(mpq(2)); + v.push_back(mpq(3)); + v.push_back(mpq(-5)); + A.push_row(v); + v.clear(); + v.push_back(mpq(0)); + v.push_back(mpq(6)); + v.push_back(mpq(-5)); + v.push_back(mpq(-6)); + v.push_back(mpq(-2)); + A.push_row(v); + v.clear(); + v.push_back(mpq(1)); + v.push_back(mpq(0)); + v.push_back(mpq(-4)); + v.push_back(mpq(-4)); + v.push_back(mpq(4)); + A.push_row(v); + v.clear(); + v.push_back(mpq(-2)); + v.push_back(mpq(3)); + v.push_back(mpq(6)); + v.push_back(mpq(-5)); + v.push_back(mpq(-1)); + A.push_row(v); + call_hnf(A); + std::cout << "test_hnf_5_5 passed" << std::endl; +} + +void test_hnf() { + test_determinant(); + test_hnf_1_2(); + test_hnf_3_3(); + test_hnf_4_4(); + test_hnf_5_5(); + test_hnf_2_2(); + for (unsigned k=1000; k>0; k--) + for (int i = 1; i < 8; i++) + test_hnf_for_dim(i); + cutting_the_mix_example_1(); + // test_hnf_m_less_than_n(); + // test_hnf_m_greater_than_n(); +} +#endif void test_gomory_cut() { test_gomory_cut_0(); test_gomory_cut_1(); @@ -3420,6 +3711,14 @@ void test_lp_local(int argn, char**argv) { } args_parser.print(); + + if (args_parser.option_is_used("-hnf")) { +#ifdef Z3DEBUG + test_hnf(); +#endif + return finalize(0); + } + if (args_parser.option_is_used("-gomory")) { test_gomory_cut(); return finalize(0); @@ -3464,7 +3763,7 @@ void test_lp_local(int argn, char**argv) { #ifdef Z3DEBUG if (args_parser.option_is_used("--test_swaps")) { - sparse_matrix m(10); + square_sparse_matrix m(10, 0); fill_matrix(m); test_swap_rows_with_permutation(m); test_swap_cols_with_permutation(m); @@ -3574,7 +3873,7 @@ void test_lp_local(int argn, char**argv) { test_init_U(); test_replace_column(); #ifdef Z3DEBUG - sparse_matrix_with_permutaions_test(); + square_sparse_matrix_with_permutaions_test(); test_dense_matrix(); test_swap_operations(); test_permutations(); @@ -3588,3 +3887,8 @@ void test_lp_local(int argn, char**argv) { void tst_lp(char ** argv, int argc, int& i) { lp::test_lp_local(argc - 2, argv + 2); } +#ifdef Z3DEBUG +namespace lp { +template void print_matrix(general_matrix&, std::ostream&); +} +#endif diff --git a/src/util/lp/CMakeLists.txt b/src/util/lp/CMakeLists.txt index 2acf08ca8..a946328f6 100644 --- a/src/util/lp/CMakeLists.txt +++ b/src/util/lp/CMakeLists.txt @@ -4,7 +4,7 @@ z3_add_component(lp binary_heap_priority_queue.cpp binary_heap_upair_queue.cpp bound_propagator.cpp - cut_solver.cpp + chase_cut_solver.cpp core_solver_pretty_printer.cpp dense_matrix.cpp eta_matrix.cpp @@ -25,7 +25,7 @@ z3_add_component(lp permutation_matrix.cpp row_eta_matrix.cpp scaler.cpp - sparse_matrix.cpp + square_sparse_matrix.cpp square_dense_submatrix.cpp static_matrix.cpp random_updater.cpp diff --git a/src/util/lp/cut_solver.cpp b/src/util/lp/chase_cut_solver.cpp similarity index 94% rename from src/util/lp/cut_solver.cpp rename to src/util/lp/chase_cut_solver.cpp index 2bfc7c229..1369e0710 100644 --- a/src/util/lp/cut_solver.cpp +++ b/src/util/lp/chase_cut_solver.cpp @@ -2,7 +2,7 @@ Copyright (c) 2017 Microsoft Corporation Author: Nikolaj Bjorner, Lev Nachmanson */ -#include "util/lp/cut_solver.h" +#include "util/lp/chase_cut_solver.h" namespace lp { mpq polynomial::m_local_zero = zero_of_type(); diff --git a/src/util/lp/cut_solver.h b/src/util/lp/chase_cut_solver.h similarity index 99% rename from src/util/lp/cut_solver.h rename to src/util/lp/chase_cut_solver.h index aaa074f23..4b41468fe 100644 --- a/src/util/lp/cut_solver.h +++ b/src/util/lp/chase_cut_solver.h @@ -35,24 +35,24 @@ #include "util/lp/indexer_of_constraints.h" #include "util/lp/lar_term.h" namespace lp { -class cut_solver; //forward definition +class chase_cut_solver; //forward definition struct pp_poly { - cut_solver const& s; + chase_cut_solver const& s; polynomial const& p; - pp_poly(cut_solver const& s, polynomial const& p): s(s), p(p) {} + pp_poly(chase_cut_solver const& s, polynomial const& p): s(s), p(p) {} }; struct pp_constraint { - cut_solver const& s; + chase_cut_solver const& s; constraint const& c; - pp_constraint(cut_solver const& s, constraint const& c): s(s), c(c) {} + pp_constraint(chase_cut_solver const& s, constraint const& c): s(s), c(c) {} }; std::ostream& operator<<(std::ostream& out, pp_poly const& p); std::ostream& operator<<(std::ostream& out, pp_constraint const& p); -class cut_solver : public column_namer { +class chase_cut_solver : public column_namer { public: enum class lbool { l_false, l_true, l_undef }; inline std::string lbool_to_string(lbool l) { @@ -356,7 +356,7 @@ public: } - ~cut_solver() { + ~chase_cut_solver() { for (constraint * c : m_asserts) delete c; for (constraint * c : m_lemmas_container.m_lemmas) @@ -386,7 +386,7 @@ public: unsigned m_bounded_search_calls; unsigned m_number_of_conflicts; vector m_scopes; - std::unordered_map m_user_vars_to_cut_solver_vars; + std::unordered_map m_user_vars_to_chase_cut_solver_vars; std::unordered_set m_explanation; // if this collection is not empty we have a conflict // the number of decisions in the current trail unsigned m_decision_level; @@ -846,7 +846,7 @@ public: } bool too_many_propagations_for_var(const var_info& vi) const { - return vi.number_of_bound_propagations() > m_settings.m_cut_solver_cycle_on_var * vi.number_of_asserts(); + return vi.number_of_bound_propagations() > m_settings.m_chase_cut_solver_cycle_on_var * vi.number_of_asserts(); } bool new_upper_bound_is_relevant(unsigned j, const mpq & v) const { @@ -1591,7 +1591,7 @@ public: unsigned find_unused_index() const { for (unsigned j = m_var_infos.size(); ; j++) - if (m_user_vars_to_cut_solver_vars.find(j) == m_user_vars_to_cut_solver_vars.end()) + if (m_user_vars_to_chase_cut_solver_vars.find(j) == m_user_vars_to_chase_cut_solver_vars.end()) return j; } @@ -1715,7 +1715,7 @@ public: m_cancelled = true; return true; } - unsigned bound = m_asserts.size() * 200 / (1 + m_settings.m_int_cut_solver_period); + unsigned bound = m_asserts.size() * 200 / (1 + m_settings.m_int_chase_cut_solver_period); if (m_trail.size() > bound || m_number_of_conflicts > bound) { m_cancelled = true; return true; @@ -1745,7 +1745,7 @@ public: if (c != nullptr) { m_number_of_conflicts++; - TRACE("propagate_and_backjump_step_int", tout << "incostistent_constraint "; trace_print_constraint(tout, c); tout << "m_number_of_conflicts = " << m_number_of_conflicts << std::endl; tout << "cut_solver_calls " << m_settings.st().m_cut_solver_calls << std::endl;); + TRACE("propagate_and_backjump_step_int", tout << "incostistent_constraint "; trace_print_constraint(tout, c); tout << "m_number_of_conflicts = " << m_number_of_conflicts << std::endl; tout << "chase_cut_solver_calls " << m_settings.st().m_chase_cut_solver_calls << std::endl;); if (at_base_lvl()) { fill_conflict_explanation(c); return lbool::l_false; @@ -2223,7 +2223,7 @@ public: public: void pop() { pop(1); } - cut_solver(std::function var_name_function, + chase_cut_solver(std::function var_name_function, std::function print_constraint_function, std::function number_of_variables_function, std::function var_value_function, @@ -2343,7 +2343,7 @@ public: bool consistent(const_constr * i) const { // an option could be to check that upper(i.poly()) <= 0 bool ret = value(i->poly()) <= zero_of_type(); - CTRACE("cut_solver_state_inconsistent", !ret, + CTRACE("chase_cut_solver_state_inconsistent", !ret, tout << "inconsistent constraint " << pp_constraint(*this, *i) << "\n"; tout << "value = " << value(i->poly()) << '\n';); return ret; @@ -2517,7 +2517,7 @@ public: constraint *c = short_lemmas[m_settings.random_next() % n]; k = - c->poly().m_a; copy_poly_coeffs_to_term(c->poly(), t); - TRACE("cut_solver_cuts", trace_print_constraint(tout, *c);); + TRACE("chase_cut_solver_cuts", trace_print_constraint(tout, *c);); return true; } }; diff --git a/src/util/lp/general_matrix.h b/src/util/lp/general_matrix.h new file mode 100644 index 000000000..2994d6101 --- /dev/null +++ b/src/util/lp/general_matrix.h @@ -0,0 +1,195 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#pragma once +namespace lp { +class general_matrix { + // fields + permutation_matrix m_row_permutation; + permutation_matrix m_column_permutation; + vector> m_data; + +public: + unsigned adjust_row(unsigned row) const{ + return m_row_permutation[row]; + } + + void push_row(vector & v) { + m_data.push_back(v); + m_row_permutation.resize(m_data.size()); + m_column_permutation.resize(v.size()); + } + + 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]; + } + + + unsigned row_count() const { return m_data.size(); } + unsigned column_count() const { return m_data[0].size(); } + + class ref_row { + general_matrix& m_matrix; + vector& m_row_data; + public: + ref_row(general_matrix& m, vector& row_data) : m_matrix(m), m_row_data(row_data) {} + mpq & operator[](unsigned col) { return m_row_data[m_matrix.adjust_column(col)]; } + }; + class ref_row_const { + const general_matrix& m_matrix; + const vector& m_row_data; + public: + ref_row_const(const general_matrix& m, const vector& row_data) : m_matrix(m), m_row_data(row_data) {} + const mpq& operator[](unsigned col) const { return m_row_data[m_matrix.adjust_column(col)]; } + }; + + ref_row operator[](unsigned i) { return ref_row(*this, m_data[adjust_row(i)]); } + ref_row_const operator[](unsigned i) const { return ref_row_const(*this, m_data[adjust_row(i)]); } + + void print(std::ostream & out, unsigned blanks = 0) const { + print_matrix(m_data, out, blanks); + } + + void clear() { m_data.clear(); } + + void print_submatrix(std::ostream & out, unsigned k, unsigned blanks = 0) const { + general_matrix m(row_count() - k, column_count() - k); + for (unsigned i = k; i < row_count(); i++) { + for (unsigned j = k; j < column_count(); j++) + m[i-k][j-k] = (*this)[i][j]; + } + print_matrix(m.m_data, out, blanks); + } + + + void copy_column_to_indexed_vector(unsigned entering, indexed_vector &w ) const { + lp_assert(false); // + } + general_matrix operator*(const general_matrix & m) const { + lp_assert(m.row_count() == column_count()); + general_matrix ret(row_count(), m.column_count()); + for (unsigned i = 0; i < row_count(); i ++) { + for (unsigned j = 0; j < m.column_count(); j++) { + mpq a(0); + for (unsigned k = 0; k < column_count(); k++) + a += ((*this)[i][k])*m[k][j]; + ret[i][j] = a; + } + } + return ret; + } + + bool elements_are_equal(const general_matrix& m) const { + for (unsigned i = 0; i < row_count(); i++) + for (unsigned j = 0; j < column_count(); j++) + if ( (*this)[i][j] != m[i][j]) + return false; + return true; + } + + bool elements_are_equal_modulo(const general_matrix& m, const mpq & d) const { + for (unsigned i = 0; i < row_count(); i++) + for (unsigned j = 0; j < column_count(); j++) + if (!is_zero(((*this)[i][j] - m[i][j]) % d)) + return false; + return true; + } + bool operator==(const general_matrix& m) const { + return row_count() == m.row_count() && column_count() == m.column_count() && elements_are_equal(m); + } + + bool operator!=(const general_matrix& m) const { + return !(*this == m); + } + + bool equal_modulo(const general_matrix& m, const mpq & d) const { + return row_count() == m.row_count() && column_count() == m.column_count() && elements_are_equal_modulo(m, d); + } + + + vector operator*(const vector & x) const { + vector r; + lp_assert(x.size() == column_count()); + for (unsigned i = 0; i < row_count(); i++) { + mpq v(0); + for (unsigned j = 0; j < column_count(); j++) { + v += (*this)[i][j] * x[j]; + } + r.push_back(v); + } + return r; + } + + // bool create_upper_triangle(general_matrix& m, vector& x) { + // for (unsigned i = 1; i < m.row_count(); i++) { + // lp_assert(false); // to be continued + // } + // } + + // bool solve_A_x_equal_b(const general_matrix& m, vector& x, const vector& b) const { + // auto m_copy = m; + // // for square matrices + // lp_assert(row_count() == b.size()); + // lp_assert(x.size() == column_count()); + // lp_assert(row_count() == column_count()); + // x = b; + // create_upper_triangle(copy_of_m, x); + // solve_on_triangle(copy_of_m, x); + // } + // + + void transpose_rows(unsigned i, unsigned l) { + lp_assert(i != l); + m_row_permutation.transpose_from_right(i, l); + } + + void transpose_columns(unsigned j, unsigned k) { + lp_assert(j != k); + m_column_permutation.transpose_from_left(j, k); + } + + general_matrix(){} + general_matrix(unsigned n) : + m_row_permutation(n), + m_column_permutation(n), + m_data(n) + { + for (auto& v : m_data){ + v.resize(n); + } + } + + general_matrix(unsigned m, unsigned n) : + m_row_permutation(m), + m_column_permutation(n), + m_data(m) { + for (auto& v : m_data){ + v.resize(n); + } + } +}; +} diff --git a/src/util/lp/hnf.h b/src/util/lp/hnf.h new file mode 100644 index 000000000..8227deb46 --- /dev/null +++ b/src/util/lp/hnf.h @@ -0,0 +1,634 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + Creates the Hermite Normal Form of a matrix in place. + We suppose that $A$ is an integral $m$ by $n$ matrix or rank $m$, where $n >= m$. + The paragraph below is applicable to the usage of HNF. +We have $H = AU$ where $H$ is in Hermite Normal Form +and $U$ is a unimodular matrix. We do not have an explicit + representation of $U$. For a given $i$ we need to find the $i$-th + row of $U^{-1}$. +Let $e_i$ be a vector of length $m$ with all elements equal to $0$ and +$1$ at $i$-th position. Then we need to find the row vector $e_iU^{-1}=t$. Noticing that $U^{-1} = H^{-1}A$, we have $e_iH^{-1}A=t$. +We find $e_iH^{-1} = f$ by solving $e_i = fH$ and then $fA$ gives us $t$. + +Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#pragma once +#include "util/lp/numeric_pair.h" +#include "util/ext_gcd.h" +namespace lp { +namespace hnf_calc { + + // d = u * a + b * b and the sum of abs(u) + abs(v) is minimal, d is positive +inline +void extended_gcd_minimal_uv(const mpq & a, const mpq & b, mpq & d, mpq & u, mpq & v) { + if (is_zero(a)) { + u = zero_of_type(); + v = one_of_type(); + d = b; + return; + } + if (is_zero(b)) { + u = one_of_type(); + v = zero_of_type(); + d = a; + return; + } + extended_gcd(a, b, d, u, v); + if (is_neg(d)) { + d = -d; + u = -u; + v = -v; + } + + if (d == a) { + u = one_of_type(); + v = zero_of_type(); + return; + } + if (d == -a) { + u = - one_of_type(); + v = zero_of_type(); + return; + } + + mpq a_over_d = abs(a) / d; + mpq r; + + mpq k = machine_div_rem(v, a_over_d, r); + if (is_neg(r)) { + r += a_over_d; + k -= one_of_type(); + } + + lp_assert(v == k * a_over_d + r); + + if (is_pos(b)) { + v = r - a_over_d; // v -= (k + 1) * a_over_d; + lp_assert(- a_over_d < v && v <= zero_of_type()); + + if (is_pos(a)) { + u += (k + 1) * (b / d); + lp_assert( one_of_type() <= u && u <= abs(b)/d); + } else { + u -= (k + 1) * (b / d); + lp_assert( one_of_type() <= -u && -u <= abs(b)/d); + } + } else { + v = r; // v -= k * a_over_d; + lp_assert(- a_over_d < -v && -v <= zero_of_type()); + if (is_pos(a)) { + u += k * (b / d); + lp_assert( one_of_type() <= u && u <= abs(b)/d); + } else { + u -= k * (b / d); + lp_assert( one_of_type() <= -u && -u <= abs(b)/d); + } + } + lp_assert(d == u * a + v * b); +} + + + +template bool prepare_pivot_for_lower_triangle(M &m, unsigned r) { + lp_assert(m.row_count() <= m.column_count()); + for (unsigned i = r; i < m.row_count(); i++) { + for (unsigned j = r; j < m.column_count(); j++) { + if (!is_zero(m[i][j])) { + if (i != r) { + m.transpose_rows(i, r); + } + if (j != r) { + m.transpose_columns(j, r); + } + return true; + } + } + } + return false; +} + +template void pivot_column_non_fractional(M &m, unsigned & r) { + lp_assert(m.row_count() <= m.column_count()); + for (unsigned j = r + 1; j < m.column_count(); j++) { + for (unsigned i = r + 1; i < m.row_count(); i++) { + m[i][j] = + (r > 0 )? + (m[r][r]*m[i][j] - m[i][r]*m[r][j]) / m[r-1][r-1] : + (m[r][r]*m[i][j] - m[i][r]*m[r][j]); + lp_assert(is_int(m[i][j])); + } + } + +} + +// returns the rank of the matrix +template unsigned to_lower_triangle_non_fractional(M &m) { + lp_assert(m.row_count() <= m.column_count()); + unsigned i = 0; + for (; i < m.row_count() - 1; i++) { + if (!prepare_pivot_for_lower_triangle(m, i)) { + return i; + } + pivot_column_non_fractional(m, i); + } + lp_assert(i == m.row_count() - 1); + // go over the last row and try to find a non-zero in the row to the right of diagonal + for (unsigned j = i; j < m.column_count(); j++) { + if (!is_zero(m[i][j])) { + return m.row_count(); + } + } + return i; +} +template +mpq gcd_of_row_starting_from_diagonal(const M& m, unsigned i) { + mpq g = zero_of_type(); + unsigned j = i; + for (; j < m.column_count() && is_zero(j); j++) { + const auto & t = m[i][j]; + if (!is_zero(t)) + g = t; + } + for (; j < m.column_count(); j++) { + const auto & t = m[i][j]; + if (!is_zero(t)) + g = gcd(g, t); + } + return g; +} + + +// it returns "r" - the rank of the matrix and the gcd of r-minors +template mpq determinant_of_rectangular_matrix(const M& m, unsigned & r) { + if (m.column_count() < m.row_count()) + throw "not implemented"; // consider working with the transposed m, or create a "transposed" version if needed + // The plan is to transform m to the lower triangular form by using non-fractional Gaussian Elimination by columns. + // Then the elements of the following elements of the last non-zero row of the matrix + // m[r-1][r-1], m[r-1][r], ..., m[r-1]m[m.column_count() - 1] give the determinants of all minors of rank r. + // The gcd of these minors is the return value + auto mc = m; + r = to_lower_triangle_non_fractional(mc); + if (r == 0) + return one_of_type(); + lp_assert(!is_zero(gcd_of_row_starting_from_diagonal(mc, r - 1))); + return gcd_of_row_starting_from_diagonal(mc, r - 1); +} + +template mpq determinant(const M& m) { + lp_assert(m.row_count() == m.column_count()); + auto mc = m; + unsigned r; + mpq d = determinant_of_rectangular_matrix(mc, r); + return r < m.row_count() ? zero_of_type() : d; +} + +} // end of namespace hnf_calc + +template // M is the matrix type +class hnf { + // fields + +#ifdef Z3DEBUG + M & m_H; + M m_U; + M m_U_reverse; +#endif + + M m_A_orig; + M m_W; + vector m_buffer; + unsigned m_m; + unsigned m_n; + mpq m_d; // it is a positive number and a multiple of gcd of r-minors of m_A_orig, where r is the rank of m_A_orig + mpq m_r; // the rank of m_A + unsigned m_i; + unsigned m_j; + mpq m_R; + mpq m_half_R; + mpq mod_R_balanced(const mpq & a) const { + mpq t = a % m_R; + return t > m_half_R? t - m_R : (t < - m_half_R? t + m_R : t); + } + + mpq mod_R(const mpq & a) const { + mpq t = a % m_R; + t = is_neg(t) ? t + m_R : t; + if(is_neg(t)){ + std::cout << "a=" << a << ", m_R= " << m_R << std::endl; + } + return t; + + } + +#ifdef Z3DEBUG + void buffer_p_col_i_plus_q_col_j_H(const mpq & p, unsigned i, const mpq & q, unsigned j) { + for (unsigned k = i; k < m_m; k++) { + m_buffer[k] = p * m_H[k][i] + q * m_H[k][j]; + } + } +#endif + bool zeros_in_column_W_above(unsigned i) { + for (unsigned k = 0; k < i; k++) + if (!is_zero(m_W[k][i])) + return false; + return true; + } + + void buffer_p_col_i_plus_q_col_j_W_modulo(const mpq & p, const mpq & q) { + lp_assert(zeros_in_column_W_above(m_i)); + for (unsigned k = m_i; k < m_m; k++) { + m_buffer[k] = mod_R_balanced(mod_R_balanced(p * m_W[k][m_i]) + mod_R_balanced(q * m_W[k][m_j])); + } + } +#ifdef Z3DEBUG + void buffer_p_col_i_plus_q_col_j_U(const mpq & p, unsigned i, const mpq & q, unsigned j) { + for (unsigned k = 0; k < m_n; k++) { + m_buffer[k] = p * m_U[k][i] + q * m_U[k][j]; + } + } + + void pivot_column_i_to_column_j_H(mpq u, unsigned i, mpq v, unsigned j) { + lp_assert(is_zero(u * m_H[i][i] + v * m_H[i][j])); + m_H[i][j] = zero_of_type(); + for (unsigned k = i + 1; k < m_m; k ++) + m_H[k][j] = u * m_H[k][i] + v * m_H[k][j]; + + } +#endif + void pivot_column_i_to_column_j_W_modulo(mpq u, mpq v) { + lp_assert(is_zero((u * m_W[m_i][m_i] + v * m_W[m_i][m_j]) % m_R)); + m_W[m_i][m_j] = zero_of_type(); + for (unsigned k = m_i + 1; k < m_m; k ++) + m_W[k][m_j] = mod_R_balanced(mod_R_balanced(u * m_W[k][m_i]) + mod_R_balanced(v * m_W[k][m_j])); + } + +#ifdef Z3DEBUG + void pivot_column_i_to_column_j_U(mpq u, unsigned i, mpq v, unsigned j) { + for (unsigned k = 0; k < m_n; k ++) + m_U[k][j] = u * m_U[k][i] + v * m_U[k][j]; + + } + + void copy_buffer_to_col_i_H(unsigned i) { + for (unsigned k = i; k < m_m; k++) { + m_H[k][i] = m_buffer[k]; + } + } + void copy_buffer_to_col_i_U(unsigned i) { + for (unsigned k = 0; k < m_n; k++) + m_U[k][i] = m_buffer[k]; + } + + // multiply by (a, b) + // (c, d) + // from the left where i and j are the modified columns + // the [i][i] = a, and [i][j] = b for the matrix we multiply by + + + void multiply_U_reverse_from_left_by(unsigned i, unsigned j, const mpq & a, const mpq & b, const mpq & c, const mpq d) { + // the new i-th row goes to the buffer + for (unsigned k = 0; k < m_n; k++) { + m_buffer[k] = a * m_U_reverse[i][k] + b * m_U_reverse[j][k]; + } + + // calculate the new j-th row in place + for (unsigned k = 0; k < m_n; k++) { + m_U_reverse[j][k] = c * m_U_reverse[i][k] + d * m_U_reverse[j][k]; + } + + // copy the buffer into i-th row + for (unsigned k = 0; k < m_n; k++) { + m_U_reverse[i][k] = m_buffer[k]; + } + } + + void handle_column_ij_in_row_i(unsigned i, unsigned j) { + lp_assert(is_correct_modulo()); + const mpq& aii = m_H[i][i]; + const mpq& aij = m_H[i][j]; + mpq p,q,r; + extended_gcd(aii, aij, r, p, q); + mpq aii_over_r = aii / r; + mpq aij_over_r = aij / r; + + + buffer_p_col_i_plus_q_col_j_H(p, i, q, j); + pivot_column_i_to_column_j_H(- aij_over_r, i, aii_over_r, j); + copy_buffer_to_col_i_H(i); + + + buffer_p_col_i_plus_q_col_j_U(p, i, q, j); + pivot_column_i_to_column_j_U(- aij_over_r, i, aii_over_r, j); + copy_buffer_to_col_i_U(i); + + // U was multiplied from the right by (p, - aij_over_r) + // (q, aii_over_r ) + // We need to multiply U_reverse by (aii_over_r, aij_over_r) + // (-q , p) + // from the left + + multiply_U_reverse_from_left_by(i, j, aii_over_r, aij_over_r, -q, p); + lp_assert(is_correct_modulo()); + } + + + void switch_sign_for_column(unsigned i) { + for (unsigned k = i; k < m_m; k++) + m_H[k][i].neg(); + for (unsigned k = 0; k < m_n; k++) + m_U[k][i].neg(); + + // switch sign for the i-th row in the reverse m_U_reverse + for (unsigned k = 0; k < m_n; k++) + m_U_reverse[i][k].neg(); + + } + + void process_row_column(unsigned i, unsigned j){ + if (is_zero(m_H[i][j])) + return; + handle_column_ij_in_row_i(i, j); + } + + void replace_column_j_by_j_minus_u_col_i_H(unsigned i, unsigned j, const mpq & u) { + lp_assert(j < i); + for (unsigned k = i; k < m_m; k++) { + m_H[k][j] -= u * m_H[k][i]; + } + } + void replace_column_j_by_j_minus_u_col_i_U(unsigned i, unsigned j, const mpq & u) { + + lp_assert(j < i); + for (unsigned k = 0; k < m_n; k++) { + m_U[k][j] -= u * m_U[k][i]; + } + // Here we multiply from m_U from the right by the matrix ( 1, 0) + // ( -u, 1). + // To adjust the reverse we multiply it from the left by (1, 0) + // (u, 1) + + for (unsigned k = 0; k < m_n; k++) { + m_U_reverse[i][k] += u * m_U_reverse[j][k]; + } + + + } + + void work_on_columns_less_than_i_in_the_triangle(unsigned i) { + const mpq & mii = m_H[i][i]; + lp_assert(is_pos(mii)); + for (unsigned j = 0; j < i; j++) { + const mpq & mij = m_H[i][j]; + if (!is_pos(mij) && - mij < mii) + continue; + mpq u = ceil(mij / mii); + replace_column_j_by_j_minus_u_col_i_H(i, j, u); + replace_column_j_by_j_minus_u_col_i_U(i, j, u); + } + } + + void process_row(unsigned i) { + + lp_assert(is_correct_modulo()); + for (unsigned j = i + 1; j < m_n; j++) { + process_row_column(i, j); + } + if (i >= m_n) { + lp_assert(m_H == m_A_orig * m_U); + return; + } + if (is_neg(m_H[i][i])) + switch_sign_for_column(i); + work_on_columns_less_than_i_in_the_triangle(i); + lp_assert(is_correct_modulo()); + } + + void calculate() { + for (unsigned i = 0; i < m_m; i++) { + process_row(i); + } + } + + void prepare_U_and_U_reverse() { + m_U = M(m_H.column_count()); + for (unsigned i = 0; i < m_U.column_count(); i++) + m_U[i][i] = 1; + + m_U_reverse = m_U; + + lp_assert(m_H == m_A_orig * m_U); + } + + bool row_is_correct_form(unsigned i) const { + if (i >= m_n) + return true; + const mpq& hii = m_H[i][i]; + if (is_neg(hii)) + return false; + for (unsigned j = 0; j < i; j++) { + const mpq & hij = m_H[i][j]; + if (is_pos(hij)) + return false; + if (- hij >= hii) + return false; + } + + return true; + } + + bool is_correct_form() const { + for (unsigned i = 0; i < m_m; i++) + if (!row_is_correct_form(i)) + return false; + return true; + } + + bool is_correct() const { + return m_H == m_A_orig * m_U && is_unit_matrix(m_U * m_U_reverse); + } + + bool is_correct_modulo() const { + return m_H.equal_modulo(m_A_orig * m_U, m_d) && is_unit_matrix(m_U * m_U_reverse); + } + + bool is_correct_final() const { + if (!is_correct()) { + std::cout << "m_H = "; m_H.print(std::cout, 17); + + std::cout << "\nm_A_orig * m_U = "; (m_A_orig * m_U).print(std::cout, 17); + + std::cout << "is_correct() does not hold" << std::endl; + return false; + } + if (!is_correct_form()) { + std::cout << "is_correct_form() does not hold" << std::endl; + return false; + } + return true; + } +public: + const M& H() const { return m_H;} + const M& U() const { return m_U;} + const M& U_reverse() const { return m_U_reverse; } +private: +#endif + void copy_buffer_to_col_i_W_modulo() { + for (unsigned k = m_i; k < m_m; k++) { + m_W[k][m_i] = m_buffer[k]; + } + } + + void replace_column_j_by_j_minus_u_col_i_W(unsigned j, const mpq & u) { + lp_assert(j < m_i); + for (unsigned k = m_i; k < m_m; k++) { + m_W[k][j] -= u * m_W[k][m_i]; + // m_W[k][j] = mod_R_balanced(m_W[k][j]); + } + } + + + + bool is_unit_matrix(const M& u) const { + unsigned m = u.row_count(); + unsigned n = u.column_count(); + if (m != n) return false; + for (unsigned i = 0; i < m; i ++) + for (unsigned j = 0; j < n; j++) { + if (i == j) { + if (one_of_type() != u[i][j]) + return false; + } else { + if (!is_zero(u[i][j])) + return false; + } + } + return true; + } + + + // follows Algorithm 2.4.8 of Henri Cohen's "A course on computational algebraic number theory", + // with some changes related to that we create a low triangle matrix + // with non-positive elements under the diagonal + void process_column_in_row_modulo() { + mpq& aii = m_W[m_i][m_i]; + const mpq& aij = m_W[m_i][m_j]; + mpq d, p,q; + hnf_calc::extended_gcd_minimal_uv(aii, aij, d, p, q); + if (is_zero(d)) + return; + mpq aii_over_d = aii / d; + mpq aij_over_d = aij / d; + buffer_p_col_i_plus_q_col_j_W_modulo(p, q); + pivot_column_i_to_column_j_W_modulo(- aij_over_d, aii_over_d); + copy_buffer_to_col_i_W_modulo(); + } + + void endl() const { + std::cout << std::endl; + } + + void fix_row_under_diagonal_W_modulo() { + mpq d, u, v; + if (is_zero(m_W[m_i][m_i])) { + m_W[m_i][m_i] = m_R; + u = one_of_type(); + d = m_R; + } else { + hnf_calc::extended_gcd_minimal_uv(m_W[m_i][m_i], m_R, d, u, v); + } + auto & mii = m_W[m_i][m_i]; + mii *= u; + mii = mod_R(mii); + if (is_zero(mii)) + mii = d; + + lp_assert(is_pos(mii)); + + // adjust column m_i + for (unsigned k = m_i + 1; k < m_m; k++) { + m_W[k][m_i] *= u; + m_W[k][m_i] = mod_R_balanced(m_W[k][m_i]); + } + + lp_assert(is_pos(mii)); + for (unsigned j = 0; j < m_i; j++) { + const mpq & mij = m_W[m_i][j]; + if (!is_pos(mij) && - mij < mii) + continue; + mpq q = ceil(mij / mii); + replace_column_j_by_j_minus_u_col_i_W(j, q); + } + } + + + void process_row_modulo() { + for (m_j = m_i + 1; m_j < m_n; m_j++) { + process_column_in_row_modulo(); + } + fix_row_under_diagonal_W_modulo(); + } + + void calculate_by_modulo() { + for (m_i = 0; m_i < m_m; m_i ++) { + process_row_modulo(); + lp_assert(is_pos(m_W[m_i][m_i])); + m_R /= m_W[m_i][m_i]; + lp_assert(is_int(m_R)); + m_half_R = floor(m_R / 2); + } + } + +public: + hnf(M & A, const mpq & d, unsigned r) : +#ifdef Z3DEBUG + m_H(A), +#endif + m_A_orig(A), + m_W(A), + m_buffer(std::max(A.row_count(), A.column_count())), + m_m(A.row_count()), + m_n(A.column_count()), + m_d(d), + m_r(r), + m_R(m_d), + m_half_R(floor(m_R / 2)) + { + lp_assert(m_m > 0 && m_n > 0); + if (is_zero(m_d)) + return; +#ifdef Z3DEBUG + prepare_U_and_U_reverse(); + calculate(); + lp_assert(is_correct_final()); +#endif + calculate_by_modulo(); +#ifdef Z3DEBUG + if (m_H != m_W) { + std::cout << "A = "; m_A_orig.print(std::cout, 4); endl(); + std::cout << "H = "; m_H.print(std::cout, 4); endl(); + std::cout << "W = "; m_W.print(std::cout, 4); endl(); + lp_assert(false); + } +#endif + } + + + +}; + +} diff --git a/src/util/lp/hnf_cutter.h b/src/util/lp/hnf_cutter.h new file mode 100644 index 000000000..c0071d8a0 --- /dev/null +++ b/src/util/lp/hnf_cutter.h @@ -0,0 +1,46 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + Lev Nachmanson (levnach) + Nikolaj Bjorner (nbjorner) + +Revision History: + + +--*/ +#pragma once +#include "util/lp/lar_term.h" +#include "util/lp/hnf.h" +#include "util/lp/general_matrix.h" +#include "util/lp/var_register.h" +namespace lp { +class hnf_cutter { + var_register m_var_register; + general_matrix m_A; + vector m_terms; +public: + void clear() { + m_A.clear(); + m_var_register.clear(); + m_terms.clear(); + } + void add_term_to_A_for_hnf(const lar_term* t, const mpq &) { + m_terms.push_back(t); + for (const auto &p : *t) { + m_var_register.register_user_var(p.var()); + } + } + void print(std::ostream & out) { + out << "terms = " << m_terms.size() << ", var = " << m_var_register.size() << std::endl; + } +}; +} diff --git a/src/util/lp/indexed_vector.cpp b/src/util/lp/indexed_vector.cpp index d198c68bd..348196195 100644 --- a/src/util/lp/indexed_vector.cpp +++ b/src/util/lp/indexed_vector.cpp @@ -50,3 +50,4 @@ template void lp::indexed_vector::resize(unsigned int); template void lp::print_vector< lp::mpq>(vector< lp::mpq> const &, std::basic_ostream > &); template void lp::print_vector >(vector> const&, std::ostream&); template void lp::indexed_vector >::erase_from_index(unsigned int); + diff --git a/src/util/lp/indexed_vector.h b/src/util/lp/indexed_vector.h index f17eb76ec..cd313f446 100644 --- a/src/util/lp/indexed_vector.h +++ b/src/util/lp/indexed_vector.h @@ -32,7 +32,7 @@ template void print_vector(const vector & t, std::ostream & out) template void print_vector(const buffer & t, std::ostream & out); template void print_sparse_vector(const vector & t, std::ostream & out); -void print_vector(const vector & t, std::ostream & out); +void print_vector_as_doubles(const vector & t, std::ostream & out); template class indexed_vector { public: diff --git a/src/util/lp/indexed_vector_def.h b/src/util/lp/indexed_vector_def.h index 796f3c956..e4a517b55 100644 --- a/src/util/lp/indexed_vector_def.h +++ b/src/util/lp/indexed_vector_def.h @@ -46,7 +46,7 @@ void print_sparse_vector(const vector & t, std::ostream & out) { out << std::endl; } -void print_vector(const vector & t, std::ostream & out) { +void print_vector_as_doubles(const vector & t, std::ostream & out) { for (unsigned i = 0; i < t.size(); i++) out << t[i].get_double() << std::setprecision(3) << " "; out << std::endl; diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index acfd45947..f5d06c816 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -5,7 +5,7 @@ #include "util/lp/int_solver.h" #include "util/lp/lar_solver.h" -#include "util/lp/cut_solver.h" +#include "util/lp/chase_cut_solver.h" #include "util/lp/lp_utils.h" #include namespace lp { @@ -414,7 +414,7 @@ unsigned int_solver::row_of_basic_column(unsigned j) const { // } -typedef cut_solver::monomial mono; +typedef chase_cut_solver::monomial mono; // it produces an inequality coeff*x <= rs template @@ -456,33 +456,33 @@ struct pivoted_rows_tracking_control { } }; -void int_solver::copy_explanations_from_cut_solver() { +void int_solver::copy_explanations_from_chase_cut_solver() { TRACE("propagate_and_backjump_step_int", - for (unsigned j: m_cut_solver.m_explanation) + for (unsigned j: m_chase_cut_solver.m_explanation) m_lar_solver->print_constraint(m_lar_solver->constraints()[j], tout);); - for (unsigned j : m_cut_solver.m_explanation) { + for (unsigned j : m_chase_cut_solver.m_explanation) { m_ex->push_justification(j); } - m_cut_solver.m_explanation.clear(); + m_chase_cut_solver.m_explanation.clear(); } -void int_solver::copy_values_from_cut_solver() { - for (unsigned j = 0; j < m_lar_solver->A_r().column_count() && j < m_cut_solver.number_of_vars(); j++) { - if (!m_cut_solver.var_is_active(j)) +void int_solver::copy_values_from_chase_cut_solver() { + for (unsigned j = 0; j < m_lar_solver->A_r().column_count() && j < m_chase_cut_solver.number_of_vars(); j++) { + if (!m_chase_cut_solver.var_is_active(j)) continue; if (!is_int(j)) { continue; } - m_lar_solver->m_mpq_lar_core_solver.m_r_x[j] = m_cut_solver.var_value(j); + m_lar_solver->m_mpq_lar_core_solver.m_r_x[j] = m_chase_cut_solver.var_value(j); lp_assert(m_lar_solver->column_value_is_int(j)); } } -void int_solver::catch_up_in_adding_constraints_to_cut_solver() { - lp_assert(m_cut_solver.number_of_asserts() <= m_lar_solver->constraints().size()); - for (unsigned j = m_cut_solver.number_of_asserts(); j < m_lar_solver->constraints().size(); j++) { - add_constraint_to_cut_solver(j, m_lar_solver->constraints()[j]); +void int_solver::catch_up_in_adding_constraints_to_chase_cut_solver() { + lp_assert(m_chase_cut_solver.number_of_asserts() <= m_lar_solver->constraints().size()); + for (unsigned j = m_chase_cut_solver.number_of_asserts(); j < m_lar_solver->constraints().size(); j++) { + add_constraint_to_chase_cut_solver(j, m_lar_solver->constraints()[j]); } } @@ -536,7 +536,7 @@ bool int_solver::tighten_terms_for_cube() { } bool int_solver::find_cube() { - if (m_branch_cut_counter % settings().m_int_find_cube_period != 0) + if (m_branch_cut_counter % settings().m_int_find_cube_period != 0) return false; settings().st().m_cube_calls++; @@ -557,13 +557,13 @@ bool int_solver::find_cube() { m_lar_solver->pop(); move_non_basic_columns_to_bounds(); find_feasible_solution(); - lp_assert(m_cut_solver.cancel() || is_feasible()); + lp_assert(m_chase_cut_solver.cancel() || is_feasible()); // it can happen that we found an integer solution here return !m_lar_solver->r_basis_has_inf_int(); } m_lar_solver->pop(); m_lar_solver->round_to_integer_solution(); - lp_assert(m_cut_solver.cancel() || is_feasible()); + lp_assert(m_chase_cut_solver.cancel() || is_feasible()); return true; } @@ -583,29 +583,29 @@ lia_move int_solver::run_gcd_test() { return lia_move::undef; } -lia_move int_solver::call_cut_solver() { - if ((m_branch_cut_counter) % settings().m_int_cut_solver_period != 0 || !all_columns_are_bounded()) +lia_move int_solver::call_chase_cut_solver() { + if ((m_branch_cut_counter) % settings().m_int_chase_cut_solver_period != 0 || !all_columns_are_bounded()) return lia_move::undef; - TRACE("check_main_int", tout<<"cut_solver";); - catch_up_in_adding_constraints_to_cut_solver(); - auto check_res = m_cut_solver.check(); - settings().st().m_cut_solver_calls++; + TRACE("check_main_int", tout<<"chase_cut_solver";); + catch_up_in_adding_constraints_to_chase_cut_solver(); + auto check_res = m_chase_cut_solver.check(); + settings().st().m_chase_cut_solver_calls++; switch (check_res) { - case cut_solver::lbool::l_false: - copy_explanations_from_cut_solver(); - settings().st().m_cut_solver_false++; + case chase_cut_solver::lbool::l_false: + copy_explanations_from_chase_cut_solver(); + settings().st().m_chase_cut_solver_false++; return lia_move::conflict; - case cut_solver::lbool::l_true: - settings().st().m_cut_solver_true++; - copy_values_from_cut_solver(); + case chase_cut_solver::lbool::l_true: + settings().st().m_chase_cut_solver_true++; + copy_values_from_chase_cut_solver(); lp_assert(m_lar_solver->all_constraints_hold()); return lia_move::sat; - case cut_solver::lbool::l_undef: - settings().st().m_cut_solver_undef++; - if (m_cut_solver.try_getting_cut(*m_t, *m_k, m_lar_solver->m_mpq_lar_core_solver.m_r_x)) { + case chase_cut_solver::lbool::l_undef: + settings().st().m_chase_cut_solver_undef++; + if (m_chase_cut_solver.try_getting_cut(*m_t, *m_k, m_lar_solver->m_mpq_lar_core_solver.m_r_x)) { m_lar_solver->subs_term_columns(*m_t); - TRACE("cut_solver_cuts", - tout<<"precut from cut_solver:"; m_lar_solver->print_term(*m_t, tout); tout << " <= " << *m_k << std::endl;); + TRACE("chase_cut_solver_cuts", + tout<<"precut from chase_cut_solver:"; m_lar_solver->print_term(*m_t, tout); tout << " <= " << *m_k << std::endl;); return lia_move::cut; } @@ -614,6 +614,63 @@ lia_move int_solver::call_cut_solver() { } } +lia_move int_solver::gomory_cut() { + TRACE("check_main_int", tout << "gomory";); + if (move_non_basic_columns_to_bounds()) { + lp_status st = m_lar_solver->find_feasible_solution(); + if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) { + TRACE("arith_int", tout << "give_up\n";); + return lia_move::undef; + } + } + int j = find_inf_int_base_column(); + if (j == -1) { + j = find_inf_int_nbasis_column(); + return j == -1? lia_move::sat : create_branch_on_column(j); + } + lia_move r = proceed_with_gomory_cut(j); + if (r != lia_move::undef) + return r; + return create_branch_on_column(j); +} + + +void int_solver::try_add_term_to_A_for_hnf(unsigned i) { + mpq rs; + const lar_term* t = m_lar_solver->terms()[i]; + for (const auto & p : *t) { + if (!is_int(p.var())) + return; // todo : the mix case! + } + if (!m_lar_solver->get_equality_for_term_on_corrent_x(i, rs)) + return; + m_hnf_cutter.add_term_to_A_for_hnf(t, rs); +} + +bool int_solver::hnf_matrix_is_empty() const { return true; } + +bool int_solver::prepare_matrix_A_for_hnf_cut() { + m_hnf_cutter.clear(); + for (unsigned i = 0; i < m_lar_solver->terms().size(); i++) + try_add_term_to_A_for_hnf(i); + m_hnf_cutter.print(std::cout); + return ! hnf_matrix_is_empty(); +} + + +lia_move int_solver::make_hnf_cut() { + if( !prepare_matrix_A_for_hnf_cut()) + return lia_move::undef; + return lia_move::undef; +} + +lia_move int_solver::hnf_cut() { + if ((m_branch_cut_counter) % settings().m_hnf_cut_period == 0) { + return make_hnf_cut(); + } + return lia_move::undef; +} + lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) { if (!has_inf_int()) return lia_move::sat; @@ -634,35 +691,17 @@ lia_move int_solver::check(lar_term& t, mpq& k, explanation& ex, bool & upper) { return lia_move::sat; } - lia_move r = call_cut_solver(); + lia_move r = call_chase_cut_solver(); + if (r != lia_move::undef) + return r; + + r = hnf_cut(); if (r != lia_move::undef) return r; if ((m_branch_cut_counter) % settings().m_int_gomory_cut_period == 0) { - TRACE("check_main_int", tout << "gomory";); - if (move_non_basic_columns_to_bounds()) { - lp_status st = m_lar_solver->find_feasible_solution(); - lp_assert(non_basic_columns_are_at_bounds()); - if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) { - TRACE("arith_int", tout << "give_up\n";); - return lia_move::undef; - } - } - int j = find_inf_int_base_column(); - if (j == -1) { - j = find_inf_int_nbasis_column(); - return j == -1? lia_move::sat : create_branch_on_column(j); - } - - TRACE("arith_int", tout << "j = " << j << " does not have an integer assignment: " << get_value(j) << "\n";); - - r = proceed_with_gomory_cut(j); - if (r != lia_move::undef) - return r; - return create_branch_on_column(j); + return gomory_cut(); } - - TRACE("check_main_int", tout << "branch"; ); int j = find_inf_int_base_column(); if (j == -1) { j = find_inf_int_nbasis_column(); @@ -969,7 +1008,7 @@ linear_combination_iterator * int_solver::get_column_iterator(unsigned j) { int_solver::int_solver(lar_solver* lar_slv) : m_lar_solver(lar_slv), m_branch_cut_counter(0), - m_cut_solver([this](unsigned j) {return m_lar_solver->get_column_name(j);}, + m_chase_cut_solver([this](unsigned j) {return m_lar_solver->get_column_name(j);}, [this](unsigned j, std::ostream &o) {m_lar_solver->print_constraint(j, o);}, [this]() {return m_lar_solver->A_r().column_count();}, [this](unsigned j) {return get_value(j);}, @@ -1328,21 +1367,21 @@ bool int_solver::is_term(unsigned j) const { return m_lar_solver->column_corresponds_to_term(j); } -void int_solver::add_constraint_to_cut_solver(unsigned ci, const lar_base_constraint * c) { +void int_solver::add_constraint_to_chase_cut_solver(unsigned ci, const lar_base_constraint * c) { vector coeffs; mpq rs; get_int_coeffs_from_constraint(c, coeffs, rs); - m_cut_solver.add_ineq(coeffs, -rs, ci); + m_chase_cut_solver.add_ineq(coeffs, -rs, ci); } void int_solver::pop(unsigned k) { - m_cut_solver.pop_trail(k); - while (m_cut_solver.number_of_asserts() > m_lar_solver->constraints().size()) - m_cut_solver.pop_last_assert(); - m_cut_solver.pop_constraints(); + m_chase_cut_solver.pop_trail(k); + while (m_chase_cut_solver.number_of_asserts() > m_lar_solver->constraints().size()) + m_chase_cut_solver.pop_last_assert(); + m_chase_cut_solver.pop_constraints(); } -void int_solver::push() { m_cut_solver.push(); } +void int_solver::push() { m_chase_cut_solver.push(); } unsigned int_solver::column_count() const { return m_lar_solver->column_count(); } diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h index dc3a64e58..7975025d8 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -22,8 +22,9 @@ Revision History: #include "util/lp/static_matrix.h" #include "util/lp/int_set.h" #include "util/lp/lar_term.h" -#include "util/lp/cut_solver.h" +#include "util/lp/chase_cut_solver.h" #include "util/lp/lar_constraints.h" +#include "util/lp/hnf_cutter.h" namespace lp { class lar_solver; @@ -54,11 +55,12 @@ public: // fields lar_solver * m_lar_solver; unsigned m_branch_cut_counter; - cut_solver m_cut_solver; + chase_cut_solver m_chase_cut_solver; lar_term* m_t; // the term to return in the cut mpq *m_k; // the right side of the cut explanation *m_ex; // the conflict explanation bool *m_upper; // we have a cut m_t*x <= k if m_upper is true nad m_t*x >= k otherwise + hnf_cutter m_hnf_cutter; // methods int_solver(lar_solver* lp); @@ -152,18 +154,18 @@ private: unsigned random(); bool has_inf_int() const; lia_move create_branch_on_column(int j); - void catch_up_in_adding_constraints_to_cut_solver(); + void catch_up_in_adding_constraints_to_chase_cut_solver(); public: template - void fill_cut_solver_vars(); + void fill_chase_cut_solver_vars(); template - void get_int_coeffs_from_constraint(const lar_base_constraint* c, vector& coeff, T & rs); + void get_int_coeffs_from_constraint(const lar_base_constraint* c, vector& coeff, T & rs); bool is_term(unsigned j) const; - void add_constraint_to_cut_solver(unsigned,const lar_base_constraint*); - void copy_explanations_from_cut_solver(); + void add_constraint_to_chase_cut_solver(unsigned,const lar_base_constraint*); + void copy_explanations_from_chase_cut_solver(); void pop(unsigned); void push(); - void copy_values_from_cut_solver(); + void copy_values_from_chase_cut_solver(); bool left_branch_is_more_narrow_than_right(unsigned); bool find_cube(); bool tighten_terms_for_cube(); @@ -174,6 +176,12 @@ public: void find_feasible_solution(); int find_inf_int_nbasis_column() const; lia_move run_gcd_test(); - lia_move call_cut_solver(); + lia_move call_chase_cut_solver(); + lia_move gomory_cut(); + lia_move hnf_cut(); + lia_move make_hnf_cut(); + bool prepare_matrix_A_for_hnf_cut(); + bool hnf_matrix_is_empty() const; + void try_add_term_to_A_for_hnf(unsigned term_index); }; } diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index b77a3dde8..c8de1807c 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -1478,7 +1478,7 @@ bool lar_solver::strategy_is_undecided() const { void lar_solver::catch_up_in_updating_int_solver() { for (unsigned i = 0; i < constraints().size(); i++) { - m_int_solver->add_constraint_to_cut_solver(i, constraints()[i]); + m_int_solver->add_constraint_to_chase_cut_solver(i, constraints()[i]); } } @@ -2217,6 +2217,40 @@ void lar_solver::round_to_integer_solution() { } } +bool lar_solver::get_equality_for_term_on_corrent_x(unsigned term_index, mpq & rs) const { + unsigned tj = term_index + m_terms_start_index; + auto it = m_ext_vars_to_columns.find(tj); + if (it == m_ext_vars_to_columns.end()) + return false; + unsigned j = it->second.internal_j(); + auto & slv = m_mpq_lar_core_solver.m_r_solver; + impq term_val; + bool term_val_ready = false; + if (slv.column_has_upper_bound(j)) { + const impq & b = slv.m_upper_bounds[j]; + lp_assert(is_zero(b.y) && is_int(b.x)); + term_val = terms()[term_index]->apply(m_mpq_lar_core_solver.m_r_x); + term_val_ready = true; + if (term_val == b) { + rs = b.x; + return true; + } + } + if (slv.column_has_lower_bound(j)) { + if (!term_val_ready) + term_val = terms()[term_index]->apply(m_mpq_lar_core_solver.m_r_x); + const impq & b = slv.m_lower_bounds[j]; + lp_assert(is_zero(b.y) && is_int(b.x)); + + if (term_val == b) { + rs = b.x; + return true; + } + } + return false; +} + + } // namespace lp diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index bb521b7af..304483817 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -579,5 +579,6 @@ public: unsigned column_count() const { return A_r().column_count(); } const vector & r_basis() const { return m_mpq_lar_core_solver.r_basis(); } const vector & r_nbasis() const { return m_mpq_lar_core_solver.r_nbasis(); } + bool get_equality_for_term_on_corrent_x(unsigned i, mpq &rs) const; }; } diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index 8fd8dc938..c43ca7142 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -69,7 +69,7 @@ public: lp_settings & m_settings; vector m_y; // the buffer for yB = cb // a device that is able to solve Bx=c, xB=d, and change the basis - lu * m_factorization; + lu> * m_factorization; const column_namer & m_column_names; indexed_vector m_w; // the vector featuring in 24.3 of the Chvatal book vector m_d; // the vector of reduced costs @@ -607,9 +607,9 @@ public: std:: cout << "m_d = " << m_d[j] << std::endl;*/ } - bool column_is_free(unsigned j) { return this->m_column_type[j] == free; } + bool column_is_free(unsigned j) const { return this->m_column_type[j] == free; } - bool column_has_upper_bound(unsigned j) { + bool column_has_upper_bound(unsigned j) const { switch(m_column_types[j]) { case column_type::free_column: case column_type::lower_bound: @@ -628,7 +628,7 @@ public: return true; } - bool column_has_lower_bound(unsigned j) { + bool column_has_lower_bound(unsigned j) const { switch(m_column_types[j]) { case column_type::free_column: case column_type::upper_bound: diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h index 0af8f0af6..26fc550b3 100644 --- a/src/util/lp/lp_primal_core_solver.h +++ b/src/util/lp/lp_primal_core_solver.h @@ -419,7 +419,7 @@ public: // returns the number of iterations unsigned solve(); - lu * factorization() {return this->m_factorization;} + lu> * factorization() {return this->m_factorization;} void delete_factorization(); diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h index 7103c5c20..33df5a180 100644 --- a/src/util/lp/lp_settings.h +++ b/src/util/lp/lp_settings.h @@ -102,10 +102,10 @@ struct stats { unsigned m_need_to_solve_inf; unsigned m_max_cols; unsigned m_max_rows; - unsigned m_cut_solver_calls; - unsigned m_cut_solver_true; - unsigned m_cut_solver_false; - unsigned m_cut_solver_undef; + unsigned m_chase_cut_solver_calls; + unsigned m_chase_cut_solver_true; + unsigned m_chase_cut_solver_false; + unsigned m_chase_cut_solver_undef; unsigned m_gcd_calls; unsigned m_gcd_conflicts; unsigned m_cube_calls; @@ -232,11 +232,12 @@ public: backup_costs(true), column_number_threshold_for_using_lu_in_lar_solver(4000), m_int_gomory_cut_period(4), - m_int_cut_solver_period(8), + m_int_chase_cut_solver_period(8), m_int_find_cube_period(4), m_int_cuts_etc_period(4), + m_hnf_cut_period(4), m_int_run_gcd_test(true), - m_cut_solver_cycle_on_var(10), + m_chase_cut_solver_cycle_on_var(10), m_int_pivot_fixed_vars_from_basis(false), m_int_patch_only_integer_values(true) {} @@ -346,11 +347,12 @@ public: bool backup_costs; unsigned column_number_threshold_for_using_lu_in_lar_solver; unsigned m_int_gomory_cut_period; - unsigned m_int_cut_solver_period; + unsigned m_int_chase_cut_solver_period; unsigned m_int_find_cube_period; unsigned m_int_cuts_etc_period; + unsigned m_hnf_cut_period; bool m_int_run_gcd_test; - unsigned m_cut_solver_cycle_on_var; + unsigned m_chase_cut_solver_cycle_on_var; bool m_int_pivot_fixed_vars_from_basis; bool m_int_patch_only_integer_values; }; // end of lp_settings class diff --git a/src/util/lp/lu.cpp b/src/util/lp/lu.cpp index 4e543a73b..e6df10908 100644 --- a/src/util/lp/lu.cpp +++ b/src/util/lp/lu.cpp @@ -23,56 +23,62 @@ Revision History: #include "util/vector.h" #include "util/debug.h" #include "util/lp/lu_def.h" -template double lp::dot_product(vector const&, vector const&); -template lp::lu::lu(lp::static_matrix const&, vector&, lp::lp_settings&); -template void lp::lu::push_matrix_to_tail(lp::tail_matrix*); -template void lp::lu::replace_column(double, lp::indexed_vector&, unsigned); -template void lp::lu::solve_Bd(unsigned int, lp::indexed_vector&, lp::indexed_vector&); -template lp::lu::~lu(); -template void lp::lu::push_matrix_to_tail(lp::tail_matrix*); -template void lp::lu::solve_Bd(unsigned int, lp::indexed_vector&, lp::indexed_vector&); -template lp::lu::~lu(); -template void lp::lu >::push_matrix_to_tail(lp::tail_matrix >*); -template void lp::lu >::solve_Bd(unsigned int, lp::indexed_vector&, lp::indexed_vector&); -template lp::lu >::~lu(); -template lp::mpq lp::dot_product(vector const&, vector const&); -template void lp::init_factorization(lp::lu*&, lp::static_matrix&, vector&, lp::lp_settings&); -template void lp::init_factorization(lp::lu*&, lp::static_matrix&, vector&, lp::lp_settings&); -template void lp::init_factorization >(lp::lu >*&, lp::static_matrix >&, vector&, lp::lp_settings&); +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 void lu>::solve_Bd(unsigned int, indexed_vector&, indexed_vector&); +template lu>::~lu(); +template void lu>::push_matrix_to_tail(tail_matrix*); +template void lu>::solve_Bd(unsigned int, indexed_vector&, indexed_vector&); +template lu>::~lu(); +template void lu>::push_matrix_to_tail(tail_matrix*); +template void lu>::solve_Bd(unsigned int, indexed_vector&, indexed_vector&); +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&); #ifdef Z3DEBUG -template void lp::print_matrix(lp::sparse_matrix&, std::ostream & out); -template void lp::print_matrix(lp::static_matrix&, std::ostream&); -template void lp::print_matrix >(lp::static_matrix >&, std::ostream&); -template void lp::print_matrix(lp::static_matrix&, std::ostream & out); -template bool lp::lu::is_correct(const vector& basis); -template bool lp::lu >::is_correct( vector const &); -template lp::dense_matrix lp::get_B(lp::lu&, const vector& basis); -template lp::dense_matrix lp::get_B(lp::lu&, vector const&); +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); +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 lp::lu::pivot_the_row(int); // NOLINT -template void lp::lu::init_vector_w(unsigned int, lp::indexed_vector&); -template void lp::lu::solve_By(vector&); -template void lp::lu::solve_By_when_y_is_ready_for_X(vector&); -template void lp::lu::solve_yB_with_error_check(vector&, const vector& basis); -template void lp::lu::solve_yB_with_error_check_indexed(lp::indexed_vector&, vector const&, const vector & basis, const lp_settings&); -template void lp::lu::replace_column(lp::mpq, lp::indexed_vector&, unsigned); -template void lp::lu::solve_By(vector&); -template void lp::lu::solve_By_when_y_is_ready_for_X(vector&); -template void lp::lu::solve_yB_with_error_check(vector&, const vector& basis); -template void lp::lu::solve_yB_with_error_check_indexed(lp::indexed_vector&, vector< int > const&, const vector & basis, const lp_settings&); -template void lp::lu >::solve_yB_with_error_check_indexed(lp::indexed_vector&, vector< int > const&, const vector & basis, const lp_settings&); -template void lp::lu >::init_vector_w(unsigned int, lp::indexed_vector&); -template void lp::lu >::replace_column(lp::mpq, lp::indexed_vector&, unsigned); -template void lp::lu >::solve_Bd_faster(unsigned int, lp::indexed_vector&); -template void lp::lu >::solve_By(vector >&); -template void lp::lu >::solve_By_when_y_is_ready_for_X(vector >&); -template void lp::lu >::solve_yB_with_error_check(vector&, const vector& basis); -template void lp::lu::solve_By(lp::indexed_vector&); -template void lp::lu::solve_By(lp::indexed_vector&); -template void lp::lu::solve_yB_indexed(lp::indexed_vector&); -template void lp::lu::solve_yB_indexed(lp::indexed_vector&); -template void lp::lu >::solve_yB_indexed(lp::indexed_vector&); -template void lp::lu::solve_By_for_T_indexed_only(lp::indexed_vector&, lp::lp_settings const&); -template void lp::lu::solve_By_for_T_indexed_only(lp::indexed_vector&, lp::lp_settings const&); +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/util/lp/lu.h b/src/util/lp/lu.h index 9713936a1..820786d87 100644 --- a/src/util/lp/lu.h +++ b/src/util/lp/lu.h @@ -8,7 +8,12 @@ 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) @@ -24,7 +29,7 @@ Revision History: #include "util/debug.h" #include #include -#include "util/lp/sparse_matrix.h" +#include "util/lp/square_sparse_matrix.h" #include "util/lp/static_matrix.h" #include #include "util/lp/numeric_pair.h" @@ -36,13 +41,11 @@ Revision History: namespace lp { #ifdef Z3DEBUG template // print the nr x nc submatrix at the top left corner -void print_submatrix(sparse_matrix & m, unsigned mr, unsigned nc); +void print_submatrix(square_sparse_matrix & m, unsigned mr, unsigned nc); -template -void print_matrix(static_matrix &m, std::ostream & out); +template +void print_matrix(M &m, std::ostream & out); -template -void print_matrix(sparse_matrix& m, std::ostream & out); #endif template @@ -123,17 +126,20 @@ 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 +template class lu { LU_status m_status; public: + typedef typename M::coefftype T; + typedef typename M::argtype X; + // the fields unsigned m_dim; - static_matrix const &m_A; + const M & m_A; permutation_matrix m_Q; permutation_matrix m_R; permutation_matrix m_r_wave; - sparse_matrix m_U; + square_sparse_matrix m_U; square_dense_submatrix* m_dense_LU; vector *> m_tail; @@ -147,10 +153,11 @@ public: // 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(static_matrix const & A, + lu(const M & A, vector& basis, lp_settings & settings); - void debug_test_of_basis(static_matrix const & A, vector & basis); + 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); @@ -222,7 +229,7 @@ public: 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, sparse_matrix& copy_of_U); + 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); @@ -245,7 +252,7 @@ public: } T B_(unsigned i, unsigned j, const vector& basis) { - return m_A.get_elem(i, basis[j]); + return m_A[i][basis[j]]; } unsigned dimension() { return m_dim; } @@ -261,11 +268,13 @@ public: void process_column(int j); bool is_correct(const vector& basis); + bool is_correct(); #ifdef Z3DEBUG dense_matrix tail_product(); dense_matrix get_left_side(const vector& basis); + dense_matrix get_left_side(); dense_matrix get_right_side(); #endif @@ -364,11 +373,14 @@ public: }; // end of lu -template -void init_factorization(lu* & factorization, static_matrix & m_A, vector & m_basis, lp_settings &m_settings); +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, const vector& basis); + +template +dense_matrix get_B(lu& f); #endif } diff --git a/src/util/lp/lu_def.h b/src/util/lp/lu_def.h index f8abdb2a5..89c959e4d 100644 --- a/src/util/lp/lu_def.h +++ b/src/util/lp/lu_def.h @@ -26,8 +26,8 @@ Revision History: #include "util/lp/lu.h" namespace lp { #ifdef Z3DEBUG -template // print the nr x nc submatrix at the top left corner -void print_submatrix(sparse_matrix & m, unsigned mr, unsigned nc, std::ostream & out) { +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++) { @@ -44,15 +44,14 @@ void print_submatrix(sparse_matrix & m, unsigned mr, unsigned nc, std::ost print_matrix_with_widths(A, widths, out); } -template -void print_matrix(static_matrix &m, std::ostream & out) { +template +void print_matrix(M &m, std::ostream & out) { vector> A; vector widths; - std::set> domain = m.get_domain(); 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(static_cast(m(i, j)))); + A[i].push_back(T_to_string(m[i][j])); } } @@ -63,23 +62,6 @@ void print_matrix(static_matrix &m, std::ostream & out) { print_matrix_with_widths(A, widths, out); } -template -void print_matrix(sparse_matrix& 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(static_cast(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); -} #endif @@ -122,10 +104,10 @@ void one_elem_on_diag::apply_from_left_to_T(indexed_vector & w, lp_sett // 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(static_matrix const & A, - vector& basis, - lp_settings & settings): +template +lu::lu(const M& A, + vector& basis, + lp_settings & settings): m_status(LU_status::OK), m_dim(A.row_count()), m_A(A), @@ -147,8 +129,30 @@ lu::lu(static_matrix const & A, // lp_assert(check_correctness()); #endif } -template -void lu::debug_test_of_basis(static_matrix const & A, vector & basis) { +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()); + std::cout << "m_U = \n"; + print_matrix(&m_U, std::cout); + 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()); @@ -157,22 +161,22 @@ void lu::debug_test_of_basis(static_matrix const & A, vector - void lu::solve_By(indexed_vector & y) { +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) { +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) { +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 @@ -189,8 +193,8 @@ void lu::solve_By_when_y_is_ready_for_X(vector & y) { } } -template -void lu::solve_By_when_y_is_ready_for_T(vector & y, vector & index) { +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 @@ -214,8 +218,8 @@ void lu::solve_By_when_y_is_ready_for_T(vector & y, vector & } } -template -void lu::solve_By_for_T_indexed_only(indexed_vector & y, const lp_settings & settings) { +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); @@ -226,8 +230,8 @@ void lu::solve_By_for_T_indexed_only(indexed_vector & y, const lp_setti m_R.apply_reverse_from_left(y); // see 24.3 from Chvatal } -template -void lu::print_matrix_compact(std::ostream & f) { +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; @@ -241,8 +245,8 @@ void lu::print_matrix_compact(std::ostream & f) { } f << "matrix_end" << std::endl; } -template -void lu::print(indexed_vector & w, const vector& basis) { +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); @@ -256,8 +260,8 @@ void lu::print(indexed_vector & w, const vector& basis) { print_indexed_vector(w, f); f.close(); } -template -void lu::solve_Bd(unsigned a_column, indexed_vector & d, indexed_vector & w) { +template +void lu< M>::solve_Bd(unsigned a_column, indexed_vector & d, indexed_vector & w) { init_vector_w(a_column, w); if (w.m_index.size() * ratio_of_index_size_to_all_size() < d.m_data.size()) { // this const might need some tuning @@ -270,14 +274,14 @@ void lu::solve_Bd(unsigned a_column, indexed_vector & d, indexed_vector } } -template -void lu::solve_Bd_faster(unsigned a_column, indexed_vector & d) { // puts the a_column into d +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::solve_yB(vector& y) { +template +void lu< M>::solve_yB(vector& y) { // first solve yU = cb*R(-1) m_R.apply_reverse_from_right_to_T(y); // got y = cb*R(-1) m_U.solve_y_U(y); // got y*U=cb*R(-1) @@ -290,8 +294,8 @@ void lu::solve_yB(vector& y) { } } -template -void lu::solve_yB_indexed(indexed_vector& y) { +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) @@ -309,15 +313,15 @@ void lu::solve_yB_indexed(indexed_vector& y) { } } -template -void lu::add_delta_to_solution(const vector& yc, vector& y){ +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::add_delta_to_solution_indexed(indexed_vector& y) { +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()); @@ -344,16 +348,16 @@ void lu::add_delta_to_solution_indexed(indexed_vector& y) { lp_assert(y.is_OK()); } -template -void lu::find_error_of_yB(vector& yc, const vector& y, const vector& m_basis) { +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::find_error_of_yB_indexed(const indexed_vector& y, const vector& heading, const lp_settings& settings) { +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; @@ -423,8 +427,8 @@ void lu::find_error_of_yB_indexed(const indexed_vector& y, const vector // solves y*B = y // y is the input -template -void lu::solve_yB_with_error_check_indexed(indexed_vector & y, const vector& heading, const vector & basis, const lp_settings & settings) { +template +void lu< M>::solve_yB_with_error_check_indexed(indexed_vector & y, const vector& heading, const vector & basis, const lp_settings & settings) { if (numeric_traits::precise()) { if (y.m_index.size() * ratio_of_index_size_to_all_size() * 3 < m_A.column_count()) { solve_yB_indexed(y); @@ -461,8 +465,8 @@ void lu::solve_yB_with_error_check_indexed(indexed_vector & y, const ve // solves y*B = y // y is the input -template -void lu::solve_yB_with_error_check(vector & y, const vector& basis) { +template +void lu< M>::solve_yB_with_error_check(vector & y, const vector& basis) { if (numeric_traits::precise()) { solve_yB(y); return; @@ -475,8 +479,8 @@ void lu::solve_yB_with_error_check(vector & y, const vector& add_delta_to_solution(yc, y); m_y_copy.clear_all(); } -template -void lu::apply_Q_R_to_U(permutation_matrix & r_wave) { +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); } @@ -488,62 +492,62 @@ void lu::apply_Q_R_to_U(permutation_matrix & r_wave) { // solving Bd = a ( to find the column d of B^{-1} A_N corresponding to the entering // variable -template -lu::~lu(){ +template +lu< M>::~lu(){ for (auto t : m_tail) { delete t; } } -template -void lu::init_vector_y(vector & y) { +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::perform_transformations_on_w(indexed_vector& w) { +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::init_vector_w(unsigned entering, indexed_vector & w) { +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::apply_lp_list_to_w(indexed_vector & 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::apply_lp_list_to_y(vector& y) { +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::swap_rows(int j, int k) { +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::swap_columns(int j, int pivot_column) { +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::pivot_the_row(int row) { +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; @@ -560,8 +564,8 @@ bool lu::pivot_the_row(int row) { return true; } // we're processing the column j now -template -eta_matrix * lu::get_eta_matrix_for_pivot(unsigned j) { +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); @@ -569,16 +573,16 @@ eta_matrix * lu::get_eta_matrix_for_pivot(unsigned j) { return ret; } // we're processing the column j now -template -eta_matrix * lu::get_eta_matrix_for_pivot(unsigned j, sparse_matrix& copy_of_U) { +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, +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); @@ -586,15 +590,15 @@ unsigned lu::transform_U_to_V_by_replacing_column(indexed_vector & w, } #ifdef Z3DEBUG -template -void lu::check_vector_w(unsigned entering) { +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) { +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); @@ -602,8 +606,8 @@ void lu::check_apply_matrix_to_vector(matrix *lp, T *w) { } } -template -void lu::check_apply_lp_lists_to_w(T * 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); } @@ -615,8 +619,8 @@ void lu::check_apply_lp_lists_to_w(T * w) { } #endif -template -void lu::process_column(int j) { +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) { @@ -637,8 +641,8 @@ void lu::process_column(int j) { m_failure = true; } } -template -bool lu::is_correct(const vector& basis) { +template +bool lu::is_correct(const vector& basis) { #ifdef Z3DEBUG if (get_status() != LU_status::OK) { return false; @@ -651,10 +655,27 @@ bool lu::is_correct(const vector& basis) { #endif } +template +bool lu::is_correct() { +#ifdef Z3DEBUG + if (get_status() != LU_status::OK) { + std::cout << " status" << std::endl; + return false; + } + dense_matrix left_side = get_left_side(); + std::cout << "ls = "; print_matrix(left_side, std::cout); + dense_matrix right_side = get_right_side(); + std::cout << "rs = "; print_matrix(right_side, std::cout); + return left_side == right_side; +#else + return true; +#endif +} + #ifdef Z3DEBUG -template -dense_matrix lu::tail_product() { +template +dense_matrix lu::tail_product() { lp_assert(tail_size() > 0); dense_matrix left_side = permutation_matrix(m_dim); for (unsigned i = 0; i < tail_size(); i++) { @@ -665,8 +686,8 @@ dense_matrix lu::tail_product() { } return left_side; } -template -dense_matrix lu::get_left_side(const vector& basis) { +template +dense_matrix lu::get_left_side(const vector& basis) { dense_matrix left_side = get_B(*this, basis); for (unsigned i = 0; i < tail_size(); i++) { matrix* lp = get_lp_matrix(i); @@ -676,8 +697,19 @@ dense_matrix lu::get_left_side(const vector& basis) { } return left_side; } -template -dense_matrix lu::get_right_side() { +template +dense_matrix lu::get_left_side() { + dense_matrix left_side = get_B(*this); + for (unsigned i = 0; i < tail_size(); i++) { + matrix* lp = get_lp_matrix(i); + lp->set_number_of_rows(m_dim); + lp->set_number_of_columns(m_dim); + left_side = ((*lp) * left_side); + } + return left_side; +} +template +dense_matrix lu::get_right_side() { auto ret = U() * R(); ret = Q() * ret; return ret; @@ -685,8 +717,8 @@ dense_matrix lu::get_right_side() { #endif // needed for debugging purposes -template -void lu::copy_w(T *buffer, indexed_vector & w) { +template +void lu::copy_w(T *buffer, indexed_vector & w) { unsigned i = m_dim; while (i--) { buffer[i] = w[i]; @@ -694,15 +726,15 @@ void lu::copy_w(T *buffer, indexed_vector & w) { } // needed for debugging purposes -template -void lu::restore_w(T *buffer, indexed_vector & w) { +template +void lu::restore_w(T *buffer, indexed_vector & w) { unsigned i = m_dim; while (i--) { w[i] = buffer[i]; } } -template -bool lu::all_columns_and_rows_are_active() { +template +bool lu::all_columns_and_rows_are_active() { unsigned i = m_dim; while (i--) { lp_assert(m_U.col_is_active(i)); @@ -710,8 +742,8 @@ bool lu::all_columns_and_rows_are_active() { } return true; } -template -bool lu::too_dense(unsigned j) const { +template +bool lu::too_dense(unsigned j) const { unsigned r = m_dim - j; if (r < 5) return false; @@ -720,8 +752,8 @@ bool lu::too_dense(unsigned j) const { // return r * r * m_settings.density_threshold <= m_U.get_number_of_nonzeroes_below_row(j); return r * r * m_settings.density_threshold <= m_U.get_n_of_active_elems(); } -template -void lu::pivot_in_dense_mode(unsigned i) { +template +void lu::pivot_in_dense_mode(unsigned i) { int j = m_dense_LU->find_pivot_column_in_row(i); if (j == -1) { m_failure = true; @@ -733,8 +765,8 @@ void lu::pivot_in_dense_mode(unsigned i) { } m_dense_LU->pivot(i, m_settings); } -template -void lu::create_initial_factorization(){ +template +void lu::create_initial_factorization(){ m_U.prepare_for_factorization(); unsigned j; for (j = 0; j < m_dim; j++) { @@ -771,8 +803,8 @@ void lu::create_initial_factorization(){ // lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings)); } -template -void lu::calculate_r_wave_and_update_U(unsigned bump_start, unsigned bump_end, permutation_matrix & r_wave) { +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; @@ -790,8 +822,8 @@ void lu::calculate_r_wave_and_update_U(unsigned bump_start, unsigned bump_ 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) { +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; @@ -805,8 +837,8 @@ void lu::scan_last_row_to_work_vector(unsigned lowest_row_of_the_bump) { } } -template -void lu::pivot_and_solve_the_system(unsigned replaced_column, unsigned lowest_row_of_the_bump) { +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++) { @@ -846,8 +878,8 @@ void lu::pivot_and_solve_the_system(unsigned replaced_column, unsigned low } // see Achim Koberstein's thesis page 58, but here we solve the system and pivot to the last // row at the same time -template -row_eta_matrix *lu::get_row_eta_matrix_and_set_row_vector(unsigned replaced_column, unsigned lowest_row_of_the_bump, const T & pivot_elem_for_checking) { +template +row_eta_matrix *lu::get_row_eta_matrix_and_set_row_vector(unsigned replaced_column, unsigned lowest_row_of_the_bump, const T & pivot_elem_for_checking) { if (replaced_column == lowest_row_of_the_bump) return nullptr; scan_last_row_to_work_vector(lowest_row_of_the_bump); pivot_and_solve_the_system(replaced_column, lowest_row_of_the_bump); @@ -861,9 +893,9 @@ row_eta_matrix *lu::get_row_eta_matrix_and_set_row_vector(unsigned r } } #ifdef Z3DEBUG - auto ret = new row_eta_matrix(replaced_column, lowest_row_of_the_bump, m_dim); + auto ret = new row_eta_matrix(replaced_column, lowest_row_of_the_bump, m_dim); #else - auto ret = new row_eta_matrix(replaced_column, lowest_row_of_the_bump); + auto ret = new row_eta_matrix(replaced_column, lowest_row_of_the_bump); #endif for (auto j : m_row_eta_work_vector.m_index) { @@ -880,8 +912,8 @@ row_eta_matrix *lu::get_row_eta_matrix_and_set_row_vector(unsigned r return ret; } -template -void lu::replace_column(T pivot_elem_for_checking, indexed_vector & w, unsigned leaving_column_of_U){ +template +void lu::replace_column(T pivot_elem_for_checking, indexed_vector & w, unsigned leaving_column_of_U){ m_refactor_counter++; unsigned replaced_column = transform_U_to_V_by_replacing_column( w, leaving_column_of_U); unsigned lowest_row_of_the_bump = m_U.lowest_row_in_column(replaced_column); @@ -903,8 +935,8 @@ void lu::replace_column(T pivot_elem_for_checking, indexed_vector & w, // lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings)); // lp_assert(w.is_OK() && m_row_eta_work_vector.is_OK()); } -template -void lu::calculate_Lwave_Pwave_for_bump(unsigned replaced_column, unsigned lowest_row_of_the_bump){ +template +void lu::calculate_Lwave_Pwave_for_bump(unsigned replaced_column, unsigned lowest_row_of_the_bump){ T diagonal_elem; if (replaced_column < lowest_row_of_the_bump) { diagonal_elem = m_row_eta_work_vector[lowest_row_of_the_bump]; @@ -922,8 +954,8 @@ void lu::calculate_Lwave_Pwave_for_bump(unsigned replaced_column, unsigned // 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) { +template +void lu::calculate_Lwave_Pwave_for_last_row(unsigned lowest_row_of_the_bump, T diagonal_element) { auto l = new one_elem_on_diag(lowest_row_of_the_bump, diagonal_element); #ifdef Z3DEBUG l->set_number_of_columns(m_dim); @@ -933,26 +965,35 @@ void lu::calculate_Lwave_Pwave_for_last_row(unsigned lowest_row_of_the_bum l->conjugate_by_permutation(m_Q); } -template -void init_factorization(lu* & factorization, static_matrix & m_A, vector & m_basis, lp_settings &m_settings) { +template +void init_factorization(lu* & factorization, M & m_A, vector & m_basis, lp_settings &m_settings) { if (factorization != nullptr) delete factorization; - factorization = new lu(m_A, m_basis, m_settings); + factorization = new lu(m_A, m_basis, m_settings); // if (factorization->get_status() != LU_status::OK) // LP_OUT(m_settings, "failing in init_factorization" << std::endl); } #ifdef Z3DEBUG -template -dense_matrix get_B(lu& f, const vector& basis) { +template +dense_matrix get_B(lu& f, const vector& basis) { lp_assert(basis.size() == f.dimension()); lp_assert(basis.size() == f.m_U.dimension()); - dense_matrix B(f.dimension(), f.dimension()); + dense_matrix B(f.dimension(), f.dimension()); for (unsigned i = 0; i < f.dimension(); i++) for (unsigned j = 0; j < f.dimension(); j++) B.set_elem(i, j, f.B_(i, j, basis)); return B; } +template +dense_matrix get_B(lu& f) { + dense_matrix B(f.dimension(), f.dimension()); + for (unsigned i = 0; i < f.dimension(); i++) + for (unsigned j = 0; j < f.dimension(); j++) + B.set_elem(i, j, f.m_A[i][j]); + + return B; +} #endif } diff --git a/src/util/lp/matrix.h b/src/util/lp/matrix.h index f6374756f..063513287 100644 --- a/src/util/lp/matrix.h +++ b/src/util/lp/matrix.h @@ -49,11 +49,25 @@ void apply_to_vector(matrix & m, T * w); unsigned get_width_of_column(unsigned j, vector> & A); -void print_matrix_with_widths(vector> & A, vector & ws, std::ostream & out); -void print_string_matrix(vector> & A); +void print_matrix_with_widths(vector> & A, vector & ws, std::ostream & out, unsigned blanks = 0); +void print_string_matrix(vector> & A, std::ostream &, unsigned blanks_in_front = 0); template void print_matrix(matrix const * m, std::ostream & out); + +template +void print_matrix(const vector> & A, std::ostream & out, unsigned blanks_in_front = 0) { + vector> s(A.size()); + for (unsigned i = 0; i < A.size(); i++) { + for (const auto & v : A[i]) { + s[i].push_back(T_to_string(v)); + } + } + + print_string_matrix(s, out, blanks_in_front); +} + + } #endif diff --git a/src/util/lp/matrix_def.h b/src/util/lp/matrix_def.h index 6eb82a9cc..ae5f05ad1 100644 --- a/src/util/lp/matrix_def.h +++ b/src/util/lp/matrix_def.h @@ -82,9 +82,11 @@ unsigned get_width_of_column(unsigned j, vector> & A) { return r; } -void print_matrix_with_widths(vector> & A, vector & ws, std::ostream & out) { +void print_matrix_with_widths(vector> & A, vector & ws, std::ostream & out, unsigned blanks_in_front) { for (unsigned i = 0; i < A.size(); i++) { for (unsigned j = 0; j < static_cast(A[i].size()); j++) { + if (i != 0 && j == 0) + print_blanks(blanks_in_front, out); print_blanks(ws[j] - static_cast(A[i][j].size()), out); out << A[i][j] << " "; } @@ -92,7 +94,7 @@ void print_matrix_with_widths(vector> & A, vector } } -void print_string_matrix(vector> & A, std::ostream & out) { +void print_string_matrix(vector> & A, std::ostream & out, unsigned blanks_in_front) { vector widths; if (A.size() > 0) @@ -100,10 +102,23 @@ void print_string_matrix(vector> & A, std::ostream & out) { widths.push_back(get_width_of_column(j, A)); } - print_matrix_with_widths(A, widths, out); + print_matrix_with_widths(A, widths, out, blanks_in_front); out << std::endl; } +template +void print_matrix(vector> & A, std::ostream & out, unsigned blanks_in_front = 0) { + vector> s(A.size()); + for (unsigned i = 0; i < A.size(); i++) { + for (const auto & v : A[i]) { + s[i].push_back(T_to_string(v)); + } + } + + print_string_matrix(s, out, blanks_in_front); +} + + template void print_matrix(matrix const * m, std::ostream & out) { vector> A(m->row_count()); diff --git a/src/util/lp/permutation_matrix.h b/src/util/lp/permutation_matrix.h index 77fb14cc3..d1422ef86 100644 --- a/src/util/lp/permutation_matrix.h +++ b/src/util/lp/permutation_matrix.h @@ -132,8 +132,6 @@ class permutation_matrix : public tail_matrix { unsigned size() const { return static_cast(m_rev.size()); } - unsigned * values() const { return m_permutation.c_ptr(); } - void resize(unsigned size) { unsigned old_size = m_permutation.size(); m_permutation.resize(size); diff --git a/src/util/lp/sparse_matrix.cpp b/src/util/lp/sparse_matrix.cpp deleted file mode 100644 index 4ce8e64ef..000000000 --- a/src/util/lp/sparse_matrix.cpp +++ /dev/null @@ -1,116 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#include -#include "util/vector.h" -#include "util/lp/lp_settings.h" -#include "util/lp/lu.h" -#include "util/lp/sparse_matrix_def.h" -#include "util/lp/dense_matrix.h" -namespace lp { -template double sparse_matrix::dot_product_with_row(unsigned int, vector const&) const; -template void sparse_matrix::add_new_element(unsigned int, unsigned int, const double&); -template void sparse_matrix::divide_row_by_constant(unsigned int, const double&, lp_settings&); -template bool sparse_matrix::fill_eta_matrix(unsigned int, eta_matrix**); -template const double & sparse_matrix::get(unsigned int, unsigned int) const; -template unsigned sparse_matrix::get_number_of_nonzeroes() const; -template bool sparse_matrix::get_pivot_for_column(unsigned int&, unsigned int&, int, unsigned int); -template unsigned sparse_matrix::lowest_row_in_column(unsigned int); -template bool sparse_matrix::pivot_row_to_row(unsigned int, const double&, unsigned int, lp_settings&); -template bool sparse_matrix::pivot_with_eta(unsigned int, eta_matrix*, lp_settings&); -template void sparse_matrix::prepare_for_factorization(); -template void sparse_matrix::remove_element(vector >&, indexed_value&); -template void sparse_matrix::replace_column(unsigned int, indexed_vector&, lp_settings&); -template void sparse_matrix::set(unsigned int, unsigned int, double); -template void sparse_matrix::set_max_in_row(vector >&); -template bool sparse_matrix::set_row_from_work_vector_and_clean_work_vector_not_adjusted(unsigned int, indexed_vector&, lp_settings&); -template bool sparse_matrix::shorten_active_matrix(unsigned int, eta_matrix*); -template void sparse_matrix::solve_y_U(vector&) const; -template sparse_matrix::sparse_matrix(static_matrix const&, vector&); -template sparse_matrix::sparse_matrix(unsigned int); -template void sparse_matrix::add_new_element(unsigned int, unsigned int, const mpq&); -template void sparse_matrix::divide_row_by_constant(unsigned int, const mpq&, lp_settings&); -template bool sparse_matrix::fill_eta_matrix(unsigned int, eta_matrix**); -template mpq const & sparse_matrix::get(unsigned int, unsigned int) const; -template unsigned sparse_matrix::get_number_of_nonzeroes() const; -template bool sparse_matrix::get_pivot_for_column(unsigned int&, unsigned int&, int, unsigned int); -template unsigned sparse_matrix::lowest_row_in_column(unsigned int); -template bool sparse_matrix::pivot_with_eta(unsigned int, eta_matrix*, lp_settings&); -template void sparse_matrix::prepare_for_factorization(); -template void sparse_matrix::remove_element(vector> &, indexed_value&); -template void sparse_matrix::replace_column(unsigned int, indexed_vector&, lp_settings&); -template void sparse_matrix::set_max_in_row(vector>&); -template bool sparse_matrix::set_row_from_work_vector_and_clean_work_vector_not_adjusted(unsigned int, indexed_vector&, lp_settings&); -template bool sparse_matrix::shorten_active_matrix(unsigned int, eta_matrix*); -template void sparse_matrix::solve_y_U(vector&) const; -template sparse_matrix::sparse_matrix(static_matrix const&, vector&); -template void sparse_matrix>::add_new_element(unsigned int, unsigned int, const mpq&); -template void sparse_matrix>::divide_row_by_constant(unsigned int, const mpq&, lp_settings&); -template bool sparse_matrix>::fill_eta_matrix(unsigned int, eta_matrix >**); -template const mpq & sparse_matrix>::get(unsigned int, unsigned int) const; -template unsigned sparse_matrix>::get_number_of_nonzeroes() const; -template bool sparse_matrix>::get_pivot_for_column(unsigned int&, unsigned int&, int, unsigned int); -template unsigned sparse_matrix>::lowest_row_in_column(unsigned int); -template bool sparse_matrix>::pivot_with_eta(unsigned int, eta_matrix >*, lp_settings&); -template void sparse_matrix>::prepare_for_factorization(); -template void sparse_matrix>::remove_element(vector>&, indexed_value&); -template void sparse_matrix>::replace_column(unsigned int, indexed_vector&, lp_settings&); -template void sparse_matrix>::set_max_in_row(vector>&); -template bool sparse_matrix>::set_row_from_work_vector_and_clean_work_vector_not_adjusted(unsigned int, indexed_vector&, lp_settings&); -template bool sparse_matrix>::shorten_active_matrix(unsigned int, eta_matrix >*); -template void sparse_matrix>::solve_y_U(vector&) const; -template sparse_matrix>::sparse_matrix(static_matrix > const&, vector&); -template void sparse_matrix::double_solve_U_y(indexed_vector&, const lp_settings &); -template void sparse_matrix::double_solve_U_y(indexed_vector&, const lp_settings&); -template void sparse_matrix>::double_solve_U_y(indexed_vector&, const lp_settings&); -template void sparse_matrix >::double_solve_U_y >(indexed_vector>&, const lp_settings&); -template void lp::sparse_matrix::solve_U_y_indexed_only(lp::indexed_vector&, const lp_settings&, vector &); -template void lp::sparse_matrix::solve_U_y_indexed_only(lp::indexed_vector&, const lp_settings &, vector &); -#ifdef Z3DEBUG -template bool sparse_matrix::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const; -template bool sparse_matrix::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const; -template bool sparse_matrix >::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const; -#endif -} -template void lp::sparse_matrix >::solve_U_y_indexed_only(lp::indexed_vector&, const lp_settings &, vector &); -template void lp::sparse_matrix::solve_U_y(vector&); -template void lp::sparse_matrix::double_solve_U_y(vector&); -template void lp::sparse_matrix::solve_U_y(vector&); -template void lp::sparse_matrix::double_solve_U_y(vector&); -template void lp::sparse_matrix >::solve_U_y >(vector >&); -template void lp::sparse_matrix >::double_solve_U_y >(vector >&); -template void lp::sparse_matrix::find_error_in_solution_U_y_indexed(lp::indexed_vector&, lp::indexed_vector&, const vector &); -template double lp::sparse_matrix::dot_product_with_row(unsigned int, lp::indexed_vector const&) const; -template void lp::sparse_matrix::find_error_in_solution_U_y_indexed(lp::indexed_vector&, lp::indexed_vector&, const vector &); -template lp::mpq lp::sparse_matrix::dot_product_with_row(unsigned int, lp::indexed_vector const&) const; -template void lp::sparse_matrix >::find_error_in_solution_U_y_indexed(lp::indexed_vector&, lp::indexed_vector&, const vector &); -template lp::mpq lp::sparse_matrix >::dot_product_with_row(unsigned int, lp::indexed_vector const&) const; -template void lp::sparse_matrix >::find_error_in_solution_U_y_indexed >(lp::indexed_vector >&, lp::indexed_vector >&, const vector &); -template lp::numeric_pair lp::sparse_matrix >::dot_product_with_row >(unsigned int, lp::indexed_vector > const&) const; -template void lp::sparse_matrix::extend_and_sort_active_rows(vector const&, vector&); - -template void lp::sparse_matrix >::extend_and_sort_active_rows(vector const&, vector&); - -template void lp::sparse_matrix >::solve_U_y(vector&); -template void lp::sparse_matrix >::double_solve_U_y(vector&); -template void lp::sparse_matrix< lp::mpq,lp::numeric_pair< lp::mpq> >::set(unsigned int,unsigned int, lp::mpq); -template void lp::sparse_matrix::solve_y_U_indexed(lp::indexed_vector&, const lp_settings & ); -template void lp::sparse_matrix::solve_y_U_indexed(lp::indexed_vector&, const lp_settings &); -template void lp::sparse_matrix >::solve_y_U_indexed(lp::indexed_vector&, const lp_settings &); - diff --git a/src/util/lp/square_dense_submatrix.cpp b/src/util/lp/square_dense_submatrix.cpp index a4592c72d..802f3d74a 100644 --- a/src/util/lp/square_dense_submatrix.cpp +++ b/src/util/lp/square_dense_submatrix.cpp @@ -20,14 +20,14 @@ Revision History: #include #include "util/vector.h" #include "util/lp/square_dense_submatrix_def.h" -template void lp::square_dense_submatrix::init(lp::sparse_matrix*, unsigned int); -template lp::square_dense_submatrix::square_dense_submatrix(lp::sparse_matrix*, unsigned int); +template void lp::square_dense_submatrix::init(lp::square_sparse_matrix*, unsigned int); +template lp::square_dense_submatrix::square_dense_submatrix(lp::square_sparse_matrix*, unsigned int); template void lp::square_dense_submatrix::update_parent_matrix(lp::lp_settings&); template bool lp::square_dense_submatrix::is_L_matrix() const; template void lp::square_dense_submatrix::conjugate_by_permutation(lp::permutation_matrix&); template int lp::square_dense_submatrix::find_pivot_column_in_row(unsigned int) const; template void lp::square_dense_submatrix::pivot(unsigned int, lp::lp_settings&); -template lp::square_dense_submatrix >::square_dense_submatrix(lp::sparse_matrix >*, unsigned int); +template lp::square_dense_submatrix >::square_dense_submatrix(lp::square_sparse_matrix >*, unsigned int); template void lp::square_dense_submatrix >::update_parent_matrix(lp::lp_settings&); template bool lp::square_dense_submatrix >::is_L_matrix() const; template void lp::square_dense_submatrix >::conjugate_by_permutation(lp::permutation_matrix >&); @@ -40,7 +40,7 @@ template void lp::square_dense_submatrix::apply_from_right(vecto template void lp::square_dense_submatrix::apply_from_left_local(lp::indexed_vector&, lp::lp_settings&); template void lp::square_dense_submatrix::apply_from_left_to_vector(vector&); -template lp::square_dense_submatrix::square_dense_submatrix(lp::sparse_matrix*, unsigned int); +template lp::square_dense_submatrix::square_dense_submatrix(lp::square_sparse_matrix*, unsigned int); template void lp::square_dense_submatrix::update_parent_matrix(lp::lp_settings&); template bool lp::square_dense_submatrix::is_L_matrix() const; template void lp::square_dense_submatrix::conjugate_by_permutation(lp::permutation_matrix&); diff --git a/src/util/lp/square_dense_submatrix.h b/src/util/lp/square_dense_submatrix.h index 266610516..0c6bbdec3 100644 --- a/src/util/lp/square_dense_submatrix.h +++ b/src/util/lp/square_dense_submatrix.h @@ -34,7 +34,7 @@ Revision History: #include "util/lp/lp_settings.h" #include "util/lp/eta_matrix.h" #include "util/lp/binary_heap_upair_queue.h" -#include "util/lp/sparse_matrix.h" +#include "util/lp/square_sparse_matrix.h" namespace lp { template class square_dense_submatrix : public tail_matrix { @@ -57,7 +57,7 @@ public: unsigned m_index_start; unsigned m_dim; vector m_v; - sparse_matrix * m_parent; + square_sparse_matrix * m_parent; permutation_matrix m_row_permutation; indexed_vector m_work_vector; public: @@ -66,9 +66,9 @@ public: square_dense_submatrix() {} - square_dense_submatrix (sparse_matrix *parent_matrix, unsigned index_start); + square_dense_submatrix (square_sparse_matrix *parent_matrix, unsigned index_start); - void init(sparse_matrix *parent_matrix, unsigned index_start); + void init(square_sparse_matrix *parent_matrix, unsigned index_start); bool is_dense() const override { return true; } diff --git a/src/util/lp/square_dense_submatrix_def.h b/src/util/lp/square_dense_submatrix_def.h index ca00285bb..54f6e5327 100644 --- a/src/util/lp/square_dense_submatrix_def.h +++ b/src/util/lp/square_dense_submatrix_def.h @@ -21,7 +21,7 @@ Revision History: #include "util/lp/square_dense_submatrix.h" namespace lp { template -square_dense_submatrix::square_dense_submatrix (sparse_matrix *parent_matrix, unsigned index_start) : +square_dense_submatrix::square_dense_submatrix (square_sparse_matrix *parent_matrix, unsigned index_start) : m_index_start(index_start), m_dim(parent_matrix->dimension() - index_start), m_v(m_dim * m_dim), @@ -40,7 +40,7 @@ square_dense_submatrix::square_dense_submatrix (sparse_matrix *paren } } -template void square_dense_submatrix::init(sparse_matrix *parent_matrix, unsigned index_start) { +template void square_dense_submatrix::init(square_sparse_matrix *parent_matrix, unsigned index_start) { m_index_start = index_start; m_dim = parent_matrix->dimension() - index_start; m_v.resize(m_dim * m_dim); diff --git a/src/util/lp/square_sparse_matrix.cpp b/src/util/lp/square_sparse_matrix.cpp new file mode 100644 index 000000000..6c65cd863 --- /dev/null +++ b/src/util/lp/square_sparse_matrix.cpp @@ -0,0 +1,119 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#include +#include "util/vector.h" +#include "util/lp/lp_settings.h" +#include "util/lp/lu.h" +#include "util/lp/square_sparse_matrix_def.h" +#include "util/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/util/lp/sparse_matrix.h b/src/util/lp/square_sparse_matrix.h similarity index 93% rename from src/util/lp/sparse_matrix.h rename to src/util/lp/square_sparse_matrix.h index 3ca6d3fa8..7c3206400 100644 --- a/src/util/lp/sparse_matrix.h +++ b/src/util/lp/square_sparse_matrix.h @@ -39,7 +39,7 @@ Revision History: namespace lp { // it is a square matrix template -class sparse_matrix +class square_sparse_matrix #ifdef Z3DEBUG : public matrix #endif @@ -96,29 +96,46 @@ public: return m_column_permutation[col]; } - void copy_column_from_static_matrix(unsigned col, static_matrix const &A, unsigned col_index_in_the_new_matrix); - void copy_B(static_matrix const &A, vector & basis); + 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: - // constructor that copies columns of the basis from A - sparse_matrix(static_matrix const &A, vector & basis); + + // 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 { - sparse_matrix & m_matrix; + square_sparse_matrix & m_matrix; unsigned m_row; unsigned m_col; public: - ref_matrix_element(sparse_matrix & m, unsigned row, unsigned col):m_matrix(m), m_row(row), m_col(col) {} + 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 { - sparse_matrix & m_matrix; + square_sparse_matrix & m_matrix; unsigned m_row; public: - ref_row(sparse_matrix & m, unsigned row) : m_matrix(m), m_row(row) {} + 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); } }; @@ -154,11 +171,6 @@ public: return m_columns[col].m_values; } - // constructor creating a zero matrix of dim*dim - sparse_matrix(unsigned dim); - - - unsigned dimension() const {return static_cast(m_row_permutation.size());} #ifdef Z3DEBUG diff --git a/src/util/lp/sparse_matrix_def.h b/src/util/lp/square_sparse_matrix_def.h similarity index 78% rename from src/util/lp/sparse_matrix_def.h rename to src/util/lp/square_sparse_matrix_def.h index 39159243c..4e9600ffc 100644 --- a/src/util/lp/sparse_matrix_def.h +++ b/src/util/lp/square_sparse_matrix_def.h @@ -19,38 +19,61 @@ Revision History: --*/ #include "util/vector.h" -#include "util/lp/sparse_matrix.h" +#include "util/lp/square_sparse_matrix.h" #include #include namespace lp { template -void sparse_matrix::copy_column_from_static_matrix(unsigned col, static_matrix const &A, unsigned col_index_in_the_new_matrix) { - vector const & A_col_vector = A.m_columns[col]; - unsigned size = static_cast(A_col_vector.size()); - vector> & new_column_vector = m_columns[col_index_in_the_new_matrix].m_values; - for (unsigned l = 0; l < size; l++) { - column_cell const & col_cell = A_col_vector[l]; +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 (const auto & c : A.column(input_column)) { unsigned col_offset = static_cast(new_column_vector.size()); - vector> & row_vector = m_rows[col_cell.m_i]; + vector> & row_vector = m_rows[c.var()]; unsigned row_offset = static_cast(row_vector.size()); - const T & val = A.get_val(col_cell); - new_column_vector.push_back(indexed_value(val, col_cell.m_i, row_offset)); - row_vector.push_back(indexed_value(val, col_index_in_the_new_matrix, col_offset)); + 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 -void sparse_matrix::copy_B(static_matrix const &A, vector & basis) { - unsigned m = A.row_count(); // this should be the size of basis +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 (const 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_static_matrix(basis[j], A, j); + 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 -sparse_matrix::sparse_matrix(static_matrix const &A, vector & basis) : +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()), @@ -59,11 +82,26 @@ sparse_matrix::sparse_matrix(static_matrix const &A, vector -void sparse_matrix::set_with_no_adjusting_for_row(unsigned row, unsigned col, T val) { // should not be used in efficient code +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) { @@ -76,7 +114,7 @@ void sparse_matrix::set_with_no_adjusting_for_row(unsigned row, unsigned c } template -void sparse_matrix::set_with_no_adjusting_for_col(unsigned row, unsigned col, T val) { // should not be used in efficient code +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) { @@ -90,13 +128,13 @@ void sparse_matrix::set_with_no_adjusting_for_col(unsigned row, unsigned c template -void sparse_matrix::set_with_no_adjusting(unsigned row, unsigned col, T val) { // should not be used in efficient code +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 sparse_matrix::set(unsigned row, unsigned col, T val) { // should not be used in efficient code +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); @@ -106,7 +144,7 @@ void sparse_matrix::set(unsigned row, unsigned col, T val) { // should not } template -T const & sparse_matrix::get_not_adjusted(unsigned row, unsigned col) const { +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; @@ -116,7 +154,7 @@ T const & sparse_matrix::get_not_adjusted(unsigned row, unsigned col) cons } template -T const & sparse_matrix::get(unsigned row, unsigned col) const { // should not be used in efficient code +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); @@ -130,7 +168,7 @@ T const & sparse_matrix::get(unsigned row, unsigned col) const { // should // constructor creating a zero matrix of dim*dim template -sparse_matrix::sparse_matrix(unsigned dim) : +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), @@ -141,21 +179,21 @@ sparse_matrix::sparse_matrix(unsigned dim) : } template -void sparse_matrix::init_row_headers() { +void square_sparse_matrix::init_row_headers() { for (unsigned l = 0; l < m_row_permutation.size(); l++) { m_rows.push_back(vector>()); } } template -void sparse_matrix::init_column_headers() { // we alway have only square sparse_matrix +void square_sparse_matrix::init_column_headers() { // we alway have only square square_sparse_matrix for (unsigned l = 0; l < m_row_permutation.size(); l++) { m_columns.push_back(col_header()); } } template -unsigned sparse_matrix::lowest_row_in_column(unsigned j) { +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) { @@ -168,7 +206,7 @@ unsigned sparse_matrix::lowest_row_in_column(unsigned j) { } template -void sparse_matrix::remove_element(vector> & row_vals, unsigned row_offset, vector> & column_vals, unsigned column_offset) { +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; @@ -187,14 +225,14 @@ void sparse_matrix::remove_element(vector> & row_vals, un } template -void sparse_matrix::remove_element(vector> & row_chunk, indexed_value & row_el_iv) { +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 sparse_matrix::put_max_index_to_0(vector> & row_vals, unsigned max_index) { +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]; @@ -209,7 +247,7 @@ void sparse_matrix::put_max_index_to_0(vector> & row_vals } template -void sparse_matrix::set_max_in_row(vector> & row_vals) { +void square_sparse_matrix::set_max_in_row(vector> & row_vals) { if (row_vals.size() == 0) return; T max_val = abs(row_vals[0].m_value); @@ -225,7 +263,7 @@ void sparse_matrix::set_max_in_row(vector> & row_vals) { } template -bool sparse_matrix::pivot_with_eta(unsigned i, eta_matrix *eta_matrix, lp_settings & settings) { +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)) { @@ -243,7 +281,7 @@ bool sparse_matrix::pivot_with_eta(unsigned i, eta_matrix *eta_matri // returns the offset of the pivot column in the row template -void sparse_matrix::scan_row_to_work_vector_and_remove_pivot_column(unsigned row, unsigned pivot_column) { +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++) { @@ -260,7 +298,7 @@ void sparse_matrix::scan_row_to_work_vector_and_remove_pivot_column(unsign #ifdef Z3DEBUG template -vector sparse_matrix::get_full_row(unsigned i) const { +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)); @@ -275,7 +313,7 @@ vector sparse_matrix::get_full_row(unsigned i) const { // 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 sparse_matrix::pivot_row_to_row(unsigned i, const T& alpha, unsigned i0, lp_settings & settings ) { +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); @@ -334,7 +372,7 @@ bool sparse_matrix::pivot_row_to_row(unsigned i, const T& alpha, unsigned // set the max val as well // returns false if the resulting row is all zeroes, and true otherwise template -bool sparse_matrix::set_row_from_work_vector_and_clean_work_vector_not_adjusted(unsigned i0, indexed_vector & work_vec, +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 @@ -358,7 +396,7 @@ bool sparse_matrix::set_row_from_work_vector_and_clean_work_vector_not_adj template -void sparse_matrix::remove_zero_elements_and_set_data_on_existing_elements(unsigned row) { +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 @@ -377,7 +415,7 @@ void sparse_matrix::remove_zero_elements_and_set_data_on_existing_elements // work_vec here has not adjusted column indices template -void sparse_matrix::remove_zero_elements_and_set_data_on_existing_elements_not_adjusted(unsigned row, indexed_vector & work_vec, lp_settings & settings) { +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 @@ -398,7 +436,7 @@ void sparse_matrix::remove_zero_elements_and_set_data_on_existing_elements // adding delta columns at the end of the matrix template -void sparse_matrix::add_columns_at_the_end(unsigned delta) { +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); @@ -407,7 +445,7 @@ void sparse_matrix::add_columns_at_the_end(unsigned delta) { } template -void sparse_matrix::delete_column(int i) { +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; @@ -417,7 +455,7 @@ void sparse_matrix::delete_column(int i) { } template -void sparse_matrix::divide_row_by_constant(unsigned i, const T & t, lp_settings & settings) { +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]) { @@ -434,7 +472,7 @@ void sparse_matrix::divide_row_by_constant(unsigned i, const T & t, lp_set // solving x * this = y, and putting the answer into y // the matrix here has to be upper triangular template -void sparse_matrix::solve_y_U(vector & y) const { // works by rows +void square_sparse_matrix::solve_y_U(vector & y) const { // works by rows #ifdef Z3DEBUG // T * rs = clone_vector(y, dimension()); #endif @@ -464,7 +502,7 @@ void sparse_matrix::solve_y_U(vector & y) const { // works by rows // solving x * this = y, and putting the answer into y // the matrix here has to be upper triangular template -void sparse_matrix::solve_y_U_indexed(indexed_vector & y, const lp_settings & settings) { +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) @@ -525,7 +563,7 @@ void sparse_matrix::solve_y_U_indexed(indexed_vector & y, const lp_sett template template -void sparse_matrix::find_error_in_solution_U_y(vector& y_orig, vector & y) { +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); @@ -534,7 +572,7 @@ void sparse_matrix::find_error_in_solution_U_y(vector& y_orig, vector template -void sparse_matrix::find_error_in_solution_U_y_indexed(indexed_vector& y_orig, indexed_vector & y, const vector& sorted_active_rows) { +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 @@ -543,7 +581,7 @@ void sparse_matrix::find_error_in_solution_U_y_indexed(indexed_vector& template template -void sparse_matrix::add_delta_to_solution(const vector& del, vector & y) { +void square_sparse_matrix::add_delta_to_solution(const vector& del, vector & y) { unsigned i = dimension(); while (i--) { y[i] += del[i]; @@ -551,7 +589,7 @@ void sparse_matrix::add_delta_to_solution(const vector& del, vector } template template -void sparse_matrix::add_delta_to_solution(const indexed_vector& del, indexed_vector & y) { +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) { @@ -560,7 +598,7 @@ void sparse_matrix::add_delta_to_solution(const indexed_vector& del, in } template template -void sparse_matrix::double_solve_U_y(indexed_vector& y, const lp_settings & settings){ +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; @@ -582,7 +620,7 @@ void sparse_matrix::double_solve_U_y(indexed_vector& y, const lp_settin } template template -void sparse_matrix::double_solve_U_y(vector& y){ +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); @@ -595,7 +633,7 @@ void sparse_matrix::double_solve_U_y(vector& y){ // the matrix here has to be upper triangular template template -void sparse_matrix::solve_U_y(vector & y) { // it is a column wise version +void square_sparse_matrix::solve_U_y(vector & y) { // it is a column wise version #ifdef Z3DEBUG // T * rs = clone_vector(y, dimension()); #endif @@ -618,7 +656,7 @@ void sparse_matrix::solve_U_y(vector & y) { // it is a column wise vers #endif } template -void sparse_matrix::process_index_recursively_for_y_U(unsigned j, vector & sorted_active_rows) { +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)]; @@ -633,7 +671,7 @@ void sparse_matrix::process_index_recursively_for_y_U(unsigned j, vector -void sparse_matrix::process_column_recursively(unsigned j, vector & sorted_active_rows) { +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) { @@ -649,7 +687,7 @@ void sparse_matrix::process_column_recursively(unsigned j, vector -void sparse_matrix::create_graph_G(const vector & index_or_right_side, vector & sorted_active_rows) { +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); @@ -662,7 +700,7 @@ void sparse_matrix::create_graph_G(const vector & index_or_right template -void sparse_matrix::extend_and_sort_active_rows(const vector & index_or_right_side, vector & sorted_active_rows) { +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); @@ -676,7 +714,7 @@ void sparse_matrix::extend_and_sort_active_rows(const vector & i template template -void sparse_matrix::solve_U_y_indexed_only(indexed_vector & y, const lp_settings & settings, vector & sorted_active_rows) { // it is a column wise version +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;) { @@ -710,7 +748,7 @@ void sparse_matrix::solve_U_y_indexed_only(indexed_vector & y, const lp template template -L sparse_matrix::dot_product_with_row (unsigned row, const vector & y) const { +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) { @@ -722,7 +760,7 @@ L sparse_matrix::dot_product_with_row (unsigned row, const vector & y) template template -L sparse_matrix::dot_product_with_row (unsigned row, const indexed_vector & y) const { +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) { @@ -734,7 +772,7 @@ L sparse_matrix::dot_product_with_row (unsigned row, const indexed_vector< template -unsigned sparse_matrix::get_number_of_nonzeroes() const { +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); @@ -743,7 +781,7 @@ unsigned sparse_matrix::get_number_of_nonzeroes() const { } template -bool sparse_matrix::get_non_zero_column_in_row(unsigned i, unsigned *j) const { +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) { @@ -754,7 +792,7 @@ bool sparse_matrix::get_non_zero_column_in_row(unsigned i, unsigned *j) co } template -void sparse_matrix::remove_element_that_is_not_in_w(vector> & column_vals, indexed_value & col_el_iv) { +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; @@ -767,7 +805,7 @@ void sparse_matrix::remove_element_that_is_not_in_w(vector -void sparse_matrix::remove_elements_that_are_not_in_w_and_update_common_elements(unsigned column_to_replace, indexed_vector & w) { +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; @@ -796,7 +834,7 @@ void sparse_matrix::remove_elements_that_are_not_in_w_and_update_common_el } template -void sparse_matrix::add_new_element(unsigned row, unsigned col, const T& val) { +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()); @@ -809,7 +847,7 @@ void sparse_matrix::add_new_element(unsigned row, unsigned col, const T& v // 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 sparse_matrix::add_new_elements_of_w_and_clear_w(unsigned column_to_replace, indexed_vector & w, lp_settings & settings) { +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 @@ -827,14 +865,14 @@ void sparse_matrix::add_new_elements_of_w_and_clear_w(unsigned column_to_r } template -void sparse_matrix::replace_column(unsigned column_to_replace, indexed_vector & w, lp_settings &settings) { +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 sparse_matrix::pivot_score(unsigned i, unsigned j) { +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, @@ -847,7 +885,7 @@ unsigned sparse_matrix::pivot_score(unsigned i, unsigned j) { } template -void sparse_matrix::enqueue_domain_into_pivot_queue() { +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]; @@ -860,7 +898,7 @@ void sparse_matrix::enqueue_domain_into_pivot_queue() { } template -void sparse_matrix::set_max_in_rows() { +void square_sparse_matrix::set_max_in_rows() { unsigned i = dimension(); while (i--) set_max_in_row(i); @@ -868,27 +906,27 @@ void sparse_matrix::set_max_in_rows() { template -void sparse_matrix::zero_shortened_markovitz_numbers() { +void square_sparse_matrix::zero_shortened_markovitz_numbers() { for (auto & ch : m_columns) ch.zero_shortened_markovitz(); } template -void sparse_matrix::prepare_for_factorization() { +void square_sparse_matrix::prepare_for_factorization() { zero_shortened_markovitz_numbers(); set_max_in_rows(); enqueue_domain_into_pivot_queue(); } template -void sparse_matrix::recover_pivot_queue(vector & rejected_pivots) { +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 sparse_matrix::elem_is_too_small(unsigned i, unsigned j, int c_partial_pivoting) { +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) { @@ -904,7 +942,7 @@ int sparse_matrix::elem_is_too_small(unsigned i, unsigned j, int c_partial } template -bool sparse_matrix::remove_row_from_active_pivots_and_shorten_columns(unsigned row) { +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); @@ -921,7 +959,7 @@ bool sparse_matrix::remove_row_from_active_pivots_and_shorten_columns(unsi } template -void sparse_matrix::remove_pivot_column(unsigned row) { +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) @@ -929,7 +967,7 @@ void sparse_matrix::remove_pivot_column(unsigned row) { } template -void sparse_matrix::update_active_pivots(unsigned row) { +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]; @@ -944,7 +982,7 @@ void sparse_matrix::update_active_pivots(unsigned row) { } template -bool sparse_matrix::shorten_active_matrix(unsigned row, eta_matrix *eta_matrix) { +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); @@ -969,7 +1007,7 @@ bool sparse_matrix::shorten_active_matrix(unsigned row, eta_matrix * } template -unsigned sparse_matrix::pivot_score_without_shortened_counters(unsigned i, unsigned j, unsigned k) { +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) { @@ -981,7 +1019,7 @@ unsigned sparse_matrix::pivot_score_without_shortened_counters(unsigned i, } #ifdef Z3DEBUG template -bool sparse_matrix::can_improve_score_for_row(unsigned row, unsigned score, T const & c_partial_pivoting, unsigned k) { +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]; @@ -1005,7 +1043,7 @@ bool sparse_matrix::can_improve_score_for_row(unsigned row, unsigned score } template -bool sparse_matrix::really_best_pivot(unsigned i, unsigned j, T const & c_partial_pivoting, unsigned k) { +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)); @@ -1013,7 +1051,7 @@ bool sparse_matrix::really_best_pivot(unsigned i, unsigned j, T const & c_ return true; } template -void sparse_matrix::print_active_matrix(unsigned k, std::ostream & out) { +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; @@ -1038,7 +1076,7 @@ void sparse_matrix::print_active_matrix(unsigned k, std::ostream & out) { } template -bool sparse_matrix::pivot_queue_is_correct_for_row(unsigned i, unsigned k) { +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) == @@ -1048,7 +1086,7 @@ bool sparse_matrix::pivot_queue_is_correct_for_row(unsigned i, unsigned k) } template -bool sparse_matrix::pivot_queue_is_correct_after_pivoting(int k) { +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()); @@ -1057,7 +1095,7 @@ bool sparse_matrix::pivot_queue_is_correct_after_pivoting(int k) { #endif template -bool sparse_matrix::get_pivot_for_column(unsigned &i, unsigned &j, int c_partial_pivoting, unsigned k) { +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); @@ -1087,7 +1125,7 @@ bool sparse_matrix::get_pivot_for_column(unsigned &i, unsigned &j, int c_p } template -bool sparse_matrix::elem_is_too_small(vector> & row_chunk, indexed_value & iv, int c_partial_pivoting) { +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 } @@ -1097,7 +1135,7 @@ bool sparse_matrix::elem_is_too_small(vector> & row_chunk } template -bool sparse_matrix::shorten_columns_by_pivot_row(unsigned i, unsigned pivot_column) { +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) { @@ -1116,7 +1154,7 @@ bool sparse_matrix::shorten_columns_by_pivot_row(unsigned i, unsigned pivo } template -bool sparse_matrix::fill_eta_matrix(unsigned j, eta_matrix ** eta) { +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) { @@ -1163,7 +1201,7 @@ bool sparse_matrix::fill_eta_matrix(unsigned j, eta_matrix ** eta) { } #ifdef Z3DEBUG template -bool sparse_matrix::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings & settings) const { +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()); @@ -1188,7 +1226,7 @@ bool sparse_matrix::is_upper_triangular_and_maximums_are_set_correctly_in_ } template -bool sparse_matrix::is_upper_triangular_until(unsigned k) const { +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); @@ -1202,7 +1240,7 @@ bool sparse_matrix::is_upper_triangular_until(unsigned k) const { } template -void sparse_matrix::check_column_vs_rows(unsigned col) { +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); @@ -1225,7 +1263,7 @@ void sparse_matrix::check_column_vs_rows(unsigned col) { } template -void sparse_matrix::check_row_vs_columns(unsigned row) { +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); @@ -1249,20 +1287,20 @@ void sparse_matrix::check_row_vs_columns(unsigned row) { } template -void sparse_matrix::check_rows_vs_columns() { +void square_sparse_matrix::check_rows_vs_columns() { for (unsigned i = 0; i < dimension(); i++) { check_row_vs_columns(i); } } template -void sparse_matrix::check_columns_vs_rows() { +void square_sparse_matrix::check_columns_vs_rows() { for (unsigned i = 0; i < dimension(); i++) { check_column_vs_rows(i); } } template -void sparse_matrix::check_matrix() { +void square_sparse_matrix::check_matrix() { check_rows_vs_columns(); check_columns_vs_rows(); } diff --git a/src/util/lp/static_matrix.h b/src/util/lp/static_matrix.h index 2f127b84a..e4bf257f4 100644 --- a/src/util/lp/static_matrix.h +++ b/src/util/lp/static_matrix.h @@ -82,17 +82,17 @@ public: ref(static_matrix & m, unsigned row, unsigned col):m_matrix(m), m_row(row), m_col(col) {} ref & operator=(T const & v) { m_matrix.set( m_row, m_col, v); return *this; } - ref & operator=(ref const & v) { m_matrix.set(m_row, m_col, v.m_matrix.get(v.m_row, v.m_col)); return *this; } + ref operator=(ref & 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_elem(m_row, m_col); } }; class ref_row { - static_matrix & m_matrix; + const static_matrix & m_matrix; unsigned m_row; public: - ref_row(static_matrix & m, unsigned row):m_matrix(m), m_row(row) {} - ref operator[](unsigned col) const { return ref(m_matrix, m_row, col); } + ref_row(const static_matrix & m, unsigned row): m_matrix(m), m_row(row) {} + const T operator[](unsigned col) const { return m_matrix.get_elem(m_row, col); } }; public: @@ -450,6 +450,10 @@ public: column_container column(unsigned j) const { return column_container(j, *this); } - + + ref_row operator[](unsigned i) const { return ref_row(*this, i);} + typedef T coefftype; + typedef X argtype; }; + } diff --git a/src/util/lp/tail_matrix.h b/src/util/lp/tail_matrix.h index 37b217205..2a851bc16 100644 --- a/src/util/lp/tail_matrix.h +++ b/src/util/lp/tail_matrix.h @@ -39,5 +39,12 @@ public: virtual void apply_from_right(indexed_vector & w) = 0; virtual ~tail_matrix() {} virtual bool is_dense() const = 0; + struct ref_row { + const tail_matrix & m_A; + unsigned m_row; + ref_row(const tail_matrix& m, unsigned row): m_A(m), m_row(row) {} + T operator[](unsigned j) const { return m_A.get_elem(m_row, j);} + }; + ref_row operator[](unsigned i) const { return ref_row(*this, i);} }; } diff --git a/src/util/mpq.h b/src/util/mpq.h index 010bb2c8a..01ad06db4 100644 --- a/src/util/mpq.h +++ b/src/util/mpq.h @@ -501,7 +501,7 @@ public: void machine_div(mpz const & a, mpz const & b, mpz & c) { mpz_manager::machine_div(a, b, c); } - void machine_div_rem(mpz const & a, mpz const & b, mpz & c, mpz& d) { mpz_manager::machine_div_rem(a, b, c, d); } + void machine_div_rem(mpz const & a, mpz const & b, mpz & c, mpz & d) { mpz_manager::machine_div_rem(a, b, c, d); } void div(mpz const & a, mpz const & b, mpz & c) { mpz_manager::div(a, b, c); } diff --git a/src/util/rational.h b/src/util/rational.h index 8047696be..3403848c4 100644 --- a/src/util/rational.h +++ b/src/util/rational.h @@ -188,6 +188,12 @@ public: return r; } + friend inline rational machine_div_rem(rational const & r1, rational const & r2, rational & rem) { + rational r; + rational::m().machine_idiv_rem(r1.m_val, r2.m_val, r.m_val, rem.m_val); + return r; + } + friend inline rational mod(rational const & r1, rational const & r2) { rational r; rational::m().mod(r1.m_val, r2.m_val, r.m_val); @@ -409,8 +415,6 @@ public: } return num_bits; } - - }; inline bool operator!=(rational const & r1, rational const & r2) {