From 9a7c99da332fae795dc52ea82e4ccba000e2b778 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Sat, 4 Mar 2023 15:26:55 -0800 Subject: [PATCH] rm lu --- src/math/lp/lu.h | 10 +- src/math/lp/lu_def.h | 235 +++---------------------------------------- src/test/lp/lp.cpp | 136 ------------------------- 3 files changed, 16 insertions(+), 365 deletions(-) diff --git a/src/math/lp/lu.h b/src/math/lp/lu.h index bcbb850434d..b3e2f0c2221 100644 --- a/src/math/lp/lu.h +++ b/src/math/lp/lu.h @@ -266,13 +266,6 @@ class lu { 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 // needed for debugging purposes void copy_w(T *buffer, indexed_vector & w); @@ -296,8 +289,7 @@ class lu { void pivot_and_solve_the_system(unsigned replaced_column, unsigned lowest_row_of_the_bump); // see Achim Koberstein's thesis page 58, but here we solve the system and pivot to the last // row at the same time - row_eta_matrix *get_row_eta_matrix_and_set_row_vector(unsigned replaced_column, unsigned lowest_row_of_the_bump, const T & pivot_elem_for_checking); - + void replace_column(T pivot_elem, indexed_vector & w, unsigned leaving_column_of_U); void calculate_Lwave_Pwave_for_bump(unsigned replaced_column, unsigned lowest_row_of_the_bump); diff --git a/src/math/lp/lu_def.h b/src/math/lp/lu_def.h index 059a430124c..d7c0c0c3a07 100644 --- a/src/math/lp/lu_def.h +++ b/src/math/lp/lu_def.h @@ -626,161 +626,41 @@ void lu::process_column(int j) { } template bool lu::is_correct(const vector& basis) { -#ifdef Z3DEBUG - if (get_status() != LU_status::OK) { - return false; - } - dense_matrix left_side = get_left_side(basis); - dense_matrix right_side = get_right_side(); - return left_side == right_side; -#else - return true; -#endif +return true; } template bool lu::is_correct() { -#ifdef Z3DEBUG - if (get_status() != LU_status::OK) { - return false; - } - dense_matrix left_side = get_left_side(); - dense_matrix right_side = get_right_side(); - return left_side == right_side; -#else return true; -#endif } -#ifdef Z3DEBUG -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++) { - 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_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); - 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_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; -} -#endif - // needed for debugging purposes template void lu::copy_w(T *buffer, indexed_vector & w) { - unsigned i = m_dim; - while (i--) { - buffer[i] = w[i]; - } + } // needed for debugging purposes 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() { - unsigned i = m_dim; - while (i--) { - lp_assert(m_U.col_is_active(i)); - lp_assert(m_U.row_is_active(i)); - } return true; } template bool lu::too_dense(unsigned j) const { - unsigned r = m_dim - j; - if (r < 5) - return false; - // if (j * 5 < m_dim * 4) // start looking for dense only at the bottom of the rows - // return false; - // 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(); + return false; } 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; - return; - } - if (i != static_cast(j)) { - swap_columns(i, j); - m_dense_LU->swap_columns(i, j); - } - m_dense_LU->pivot(i, m_settings); + } template void lu::create_initial_factorization(){ - m_U.prepare_for_factorization(); - unsigned j; - for (j = 0; j < m_dim; j++) { - process_column(j); - if (m_failure) { - set_status(LU_status::Degenerated); - return; - } - if (too_dense(j)) { - break; - } - } - if (j == m_dim) { - // TBD does not compile: lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings)); - // lp_assert(is_correct()); - // lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings)); - return; - } - j++; - m_dense_LU = new square_dense_submatrix(&m_U, j); - for (; j < m_dim; j++) { - pivot_in_dense_mode(j); - if (m_failure) { - set_status(LU_status::Degenerated); - return; - } - } - m_dense_LU->update_parent_matrix(m_settings); - lp_assert(m_dense_LU->is_L_matrix()); - m_dense_LU->conjugate_by_permutation(m_Q); - push_matrix_to_tail(m_dense_LU); - m_refactor_counter = 0; - // lp_assert(is_correct()); - // lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings)); + } template @@ -856,123 +736,38 @@ void lu::pivot_and_solve_the_system(unsigned replaced_column, unsigned lowest } } } -// 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) { - 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); - if (numeric_traits::precise() == false && !is_zero(pivot_elem_for_checking)) { - T denom = std::max(T(1), abs(pivot_elem_for_checking)); - if ( - !m_settings.abs_val_is_smaller_than_pivot_tolerance((m_row_eta_work_vector[lowest_row_of_the_bump] - pivot_elem_for_checking) / denom)) { - set_status(LU_status::Degenerated); - // LP_OUT(m_settings, "diagonal element is off" << std::endl); - return nullptr; - } - } -#ifdef Z3DEBUG - 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); -#endif - - for (auto j : m_row_eta_work_vector.m_index) { - if (j < lowest_row_of_the_bump) { - auto & v = m_row_eta_work_vector[j]; - if (!is_zero(v)) { - if (!m_settings.abs_val_is_smaller_than_drop_tolerance(v)){ - ret->push_back(j, v); - } - v = numeric_traits::zero(); - } - } - } // now the lowest_row_of_the_bump contains the rest of the row to the right of the bump with correct values - return ret; -} 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); - m_r_wave.init(m_dim); - calculate_r_wave_and_update_U(replaced_column, lowest_row_of_the_bump, m_r_wave); - auto row_eta = get_row_eta_matrix_and_set_row_vector(replaced_column, lowest_row_of_the_bump, pivot_elem_for_checking); - - if (get_status() == LU_status::Degenerated) { - m_row_eta_work_vector.clear_all(); - return; - } - m_Q.multiply_by_permutation_from_right(m_r_wave); - m_R.multiply_by_permutation_reverse_from_left(m_r_wave); - if (row_eta != nullptr) { - row_eta->conjugate_by_permutation(m_Q); - push_matrix_to_tail(row_eta); - } - calculate_Lwave_Pwave_for_bump(replaced_column, lowest_row_of_the_bump); - // 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()); + lp_assert(false); } 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]; - // lp_assert(m_row_eta_work_vector.is_OK()); - m_U.set_row_from_work_vector_and_clean_work_vector_not_adjusted(m_U.adjust_row(lowest_row_of_the_bump), m_row_eta_work_vector, m_settings); - } else { - diagonal_elem = m_U(lowest_row_of_the_bump, lowest_row_of_the_bump); // todo - get it more efficiently - } - if (m_settings.abs_val_is_smaller_than_pivot_tolerance(diagonal_elem)) { - set_status(LU_status::Degenerated); - return; - } - - calculate_Lwave_Pwave_for_last_row(lowest_row_of_the_bump, diagonal_elem); - // lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings)); + lp_assert(false);// lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings)); } template void lu::calculate_Lwave_Pwave_for_last_row(unsigned lowest_row_of_the_bump, T diagonal_element) { - auto l = new one_elem_on_diag(lowest_row_of_the_bump, diagonal_element); -#ifdef Z3DEBUG - l->set_number_of_columns(m_dim); -#endif - push_matrix_to_tail(l); - m_U.divide_row_by_constant(lowest_row_of_the_bump, diagonal_element, m_settings); - l->conjugate_by_permutation(m_Q); + lp_assert(false); } 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); - // if (factorization->get_status() != LU_status::OK) - // LP_OUT(m_settings, "failing in init_factorization" << std::endl); + lp_assert(false); } #ifdef Z3DEBUG 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()); - 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)); - + lp_assert(false); + + dense_matrix B(0, 0); 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]); - + lp_assert(false); + dense_matrix B(0,0); return B; } #endif diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index ceca92b4ecc..1a2ee2338b6 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -611,142 +611,6 @@ void fill_larger_square_sparse_matrix(static_matrix & m){ int perm_id = 0; -#ifdef Z3DEBUG -void test_larger_lu_exp(lp_settings & settings) { - std::cout << " test_larger_lu_exp" << std::endl; - static_matrix m(6, 12); - vector basis(6); - basis[0] = 1; - basis[1] = 3; - basis[2] = 0; - basis[3] = 4; - basis[4] = 5; - basis[5] = 6; - - - 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); - - dense_matrix left_side = l.get_left_side(basis); - dense_matrix right_side = l.get_right_side(); - lp_assert(left_side == right_side); - int leaving = 3; - int entering = 8; - for (unsigned i = 0; i < m.row_count(); i++) { - std::cout << static_cast(m(i, entering)) << std::endl; - } - - indexed_vector w(m.row_count()); - - l.prepare_entering(entering, w); - l.replace_column(0, w, heading[leaving]); - change_basis(entering, leaving, basis, non_basic_columns, heading); - lp_assert(l.is_correct(basis)); - - l.prepare_entering(11, w); // to init vector w - l.replace_column(0, w, heading[0]); - change_basis(11, 0, basis, non_basic_columns, heading); - lp_assert(l.is_correct(basis)); -} - -void test_larger_lu_with_holes(lp_settings & settings) { - std::cout << " test_larger_lu_with_holes" << std::endl; - static_matrix m(8, 9); - vector basis(8); - for (unsigned i = 0; i < m.row_count(); i++) { - basis[i] = i; - } - m(0, 0) = 1; m(0, 1) = 2; m(0, 2) = 3; m(0, 3) = 4; m(0, 4) = 5; m(0, 8) = 99; - /* */ m(1, 1) =- 6; m(1, 2) = 7; m(1, 3) = 8; m(1, 4) = 9; - /* */ m(2, 2) = 10; - /* */ m(3, 2) = 11; m(3, 3) = -12; - /* */ m(4, 2) = 13; m(4, 3) = 14; m(4, 4) = 15; - // the rest of the matrix is denser - m(5, 4) = 28; m(5, 5) = -18; m(5, 6) = 19; m(5, 7) = 25; - /* */ m(6, 5) = 20; m(6, 6) = -21; - /* */ m(7, 5) = 22; m(7, 6) = 23; m(7, 7) = 24; m(7, 8) = 88; - 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); - 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); - } - - dense_matrix left_side = l.get_left_side(basis); - dense_matrix right_side = l.get_right_side(); - if (!(left_side == right_side)) { - std::cout << "different sides" << std::endl; - } - - indexed_vector w(m.row_count()); - l.prepare_entering(8, w); // to init vector w - l.replace_column(0, w, heading[0]); - change_basis(8, 0, basis, non_basic_columns, heading); - lp_assert(l.is_correct(basis)); -} - - -void test_larger_lu(lp_settings& settings) { - std::cout << " test_larger_lu" << std::endl; - static_matrix m(6, 12); - vector basis(6); - basis[0] = 1; - basis[1] = 3; - basis[2] = 0; - basis[3] = 4; - basis[4] = 5; - basis[5] = 6; - - - 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); - // std::cout << "printing factorization" << std::endl; - // for (int i = lu.tail_size() - 1; i >=0; i--) { - // auto lp = lu.get_lp_matrix(i); - // lp->set_number_of_columns(m.row_count()); - // lp->set_number_of_rows(m.row_count()); - // print_matrix(* lp); - // } - - dense_matrix left_side = l.get_left_side(basis); - dense_matrix right_side = l.get_right_side(); - if (!(left_side == right_side)) { - std::cout << "left side" << std::endl; - print_matrix(&left_side, std::cout); - std::cout << "right side" << std::endl; - print_matrix(&right_side, std::cout); - - std::cout << "different sides" << std::endl; - std::cout << "initial factorization is incorrect" << std::endl; - exit(1); - } - indexed_vector w(m.row_count()); - l.prepare_entering(9, w); // to init vector w - l.replace_column(0, w, heading[0]); - change_basis(9, 0, basis, non_basic_columns, heading); - lp_assert(l.is_correct(basis)); -} - - -void test_lu(lp_settings & settings) { - -} -#endif