mirror of
https://github.com/Z3Prover/z3
synced 2025-04-13 12:28:44 +00:00
rm lu
This commit is contained in:
Lev Nachmanson
parent
6132bf93f7
commit
377ceba6d5
|
|
@ -21,7 +21,6 @@ z3_add_component(lp
|
|||
lp_core_solver_base.cpp
|
||||
lp_primal_core_solver.cpp
|
||||
lp_settings.cpp
|
||||
lu.cpp
|
||||
lp_utils.cpp
|
||||
matrix.cpp
|
||||
mon_eq.cpp
|
||||
|
|
|
|||
|
|
@ -25,11 +25,21 @@ Revision History:
|
|||
#include "math/lp/core_solver_pretty_printer.h"
|
||||
#include "math/lp/numeric_pair.h"
|
||||
#include "math/lp/static_matrix.h"
|
||||
#include "math/lp/lu.h"
|
||||
#include "math/lp/permutation_matrix.h"
|
||||
#include "math/lp/column_namer.h"
|
||||
#include "math/lp/u_set.h"
|
||||
|
||||
|
||||
namespace lp {
|
||||
template <typename T, typename X>
|
||||
X dot_product(const vector<T> & a, const vector<X> & b) {
|
||||
lp_assert(a.size() == b.size());
|
||||
auto r = zero_of_type<X>();
|
||||
for (unsigned i = 0; i < a.size(); i++) {
|
||||
r += a[i] * b[i];
|
||||
}
|
||||
return r;
|
||||
}
|
||||
|
||||
template <typename T, typename X> // X represents the type of the x variable and the bounds
|
||||
class lp_core_solver_base {
|
||||
|
|
@ -156,6 +166,8 @@ public:
|
|||
|
||||
void pretty_print(std::ostream & out);
|
||||
|
||||
|
||||
|
||||
X get_cost() const {
|
||||
return dot_product(m_costs, m_x);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -29,7 +29,6 @@ Revision History:
|
|||
#include <math.h>
|
||||
#include <cstdlib>
|
||||
#include <algorithm>
|
||||
#include "math/lp/lu.h"
|
||||
#include "math/lp/static_matrix.h"
|
||||
#include "math/lp/core_solver_pretty_printer.h"
|
||||
#include "math/lp/lp_core_solver_base.h"
|
||||
|
|
@ -418,9 +417,6 @@ public:
|
|||
// returns the number of iterations
|
||||
unsigned solve();
|
||||
|
||||
lu<static_matrix<T, X>> * factorization() {return nullptr;}
|
||||
|
||||
void delete_factorization();
|
||||
|
||||
// according to Swietanowski, " A new steepest edge approximation for the simplex method for linear programming"
|
||||
void init_column_norms();
|
||||
|
|
|
|||
|
|
@ -26,6 +26,7 @@ Revision History:
|
|||
#include <set>
|
||||
#include <string>
|
||||
#include "math/lp/lp_primal_core_solver.h"
|
||||
#include "math/lp/dense_matrix.h"
|
||||
namespace lp {
|
||||
// This core solver solves (Ax=b, lower_bound_values \leq x \leq upper_bound_values, maximize costs*x )
|
||||
// The right side b is given implicitly by x and the basis
|
||||
|
|
@ -733,8 +734,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
|
|||
default:
|
||||
break; // do nothing
|
||||
}
|
||||
} while (this->get_status() != lp_status::FLOATING_POINT_ERROR
|
||||
&&
|
||||
} while (
|
||||
this->get_status() != lp_status::UNBOUNDED
|
||||
&&
|
||||
this->get_status() != lp_status::OPTIMAL
|
||||
|
|
@ -745,18 +745,13 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
|
|||
&&
|
||||
!(this->current_x_is_feasible() && this->m_look_for_feasible_solution_only));
|
||||
|
||||
lp_assert(this->get_status() == lp_status::FLOATING_POINT_ERROR
|
||||
||
|
||||
lp_assert(
|
||||
this->current_x_is_feasible() == false
|
||||
||
|
||||
this->calc_current_x_is_feasible_include_non_basis());
|
||||
return this->total_iterations();
|
||||
}
|
||||
|
||||
template <typename T, typename X> void lp_primal_core_solver<T, X>::delete_factorization() {
|
||||
lp_assert(false);
|
||||
}
|
||||
|
||||
// according to Swietanowski, " A new steepest edge approximation for the simplex method for linear programming"
|
||||
template <typename T, typename X> void lp_primal_core_solver<T, X>::init_column_norms() {
|
||||
lp_assert(numeric_traits<T>::precise() == false);
|
||||
|
|
@ -860,7 +855,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::fill_breakpo
|
|||
|
||||
|
||||
template <typename T, typename X> bool lp_primal_core_solver<T, X>::done() {
|
||||
if (this->get_status() == lp_status::OPTIMAL || this->get_status() == lp_status::FLOATING_POINT_ERROR) return true;
|
||||
if (this->get_status() == lp_status::OPTIMAL) return true;
|
||||
if (this->get_status() == lp_status::INFEASIBLE) {
|
||||
return true;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -157,8 +157,7 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
|
|||
this->set_status(lp_status::CANCELLED);
|
||||
break; // from the loop
|
||||
}
|
||||
} while (this->get_status() != lp_status::FLOATING_POINT_ERROR
|
||||
&&
|
||||
} while (
|
||||
this->get_status() != lp_status::UNBOUNDED
|
||||
&&
|
||||
this->get_status() != lp_status::OPTIMAL
|
||||
|
|
@ -168,8 +167,7 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
|
|||
!(this->current_x_is_feasible() && this->m_look_for_feasible_solution_only)
|
||||
);
|
||||
|
||||
lp_assert(this->get_status() == lp_status::FLOATING_POINT_ERROR
|
||||
||
|
||||
lp_assert(
|
||||
this->get_status() == lp_status::CANCELLED
|
||||
||
|
||||
this->current_x_is_feasible() == false
|
||||
|
|
|
|||
|
|
@ -69,7 +69,6 @@ enum class lp_status {
|
|||
DUAL_UNBOUNDED,
|
||||
OPTIMAL,
|
||||
FEASIBLE,
|
||||
FLOATING_POINT_ERROR,
|
||||
TIME_EXHAUSTED,
|
||||
EMPTY,
|
||||
UNSTABLE,
|
||||
|
|
|
|||
|
|
@ -45,7 +45,6 @@ const char* lp_status_to_string(lp_status status) {
|
|||
case lp_status::DUAL_UNBOUNDED: return "DUAL_UNBOUNDED";
|
||||
case lp_status::OPTIMAL: return "OPTIMAL";
|
||||
case lp_status::FEASIBLE: return "FEASIBLE";
|
||||
case lp_status::FLOATING_POINT_ERROR: return "FLOATING_POINT_ERROR";
|
||||
case lp_status::TIME_EXHAUSTED: return "TIME_EXHAUSTED";
|
||||
case lp_status::EMPTY: return "EMPTY";
|
||||
case lp_status::UNSTABLE: return "UNSTABLE";
|
||||
|
|
@ -62,7 +61,6 @@ lp_status lp_status_from_string(std::string status) {
|
|||
if (status == "UNBOUNDED") return lp_status::UNBOUNDED;
|
||||
if (status == "OPTIMAL") return lp_status::OPTIMAL;
|
||||
if (status == "FEASIBLE") return lp_status::FEASIBLE;
|
||||
if (status == "FLOATING_POINT_ERROR") return lp_status::FLOATING_POINT_ERROR;
|
||||
if (status == "TIME_EXHAUSTED") return lp_status::TIME_EXHAUSTED;
|
||||
if (status == "EMPTY") return lp_status::EMPTY;
|
||||
lp_unreachable();
|
||||
|
|
|
|||
|
|
@ -1,81 +0,0 @@
|
|||
/*++
|
||||
Copyright (c) 2017 Microsoft Corporation
|
||||
|
||||
Module Name:
|
||||
|
||||
<name>
|
||||
|
||||
Abstract:
|
||||
|
||||
<abstract>
|
||||
|
||||
Author:
|
||||
|
||||
Lev Nachmanson (levnach)
|
||||
|
||||
Revision History:
|
||||
|
||||
|
||||
--*/
|
||||
#include <utility>
|
||||
#include <memory>
|
||||
#include <string>
|
||||
#include "util/vector.h"
|
||||
#include "util/debug.h"
|
||||
#include "math/lp/lu_def.h"
|
||||
namespace lp {
|
||||
template double dot_product<double, double>(vector<double> const&, vector<double> const&);
|
||||
template lu<static_matrix<double, double>>::lu(static_matrix<double, double> const&, vector<unsigned int>&, lp_settings&);
|
||||
template void lu<static_matrix<double, double>>::push_matrix_to_tail(tail_matrix<double, double>*);
|
||||
template void lu<static_matrix<double, double>>::replace_column(double, indexed_vector<double>&, unsigned);
|
||||
template lu<static_matrix<double, double>>::~lu();
|
||||
template void lu<static_matrix<mpq, mpq>>::push_matrix_to_tail(tail_matrix<mpq, mpq>*);
|
||||
template lu<static_matrix<mpq, mpq>>::~lu();
|
||||
template void lu<static_matrix<mpq, impq>>::push_matrix_to_tail(tail_matrix<mpq, impq >*);
|
||||
template lu<static_matrix<mpq, impq>>::~lu();
|
||||
template mpq dot_product<mpq, mpq>(vector<mpq > const&, vector<mpq > const&);
|
||||
template void init_factorization<static_matrix<double, double>>
|
||||
(lu<static_matrix<double, double>>*&, static_matrix<double, double>&, vector<unsigned int>&, lp_settings&);
|
||||
template void init_factorization<static_matrix<mpq, mpq>>
|
||||
(lu<static_matrix<mpq,mpq>>*&, static_matrix<mpq, mpq>&, vector<unsigned int>&, lp_settings&);
|
||||
template void init_factorization<static_matrix<mpq, impq>>(lu<static_matrix<mpq, impq> >*&, static_matrix<mpq, impq >&, vector<unsigned int>&, lp_settings&);
|
||||
template void print_matrix<square_sparse_matrix<double, double>>(square_sparse_matrix<double, double>&, std::ostream & out);
|
||||
template void print_matrix<static_matrix<mpq,mpq>>(static_matrix<mpq, mpq>&, std::ostream&);
|
||||
template void print_matrix<static_matrix<mpq, impq> >(static_matrix<mpq, impq >&, std::ostream&);
|
||||
template void print_matrix<static_matrix<double, double>>(static_matrix<double, double>&, std::ostream & out);
|
||||
#ifdef Z3DEBUG
|
||||
template bool lu<static_matrix<double, double>>::is_correct(const vector<unsigned>& basis);
|
||||
template bool lu<static_matrix<mpq, impq>>::is_correct( vector<unsigned int> const &);
|
||||
template dense_matrix<double, double> get_B<static_matrix<double, double>>(lu<static_matrix<double, double>>&, const vector<unsigned>& basis);
|
||||
template dense_matrix<mpq, mpq> get_B<static_matrix<mpq, mpq>>(lu<static_matrix<mpq, mpq>>&, vector<unsigned int> const&);
|
||||
|
||||
#endif
|
||||
|
||||
template bool lu<static_matrix<double, double>>::pivot_the_row(int); // NOLINT
|
||||
template void lu<static_matrix<double, double>>::init_vector_w(unsigned int, indexed_vector<double>&);
|
||||
template void lu<static_matrix<double, double>>::solve_By(vector<double>&);
|
||||
template void lu<static_matrix<double, double>>::solve_By_when_y_is_ready_for_X(vector<double>&);
|
||||
template void lu<static_matrix<double, double>>::solve_yB_with_error_check(vector<double>&, const vector<unsigned>& basis);
|
||||
template void lu<static_matrix<double, double>>::solve_yB_with_error_check_indexed(indexed_vector<double>&, vector<int> const&, const vector<unsigned> & basis, const lp_settings&);
|
||||
template void lu<static_matrix<mpq, mpq>>::replace_column(mpq, indexed_vector<mpq>&, unsigned);
|
||||
template void lu<static_matrix<mpq, mpq>>::solve_By(vector<mpq >&);
|
||||
template void lu<static_matrix<mpq, mpq>>::solve_By_when_y_is_ready_for_X(vector<mpq >&);
|
||||
template void lu<static_matrix<mpq, mpq>>::solve_yB_with_error_check(vector<mpq >&, const vector<unsigned>& basis);
|
||||
template void lu<static_matrix<mpq, mpq>>::solve_yB_with_error_check_indexed(indexed_vector<mpq>&, vector< int > const&, const vector<unsigned> & basis, const lp_settings&);
|
||||
template void lu<static_matrix<mpq, impq> >::solve_yB_with_error_check_indexed(indexed_vector<mpq>&, vector< int > const&, const vector<unsigned> & basis, const lp_settings&);
|
||||
template void lu<static_matrix<mpq, impq> >::init_vector_w(unsigned int, indexed_vector<mpq>&);
|
||||
template void lu<static_matrix<mpq, impq> >::replace_column(mpq, indexed_vector<mpq>&, unsigned);
|
||||
template void lu<static_matrix<mpq, impq> >::solve_Bd_faster(unsigned int, indexed_vector<mpq>&);
|
||||
template void lu<static_matrix<mpq, impq> >::solve_By(vector<impq >&);
|
||||
template void lu<static_matrix<mpq, impq> >::solve_By_when_y_is_ready_for_X(vector<impq >&);
|
||||
template void lu<static_matrix<mpq, impq> >::solve_yB_with_error_check(vector<mpq >&, const vector<unsigned>& basis);
|
||||
template void lu<static_matrix<mpq, mpq>>::solve_By(indexed_vector<mpq>&);
|
||||
template void lu<static_matrix<double, double>>::solve_By(indexed_vector<double>&);
|
||||
template void lu<static_matrix<double, double>>::solve_yB_indexed(indexed_vector<double>&);
|
||||
template void lu<static_matrix<mpq, impq> >::solve_yB_indexed(indexed_vector<mpq>&);
|
||||
template void lu<static_matrix<mpq, mpq>>::solve_By_for_T_indexed_only(indexed_vector<mpq>&, lp_settings const&);
|
||||
template void lu<static_matrix<double, double>>::solve_By_for_T_indexed_only(indexed_vector<double>&, lp_settings const&);
|
||||
#ifdef Z3DEBUG
|
||||
template void print_matrix<tail_matrix<double, double>>(tail_matrix<double, double>&, std::ostream&);
|
||||
#endif
|
||||
}
|
||||
325
src/math/lp/lu.h
325
src/math/lp/lu.h
|
|
@ -1,325 +0,0 @@
|
|||
/*++
|
||||
Copyright (c) 2017 Microsoft Corporation
|
||||
|
||||
Module Name:
|
||||
|
||||
<name>
|
||||
|
||||
Abstract:
|
||||
|
||||
<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 <algorithm>
|
||||
#include <set>
|
||||
#include "math/lp/square_sparse_matrix.h"
|
||||
#include "math/lp/static_matrix.h"
|
||||
#include <string>
|
||||
#include "math/lp/numeric_pair.h"
|
||||
#include <ostream>
|
||||
#include <fstream>
|
||||
#include "math/lp/row_eta_matrix.h"
|
||||
#include "math/lp/square_dense_submatrix.h"
|
||||
#include "math/lp/dense_matrix.h"
|
||||
namespace lp {
|
||||
template <typename T, typename X> // print the nr x nc submatrix at the top left corner
|
||||
void print_submatrix(square_sparse_matrix<T, X> & m, unsigned mr, unsigned nc);
|
||||
|
||||
template <typename M>
|
||||
void print_matrix(M &m, std::ostream & out);
|
||||
|
||||
template <typename T, typename X>
|
||||
X dot_product(const vector<T> & a, const vector<X> & b) {
|
||||
lp_assert(a.size() == b.size());
|
||||
auto r = zero_of_type<X>();
|
||||
for (unsigned i = 0; i < a.size(); i++) {
|
||||
r += a[i] * b[i];
|
||||
}
|
||||
return r;
|
||||
}
|
||||
|
||||
|
||||
template <typename T, typename X>
|
||||
class one_elem_on_diag: public tail_matrix<T, X> {
|
||||
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<T>::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<X> & w, lp_settings &) override {
|
||||
w[m_i] /= m_val;
|
||||
}
|
||||
|
||||
void apply_from_right(vector<T> & w) override {
|
||||
w[m_i] /= m_val;
|
||||
}
|
||||
|
||||
void apply_from_right(indexed_vector<T> & 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<T>();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void apply_from_left_to_T(indexed_vector<T> & w, lp_settings & settings) override;
|
||||
|
||||
void conjugate_by_permutation(permutation_matrix<T, X> & 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 <typename M>
|
||||
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<T, X> m_Q;
|
||||
permutation_matrix<T, X> m_R;
|
||||
permutation_matrix<T, X> m_r_wave;
|
||||
square_sparse_matrix<T, X> m_U;
|
||||
square_dense_submatrix<T, X>* m_dense_LU;
|
||||
|
||||
vector<tail_matrix<T, X> *> m_tail;
|
||||
lp_settings & m_settings;
|
||||
bool m_failure;
|
||||
indexed_vector<T> m_row_eta_work_vector;
|
||||
indexed_vector<T> m_w_for_extension;
|
||||
indexed_vector<T> m_y_copy;
|
||||
indexed_vector<unsigned> 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<unsigned>& basis,
|
||||
lp_settings & settings);
|
||||
lu(const M & A, lp_settings&);
|
||||
void debug_test_of_basis(const M & A, vector<unsigned> & basis);
|
||||
void solve_Bd_when_w_is_ready(vector<T> & d, indexed_vector<T>& w );
|
||||
void solve_By(indexed_vector<X> & y);
|
||||
|
||||
void solve_By(vector<X> & y);
|
||||
|
||||
void solve_By_for_T_indexed_only(indexed_vector<T>& y, const lp_settings &);
|
||||
|
||||
template <typename L>
|
||||
void solve_By_when_y_is_ready(indexed_vector<L> & y);
|
||||
void solve_By_when_y_is_ready_for_X(vector<X> & y);
|
||||
void solve_By_when_y_is_ready_for_T(vector<T> & y, vector<unsigned> & index);
|
||||
void print_indexed_vector(indexed_vector<T> & w, std::ofstream & f);
|
||||
|
||||
void print_matrix_compact(std::ostream & f);
|
||||
|
||||
void print(indexed_vector<T> & w, const vector<unsigned>& basis);
|
||||
void solve_Bd_faster(unsigned a_column, indexed_vector<T> & d); // d is the right side on the input and the solution at the exit
|
||||
|
||||
|
||||
void solve_yB_indexed(indexed_vector<T>& y);
|
||||
|
||||
void add_delta_to_solution_indexed(indexed_vector<T>& y);
|
||||
|
||||
void add_delta_to_solution(const vector<T>& yc, vector<T>& y);
|
||||
|
||||
|
||||
void find_error_of_yB(vector<T>& yc, const vector<T>& y,
|
||||
const vector<unsigned>& basis);
|
||||
|
||||
void find_error_of_yB_indexed(const indexed_vector<T>& y,
|
||||
const vector<int>& heading, const lp_settings& settings);
|
||||
|
||||
|
||||
void solve_yB_with_error_check(vector<T> & y, const vector<unsigned>& basis);
|
||||
|
||||
void solve_yB_with_error_check_indexed(indexed_vector<T> & y, const vector<int>& heading, const vector<unsigned> & basis, const lp_settings &);
|
||||
|
||||
void apply_Q_R_to_U(permutation_matrix<T, X> & r_wave);
|
||||
|
||||
|
||||
LU_status get_status() { return m_status; }
|
||||
|
||||
void set_status(LU_status status) {
|
||||
m_status = status;
|
||||
}
|
||||
|
||||
~lu();
|
||||
|
||||
void init_vector_y(vector<X> & y);
|
||||
|
||||
void perform_transformations_on_w(indexed_vector<T>& w);
|
||||
|
||||
void init_vector_w(unsigned entering, indexed_vector<T> & w);
|
||||
void apply_lp_list_to_w(indexed_vector<T> & w);
|
||||
void apply_lp_list_to_y(vector<X>& y);
|
||||
|
||||
void swap_rows(int j, int k);
|
||||
|
||||
void swap_columns(int j, int pivot_column);
|
||||
|
||||
void push_matrix_to_tail(tail_matrix<T, X>* tm) {
|
||||
m_tail.push_back(tm);
|
||||
}
|
||||
|
||||
bool pivot_the_row(int row);
|
||||
|
||||
eta_matrix<T, X> * get_eta_matrix_for_pivot(unsigned j);
|
||||
// we're processing the column j now
|
||||
eta_matrix<T, X> * get_eta_matrix_for_pivot(unsigned j, square_sparse_matrix<T, X>& copy_of_U);
|
||||
|
||||
// see page 407 of Chvatal
|
||||
unsigned transform_U_to_V_by_replacing_column(indexed_vector<T> & w, unsigned leaving_column_of_U);
|
||||
|
||||
#ifdef Z3DEBUG
|
||||
void check_vector_w(unsigned entering);
|
||||
|
||||
void check_apply_matrix_to_vector(matrix<T, X> *lp, T *w);
|
||||
|
||||
void check_apply_lp_lists_to_w(T * w);
|
||||
|
||||
// provide some access operators for testing
|
||||
permutation_matrix<T, X> & Q() { return m_Q; }
|
||||
permutation_matrix<T, X> & R() { return m_R; }
|
||||
matrix<T, X> & U() { return m_U; }
|
||||
unsigned tail_size() { return m_tail.size(); }
|
||||
|
||||
tail_matrix<T, X> * get_lp_matrix(unsigned i) {
|
||||
return m_tail[i];
|
||||
}
|
||||
|
||||
T B_(unsigned i, unsigned j, const vector<unsigned>& 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<unsigned>& basis);
|
||||
bool is_correct();
|
||||
|
||||
|
||||
|
||||
// needed for debugging purposes
|
||||
void copy_w(T *buffer, indexed_vector<T> & w);
|
||||
|
||||
// needed for debugging purposes
|
||||
void restore_w(T *buffer, indexed_vector<T> & 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<T, X> & 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<T> & 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<T> & 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 <typename M>
|
||||
void init_factorization(lu<M>* & factorization, M & m_A, vector<unsigned> & m_basis, lp_settings &m_settings);
|
||||
|
||||
#ifdef Z3DEBUG
|
||||
template <typename T, typename X, typename M>
|
||||
dense_matrix<T, X> get_B(lu<M>& f, const vector<unsigned>& basis);
|
||||
|
||||
template <typename T, typename X, typename M>
|
||||
dense_matrix<T, X> get_B(lu<M>& f);
|
||||
#endif
|
||||
}
|
||||
|
|
@ -1,720 +0,0 @@
|
|||
/*++
|
||||
Copyright (c) 2017 Microsoft Corporation
|
||||
|
||||
Module Name:
|
||||
|
||||
<name>
|
||||
|
||||
Abstract:
|
||||
|
||||
<abstract>
|
||||
|
||||
Author:
|
||||
|
||||
Lev Nachmanson (levnach)
|
||||
|
||||
Revision History:
|
||||
|
||||
|
||||
--*/
|
||||
#pragma once
|
||||
|
||||
#include <string>
|
||||
#include <algorithm>
|
||||
#include <set>
|
||||
#include "util/vector.h"
|
||||
#include <utility>
|
||||
#include "util/debug.h"
|
||||
#include "math/lp/lu.h"
|
||||
namespace lp {
|
||||
template <typename T, typename X, typename M> // print the nr x nc submatrix at the top left corner
|
||||
void print_submatrix(square_sparse_matrix<T, X> & m, unsigned mr, unsigned nc, std::ostream & out) {
|
||||
vector<vector<std::string>> A;
|
||||
vector<unsigned> widths;
|
||||
for (unsigned i = 0; i < m.row_count() && i < mr ; i++) {
|
||||
A.push_back(vector<std::string>());
|
||||
for (unsigned j = 0; j < m.column_count() && j < nc; j++) {
|
||||
A[i].push_back(T_to_string(static_cast<T>(m(i, j))));
|
||||
}
|
||||
}
|
||||
|
||||
for (unsigned j = 0; j < m.column_count() && j < nc; j++) {
|
||||
widths.push_back(get_width_of_column(j, A));
|
||||
}
|
||||
|
||||
print_matrix_with_widths(A, widths, out);
|
||||
}
|
||||
|
||||
template<typename M>
|
||||
void print_matrix(M &m, std::ostream & out) {
|
||||
vector<vector<std::string>> A;
|
||||
vector<unsigned> widths;
|
||||
for (unsigned i = 0; i < m.row_count(); i++) {
|
||||
A.push_back(vector<std::string>());
|
||||
for (unsigned j = 0; j < m.column_count(); j++) {
|
||||
A[i].push_back(T_to_string(m[i][j]));
|
||||
}
|
||||
}
|
||||
|
||||
for (unsigned j = 0; j < m.column_count(); j++) {
|
||||
widths.push_back(get_width_of_column(j, A));
|
||||
}
|
||||
|
||||
print_matrix_with_widths(A, widths, out);
|
||||
}
|
||||
|
||||
template <typename T, typename X>
|
||||
one_elem_on_diag<T, X>::one_elem_on_diag(const one_elem_on_diag & o) {
|
||||
m_i = o.m_i;
|
||||
m_val = o.m_val;
|
||||
#ifdef Z3DEBUG
|
||||
m_m = m_n = o.m_m;
|
||||
m_one_over_val = numeric_traits<T>::one() / o.m_val;
|
||||
#endif
|
||||
}
|
||||
|
||||
#ifdef Z3DEBUG
|
||||
template <typename T, typename X>
|
||||
T one_elem_on_diag<T, X>::get_elem(unsigned i, unsigned j) const {
|
||||
if (i == j){
|
||||
if (j == m_i) {
|
||||
return m_one_over_val;
|
||||
}
|
||||
return numeric_traits<T>::one();
|
||||
}
|
||||
|
||||
return numeric_traits<T>::zero();
|
||||
}
|
||||
#endif
|
||||
template <typename T, typename X>
|
||||
void one_elem_on_diag<T, X>::apply_from_left_to_T(indexed_vector<T> & w, lp_settings & settings) {
|
||||
T & t = w[m_i];
|
||||
if (numeric_traits<T>::is_zero(t)) {
|
||||
return;
|
||||
}
|
||||
t /= m_val;
|
||||
if (numeric_traits<T>::precise()) return;
|
||||
if (settings.abs_val_is_smaller_than_drop_tolerance(t)) {
|
||||
w.erase_from_index(m_i);
|
||||
t = numeric_traits<T>::zero();
|
||||
}
|
||||
}
|
||||
|
||||
// This class supports updates of the columns of B, and solves systems Bx=b,and yB=c
|
||||
// Using Suhl-Suhl method described in the dissertation of Achim Koberstein, Chapter 5
|
||||
template <typename M>
|
||||
lu<M>::lu(const M& A,
|
||||
vector<unsigned>& basis,
|
||||
lp_settings & settings):
|
||||
m_status(LU_status::OK),
|
||||
m_dim(A.row_count()),
|
||||
m_A(A),
|
||||
m_Q(m_dim),
|
||||
m_R(m_dim),
|
||||
m_r_wave(m_dim),
|
||||
m_U(A, basis), // create the square matrix that eventually will be factorized
|
||||
m_settings(settings),
|
||||
m_failure(false),
|
||||
m_row_eta_work_vector(A.row_count()),
|
||||
m_refactor_counter(0) {
|
||||
lp_assert(!(numeric_traits<T>::precise() ));
|
||||
#ifdef Z3DEBUG
|
||||
debug_test_of_basis(A, basis);
|
||||
#endif
|
||||
++m_settings.stats().m_num_factorizations;
|
||||
create_initial_factorization();
|
||||
#ifdef Z3DEBUG
|
||||
// lp_assert(check_correctness());
|
||||
#endif
|
||||
}
|
||||
template <typename M>
|
||||
lu<M>::lu(const M& A,
|
||||
lp_settings & settings):
|
||||
m_status(LU_status::OK),
|
||||
m_dim(A.row_count()),
|
||||
m_A(A),
|
||||
m_Q(m_dim),
|
||||
m_R(m_dim),
|
||||
m_r_wave(m_dim),
|
||||
m_U(A), // create the square matrix that eventually will be factorized
|
||||
m_settings(settings),
|
||||
m_failure(false),
|
||||
m_row_eta_work_vector(A.row_count()),
|
||||
m_refactor_counter(0) {
|
||||
lp_assert(A.row_count() == A.column_count());
|
||||
create_initial_factorization();
|
||||
#ifdef Z3DEBUG
|
||||
lp_assert(is_correct());
|
||||
#endif
|
||||
}
|
||||
template <typename M>
|
||||
void lu<M>::debug_test_of_basis( M const & A, vector<unsigned> & basis) {
|
||||
std::set<unsigned> set;
|
||||
for (unsigned i = 0; i < A.row_count(); i++) {
|
||||
lp_assert(basis[i]< A.column_count());
|
||||
set.insert(basis[i]);
|
||||
}
|
||||
lp_assert(set.size() == A.row_count());
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu<M>::solve_By(indexed_vector<X> & y) {
|
||||
lp_assert(false); // not implemented
|
||||
// init_vector_y(y);
|
||||
// solve_By_when_y_is_ready(y);
|
||||
}
|
||||
|
||||
|
||||
template <typename M>
|
||||
void lu<M>::solve_By(vector<X> & y) {
|
||||
init_vector_y(y);
|
||||
solve_By_when_y_is_ready_for_X(y);
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu<M>::solve_By_when_y_is_ready_for_X(vector<X> & y) {
|
||||
if (numeric_traits<T>::precise()) {
|
||||
m_U.solve_U_y(y);
|
||||
m_R.apply_reverse_from_left_to_X(y); // see 24.3 from Chvatal
|
||||
return;
|
||||
}
|
||||
m_U.double_solve_U_y(y);
|
||||
m_R.apply_reverse_from_left_to_X(y); // see 24.3 from Chvatal
|
||||
unsigned i = m_dim;
|
||||
while (i--) {
|
||||
if (is_zero(y[i])) continue;
|
||||
if (m_settings.abs_val_is_smaller_than_drop_tolerance(y[i])){
|
||||
y[i] = zero_of_type<X>();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu<M>::solve_By_when_y_is_ready_for_T(vector<T> & y, vector<unsigned> & index) {
|
||||
if (numeric_traits<T>::precise()) {
|
||||
m_U.solve_U_y(y);
|
||||
m_R.apply_reverse_from_left_to_T(y); // see 24.3 from Chvatal
|
||||
unsigned j = m_dim;
|
||||
while (j--) {
|
||||
if (!is_zero(y[j]))
|
||||
index.push_back(j);
|
||||
}
|
||||
return;
|
||||
}
|
||||
m_U.double_solve_U_y(y);
|
||||
m_R.apply_reverse_from_left_to_T(y); // see 24.3 from Chvatal
|
||||
unsigned i = m_dim;
|
||||
while (i--) {
|
||||
if (is_zero(y[i])) continue;
|
||||
if (m_settings.abs_val_is_smaller_than_drop_tolerance(y[i])){
|
||||
y[i] = zero_of_type<T>();
|
||||
} else {
|
||||
index.push_back(i);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu<M>::solve_By_for_T_indexed_only(indexed_vector<T> & y, const lp_settings & settings) {
|
||||
if (numeric_traits<T>::precise()) {
|
||||
vector<unsigned> active_rows;
|
||||
m_U.solve_U_y_indexed_only(y, settings, active_rows);
|
||||
m_R.apply_reverse_from_left(y); // see 24.3 from Chvatal
|
||||
return;
|
||||
}
|
||||
m_U.double_solve_U_y(y, m_settings);
|
||||
m_R.apply_reverse_from_left(y); // see 24.3 from Chvatal
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu<M>::print_matrix_compact(std::ostream & f) {
|
||||
f << "matrix_start" << std::endl;
|
||||
f << "nrows " << m_A.row_count() << std::endl;
|
||||
f << "ncolumns " << m_A.column_count() << std::endl;
|
||||
for (unsigned i = 0; i < m_A.row_count(); i++) {
|
||||
auto & row = m_A.m_rows[i];
|
||||
f << "row " << i << std::endl;
|
||||
for (auto & t : row) {
|
||||
f << "column " << t.m_j << " value " << t.m_value << std::endl;
|
||||
}
|
||||
f << "row_end" << std::endl;
|
||||
}
|
||||
f << "matrix_end" << std::endl;
|
||||
}
|
||||
template <typename M>
|
||||
void lu< M>::print(indexed_vector<T> & w, const vector<unsigned>& basis) {
|
||||
std::string dump_file_name("/tmp/lu");
|
||||
remove(dump_file_name.c_str());
|
||||
std::ofstream f(dump_file_name);
|
||||
if (!f.is_open()) {
|
||||
LP_OUT(m_settings, "cannot open file " << dump_file_name << std::endl);
|
||||
return;
|
||||
}
|
||||
LP_OUT(m_settings, "writing lu dump to " << dump_file_name << std::endl);
|
||||
print_matrix_compact(f);
|
||||
print_vector(basis, f);
|
||||
print_indexed_vector(w, f);
|
||||
f.close();
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu< M>::solve_Bd_faster(unsigned a_column, indexed_vector<T> & d) { // puts the a_column into d
|
||||
init_vector_w(a_column, d);
|
||||
solve_By_for_T_indexed_only(d, m_settings);
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu< M>::solve_yB_indexed(indexed_vector<T>& y) {
|
||||
lp_assert(y.is_OK());
|
||||
// first solve yU = cb*R(-1)
|
||||
m_R.apply_reverse_from_right_to_T(y); // got y = cb*R(-1)
|
||||
lp_assert(y.is_OK());
|
||||
m_U.solve_y_U_indexed(y, m_settings); // got y*U=cb*R(-1)
|
||||
lp_assert(y.is_OK());
|
||||
m_Q.apply_reverse_from_right_to_T(y);
|
||||
lp_assert(y.is_OK());
|
||||
for (auto e = m_tail.rbegin(); e != m_tail.rend(); ++e) {
|
||||
#ifdef Z3DEBUG
|
||||
(*e)->set_number_of_columns(m_dim);
|
||||
#endif
|
||||
(*e)->apply_from_right(y);
|
||||
lp_assert(y.is_OK());
|
||||
}
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu< M>::add_delta_to_solution(const vector<T>& yc, vector<T>& y){
|
||||
unsigned i = static_cast<unsigned>(y.size());
|
||||
while (i--)
|
||||
y[i]+=yc[i];
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu< M>::add_delta_to_solution_indexed(indexed_vector<T>& y) {
|
||||
// the delta sits in m_y_copy, put result into y
|
||||
lp_assert(y.is_OK());
|
||||
lp_assert(m_y_copy.is_OK());
|
||||
m_ii.clear();
|
||||
m_ii.resize(y.data_size());
|
||||
for (unsigned i : y.m_index)
|
||||
m_ii.set_value(1, i);
|
||||
for (unsigned i : m_y_copy.m_index) {
|
||||
y.m_data[i] += m_y_copy[i];
|
||||
if (m_ii[i] == 0)
|
||||
m_ii.set_value(1, i);
|
||||
}
|
||||
lp_assert(m_ii.is_OK());
|
||||
y.m_index.clear();
|
||||
|
||||
for (unsigned i : m_ii.m_index) {
|
||||
T & v = y.m_data[i];
|
||||
if (!lp_settings::is_eps_small_general(v, 1e-14))
|
||||
y.m_index.push_back(i);
|
||||
else if (!numeric_traits<T>::is_zero(v))
|
||||
v = zero_of_type<T>();
|
||||
}
|
||||
|
||||
lp_assert(y.is_OK());
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu< M>::find_error_of_yB(vector<T>& yc, const vector<T>& y, const vector<unsigned>& m_basis) {
|
||||
unsigned i = m_dim;
|
||||
while (i--) {
|
||||
yc[i] -= m_A.dot_product_with_column(y, m_basis[i]);
|
||||
}
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu< M>::find_error_of_yB_indexed(const indexed_vector<T>& y, const vector<int>& heading, const lp_settings& settings) {
|
||||
#if 0 == 1
|
||||
// it is a non efficient version
|
||||
indexed_vector<T> yc = m_y_copy;
|
||||
yc.m_index.clear();
|
||||
lp_assert(!numeric_traits<T>::precise());
|
||||
{
|
||||
|
||||
vector<unsigned> d_basis(y.m_data.size());
|
||||
for (unsigned j = 0; j < heading.size(); j++) {
|
||||
if (heading[j] >= 0) {
|
||||
d_basis[heading[j]] = j;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
unsigned i = m_dim;
|
||||
while (i--) {
|
||||
T & v = yc.m_data[i] -= m_A.dot_product_with_column(y.m_data, d_basis[i]);
|
||||
if (settings.abs_val_is_smaller_than_drop_tolerance(v))
|
||||
v = zero_of_type<T>();
|
||||
else
|
||||
yc.m_index.push_back(i);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
lp_assert(m_ii.is_OK());
|
||||
m_ii.clear();
|
||||
m_ii.resize(y.data_size());
|
||||
lp_assert(m_y_copy.is_OK());
|
||||
// put the error into m_y_copy
|
||||
for (auto k : y.m_index) {
|
||||
auto & row = m_A.m_rows[k];
|
||||
const T & y_k = y.m_data[k];
|
||||
for (auto & c : row) {
|
||||
unsigned j = c.var();
|
||||
int hj = heading[j];
|
||||
if (hj < 0) continue;
|
||||
if (m_ii.m_data[hj] == 0)
|
||||
m_ii.set_value(1, hj);
|
||||
m_y_copy.m_data[hj] -= c.coeff() * y_k;
|
||||
}
|
||||
}
|
||||
// add the index of m_y_copy to m_ii
|
||||
for (unsigned i : m_y_copy.m_index) {
|
||||
if (m_ii.m_data[i] == 0)
|
||||
m_ii.set_value(1, i);
|
||||
}
|
||||
|
||||
// there is no guarantee that m_y_copy is OK here, but its index
|
||||
// is contained in m_ii index
|
||||
m_y_copy.m_index.clear();
|
||||
// setup the index of m_y_copy
|
||||
for (auto k : m_ii.m_index) {
|
||||
T& v = m_y_copy.m_data[k];
|
||||
if (settings.abs_val_is_smaller_than_drop_tolerance(v))
|
||||
v = zero_of_type<T>();
|
||||
else {
|
||||
m_y_copy.set_value(v, k);
|
||||
}
|
||||
}
|
||||
lp_assert(m_y_copy.is_OK());
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
// solves y*B = y
|
||||
// y is the input
|
||||
template <typename M>
|
||||
void lu< M>::solve_yB_with_error_check_indexed(indexed_vector<T> & y, const vector<int>& heading, const vector<unsigned> & basis, const lp_settings & settings) {
|
||||
lp_assert(false);
|
||||
}
|
||||
|
||||
|
||||
// solves y*B = y
|
||||
// y is the input
|
||||
template <typename M>
|
||||
void lu< M>::solve_yB_with_error_check(vector<T> & y, const vector<unsigned>& basis) {
|
||||
lp_assert(false);
|
||||
}
|
||||
template <typename M>
|
||||
void lu< M>::apply_Q_R_to_U(permutation_matrix<T, X> & r_wave) {
|
||||
m_U.multiply_from_right(r_wave);
|
||||
m_U.multiply_from_left_with_reverse(r_wave);
|
||||
}
|
||||
|
||||
|
||||
// Solving yB = cb to find the entering variable,
|
||||
// where cb is the cost vector projected to B.
|
||||
// The result is stored in cb.
|
||||
|
||||
// solving Bd = a ( to find the column d of B^{-1} A_N corresponding to the entering
|
||||
// variable
|
||||
template <typename M>
|
||||
lu< M>::~lu(){
|
||||
for (auto t : m_tail) {
|
||||
delete t;
|
||||
}
|
||||
}
|
||||
template <typename M>
|
||||
void lu< M>::init_vector_y(vector<X> & y) {
|
||||
apply_lp_list_to_y(y);
|
||||
m_Q.apply_reverse_from_left_to_X(y);
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu< M>::perform_transformations_on_w(indexed_vector<T>& w) {
|
||||
apply_lp_list_to_w(w);
|
||||
m_Q.apply_reverse_from_left(w);
|
||||
// TBD does not compile: lp_assert(numeric_traits<T>::precise() || check_vector_for_small_values(w, m_settings));
|
||||
}
|
||||
|
||||
// see Chvatal 24.3
|
||||
template <typename M>
|
||||
void lu< M>::init_vector_w(unsigned entering, indexed_vector<T> & w) {
|
||||
w.clear();
|
||||
m_A.copy_column_to_indexed_vector(entering, w); // w = a, the column
|
||||
perform_transformations_on_w(w);
|
||||
}
|
||||
template <typename M>
|
||||
void lu< M>::apply_lp_list_to_w(indexed_vector<T> & 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 <typename M>
|
||||
void lu< M>::apply_lp_list_to_y(vector<X>& y) {
|
||||
for (unsigned i = 0; i < m_tail.size(); i++) {
|
||||
m_tail[i]->apply_from_left(y, m_settings);
|
||||
}
|
||||
}
|
||||
template <typename M>
|
||||
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 <typename M>
|
||||
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 <typename M>
|
||||
bool lu< M>::pivot_the_row(int row) {
|
||||
eta_matrix<T, X> * eta_matrix = get_eta_matrix_for_pivot(row);
|
||||
if (get_status() != LU_status::OK) {
|
||||
return false;
|
||||
}
|
||||
|
||||
if (eta_matrix == nullptr) {
|
||||
m_U.shorten_active_matrix(row, nullptr);
|
||||
return true;
|
||||
}
|
||||
if (!m_U.pivot_with_eta(row, eta_matrix, m_settings))
|
||||
return false;
|
||||
eta_matrix->conjugate_by_permutation(m_Q);
|
||||
push_matrix_to_tail(eta_matrix);
|
||||
return true;
|
||||
}
|
||||
// we're processing the column j now
|
||||
template <typename M>
|
||||
eta_matrix<typename M::coefftype, typename M::argtype> * lu< M>::get_eta_matrix_for_pivot(unsigned j) {
|
||||
eta_matrix<T, X> *ret;
|
||||
if(!m_U.fill_eta_matrix(j, &ret)) {
|
||||
set_status(LU_status::Degenerated);
|
||||
}
|
||||
return ret;
|
||||
}
|
||||
// we're processing the column j now
|
||||
template <typename M>
|
||||
eta_matrix<typename M::coefftype, typename M::argtype> * lu<M>::get_eta_matrix_for_pivot(unsigned j, square_sparse_matrix<T, X>& copy_of_U) {
|
||||
eta_matrix<T, X> *ret;
|
||||
copy_of_U.fill_eta_matrix(j, &ret);
|
||||
return ret;
|
||||
}
|
||||
|
||||
// see page 407 of Chvatal
|
||||
template <typename M>
|
||||
unsigned lu<M>::transform_U_to_V_by_replacing_column(indexed_vector<T> & w,
|
||||
unsigned leaving_column) {
|
||||
unsigned column_to_replace = m_R.apply_reverse(leaving_column);
|
||||
m_U.replace_column(column_to_replace, w, m_settings);
|
||||
return column_to_replace;
|
||||
}
|
||||
|
||||
#ifdef Z3DEBUG
|
||||
template <typename M>
|
||||
void lu<M>::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 <typename M>
|
||||
void lu<M>::check_apply_matrix_to_vector(matrix<T, X> *lp, T *w) {
|
||||
if (lp != nullptr) {
|
||||
lp -> set_number_of_rows(m_dim);
|
||||
lp -> set_number_of_columns(m_dim);
|
||||
apply_to_vector(*lp, w);
|
||||
}
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu<M>::check_apply_lp_lists_to_w(T * w) {
|
||||
for (unsigned i = 0; i < m_tail.size(); i++) {
|
||||
check_apply_matrix_to_vector(m_tail[i], w);
|
||||
}
|
||||
permutation_matrix<T, X> qr = m_Q.get_reverse();
|
||||
apply_to_vector(qr, w);
|
||||
for (int i = m_dim - 1; i >= 0; i--) {
|
||||
lp_assert(abs(w[i] - w[i]) < 0.0000001);
|
||||
}
|
||||
}
|
||||
|
||||
#endif
|
||||
template <typename M>
|
||||
void lu<M>::process_column(int j) {
|
||||
unsigned pi, pj;
|
||||
bool success = m_U.get_pivot_for_column(pi, pj, m_settings.c_partial_pivoting, j);
|
||||
if (!success) {
|
||||
// LP_OUT(m_settings, "get_pivot returned false: cannot find the pivot for column " << j << std::endl);
|
||||
m_failure = true;
|
||||
return;
|
||||
}
|
||||
|
||||
if (static_cast<int>(pi) == -1) {
|
||||
// LP_OUT(m_settings, "cannot find the pivot for column " << j << std::endl);
|
||||
m_failure = true;
|
||||
return;
|
||||
}
|
||||
swap_columns(j, pj);
|
||||
swap_rows(j, pi);
|
||||
if (!pivot_the_row(j)) {
|
||||
// LP_OUT(m_settings, "pivot_the_row(" << j << ") failed" << std::endl);
|
||||
m_failure = true;
|
||||
}
|
||||
}
|
||||
template <typename M>
|
||||
bool lu<M>::is_correct(const vector<unsigned>& basis) {
|
||||
return true;
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
bool lu<M>::is_correct() {
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
// needed for debugging purposes
|
||||
template <typename M>
|
||||
void lu<M>::copy_w(T *buffer, indexed_vector<T> & w) {
|
||||
|
||||
}
|
||||
|
||||
// needed for debugging purposes
|
||||
template <typename M>
|
||||
void lu<M>::restore_w(T *buffer, indexed_vector<T> & w) {
|
||||
|
||||
}
|
||||
template <typename M>
|
||||
bool lu<M>::all_columns_and_rows_are_active() {
|
||||
return true;
|
||||
}
|
||||
template <typename M>
|
||||
bool lu<M>::too_dense(unsigned j) const {
|
||||
return false;
|
||||
}
|
||||
template <typename M>
|
||||
void lu<M>::pivot_in_dense_mode(unsigned i) {
|
||||
|
||||
}
|
||||
template <typename M>
|
||||
void lu<M>::create_initial_factorization(){
|
||||
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu<M>::calculate_r_wave_and_update_U(unsigned bump_start, unsigned bump_end, permutation_matrix<T, X> & r_wave) {
|
||||
if (bump_start > bump_end) {
|
||||
set_status(LU_status::Degenerated);
|
||||
return;
|
||||
}
|
||||
if (bump_start == bump_end) {
|
||||
return;
|
||||
}
|
||||
|
||||
r_wave[bump_start] = bump_end; // sending the offensive column to the end of the bump
|
||||
|
||||
for ( unsigned i = bump_start + 1 ; i <= bump_end; i++ ) {
|
||||
r_wave[i] = i - 1;
|
||||
}
|
||||
|
||||
m_U.multiply_from_right(r_wave);
|
||||
m_U.multiply_from_left_with_reverse(r_wave);
|
||||
}
|
||||
template <typename M>
|
||||
void lu<M>::scan_last_row_to_work_vector(unsigned lowest_row_of_the_bump) {
|
||||
vector<indexed_value<T>> & last_row_vec = m_U.get_row_values(m_U.adjust_row(lowest_row_of_the_bump));
|
||||
for (auto & iv : last_row_vec) {
|
||||
if (is_zero(iv.m_value)) continue;
|
||||
lp_assert(!m_settings.abs_val_is_smaller_than_drop_tolerance(iv.m_value));
|
||||
unsigned adjusted_col = m_U.adjust_column_inverse(iv.m_index);
|
||||
if (adjusted_col < lowest_row_of_the_bump) {
|
||||
m_row_eta_work_vector.set_value(-iv.m_value, adjusted_col);
|
||||
} else {
|
||||
m_row_eta_work_vector.set_value(iv.m_value, adjusted_col); // preparing to calculate the real value in the matrix
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu<M>::pivot_and_solve_the_system(unsigned replaced_column, unsigned lowest_row_of_the_bump) {
|
||||
// we have the system right side at m_row_eta_work_vector now
|
||||
// solve the system column wise
|
||||
for (unsigned j = replaced_column; j < lowest_row_of_the_bump; j++) {
|
||||
T v = m_row_eta_work_vector[j];
|
||||
if (numeric_traits<T>::is_zero(v)) continue; // this column does not contribute to the solution
|
||||
unsigned aj = m_U.adjust_row(j);
|
||||
vector<indexed_value<T>> & row = m_U.get_row_values(aj);
|
||||
for (auto & iv : row) {
|
||||
unsigned col = m_U.adjust_column_inverse(iv.m_index);
|
||||
lp_assert(col >= j || numeric_traits<T>::is_zero(iv.m_value));
|
||||
if (col == j) continue;
|
||||
if (numeric_traits<T>::is_zero(iv.m_value)) {
|
||||
continue;
|
||||
}
|
||||
// the -v is for solving the system ( to zero the last row), and +v is for pivoting
|
||||
T delta = col < lowest_row_of_the_bump? -v * iv.m_value: v * iv.m_value;
|
||||
lp_assert(numeric_traits<T>::is_zero(delta) == false);
|
||||
|
||||
|
||||
|
||||
// m_row_eta_work_vector.add_value_at_index_with_drop_tolerance(col, delta);
|
||||
if (numeric_traits<T>::is_zero(m_row_eta_work_vector[col])) {
|
||||
if (!m_settings.abs_val_is_smaller_than_drop_tolerance(delta)){
|
||||
m_row_eta_work_vector.set_value(delta, col);
|
||||
}
|
||||
} else {
|
||||
T t = (m_row_eta_work_vector[col] += delta);
|
||||
if (m_settings.abs_val_is_smaller_than_drop_tolerance(t)){
|
||||
m_row_eta_work_vector[col] = numeric_traits<T>::zero();
|
||||
auto it = std::find(m_row_eta_work_vector.m_index.begin(), m_row_eta_work_vector.m_index.end(), col);
|
||||
if (it != m_row_eta_work_vector.m_index.end())
|
||||
m_row_eta_work_vector.m_index.erase(it);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu<M>::replace_column(T pivot_elem_for_checking, indexed_vector<T> & w, unsigned leaving_column_of_U){
|
||||
lp_assert(false);
|
||||
}
|
||||
template <typename M>
|
||||
void lu<M>::calculate_Lwave_Pwave_for_bump(unsigned replaced_column, unsigned lowest_row_of_the_bump){
|
||||
lp_assert(false);// lp_assert(m_U.is_upper_triangular_and_maximums_are_set_correctly_in_rows(m_settings));
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void lu<M>::calculate_Lwave_Pwave_for_last_row(unsigned lowest_row_of_the_bump, T diagonal_element) {
|
||||
lp_assert(false);
|
||||
}
|
||||
|
||||
template <typename M>
|
||||
void init_factorization(lu<M>* & factorization, M & m_A, vector<unsigned> & m_basis, lp_settings &m_settings) {
|
||||
lp_assert(false);
|
||||
}
|
||||
|
||||
#ifdef Z3DEBUG
|
||||
template <typename M>
|
||||
dense_matrix<typename M::coefftype, typename M::argtype> get_B(lu<M>& f, const vector<unsigned>& basis) {
|
||||
lp_assert(false);
|
||||
|
||||
dense_matrix<typename M::coefftype, typename M::argtype> B(0, 0);
|
||||
return B;
|
||||
}
|
||||
template <typename M>
|
||||
dense_matrix<typename M::coefftype, typename M::argtype> get_B(lu<M>& f) {
|
||||
lp_assert(false);
|
||||
dense_matrix<typename M::coefftype, typename M::argtype> B(0,0);
|
||||
return B;
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
|
@ -20,7 +20,6 @@ Revision History:
|
|||
#include <memory>
|
||||
#include "util/vector.h"
|
||||
#include "math/lp/row_eta_matrix_def.h"
|
||||
#include "math/lp/lu.h"
|
||||
namespace lp {
|
||||
template void row_eta_matrix<double, double>::conjugate_by_permutation(permutation_matrix<double, double>&);
|
||||
template void row_eta_matrix<mpq, numeric_pair<mpq> >::conjugate_by_permutation(permutation_matrix<mpq, numeric_pair<mpq> >&);
|
||||
|
|
|
|||
|
|
@ -20,7 +20,6 @@ Revision History:
|
|||
#include <memory>
|
||||
#include "util/vector.h"
|
||||
#include "math/lp/lp_settings.h"
|
||||
#include "math/lp/lu.h"
|
||||
#include "math/lp/square_sparse_matrix_def.h"
|
||||
#include "math/lp/dense_matrix.h"
|
||||
namespace lp {
|
||||
|
|
|
|||
|
|
@ -21,6 +21,8 @@ Revision History:
|
|||
|
||||
#include "util/vector.h"
|
||||
#include "math/lp/square_sparse_matrix.h"
|
||||
#include "math/lp/dense_matrix.h"
|
||||
|
||||
#include <set>
|
||||
#include <queue>
|
||||
namespace lp {
|
||||
|
|
|
|||
|
|
@ -3126,7 +3126,7 @@ public:
|
|||
return l_false;
|
||||
TRACE("arith", tout << "status treated as inconclusive: " << status << "\n";);
|
||||
// TENTATIVE_UNBOUNDED, UNBOUNDED, TENTATIVE_DUAL_UNBOUNDED, DUAL_UNBOUNDED,
|
||||
// FLOATING_POINT_ERROR, TIME_EXAUSTED, EMPTY, UNSTABLE
|
||||
// TIME_EXAUSTED, EMPTY, UNSTABLE
|
||||
return l_undef;
|
||||
}
|
||||
|
||||
|
|
|
|||
Loading…
Reference in a new issue