/*++ Copyright (c) 2017 Microsoft Corporation Module Name: Abstract: for matrix B we have t0*...*tn-1*B = Q*U*R here ti are matrices corresponding to pivot operations, including columns and rows swaps, or a multiplication matrix row by a number Q, R - permutations and U is an upper triangular matrix Author: Lev Nachmanson (levnach) Revision History: --*/ #pragma once #include "util/vector.h" #include "util/debug.h" #include #include #include "math/lp/square_sparse_matrix.h" #include "math/lp/static_matrix.h" #include #include "math/lp/numeric_pair.h" #include #include #include "math/lp/row_eta_matrix.h" #include "math/lp/square_dense_submatrix.h" #include "math/lp/dense_matrix.h" namespace lp { template // print the nr x nc submatrix at the top left corner void print_submatrix(square_sparse_matrix & m, unsigned mr, unsigned nc); template void print_matrix(M &m, std::ostream & out); template X dot_product(const vector & a, const vector & b) { lp_assert(a.size() == b.size()); auto r = zero_of_type(); for (unsigned i = 0; i < a.size(); i++) { r += a[i] * b[i]; } return r; } template class one_elem_on_diag: public tail_matrix { unsigned m_i; T m_val; public: one_elem_on_diag(unsigned i, T val) : m_i(i), m_val(val) { #ifdef Z3DEBUG m_one_over_val = numeric_traits::one() / m_val; #endif } bool is_dense() const override { return false; } one_elem_on_diag(const one_elem_on_diag & o); #ifdef Z3DEBUG unsigned m_m; unsigned m_n; void set_number_of_rows(unsigned m) override { m_m = m; m_n = m; } void set_number_of_columns(unsigned n) override { m_m = n; m_n = n; } T m_one_over_val; T get_elem (unsigned i, unsigned j) const override; unsigned row_count() const override { return m_m; } // not defined } unsigned column_count() const override { return m_m; } // not defined } #endif void apply_from_left(vector & w, lp_settings &) override { w[m_i] /= m_val; } void apply_from_right(vector & w) override { w[m_i] /= m_val; } void apply_from_right(indexed_vector & w) override { if (is_zero(w.m_data[m_i])) return; auto & v = w.m_data[m_i] /= m_val; if (lp_settings::is_eps_small_general(v, 1e-14)) { w.erase_from_index(m_i); v = zero_of_type(); } } void apply_from_left_to_T(indexed_vector & w, lp_settings & settings) override; void conjugate_by_permutation(permutation_matrix & p) { // this = p * this * p(-1) #ifdef Z3DEBUG // auto rev = p.get_reverse(); // auto deb = ((*this) * rev); // deb = p * deb; #endif m_i = p.apply_reverse(m_i); #ifdef Z3DEBUG // lp_assert(*this == deb); #endif } }; // end of one_elem_on_diag enum class LU_status { OK, Degenerated}; // This class supports updates of the columns of B, and solves systems Bx=b,and yB=c // Using Suhl-Suhl method described in the dissertation of Achim Koberstein, Chapter 5 template class lu { LU_status m_status; public: typedef typename M::coefftype T; typedef typename M::argtype X; // the fields unsigned m_dim; const M & m_A; permutation_matrix m_Q; permutation_matrix m_R; permutation_matrix m_r_wave; square_sparse_matrix m_U; square_dense_submatrix* m_dense_LU; vector *> m_tail; lp_settings & m_settings; bool m_failure; indexed_vector m_row_eta_work_vector; indexed_vector m_w_for_extension; indexed_vector m_y_copy; indexed_vector m_ii; //to optimize the work with the m_index fields unsigned m_refactor_counter; // constructor // if A is an m by n matrix then basis has length m and values in [0,n); the values are all different // they represent the set of m columns lu(const M & A, vector& basis, lp_settings & settings); lu(const M & A, lp_settings&); void debug_test_of_basis(const M & A, vector & basis); void solve_Bd_when_w_is_ready(vector & d, indexed_vector& w ); void solve_By(indexed_vector & y); void solve_By(vector & y); void solve_By_for_T_indexed_only(indexed_vector& y, const lp_settings &); template void solve_By_when_y_is_ready(indexed_vector & y); void solve_By_when_y_is_ready_for_X(vector & y); void solve_By_when_y_is_ready_for_T(vector & y, vector & index); void print_indexed_vector(indexed_vector & w, std::ofstream & f); void print_matrix_compact(std::ostream & f); void print(indexed_vector & w, const vector& basis); void solve_Bd_faster(unsigned a_column, indexed_vector & d); // d is the right side on the input and the solution at the exit void solve_yB_indexed(indexed_vector& y); void add_delta_to_solution_indexed(indexed_vector& y); void add_delta_to_solution(const vector& yc, vector& y); void find_error_of_yB(vector& yc, const vector& y, const vector& basis); void find_error_of_yB_indexed(const indexed_vector& y, const vector& heading, const lp_settings& settings); void solve_yB_with_error_check(vector & y, const vector& basis); void solve_yB_with_error_check_indexed(indexed_vector & y, const vector& heading, const vector & basis, const lp_settings &); void apply_Q_R_to_U(permutation_matrix & r_wave); LU_status get_status() { return m_status; } void set_status(LU_status status) { m_status = status; } ~lu(); void init_vector_y(vector & y); void perform_transformations_on_w(indexed_vector& w); void init_vector_w(unsigned entering, indexed_vector & w); void apply_lp_list_to_w(indexed_vector & w); void apply_lp_list_to_y(vector& y); void swap_rows(int j, int k); void swap_columns(int j, int pivot_column); void push_matrix_to_tail(tail_matrix* tm) { m_tail.push_back(tm); } bool pivot_the_row(int row); eta_matrix * get_eta_matrix_for_pivot(unsigned j); // we're processing the column j now eta_matrix * get_eta_matrix_for_pivot(unsigned j, square_sparse_matrix& copy_of_U); // see page 407 of Chvatal unsigned transform_U_to_V_by_replacing_column(indexed_vector & w, unsigned leaving_column_of_U); #ifdef Z3DEBUG void check_vector_w(unsigned entering); void check_apply_matrix_to_vector(matrix *lp, T *w); void check_apply_lp_lists_to_w(T * w); // provide some access operators for testing permutation_matrix & Q() { return m_Q; } permutation_matrix & R() { return m_R; } matrix & U() { return m_U; } unsigned tail_size() { return m_tail.size(); } tail_matrix * get_lp_matrix(unsigned i) { return m_tail[i]; } T B_(unsigned i, unsigned j, const vector& basis) { return m_A[i][basis[j]]; } unsigned dimension() { return m_dim; } #endif unsigned get_number_of_nonzeroes() { return m_U.get_number_of_nonzeroes(); } void process_column(int j); bool is_correct(const vector& basis); bool is_correct(); // needed for debugging purposes void copy_w(T *buffer, indexed_vector & w); // needed for debugging purposes void restore_w(T *buffer, indexed_vector & w); bool all_columns_and_rows_are_active(); bool too_dense(unsigned j) const; void pivot_in_dense_mode(unsigned i); void create_initial_factorization(); void calculate_r_wave_and_update_U(unsigned bump_start, unsigned bump_end, permutation_matrix & r_wave); void scan_last_row_to_work_vector(unsigned lowest_row_of_the_bump); bool diagonal_element_is_off(T /* diag_element */) { return false; } void pivot_and_solve_the_system(unsigned replaced_column, unsigned lowest_row_of_the_bump); // see Achim Koberstein's thesis page 58, but here we solve the system and pivot to the last // row at the same time void replace_column(T pivot_elem, indexed_vector & w, unsigned leaving_column_of_U); void calculate_Lwave_Pwave_for_bump(unsigned replaced_column, unsigned lowest_row_of_the_bump); void calculate_Lwave_Pwave_for_last_row(unsigned lowest_row_of_the_bump, T diagonal_element); void prepare_entering(unsigned entering, indexed_vector & w) { lp_assert(false); } bool need_to_refactor() { lp_assert(false); return m_refactor_counter >= 200; } void adjust_dimension_with_matrix_A() { lp_assert(false); } }; // end of lu template void init_factorization(lu* & factorization, M & m_A, vector & m_basis, lp_settings &m_settings); #ifdef Z3DEBUG template dense_matrix get_B(lu& f, const vector& basis); template dense_matrix get_B(lu& f); #endif }