mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-13 12:28:44 +00:00
Code Activity

rm lu

Browse source
This commit is contained in:
Lev Nachmanson 2023-03-06 13:40:21 -08:00
parent e04e726f45
commit 6eedbd4f35
5 changed files with 4 additions and 1232 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
src/math/lp/CMakeLists.txt
View file

@ -42,7 +42,6 @@ z3_add_component(lp
random_updater.cpp
row_eta_matrix.cpp
scaler.cpp
square_dense_submatrix.cpp
square_sparse_matrix.cpp
static_matrix.cpp
COMPONENT_DEPENDENCIES

48
src/math/lp/square_dense_submatrix.cpp
View file

@ -1,48 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <memory>
#include "util/vector.h"
#include "math/lp/square_dense_submatrix_def.h"
template void lp::square_dense_submatrix<double, double>::init(lp::square_sparse_matrix<double, double>*, unsigned int);
template lp::square_dense_submatrix<double, double>::square_dense_submatrix(lp::square_sparse_matrix<double, double>*, unsigned int);
template void lp::square_dense_submatrix<double, double>::update_parent_matrix(lp::lp_settings&);
template bool lp::square_dense_submatrix<double, double>::is_L_matrix() const;
template void lp::square_dense_submatrix<double, double>::conjugate_by_permutation(lp::permutation_matrix<double, double>&);
template int lp::square_dense_submatrix<double, double>::find_pivot_column_in_row(unsigned int) const;
template void lp::square_dense_submatrix<double, double>::pivot(unsigned int, lp::lp_settings&);
template lp::square_dense_submatrix<lp::mpq, lp::numeric_pair<lp::mpq> >::square_dense_submatrix(lp::square_sparse_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >*, unsigned int);
template void lp::square_dense_submatrix<lp::mpq, lp::numeric_pair<lp::mpq> >::update_parent_matrix(lp::lp_settings&);
template bool lp::square_dense_submatrix<lp::mpq, lp::numeric_pair<lp::mpq> >::is_L_matrix() const;
template void lp::square_dense_submatrix<lp::mpq, lp::numeric_pair<lp::mpq> >::conjugate_by_permutation(lp::permutation_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >&);
template int lp::square_dense_submatrix<lp::mpq, lp::numeric_pair<lp::mpq> >::find_pivot_column_in_row(unsigned int) const;
template void lp::square_dense_submatrix<lp::mpq, lp::numeric_pair<lp::mpq> >::pivot(unsigned int, lp::lp_settings&);
#ifdef Z3DEBUG
template double lp::square_dense_submatrix<double, double>::get_elem(unsigned int, unsigned int) const;
#endif
template void lp::square_dense_submatrix<double, double>::apply_from_right(vector<double>&);
template void lp::square_dense_submatrix<double, double>::apply_from_left_local<double>(lp::indexed_vector<double>&, lp::lp_settings&);
template void lp::square_dense_submatrix<double, double>::apply_from_left_to_vector<double>(vector<double>&);
template lp::square_dense_submatrix<lp::mpq, lp::mpq>::square_dense_submatrix(lp::square_sparse_matrix<lp::mpq, lp::mpq>*, unsigned int);
template void lp::square_dense_submatrix<lp::mpq, lp::mpq>::update_parent_matrix(lp::lp_settings&);
template bool lp::square_dense_submatrix<lp::mpq, lp::mpq>::is_L_matrix() const;
template void lp::square_dense_submatrix<lp::mpq, lp::mpq>::conjugate_by_permutation(lp::permutation_matrix<lp::mpq, lp::mpq>&);
template int lp::square_dense_submatrix<lp::mpq, lp::mpq>::find_pivot_column_in_row(unsigned int) const;
template void lp::square_dense_submatrix<lp::mpq, lp::mpq>::pivot(unsigned int, lp::lp_settings&);

225
src/math/lp/square_dense_submatrix.h
View file

@ -1,225 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
#include "math/lp/permutation_matrix.h"
#include <unordered_map>
#include "math/lp/static_matrix.h"
#include <set>
#include <utility>
#include <string>
#include <algorithm>
#include <queue>
#include "math/lp/indexed_value.h"
#include "math/lp/indexed_vector.h"
#include <functional>
#include "math/lp/lp_settings.h"
#include "math/lp/eta_matrix.h"
#include "math/lp/binary_heap_upair_queue.h"
#include "math/lp/square_sparse_matrix.h"
namespace lp {
template <typename T, typename X>
class square_dense_submatrix : public tail_matrix<T, X> {
// the submatrix uses the permutations of the parent matrix to access the elements
struct ref {
unsigned m_i_offset;
square_dense_submatrix & m_s;
ref(unsigned i, square_dense_submatrix & s) :
m_i_offset((i - s.m_index_start) * s.m_dim), m_s(s){}
T & operator[] (unsigned j) {
lp_assert(j >= m_s.m_index_start);
return m_s.m_v[m_i_offset + m_s.adjust_column(j) - m_s.m_index_start];
}
const T & operator[] (unsigned j) const {
lp_assert(j >= m_s.m_index_start);
return m_s.m_v[m_i_offset + m_s.adjust_column(j) - m_s.m_index_start];
}
};
public:
unsigned m_index_start;
unsigned m_dim;
vector<T> m_v;
square_sparse_matrix<T, X> * m_parent;
permutation_matrix<T, X> m_row_permutation;
indexed_vector<T> m_work_vector;
public:
permutation_matrix<T, X> m_column_permutation;
bool is_active() const { return m_parent != nullptr; }
square_dense_submatrix() {}
square_dense_submatrix (square_sparse_matrix<T, X> *parent_matrix, unsigned index_start);
void init(square_sparse_matrix<T, X> *parent_matrix, unsigned index_start);
bool is_dense() const override { return true; }
ref operator[] (unsigned i) {
lp_assert(i >= m_index_start);
lp_assert(i < m_parent->dimension());
return ref(i, *this);
}
int find_pivot_column_in_row(unsigned i) const;
void swap_columns(unsigned i, unsigned j) {
if (i != j)
m_column_permutation.transpose_from_left(i, j);
}
unsigned adjust_column(unsigned col) const{
if (col >= m_column_permutation.size())
return col;
return m_column_permutation.apply_reverse(col);
}
unsigned adjust_column_inverse(unsigned col) const{
if (col >= m_column_permutation.size())
return col;
return m_column_permutation[col];
}
unsigned adjust_row(unsigned row) const{
if (row >= m_row_permutation.size())
return row;
return m_row_permutation[row];
}
unsigned adjust_row_inverse(unsigned row) const{
if (row >= m_row_permutation.size())
return row;
return m_row_permutation.apply_reverse(row);
}
void pivot(unsigned i, lp_settings & settings);
void pivot_row_to_row(unsigned i, unsigned row, lp_settings & settings);;
void divide_row_by_pivot(unsigned i);
void update_parent_matrix(lp_settings & settings);
void update_existing_or_delete_in_parent_matrix_for_row(unsigned i, lp_settings & settings);
void push_new_elements_to_parent_matrix(lp_settings & settings);
template <typename L>
L row_by_vector_product(unsigned i, const vector<L> & v);
template <typename L>
L column_by_vector_product(unsigned j, const vector<L> & v);
template <typename L>
L row_by_indexed_vector_product(unsigned i, const indexed_vector<L> & v);
template <typename L>
void apply_from_left_local(indexed_vector<L> & w, lp_settings & settings);
template <typename L>
void apply_from_left_to_vector(vector<L> & w);
bool is_L_matrix() const;
void apply_from_left_to_T(indexed_vector<T> & w, lp_settings & settings) override {
apply_from_left_local(w, settings);
}
void apply_from_right(indexed_vector<T> & w) override {
#if 1==0
indexed_vector<T> wcopy = w;
apply_from_right(wcopy.m_data);
wcopy.m_index.clear();
if (numeric_traits<T>::precise()) {
for (unsigned i = 0; i < m_parent->dimension(); i++) {
if (!is_zero(wcopy.m_data[i]))
wcopy.m_index.push_back(i);
}
} else {
for (unsigned i = 0; i < m_parent->dimension(); i++) {
T & v = wcopy.m_data[i];
if (!lp_settings::is_eps_small_general(v, 1e-14)){
wcopy.m_index.push_back(i);
} else {
v = zero_of_type<T>();
}
}
}
lp_assert(wcopy.is_OK());
apply_from_right(w.m_data);
w.m_index.clear();
if (numeric_traits<T>::precise()) {
for (unsigned i = 0; i < m_parent->dimension(); i++) {
if (!is_zero(w.m_data[i]))
w.m_index.push_back(i);
}
} else {
for (unsigned i = 0; i < m_parent->dimension(); i++) {
T & v = w.m_data[i];
if (!lp_settings::is_eps_small_general(v, 1e-14)){
w.m_index.push_back(i);
} else {
v = zero_of_type<T>();
}
}
}
#else
lp_assert(w.is_OK());
lp_assert(m_work_vector.is_OK());
m_work_vector.resize(w.data_size());
m_work_vector.clear();
lp_assert(m_work_vector.is_OK());
unsigned end = m_index_start + m_dim;
for (unsigned k : w.m_index) {
// find j such that k = adjust_row_inverse(j)
unsigned j = adjust_row(k);
if (j < m_index_start || j >= end) {
m_work_vector.set_value(w[k], adjust_column_inverse(j));
} else { // j >= m_index_start and j < end
unsigned offset = (j - m_index_start) * m_dim; // this is the row start
const T& wv = w[k];
for (unsigned col = m_index_start; col < end; col++, offset ++) {
unsigned adj_col = adjust_column_inverse(col);
m_work_vector.add_value_at_index(adj_col, m_v[offset] * wv);
}
}
}
m_work_vector.clean_up();
lp_assert(m_work_vector.is_OK());
w = m_work_vector;
#endif
}
void apply_from_left(vector<X> & w, lp_settings & /*settings*/) override {
apply_from_left_to_vector(w);// , settings);
}
void apply_from_right(vector<T> & w) override;
#ifdef Z3DEBUG
T get_elem (unsigned i, unsigned j) const override;
unsigned row_count() const override { return m_parent->row_count();}
unsigned column_count() const override { return row_count();}
void set_number_of_rows(unsigned) override {}
void set_number_of_columns(unsigned) override {}
#endif
void conjugate_by_permutation(permutation_matrix<T, X> & q);
};
}

370
src/math/lp/square_dense_submatrix_def.h
View file

@ -1,370 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "util/vector.h"
#include "math/lp/square_dense_submatrix.h"
namespace lp {
template <typename T, typename X>
square_dense_submatrix<T, X>::square_dense_submatrix (square_sparse_matrix<T, X> *parent_matrix, unsigned index_start) :
m_index_start(index_start),
m_dim(parent_matrix->dimension() - index_start),
m_v(m_dim * m_dim),
m_parent(parent_matrix),
m_row_permutation(m_parent->dimension()),
m_column_permutation(m_parent->dimension()) {
int row_offset = - static_cast<int>(m_index_start);
for (unsigned i = index_start; i < parent_matrix->dimension(); i++) {
unsigned row = parent_matrix->adjust_row(i);
for (auto & iv : parent_matrix->get_row_values(row)) {
unsigned j = parent_matrix->adjust_column_inverse(iv.m_index);
lp_assert(j>= m_index_start);
m_v[row_offset + j] = iv.m_value;
}
row_offset += m_dim;
}
}
template <typename T, typename X> void square_dense_submatrix<T, X>::init(square_sparse_matrix<T, X> *parent_matrix, unsigned index_start) {
m_index_start = index_start;
m_dim = parent_matrix->dimension() - index_start;
m_v.resize(m_dim * m_dim);
m_parent = parent_matrix;
m_column_permutation.init(m_parent->dimension());
for (unsigned i = index_start; i < parent_matrix->dimension(); i++) {
unsigned row = parent_matrix->adjust_row(i);
for (auto & iv : parent_matrix->get_row_values(row)) {
unsigned j = parent_matrix->adjust_column_inverse(iv.m_index);
(*this)[i][j] = iv.m_value;
}
}
}
template <typename T, typename X> int square_dense_submatrix<T, X>::find_pivot_column_in_row(unsigned i) const {
int j = -1;
T max = zero_of_type<T>();
lp_assert(i >= m_index_start);
unsigned row_start = (i - m_index_start) * m_dim;
for (unsigned k = i; k < m_parent->dimension(); k++) {
unsigned col = adjust_column(k); // this is where the column is in the row
unsigned offs = row_start + col - m_index_start;
T t = abs(m_v[offs]);
if (t > max) {
j = k;
max = t;
}
}
return j;
}
template <typename T, typename X> void square_dense_submatrix<T, X>::pivot(unsigned i, lp_settings & settings) {
divide_row_by_pivot(i);
for (unsigned k = i + 1; k < m_parent->dimension(); k++)
pivot_row_to_row(i, k, settings);
}
template <typename T, typename X> void square_dense_submatrix<T, X>::pivot_row_to_row(unsigned i, unsigned row, lp_settings & settings) {
lp_assert(i < row);
unsigned pj = adjust_column(i); // the pivot column
unsigned pjd = pj - m_index_start;
unsigned pivot_row_offset = (i-m_index_start)*m_dim;
T pivot = m_v[pivot_row_offset + pjd];
unsigned row_offset= (row-m_index_start)*m_dim;
T m = m_v[row_offset + pjd];
lp_assert(!is_zero(pivot));
m_v[row_offset + pjd] = -m * pivot; // creating L matrix
for (unsigned j = m_index_start; j < m_parent->dimension(); j++) {
if (j == pj) {
pivot_row_offset++;
row_offset++;
continue;
}
auto t = m_v[row_offset] - m_v[pivot_row_offset] * m;
if (settings.abs_val_is_smaller_than_drop_tolerance(t)) {
m_v[row_offset] = zero_of_type<T>();
} else {
m_v[row_offset] = t;
}
row_offset++; pivot_row_offset++;
// at the same time we pivot the L too
}
}
template <typename T, typename X> void square_dense_submatrix<T, X>::divide_row_by_pivot(unsigned i) {
unsigned pj = adjust_column(i); // the pivot column
unsigned irow_offset = (i - m_index_start) * m_dim;
T pivot = m_v[irow_offset + pj - m_index_start];
lp_assert(!is_zero(pivot));
for (unsigned k = m_index_start; k < m_parent->dimension(); k++) {
if (k == pj){
m_v[irow_offset++] = one_of_type<T>() / pivot; // creating the L matrix diagonal
continue;
}
m_v[irow_offset++] /= pivot;
}
}
template <typename T, typename X> void square_dense_submatrix<T, X>::update_parent_matrix(lp_settings & settings) {
for (unsigned i = m_index_start; i < m_parent->dimension(); i++)
update_existing_or_delete_in_parent_matrix_for_row(i, settings);
push_new_elements_to_parent_matrix(settings);
for (unsigned i = m_index_start; i < m_parent->dimension(); i++)
m_parent->set_max_in_row(m_parent->adjust_row(i));
}
template <typename T, typename X> void square_dense_submatrix<T, X>::update_existing_or_delete_in_parent_matrix_for_row(unsigned i, lp_settings & settings) {
bool diag_updated = false;
unsigned ai = m_parent->adjust_row(i);
auto & row_vals = m_parent->get_row_values(ai);
for (unsigned k = 0; k < row_vals.size(); k++) {
auto & iv = row_vals[k];
unsigned j = m_parent->adjust_column_inverse(iv.m_index);
if (j < i) {
m_parent->remove_element(row_vals, iv);
k--;
} else if (i == j) {
m_parent->m_columns[iv.m_index].m_values[iv.m_other].set_value(iv.m_value = one_of_type<T>());
diag_updated = true;
} else { // j > i
T & v = (*this)[i][j];
if (settings.abs_val_is_smaller_than_drop_tolerance(v)) {
m_parent->remove_element(row_vals, iv);
k--;
} else {
m_parent->m_columns[iv.m_index].m_values[iv.m_other].set_value(iv.m_value = v);
v = zero_of_type<T>(); // only new elements are left above the diagonal
}
}
}
if (!diag_updated) {
unsigned aj = m_parent->adjust_column(i);
m_parent->add_new_element(ai, aj, one_of_type<T>());
}
}
template <typename T, typename X> void square_dense_submatrix<T, X>::push_new_elements_to_parent_matrix(lp_settings & settings) {
for (unsigned i = m_index_start; i < m_parent->dimension() - 1; i++) {
unsigned ai = m_parent->adjust_row(i);
for (unsigned j = i + 1; j < m_parent->dimension(); j++) {
T & v = (*this)[i][j];
if (!settings.abs_val_is_smaller_than_drop_tolerance(v)) {
unsigned aj = m_parent->adjust_column(j);
m_parent->add_new_element(ai, aj, v);
}
v = zero_of_type<T>(); // leave only L elements now
}
}
}
template <typename T, typename X>
template <typename L>
L square_dense_submatrix<T, X>::row_by_vector_product(unsigned i, const vector<L> & v) {
lp_assert(i >= m_index_start);
unsigned row_in_subm = i - m_index_start;
unsigned row_offset = row_in_subm * m_dim;
L r = zero_of_type<L>();
for (unsigned j = 0; j < m_dim; j++)
r += m_v[row_offset + j] * v[adjust_column_inverse(m_index_start + j)];
return r;
}
template <typename T, typename X>
template <typename L>
L square_dense_submatrix<T, X>::column_by_vector_product(unsigned j, const vector<L> & v) {
lp_assert(j >= m_index_start);
unsigned offset = j - m_index_start;
L r = zero_of_type<L>();
for (unsigned i = 0; i < m_dim; i++, offset += m_dim)
r += m_v[offset] * v[adjust_row_inverse(m_index_start + i)];
return r;
}
template <typename T, typename X>
template <typename L>
L square_dense_submatrix<T, X>::row_by_indexed_vector_product(unsigned i, const indexed_vector<L> & v) {
lp_assert(i >= m_index_start);
unsigned row_in_subm = i - m_index_start;
unsigned row_offset = row_in_subm * m_dim;
L r = zero_of_type<L>();
for (unsigned j = 0; j < m_dim; j++)
r += m_v[row_offset + j] * v[adjust_column_inverse(m_index_start + j)];
return r;
}
template <typename T, typename X>
template <typename L>
void square_dense_submatrix<T, X>::apply_from_left_local(indexed_vector<L> & w, lp_settings & settings) {
#ifdef Z3DEBUG
// dense_matrix<T, X> deb(*this);
// vector<L> deb_w(w.m_data.size());
// for (unsigned i = 0; i < w.m_data.size(); i++)
// deb_w[i] = w[i];
// deb.apply_from_left(deb_w);
#endif // use indexed vector here
#ifndef DO_NOT_USE_INDEX
vector<L> t(m_parent->dimension(), zero_of_type<L>());
for (auto k : w.m_index) {
unsigned j = adjust_column(k); // k-th element will contribute only to column j
if (j < m_index_start || j >= this->m_index_start + this->m_dim) { // it is a unit matrix outside
t[adjust_row_inverse(j)] = w[k];
} else {
const L & v = w[k];
for (unsigned i = 0; i < m_dim; i++) {
unsigned row = adjust_row_inverse(m_index_start + i);
unsigned offs = i * m_dim + j - m_index_start;
t[row] += m_v[offs] * v;
}
}
}
w.m_index.clear();
for (unsigned i = 0; i < m_parent->dimension(); i++) {
const L & v = t[i];
if (!settings.abs_val_is_smaller_than_drop_tolerance(v)){
w.m_index.push_back(i);
w.m_data[i] = v;
} else {
w.m_data[i] = zero_of_type<L>();
}
}
#else
vector<L> t(m_parent->dimension());
for (unsigned i = 0; i < m_index_start; i++) {
t[adjust_row_inverse(i)] = w[adjust_column_inverse(i)];
}
for (unsigned i = m_index_start; i < m_parent->dimension(); i++){
t[adjust_row_inverse(i)] = row_by_indexed_vector_product(i, w);
}
for (unsigned i = 0; i < m_parent->dimension(); i++) {
w.set_value(t[i], i);
}
for (unsigned i = 0; i < m_parent->dimension(); i++) {
const L & v = t[i];
if (!is_zero(v))
w.m_index.push_back(i);
w.m_data[i] = v;
}
#endif
#ifdef Z3DEBUG
// cout << "w final" << endl;
// print_vector(w.m_data);
// lp_assert(vectors_are_equal<T>(deb_w, w.m_data));
// lp_assert(w.is_OK());
#endif
}
template <typename T, typename X>
template <typename L>
void square_dense_submatrix<T, X>::apply_from_left_to_vector(vector<L> & w) {
// lp_settings & settings) {
// dense_matrix<T, L> deb(*this);
// vector<L> deb_w(w);
// deb.apply_from_left_to_X(deb_w, settings);
// // cout << "deb" << endl;
// // print_matrix(deb);
// // cout << "w" << endl;
// // print_vector(w.m_data);
// // cout << "deb_w" << endl;
// // print_vector(deb_w);
vector<L> t(m_parent->dimension());
for (unsigned i = 0; i < m_index_start; i++) {
t[adjust_row_inverse(i)] = w[adjust_column_inverse(i)];
}
for (unsigned i = m_index_start; i < m_parent->dimension(); i++){
t[adjust_row_inverse(i)] = row_by_vector_product(i, w);
}
for (unsigned i = 0; i < m_parent->dimension(); i++) {
w[i] = t[i];
}
#ifdef Z3DEBUG
// cout << "w final" << endl;
// print_vector(w.m_data);
// lp_assert(vectors_are_equal<L>(deb_w, w));
#endif
}
template <typename T, typename X> bool square_dense_submatrix<T, X>::is_L_matrix() const {
#ifdef Z3DEBUG
lp_assert(m_row_permutation.is_identity());
for (unsigned i = 0; i < m_parent->dimension(); i++) {
if (i < m_index_start) {
lp_assert(m_column_permutation[i] == i);
continue;
}
unsigned row_offs = (i-m_index_start)*m_dim;
for (unsigned k = 0; k < m_dim; k++) {
unsigned j = m_index_start + k;
unsigned jex = adjust_column_inverse(j);
if (jex > i) {
lp_assert(is_zero(m_v[row_offs + k]));
} else if (jex == i) {
lp_assert(!is_zero(m_v[row_offs + k]));
}
}
}
#endif
return true;
}
template <typename T, typename X> void square_dense_submatrix<T, X>::apply_from_right(vector<T> & w) {
#ifdef Z3DEBUG
// dense_matrix<T, X> deb(*this);
// vector<T> deb_w(w);
// deb.apply_from_right(deb_w);
#endif
vector<T> t(w.size());
for (unsigned j = 0; j < m_index_start; j++) {
t[adjust_column_inverse(j)] = w[adjust_row_inverse(j)];
}
unsigned end = m_index_start + m_dim;
for (unsigned j = end; j < m_parent->dimension(); j++) {
t[adjust_column_inverse(j)] = w[adjust_row_inverse(j)];
}
for (unsigned j = m_index_start; j < end; j++) {
t[adjust_column_inverse(j)] = column_by_vector_product(j, w);
}
w = t;
#ifdef Z3DEBUG
// lp_assert(vector_are_equal<T>(deb_w, w));
#endif
}
#ifdef Z3DEBUG
template <typename T, typename X> T square_dense_submatrix<T, X>::get_elem (unsigned i, unsigned j) const {
i = adjust_row(i);
j = adjust_column(j);
if (i < m_index_start || j < m_index_start)
return i == j? one_of_type<T>() : zero_of_type<T>();
unsigned offs = (i - m_index_start)* m_dim + j - m_index_start;
return m_v[offs];
}
#endif
template <typename T, typename X> void square_dense_submatrix<T, X>::conjugate_by_permutation(permutation_matrix<T, X> & q) {
m_row_permutation.multiply_by_permutation_from_left(q);
m_column_permutation.multiply_by_reverse_from_right(q);
}
}

592
src/test/lp/lp.cpp
View file

@ -48,7 +48,6 @@
#include "test/lp/gomory_test.h"
#include "math/lp/matrix.h"
#include "math/lp/hnf.h"
#include "math/lp/square_sparse_matrix_def.h"
#include "math/lp/general_matrix.h"
#include "math/lp/lp_bound_propagator.h"
#include "math/lp/nla_solver.h"
@ -56,6 +55,8 @@
#include "math/lp/cross_nested.h"
#include "math/lp/int_cube.h"
#include "math/lp/emonics.h"
#include "math/lp/static_matrix.h"
#include "math/lp/dense_matrix.h"
bool my_white_space(const char & a) {
return a == ' ' || a == '\t';
@ -388,122 +389,6 @@ struct simple_column_namer:public column_namer
};
template <typename T, typename X>
void test_matrix(square_sparse_matrix<T, X> & a) {
auto m = a.dimension();
// copy a to b in the reversed order
square_sparse_matrix<T, X> 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 --) {
b(row, col) = (T const&) a(row, col);
}
std::cout << "zeroing b in the reverse order"<< std::endl;
for (int row = m - 1; row >= 0; row--)
for (int col = m - 1; col >= 0; col --)
b.set(row, col, T(0));
for (unsigned row = 0; row < m; row ++)
for (unsigned col = 0; col < m; col ++)
a.set(row, col, T(0));
unsigned i = my_random() % m;
unsigned j = my_random() % m;
auto t = T(1);
a.set(i, j, t);
lp_assert(a.get(i, j) == t);
unsigned j1;
if (j < m - 1) {
j1 = m - 1;
a.set(i, j1, T(2));
}
}
void tst1() {
std::cout << "testing the minimal matrix with 1 row and 1 column" << std::endl;
square_sparse_matrix<double, double> m0(1, 1);
m0.set(0, 0, 1);
// print_matrix(m0);
m0.set(0, 0, 0);
// print_matrix(m0);
test_matrix(m0);
unsigned rows = 2;
square_sparse_matrix<double, double> m(rows, rows);
std::cout << "setting m(0,1)=" << std::endl;
m.set(0, 1, 11);
m.set(0, 0, 12);
// print_matrix<double>(m);
test_matrix(m);
square_sparse_matrix<double, double> m1(2, 2);
m1.set(0, 0, 2);
m1.set(1, 0, 3);
// print_matrix(m1);
std::cout << " zeroing matrix 2 by 2" << std::endl;
m1.set(0, 0, 0);
m1.set(1, 0, 0);
// print_matrix(m1);
test_matrix(m1);
std::cout << "printing zero matrix 3 by 1" << std::endl;
square_sparse_matrix<double, double> m2(3, 3);
// print_matrix(m2);
m2.set(0, 0, 1);
m2.set(2, 0, 2);
std::cout << "printing matrix 3 by 1 with a gap" << std::endl;
// print_matrix(m2);
test_matrix(m2);
square_sparse_matrix<double, double> m10by9(10, 10);
m10by9.set(0, 1, 1);
m10by9(0, 1) = 4;
double test = m10by9(0, 1);
std::cout << "got " << test << std::endl;
m10by9.set(0, 8, 8);
m10by9.set(3, 4, 7);
m10by9.set(3, 2, 5);
m10by9.set(3, 8, 99);
m10by9.set(3, 2, 6);
m10by9.set(1, 8, 9);
m10by9.set(4, 0, 40);
m10by9.set(0, 0, 10);
std::cout << "printing matrix 10 by 9" << std::endl;
// print_matrix(m10by9);
test_matrix(m10by9);
std::cout <<"zeroing m10by9\n";
#ifdef Z3DEBUG
for (unsigned int i = 0; i < m10by9.dimension(); i++)
for (unsigned int j = 0; j < m10by9.column_count(); j++)
m10by9.set(i, j, 0);
#endif
// print_matrix(m10by9);
}
vector<int> allocate_basis_heading(unsigned count) { // the rest of initialization will be handled by lu_QR
vector<int> basis_heading(count, -1);
@ -555,12 +440,6 @@ void test_small_lu(lp_settings & settings) {
#endif
void fill_long_row(square_sparse_matrix<double, double> &m, int i) {
int n = m.dimension();
for (int j = 0; j < n; j ++) {
m (i, (j + i) % n) = j * j;
}
}
void fill_long_row(static_matrix<double, double> &m, int i) {
int n = m.column_count();
@ -570,13 +449,6 @@ void fill_long_row(static_matrix<double, double> &m, int i) {
}
void fill_long_row_exp(square_sparse_matrix<double, double> &m, int i) {
int n = m.dimension();
for (int j = 0; j < n; j ++) {
m(i, j) = my_random() % 20;
}
}
void fill_long_row_exp(static_matrix<double, double> &m, int i) {
int n = m.column_count();
@ -586,26 +458,9 @@ void fill_long_row_exp(static_matrix<double, double> &m, int i) {
}
}
void fill_larger_square_sparse_matrix_exp(square_sparse_matrix<double, double> & m){
for ( unsigned i = 0; i < m.dimension(); i++ )
fill_long_row_exp(m, i);
}
void fill_larger_square_sparse_matrix_exp(static_matrix<double, double> & m){
for ( unsigned i = 0; i < m.row_count(); i++ )
fill_long_row_exp(m, i);
}
void fill_larger_square_sparse_matrix(square_sparse_matrix<double, double> & m){
for ( unsigned i = 0; i < m.dimension(); i++ )
fill_long_row(m, i);
}
void fill_larger_square_sparse_matrix(static_matrix<double, double> & m){
for ( unsigned i = 0; i < m.row_count(); i++ )
fill_long_row(m, i);
}
int perm_id = 0;
@ -616,11 +471,6 @@ int perm_id = 0;
void init_b(vector<double> & b, square_sparse_matrix<double, double> & m, vector<double>& x) {
for (unsigned i = 0; i < m.dimension(); i++) {
b.push_back(m.dot_product_with_row(i, x));
}
}
void init_b(vector<double> & b, static_matrix<double, double> & m, vector<double> & x) {
for (unsigned i = 0; i < m.row_count(); i++) {
@ -671,7 +521,7 @@ void test_lp_1() {
m(2, 0) = 2; m(2, 1) = -1; m(2, 2) = 2; m(2, 5) = 1;
m(3, 0) = 2; m(3, 1) = 3; m(3, 2) = -1; m(3, 6) = 1;
#ifdef Z3DEBUG
print_matrix(m, std::cout);
//print_matrix(m, std::cout);
#endif
vector<double> x_star(7);
x_star[0] = 0; x_star[1] = 0; x_star[2] = 0;
@ -724,213 +574,9 @@ void test_lp_primal_core_solver() {
}
#ifdef Z3DEBUG
template <typename T, typename X>
void test_swap_rows_with_permutation(square_sparse_matrix<T, X>& m){
std::cout << "testing swaps" << std::endl;
unsigned dim = m.row_count();
dense_matrix<double, double> original(&m);
permutation_matrix<double, double> q(dim);
print_matrix(m, std::cout);
lp_assert(original == q * m);
for (int i = 0; i < 100; i++) {
unsigned row1 = my_random() % dim;
unsigned row2 = my_random() % dim;
if (row1 == row2) continue;
std::cout << "swap " << row1 << " " << row2 << std::endl;
m.swap_rows(row1, row2);
q.transpose_from_left(row1, row2);
lp_assert(original == q * m);
print_matrix(m, std::cout);
std::cout << std::endl;
}
}
#endif
template <typename T, typename X>
void fill_matrix(square_sparse_matrix<T, X>& m); // forward definition
#ifdef Z3DEBUG
template <typename T, typename X>
void test_swap_cols_with_permutation(square_sparse_matrix<T, X>& m){
std::cout << "testing swaps" << std::endl;
unsigned dim = m.row_count();
dense_matrix<double, double> original(&m);
permutation_matrix<double, double> q(dim);
print_matrix(m, std::cout);
lp_assert(original == q * m);
for (int i = 0; i < 100; i++) {
unsigned row1 = my_random() % dim;
unsigned row2 = my_random() % dim;
if (row1 == row2) continue;
std::cout << "swap " << row1 << " " << row2 << std::endl;
m.swap_rows(row1, row2);
q.transpose_from_right(row1, row2);
lp_assert(original == q * m);
print_matrix(m, std::cout);
std::cout << std::endl;
}
}
template <typename T, typename X>
void test_swap_rows(square_sparse_matrix<T, X>& m, unsigned i0, unsigned i1){
std::cout << "test_swap_rows(" << i0 << "," << i1 << ")" << std::endl;
square_sparse_matrix<T, X> 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);
}
std::cout << "swapping rows "<< i0 << "," << i1 << std::endl;
m.swap_rows(i0, i1);
for (unsigned j = 0; j < m.dimension(); j++) {
lp_assert(mcopy(i0, j) == m(i1, j));
lp_assert(mcopy(i1, j) == m(i0, j));
}
}
template <typename T, typename X>
void test_swap_columns(square_sparse_matrix<T, X>& m, unsigned i0, unsigned i1){
std::cout << "test_swap_columns(" << i0 << "," << i1 << ")" << std::endl;
square_sparse_matrix<T, X> 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);
}
m.swap_columns(i0, i1);
for (unsigned j = 0; j < m.dimension(); j++) {
lp_assert(mcopy(j, i0) == m(j, i1));
lp_assert(mcopy(j, i1) == m(j, i0));
}
for (unsigned i = 0; i < m.dimension(); i++) {
if (i == i0 || i == i1)
continue;
for (unsigned j = 0; j < m.dimension(); j++) {
lp_assert(mcopy(j, i)== m(j, i));
}
}
}
#endif
template <typename T, typename X>
void fill_matrix(square_sparse_matrix<T, X>& m){
int v = 0;
for (int i = m.dimension() - 1; i >= 0; i--) {
for (int j = m.dimension() - 1; j >=0; j--){
m(i, j) = v++;
}
}
}
void test_pivot_like_swaps_and_pivot(){
square_sparse_matrix<double, double> m(10, 10);
fill_matrix(m);
// print_matrix(m);
// pivot at 2,7
m.swap_columns(0, 7);
// print_matrix(m);
m.swap_rows(2, 0);
// print_matrix(m);
for (unsigned i = 1; i < m.dimension(); i++) {
m(i, 0) = 0;
}
// print_matrix(m);
// say pivot at 3,4
m.swap_columns(1, 4);
// print_matrix(m);
m.swap_rows(1, 3);
// print_matrix(m);
vector<double> row;
double alpha = 2.33;
unsigned pivot_row = 1;
unsigned target_row = 2;
unsigned pivot_row_0 = 3;
double beta = 3.1;
m(target_row, 3) = 0;
m(target_row, 5) = 0;
m(pivot_row, 6) = 0;
#ifdef Z3DEBUG
print_matrix(m, std::cout);
#endif
for (unsigned j = 0; j < m.dimension(); j++) {
row.push_back(m(target_row, j) + alpha * m(pivot_row, j) + beta * m(pivot_row_0, j));
}
for (auto & t : row) {
std::cout << t << ",";
}
std::cout << std::endl;
lp_settings settings;
m.pivot_row_to_row(pivot_row, alpha, target_row, settings);
m.pivot_row_to_row(pivot_row_0, beta, target_row, settings);
// print_matrix(m);
for (unsigned j = 0; j < m.dimension(); j++) {
lp_assert(abs(row[j] - m(target_row, j)) < 0.00000001);
}
}
#ifdef Z3DEBUG
void test_swap_rows() {
square_sparse_matrix<double, double> m(10, 10);
fill_matrix(m);
// print_matrix(m);
test_swap_rows(m, 3, 5);
test_swap_rows(m, 1, 3);
test_swap_rows(m, 1, 3);
test_swap_rows(m, 1, 7);
test_swap_rows(m, 3, 7);
test_swap_rows(m, 0, 7);
m(0, 4) = 1;
// print_matrix(m);
test_swap_rows(m, 0, 7);
// go over some corner cases
square_sparse_matrix<double, double> m0(2, 2);
test_swap_rows(m0, 0, 1);
m0(0, 0) = 3;
test_swap_rows(m0, 0, 1);
m0(1, 0) = 3;
test_swap_rows(m0, 0, 1);
square_sparse_matrix<double, double> m1(10, 10);
test_swap_rows(m1, 0, 1);
m1(0, 0) = 3;
test_swap_rows(m1, 0, 1);
m1(1, 0) = 3;
m1(0, 3) = 5;
m1(1, 3) = 4;
m1(1, 8) = 8;
m1(1, 9) = 8;
test_swap_rows(m1, 0, 1);
square_sparse_matrix<double, double> m2(3, 3);
test_swap_rows(m2, 0, 1);
m2(0, 0) = 3;
test_swap_rows(m2, 0, 1);
m2(2, 0) = 3;
test_swap_rows(m2, 0, 2);
}
void fill_uniformly(square_sparse_matrix<double, double> & m, unsigned dim) {
int v = 0;
for (unsigned i = 0; i < dim; i++) {
for (unsigned j = 0; j < dim; j++) {
m(i, j) = v++;
}
}
}
void fill_uniformly(dense_matrix<double, double> & m, unsigned dim) {
int v = 0;
@ -941,146 +587,8 @@ void fill_uniformly(dense_matrix<double, double> & m, unsigned dim) {
}
}
void square_sparse_matrix_with_permutations_test() {
unsigned dim = 4;
square_sparse_matrix<double, double> m(dim, dim);
fill_uniformly(m, dim);
dense_matrix<double, double> dm(dim, dim);
fill_uniformly(dm, dim);
dense_matrix<double, double> dm0(dim, dim);
fill_uniformly(dm0, dim);
permutation_matrix<double, double> q0(dim);
q0[0] = 1;
q0[1] = 0;
q0[2] = 3;
q0[3] = 2;
permutation_matrix<double, double> q1(dim);
q1[0] = 1;
q1[1] = 2;
q1[2] = 3;
q1[3] = 0;
permutation_matrix<double, double> p0(dim);
p0[0] = 1;
p0[1] = 0;
p0[2] = 3;
p0[3] = 2;
permutation_matrix<double, double> p1(dim);
p1[0] = 1;
p1[1] = 2;
p1[2] = 3;
p1[3] = 0;
m.multiply_from_left(q0);
for (unsigned i = 0; i < dim; i++) {
for (unsigned j = 0; j < dim; j++) {
lp_assert(m(i, j) == dm0.get_elem(q0[i], j));
}
}
auto q0_dm = q0 * dm;
lp_assert(m == q0_dm);
m.multiply_from_left(q1);
for (unsigned i = 0; i < dim; i++) {
for (unsigned j = 0; j < dim; j++) {
lp_assert(m(i, j) == dm0.get_elem(q0[q1[i]], j));
}
}
auto q1_q0_dm = q1 * q0_dm;
lp_assert(m == q1_q0_dm);
m.multiply_from_right(p0);
for (unsigned i = 0; i < dim; i++) {
for (unsigned j = 0; j < dim; j++) {
lp_assert(m(i, j) == dm0.get_elem(q0[q1[i]], p0[j]));
}
}
auto q1_q0_dm_p0 = q1_q0_dm * p0;
lp_assert(m == q1_q0_dm_p0);
m.multiply_from_right(p1);
for (unsigned i = 0; i < dim; i++) {
for (unsigned j = 0; j < dim; j++) {
lp_assert(m(i, j) == dm0.get_elem(q0[q1[i]], p1[p0[j]]));
}
}
auto q1_q0_dm_p0_p1 = q1_q0_dm_p0 * p1;
lp_assert(m == q1_q0_dm_p0_p1);
m.multiply_from_right(p1);
for (unsigned i = 0; i < dim; i++) {
for (unsigned j = 0; j < dim; j++) {
lp_assert(m(i, j) == dm0.get_elem(q0[q1[i]], p1[p1[p0[j]]]));
}
}
auto q1_q0_dm_p0_p1_p1 = q1_q0_dm_p0_p1 * p1;
lp_assert(m == q1_q0_dm_p0_p1_p1);
}
void test_swap_columns() {
square_sparse_matrix<double, double> m(10, 10);
fill_matrix(m);
// print_matrix(m);
test_swap_columns(m, 3, 5);
test_swap_columns(m, 1, 3);
test_swap_columns(m, 1, 3);
// print_matrix(m);
test_swap_columns(m, 1, 7);
test_swap_columns(m, 3, 7);
test_swap_columns(m, 0, 7);
test_swap_columns(m, 0, 7);
// go over some corner cases
square_sparse_matrix<double, double> m0(2, 2);
test_swap_columns(m0, 0, 1);
m0(0, 0) = 3;
test_swap_columns(m0, 0, 1);
m0(0, 1) = 3;
test_swap_columns(m0, 0, 1);
square_sparse_matrix<double, double> m1(10, 10);
test_swap_columns(m1, 0, 1);
m1(0, 0) = 3;
test_swap_columns(m1, 0, 1);
m1(0, 1) = 3;
m1(3, 0) = 5;
m1(3, 1) = 4;
m1(8, 1) = 8;
m1(9, 1) = 8;
test_swap_columns(m1, 0, 1);
square_sparse_matrix<double, double> m2(3, 3);
test_swap_columns(m2, 0, 1);
m2(0, 0) = 3;
test_swap_columns(m2, 0, 1);
m2(0, 2) = 3;
test_swap_columns(m2, 0, 2);
}
void test_swap_operations() {
test_swap_rows();
test_swap_columns();
}
void test_dense_matrix() {
dense_matrix<double, double> d(3, 2);
@ -1110,7 +618,7 @@ void test_dense_matrix() {
p2.set_elem(1, 0, 1);
auto c2 = d * p2;
}
#endif
@ -1471,59 +979,6 @@ bool contains(std::string const & s, char const * pattern) {
void test_init_U() {
static_matrix<double, double> m(3, 7);
m(0, 0) = 10; m(0, 1) = 11; m(0, 2) = 12; m(0, 3) = 13; m(0, 4) = 14;
m(1, 0) = 20; m(1, 1) = 21; m(1, 2) = 22; m(1, 3) = 23; m(1, 5) = 24;
m(2, 0) = 30; m(2, 1) = 31; m(2, 2) = 32; m(2, 3) = 33; m(2, 6) = 34;
#ifdef Z3DEBUG
print_matrix(m, std::cout);
#endif
vector<unsigned> basis(3);
basis[0] = 1;
basis[1] = 2;
basis[2] = 4;
square_sparse_matrix<double, double> u(m, basis);
for (unsigned i = 0; i < 3; i++) {
for (unsigned j = 0; j < 3; j ++) {
lp_assert(m(i, basis[j]) == u(i, j));
}
}
// print_matrix(m);
// print_matrix(u);
}
void test_replace_column() {
square_sparse_matrix<double, double> m(10, 10);
fill_matrix(m);
m.swap_columns(0, 7);
m.swap_columns(6, 3);
m.swap_rows(2, 0);
for (unsigned i = 1; i < m.dimension(); i++) {
m(i, 0) = 0;
}
indexed_vector<double> w(m.dimension());
for (unsigned i = 0; i < m.dimension(); i++) {
w.set_value(i % 3, i);
}
lp_settings settings;
for (unsigned column_to_replace = 0; column_to_replace < m.dimension(); column_to_replace ++) {
m.replace_column(column_to_replace, w, settings);
for (unsigned i = 0; i < m.dimension(); i++) {
lp_assert(abs(w[i] - m(i, column_to_replace)) < 0.00000001);
}
}
}
void setup_args_parser(argument_parser & parser) {
parser.add_option_with_help_string("-monics", "test emonics");
parser.add_option_with_help_string("-nex_order", "test nex order");
@ -1544,8 +999,6 @@ void setup_args_parser(argument_parser & parser) {
parser.add_option_with_help_string("-xyz_sample", "run a small interactive scenario");
parser.add_option_with_after_string_with_help("--density", "the percentage of non-zeroes in the matrix below which it is not dense");
parser.add_option_with_after_string_with_help("--harris_toler", "harris tolerance");
parser.add_option_with_help_string("--test_swaps", "test row swaps with a permutation");
parser.add_option_with_help_string("--test_perm", "test permutations");
parser.add_option_with_after_string_with_help("--checklu", "the file name for lu checking");
parser.add_option_with_after_string_with_help("--partial_pivot", "the partial pivot constant, a number somewhere between 10 and 100");
parser.add_option_with_after_string_with_help("--percent_for_enter", "which percent of columns check for entering column");
@ -1871,38 +1324,6 @@ void check_lu_from_file(std::string lufile_name) {
lp_assert(false);
}
void test_square_dense_submatrix() {
std::cout << "testing square_dense_submatrix" << std::endl;
unsigned parent_dim = 7;
square_sparse_matrix<double, double> parent(parent_dim, 0);
fill_matrix(parent);
unsigned index_start = 3;
square_dense_submatrix<double, double> d;
d.init(&parent, index_start);
for (unsigned i = index_start; i < parent_dim; i++)
for (unsigned j = index_start; j < parent_dim; j++)
d[i][j] = i*3+j*2;
#ifdef Z3DEBUG
unsigned dim = parent_dim - index_start;
dense_matrix<double, double> m(dim, dim);
for (unsigned i = index_start; i < parent_dim; i++)
for (unsigned j = index_start; j < parent_dim; j++)
m[i-index_start][j-index_start] = d[i][j];
print_matrix(&m, std::cout);
#endif
for (unsigned i = index_start; i < parent_dim; i++)
for (unsigned j = index_start; j < parent_dim; j++)
d[i][j] = d[j][i];
#ifdef Z3DEBUG
for (unsigned i = index_start; i < parent_dim; i++)
for (unsigned j = index_start; j < parent_dim; j++)
m[i-index_start][j-index_start] = d[i][j];
print_matrix(&m, std::cout);
std::cout << std::endl;
#endif
}
void print_st(lp_status status) {
@ -2958,8 +2379,3 @@ 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>(general_matrix&, std::ostream&);
}
#endif