michael/lean2
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

chore(lp): improve formatting

Browse source 
Signed-off-by: Lev Nachmanson <levnach@microsoft.com>
This commit is contained in:
Lev Nachmanson 2016-01-21 11:33:09 -08:00 committed by Leonardo de Moura
parent e9cd621855
commit fc858d98c0
15 changed files with 837 additions and 971 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

296
src/util/lp/binary_heap_priority_queue.cpp
View file

@ -6,11 +6,11 @@
*/
#include "util/lp/binary_heap_priority_queue.h"
namespace lean {
// is is the child place in heap
// is is the child place in heap
template <typename T> void binary_heap_priority_queue<T>::swap_with_parent(unsigned i) {
unsigned parent = m_heap[i >> 1];
put_at(i >> 1, m_heap[i]);
put_at(i, parent);
unsigned parent = m_heap[i >> 1];
put_at(i >> 1, m_heap[i]);
put_at(i, parent);
}
template <typename T> void binary_heap_priority_queue<T>::put_at(unsigned i, unsigned h) {
@ -19,183 +19,183 @@ template <typename T> void binary_heap_priority_queue<T>::put_at(unsigned i, uns
}
template <typename T> void binary_heap_priority_queue<T>::decrease_priority(unsigned o, T newPriority) {
m_priorities[o] = newPriority;
int i = m_heap_inverse[o];
while (i > 1) {
if (m_priorities[m_heap[i]] < m_priorities[m_heap[i >> 1]])
swap_with_parent(i);
else
break;
i >>= 1;
}
m_priorities[o] = newPriority;
int i = m_heap_inverse[o];
while (i > 1) {
if (m_priorities[m_heap[i]] < m_priorities[m_heap[i >> 1]])
swap_with_parent(i);
else
break;
i >>= 1;
}
}
#ifdef LEAN_DEBUG
template <typename T> bool binary_heap_priority_queue<T>::is_consistent() const {
for (int i = 0; i < m_heap_inverse.size(); i++) {
int i_index = m_heap_inverse[i];
lean_assert(i_index <= static_cast<int>(m_heap_size));
lean_assert(i_index == -1 || m_heap[i_index] == i);
}
for (unsigned i = 1; i < m_heap_size; i++) {
unsigned ch = i << 1;
for (int k = 0; k < 2; k++) {
if (ch > m_heap_size) break;
if (!(m_priorities[m_heap[i]] <= m_priorities[m_heap[ch]])){
std::cout << "m_heap_size = " << m_heap_size << std::endl;
std::cout << "i = " << i << std::endl;
std::cout << "m_heap[i] = " << m_heap[i] << std::endl;
std::cout << "ch = " << ch << std::endl;
std::cout << "m_heap[ch] = " << m_heap[ch] << std::endl;
std::cout << "m_priorities[m_heap[i]] = " << m_priorities[m_heap[i]] << std::endl;
std::cout << "m_priorities[m_heap[ch]] = " << m_priorities[m_heap[ch]]<< std::endl;
return false;
}
ch++;
}
}
return true;
for (int i = 0; i < m_heap_inverse.size(); i++) {
int i_index = m_heap_inverse[i];
lean_assert(i_index <= static_cast<int>(m_heap_size));
lean_assert(i_index == -1 || m_heap[i_index] == i);
}
for (unsigned i = 1; i < m_heap_size; i++) {
unsigned ch = i << 1;
for (int k = 0; k < 2; k++) {
if (ch > m_heap_size) break;
if (!(m_priorities[m_heap[i]] <= m_priorities[m_heap[ch]])){
std::cout << "m_heap_size = " << m_heap_size << std::endl;
std::cout << "i = " << i << std::endl;
std::cout << "m_heap[i] = " << m_heap[i] << std::endl;
std::cout << "ch = " << ch << std::endl;
std::cout << "m_heap[ch] = " << m_heap[ch] << std::endl;
std::cout << "m_priorities[m_heap[i]] = " << m_priorities[m_heap[i]] << std::endl;
std::cout << "m_priorities[m_heap[ch]] = " << m_priorities[m_heap[ch]]<< std::endl;
return false;
}
ch++;
}
}
return true;
}
#endif
template <typename T> void binary_heap_priority_queue<T>::remove(unsigned o) {
T priority_of_o = m_priorities[o];
int o_in_heap = m_heap_inverse[o];
if (o_in_heap == -1) {
return; // nothing to do
}
lean_assert(o_in_heap <= m_heap_size);
if (o_in_heap < m_heap_size) {
put_at(o_in_heap, m_heap[m_heap_size--]);
if (m_priorities[m_heap[o_in_heap]] > priority_of_o) {
fix_heap_under(o_in_heap);
} else { // we need to propogate the m_heap[o_in_heap] up
unsigned i = o_in_heap;
while (i > 1) {
unsigned ip = i >> 1;
if (m_priorities[m_heap[i]] < m_priorities[m_heap[ip]])
swap_with_parent(i);
else
break;
i = ip;
}
}
} else {
lean_assert(o_in_heap == m_heap_size);
m_heap_size--;
}
m_heap_inverse[o] = -1;
// lean_assert(is_consistent());
T priority_of_o = m_priorities[o];
int o_in_heap = m_heap_inverse[o];
if (o_in_heap == -1) {
return; // nothing to do
}
// n is the initial queue capacity.
// The capacity will be enlarged two times automatically if needed
lean_assert(o_in_heap <= m_heap_size);
if (o_in_heap < m_heap_size) {
put_at(o_in_heap, m_heap[m_heap_size--]);
if (m_priorities[m_heap[o_in_heap]] > priority_of_o) {
fix_heap_under(o_in_heap);
} else { // we need to propogate the m_heap[o_in_heap] up
unsigned i = o_in_heap;
while (i > 1) {
unsigned ip = i >> 1;
if (m_priorities[m_heap[i]] < m_priorities[m_heap[ip]])
swap_with_parent(i);
else
break;
i = ip;
}
}
} else {
lean_assert(o_in_heap == m_heap_size);
m_heap_size--;
}
m_heap_inverse[o] = -1;
// lean_assert(is_consistent());
}
// n is the initial queue capacity.
// The capacity will be enlarged two times automatically if needed
template <typename T> binary_heap_priority_queue<T>::binary_heap_priority_queue(unsigned n) :
m_priorities(n),
m_heap(n + 1), // because the indexing for A starts from 1
m_heap_inverse(n, -1)
{ }
m_priorities(n),
m_heap(n + 1), // because the indexing for A starts from 1
m_heap_inverse(n, -1)
{ }
template <typename T> void binary_heap_priority_queue<T>::resize(unsigned n) {
m_priorities.resize(n);
m_heap.resize(n + 1);
m_heap_inverse.resize(n, -1);
}
m_priorities.resize(n);
m_heap.resize(n + 1);
m_heap_inverse.resize(n, -1);
}
template <typename T> void binary_heap_priority_queue<T>::put_to_heap(unsigned i, unsigned o) {
m_heap[i] = o;
m_heap_inverse[o] = i;
}
m_heap[i] = o;
m_heap_inverse[o] = i;
}
template <typename T> void binary_heap_priority_queue<T>::enqueue_new(unsigned o, const T& priority) {
m_heap_size++;
int i = m_heap_size;
lean_assert(o < m_priorities.size());
m_priorities[o] = priority;
put_at(i, o);
while (i > 1 && m_priorities[m_heap[i >> 1]] > priority) {
swap_with_parent(i);
i >>= 1;
}
m_heap_size++;
int i = m_heap_size;
lean_assert(o < m_priorities.size());
m_priorities[o] = priority;
put_at(i, o);
while (i > 1 && m_priorities[m_heap[i >> 1]] > priority) {
swap_with_parent(i);
i >>= 1;
}
// This method can work with an element that is already in the queue.
// In this case the priority will be changed and the queue adjusted.
}
// This method can work with an element that is already in the queue.
// In this case the priority will be changed and the queue adjusted.
template <typename T> void binary_heap_priority_queue<T>::enqueue(unsigned o, const T & priority) {
if (o >= m_priorities.size()) {
resize(o << 1); // make the size twice larger
}
if (m_heap_inverse[o] == -1)
enqueue_new(o, priority);
else
change_priority_for_existing(o, priority);
if (o >= m_priorities.size()) {
resize(o << 1); // make the size twice larger
}
if (m_heap_inverse[o] == -1)
enqueue_new(o, priority);
else
change_priority_for_existing(o, priority);
}
template <typename T> void binary_heap_priority_queue<T>::change_priority_for_existing(unsigned o, const T & priority) {
if (m_priorities[o] > priority) {
decrease_priority(o, priority);
} else {
m_priorities[o] = priority;
fix_heap_under(m_heap_inverse[o]);
}
if (m_priorities[o] > priority) {
decrease_priority(o, priority);
} else {
m_priorities[o] = priority;
fix_heap_under(m_heap_inverse[o]);
}
}
/// return the first element of the queue and removes it from the queue
/// return the first element of the queue and removes it from the queue
template <typename T> unsigned binary_heap_priority_queue<T>::dequeue_and_get_priority(T & priority) {
lean_assert(m_heap_size != 0);
int ret = m_heap[1];
priority = m_priorities[ret];
put_the_last_at_the_top_and_fix_the_heap();
return ret;
}
lean_assert(m_heap_size != 0);
int ret = m_heap[1];
priority = m_priorities[ret];
put_the_last_at_the_top_and_fix_the_heap();
return ret;
}
template <typename T> void binary_heap_priority_queue<T>::fix_heap_under(unsigned i) {
while (true) {
int smallest = i;
int l = i << 1;
if (l <= m_heap_size && m_priorities[m_heap[l]] < m_priorities[m_heap[i]])
smallest = l;
int r = l + 1;
if (r <= m_heap_size && m_priorities[m_heap[r]] < m_priorities[m_heap[smallest]])
smallest = r;
if (smallest != i)
swap_with_parent(smallest);
else
break;
i = smallest;
}
while (true) {
int smallest = i;
int l = i << 1;
if (l <= m_heap_size && m_priorities[m_heap[l]] < m_priorities[m_heap[i]])
smallest = l;
int r = l + 1;
if (r <= m_heap_size && m_priorities[m_heap[r]] < m_priorities[m_heap[smallest]])
smallest = r;
if (smallest != i)
swap_with_parent(smallest);
else
break;
i = smallest;
}
}
template <typename T> void binary_heap_priority_queue<T>::put_the_last_at_the_top_and_fix_the_heap() {
if (m_heap_size > 1) {
put_at(1, m_heap[m_heap_size--]);
fix_heap_under(1);
} else {
m_heap_size--;
}
if (m_heap_size > 1) {
put_at(1, m_heap[m_heap_size--]);
fix_heap_under(1);
} else {
m_heap_size--;
}
/// return the first element of the queue and removes it from the queue
}
/// return the first element of the queue and removes it from the queue
template <typename T> unsigned binary_heap_priority_queue<T>::dequeue() {
lean_assert(m_heap_size);
int ret = m_heap[1];
put_the_last_at_the_top_and_fix_the_heap();
m_heap_inverse[ret] = -1;
return ret;
}
lean_assert(m_heap_size);
int ret = m_heap[1];
put_the_last_at_the_top_and_fix_the_heap();
m_heap_inverse[ret] = -1;
return ret;
}
#ifdef LEAN_DEBUG
template <typename T> void binary_heap_priority_queue<T>::print() {
std::vector<int> index;
std::vector<T> prs;
while (size()) {
T prior;
int j = dequeue_and_get_priority(prior);
index.push_back(j);
prs.push_back(prior);
std::cout << "(" << j << ", " << prior << ")";
}
std::cout << std::endl;
// restore the queue
for (int i = 0; i < index.size(); i++)
enqueue(index[i], prs[i]);
std::vector<int> index;
std::vector<T> prs;
while (size()) {
T prior;
int j = dequeue_and_get_priority(prior);
index.push_back(j);
prs.push_back(prior);
std::cout << "(" << j << ", " << prior << ")";
}
#endif
std::cout << std::endl;
// restore the queue
for (int i = 0; i < index.size(); i++)
enqueue(index[i], prs[i]);
}
#endif
}

282
src/util/lp/binary_heap_upair_queue.cpp
View file

@ -9,242 +9,108 @@
#include "util/lp/binary_heap_upair_queue.h"
namespace lean {
template <typename T> binary_heap_upair_queue<T>::binary_heap_upair_queue(unsigned size) : m_q(size), m_pairs(size) {
lean_assert(size);
for (unsigned i = 0; i < size; i++)
m_available_spots.push_back(i);
lean_assert(size);
for (unsigned i = 0; i < size; i++)
m_available_spots.push_back(i);
}
template <typename T> unsigned
binary_heap_upair_queue<T>::dequeue_available_spot() {
binary_heap_upair_queue<T>::dequeue_available_spot() {
lean_assert(m_available_spots.empty() == false);
unsigned ret = m_available_spots.back();
m_available_spots.pop_back();
return ret;
}
template <typename T> void binary_heap_upair_queue<T>::remove(unsigned i, unsigned j) {
upair p(i, j);
auto it = m_pairs_to_index.find(p);
if (it == m_pairs_to_index.end())
return; // nothing to do
m_q.remove(it->second);
m_available_spots.push_back(it->second);
m_pairs_to_index.erase(it);
}
upair p(i, j);
auto it = m_pairs_to_index.find(p);
if (it == m_pairs_to_index.end())
return; // nothing to do
m_q.remove(it->second);
m_available_spots.push_back(it->second);
m_pairs_to_index.erase(it);
}
template <typename T> bool binary_heap_upair_queue<T>::ij_index_is_new(unsigned ij_index) const {
for (auto it : m_pairs_to_index) {
if (it.second == ij_index)
return false;
}
return true;
for (auto it : m_pairs_to_index) {
if (it.second == ij_index)
return false;
}
return true;
}
template <typename T> void binary_heap_upair_queue<T>::enqueue(unsigned i, unsigned j, const T & priority) {
upair p(i, j);
auto it = m_pairs_to_index.find(p);
unsigned ij_index;
if (it == m_pairs_to_index.end()) {
// it is a new pair, let us find a spot for it
if (m_available_spots.empty()) {
// we ran out of empty spots
unsigned size_was = m_pairs.size();
unsigned new_size = size_was << 1;
for (unsigned i = size_was; i < new_size; i++)
m_available_spots.push_back(i);
m_pairs.resize(new_size);
}
ij_index = dequeue_available_spot();
// lean_assert(ij_index<m_pairs.size() && ij_index_is_new(ij_index));
m_pairs[ij_index] = p;
m_pairs_to_index[p] = ij_index;
} else {
ij_index = it->second;
upair p(i, j);
auto it = m_pairs_to_index.find(p);
unsigned ij_index;
if (it == m_pairs_to_index.end()) {
// it is a new pair, let us find a spot for it
if (m_available_spots.empty()) {
// we ran out of empty spots
unsigned size_was = m_pairs.size();
unsigned new_size = size_was << 1;
for (unsigned i = size_was; i < new_size; i++)
m_available_spots.push_back(i);
m_pairs.resize(new_size);
}
m_q.enqueue(ij_index, priority);
ij_index = dequeue_available_spot();
// lean_assert(ij_index<m_pairs.size() && ij_index_is_new(ij_index));
m_pairs[ij_index] = p;
m_pairs_to_index[p] = ij_index;
} else {
ij_index = it->second;
}
m_q.enqueue(ij_index, priority);
}
template <typename T> void binary_heap_upair_queue<T>::dequeue(unsigned & i, unsigned &j) {
lean_assert(!m_q.is_empty());
unsigned ij_index = m_q.dequeue();
upair & p = m_pairs[ij_index];
i = p.first;
j = p.second;
m_available_spots.push_back(ij_index);
m_pairs_to_index.erase(p);
}
lean_assert(!m_q.is_empty());
unsigned ij_index = m_q.dequeue();
upair & p = m_pairs[ij_index];
i = p.first;
j = p.second;
m_available_spots.push_back(ij_index);
m_pairs_to_index.erase(p);
}
template <typename T> T binary_heap_upair_queue<T>::get_priority(unsigned i, unsigned j) const {
auto it = m_pairs_to_index.find(std::make_pair(i, j));
if (it == m_pairs_to_index.end())
return T(0xFFFFFF); // big number
return m_q.get_priority(it->second);
}
auto it = m_pairs_to_index.find(std::make_pair(i, j));
if (it == m_pairs_to_index.end())
return T(0xFFFFFF); // big number
return m_q.get_priority(it->second);
}
#ifdef LEAN_DEBUG
template <typename T> bool binary_heap_upair_queue<T>::pair_to_index_is_a_bijection() const {
std::set<int> tmp;
for (auto p : m_pairs_to_index) {
unsigned j = p.second;
auto it = tmp.find(j);
if (it != tmp.end()) {
std::cout << "for pair (" << p.first.first << ", " << p.first.second << "), the index " << j
<< " is already inside " << std::endl;
lean_assert(false);
} else {
tmp.insert(j);
}
std::set<int> tmp;
for (auto p : m_pairs_to_index) {
unsigned j = p.second;
auto it = tmp.find(j);
if (it != tmp.end()) {
std::cout << "for pair (" << p.first.first << ", " << p.first.second << "), the index " << j
<< " is already inside " << std::endl;
lean_assert(false);
} else {
tmp.insert(j);
}
return true;
}
return true;
}
template <typename T> bool binary_heap_upair_queue<T>::available_spots_are_correct() const {
std::set<int> tmp;
for (auto p : m_available_spots){
tmp.insert(p);
}
if (tmp.size() != m_available_spots.size())
return false;
for (auto it : m_pairs_to_index)
if (tmp.find(it.second) != tmp.end())
return false;
return true;
std::set<int> tmp;
for (auto p : m_available_spots){
tmp.insert(p);
}
if (tmp.size() != m_available_spots.size())
return false;
for (auto it : m_pairs_to_index)
if (tmp.find(it.second) != tmp.end())
return false;
return true;
}
#endif
}
/*
#include <unordered_set>
#include <queue>
#include <vector>
#include <set>
#include <utility>
typedef std::pair<unsigned, unsigned> upair;
namespace lean {
template <typename T>
class binary_heap_upair_queue {
binary_heap_priority_queue<T> m_q;
std::unordered_map<upair, unsigned> m_pairs_to_index;
std::vector<upair> m_pairs; // inverse to index
std::vector<unsigned> m_available_spots;
public:
binary_heap_upair_queue(unsigned size) : m_q(size), m_pairs(size) {
lean_assert(size);
for (unsigned i = 0; i < size; i++)
m_available_spots.push_back(i);
}
unsigned dequeue_available_spot() {
lean_assert(m_available_spots.empty() == false);
unsigned ret = m_available_spots.back();
m_available_spots.pop_back();
return ret;
}
void remove(unsigned i, unsigned j) {
upair p(i, j);
auto it = m_pairs_to_index.find(p);
if (it == m_pairs_to_index.end())
return; // nothing to do
m_q.remove(it->second);
m_available_spots.push_back(it->second);
m_pairs_to_index.erase(it);
}
bool ij_index_is_new(unsigned ij_index) const {
for (auto it : m_pairs_to_index) {
if (it.second == ij_index)
return false;
}
return true;
}
void enqueue(unsigned i, unsigned j, const T & priority) {
upair p(i, j);
auto it = m_pairs_to_index.find(p);
unsigned ij_index;
if (it == m_pairs_to_index.end()) {
// it is a new pair, let us find a spot for it
if (m_available_spots.empty()) {
// we ran out of empty spots
unsigned size_was = m_pairs.size();
unsigned new_size = size_was << 1;
for (unsigned i = size_was; i < new_size; i++)
m_available_spots.push_back(i);
m_pairs.resize(new_size);
}
ij_index = dequeue_available_spot();
// lean_assert(ij_index<m_pairs.size() && ij_index_is_new(ij_index));
m_pairs[ij_index] = p;
m_pairs_to_index[p] = ij_index;
} else {
ij_index = it->second;
}
m_q.enqueue(ij_index, priority);
}
void dequeue(unsigned & i, unsigned &j) {
lean_assert(!m_q.is_empty());
unsigned ij_index = m_q.dequeue();
upair & p = m_pairs[ij_index];
i = p.first;
j = p.second;
m_available_spots.push_back(ij_index);
m_pairs_to_index.erase(p);
}
bool is_empty() const {
return m_q.is_empty();
}
unsigned size() const {
return m_q.size();
}
bool contains(unsigned i, unsigned j) const {
return m_pairs_to_index.find(std::make_pair(i, j)) != m_pairs_to_index.end();
}
T get_priority(unsigned i, unsigned j) const {
auto it = m_pairs_to_index.find(std::make_pair(i, j));
if (it == m_pairs_to_index.end())
return T(0xFFFFFF); // big number
return m_q.get_priority(it->second);
}
bool pair_to_index_is_a_bijection() const {
std::set<int> tmp;
for (auto p : m_pairs_to_index) {
unsigned j = p.second;
auto it = tmp.find(j);
if (it != tmp.end()) {
std::cout << "for pair (" << p.first.first << ", " << p.first.second << "), the index " << j
<< " is already inside " << std::endl;
lean_assert(false);
} else {
tmp.insert(j);
}
}
return true;
}
bool available_spots_are_correct() const {
std::set<int> tmp;
for (auto p : m_available_spots){
tmp.insert(p);
}
if (tmp.size() != m_available_spots.size())
return false;
for (auto it : m_pairs_to_index)
if (tmp.find(it.second) != tmp.end())
return false;
return true;
}
bool is_correct() const {
return m_q.is_consistent() && pair_to_index_is_a_bijection() && available_spots_are_correct();
}
};
}
*/

498
src/util/lp/core_solver_pretty_printer.cpp
View file

@ -11,329 +11,329 @@ namespace lean {
template <typename T, typename X>
core_solver_pretty_printer<T, X>::core_solver_pretty_printer(lp_core_solver_base<T, X > & core_solver): m_core_solver(core_solver),
m_column_widths(core_solver.m_A.column_count(), 0),
m_A(core_solver.m_A.row_count(), vector<string>(core_solver.m_A.column_count(), "")),
m_signs(core_solver.m_A.row_count(), vector<string>(core_solver.m_A.column_count(), " ")),
m_costs(ncols(), ""),
m_cost_signs(ncols(), " "),
m_rs(ncols(), zero_of_type<X>()) {
m_w_buff = new T[m_core_solver.m_m];
m_ed_buff = new T[m_core_solver.m_m];
m_core_solver.save_state(m_w_buff, m_ed_buff);
init_m_A_and_signs();
init_costs();
init_column_widths();
init_rs_width();
m_cost_title = "costs";
m_basis_heading_title = "heading";
m_x_title = "x*";
m_title_width = std::max(std::max(m_cost_title.size(), std::max(m_basis_heading_title.size(), m_x_title.size())), m_approx_norm_title.size());
}
m_column_widths(core_solver.m_A.column_count(), 0),
m_A(core_solver.m_A.row_count(), vector<string>(core_solver.m_A.column_count(), "")),
m_signs(core_solver.m_A.row_count(), vector<string>(core_solver.m_A.column_count(), " ")),
m_costs(ncols(), ""),
m_cost_signs(ncols(), " "),
m_rs(ncols(), zero_of_type<X>()) {
m_w_buff = new T[m_core_solver.m_m];
m_ed_buff = new T[m_core_solver.m_m];
m_core_solver.save_state(m_w_buff, m_ed_buff);
init_m_A_and_signs();
init_costs();
init_column_widths();
init_rs_width();
m_cost_title = "costs";
m_basis_heading_title = "heading";
m_x_title = "x*";
m_title_width = std::max(std::max(m_cost_title.size(), std::max(m_basis_heading_title.size(), m_x_title.size())), m_approx_norm_title.size());
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_costs() {
vector<T> local_y(m_core_solver.m_m);
m_core_solver.solve_yB(local_y);
for (unsigned i = 0; i < ncols(); i++) {
if (m_core_solver.m_basis_heading[i] < 0) {
T t = m_core_solver.m_costs[i] - m_core_solver.m_A.dot_product_with_column(local_y, i);
set_coeff(m_costs, m_cost_signs, i, t, m_core_solver.column_name(i));
}
vector<T> local_y(m_core_solver.m_m);
m_core_solver.solve_yB(local_y);
for (unsigned i = 0; i < ncols(); i++) {
if (m_core_solver.m_basis_heading[i] < 0) {
T t = m_core_solver.m_costs[i] - m_core_solver.m_A.dot_product_with_column(local_y, i);
set_coeff(m_costs, m_cost_signs, i, t, m_core_solver.column_name(i));
}
}
}
template <typename T, typename X> core_solver_pretty_printer<T, X>::~core_solver_pretty_printer() {
m_core_solver.restore_state(m_w_buff, m_ed_buff);
delete [] m_w_buff;
delete [] m_ed_buff;
}
m_core_solver.restore_state(m_w_buff, m_ed_buff);
delete [] m_w_buff;
delete [] m_ed_buff;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_rs_width() {
m_rs_width = T_to_string(m_core_solver.get_cost()).size();
for (unsigned i = 0; i < nrows(); i++) {
auto wt = T_to_string(m_rs[i]).size();
if (wt > m_rs_width) {
m_rs_width = wt;
}
m_rs_width = T_to_string(m_core_solver.get_cost()).size();
for (unsigned i = 0; i < nrows(); i++) {
auto wt = T_to_string(m_rs[i]).size();
if (wt > m_rs_width) {
m_rs_width = wt;
}
}
}
template <typename T, typename X> T core_solver_pretty_printer<T, X>::current_column_norm() {
T ret = zero_of_type<T>();
for (T & ed : m_core_solver.m_ed)
ret += ed * ed;
return ret;
}
T ret = zero_of_type<T>();
for (T & ed : m_core_solver.m_ed)
ret += ed * ed;
return ret;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_m_A_and_signs() {
for (unsigned column = 0; column < ncols(); column++) {
m_core_solver.solve_Bd(column); // puts the result into m_core_solver.m_ed
string name = m_core_solver.column_name(column);
for (unsigned row = 0; row < nrows(); row ++) {
set_coeff(
m_A[row],
m_signs[row],
column,
m_core_solver.m_ed[row],
name);
m_rs[row] += m_core_solver.m_ed[row] * m_core_solver.m_x[column];
}
m_exact_column_norms.push_back(current_column_norm() + 1);
for (unsigned column = 0; column < ncols(); column++) {
m_core_solver.solve_Bd(column); // puts the result into m_core_solver.m_ed
string name = m_core_solver.column_name(column);
for (unsigned row = 0; row < nrows(); row ++) {
set_coeff(
m_A[row],
m_signs[row],
column,
m_core_solver.m_ed[row],
name);
m_rs[row] += m_core_solver.m_ed[row] * m_core_solver.m_x[column];
}
m_exact_column_norms.push_back(current_column_norm() + 1);
}
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_column_widths() {
for (unsigned i = 0; i < ncols(); i++) {
m_column_widths[i] = get_column_width(i);
}
for (unsigned i = 0; i < ncols(); i++) {
m_column_widths[i] = get_column_width(i);
}
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::adjust_width_with_low_bound(unsigned column, unsigned & w) {
if (!m_core_solver.low_bounds_are_set()) return;
w = std::max(w, (unsigned)T_to_string(m_core_solver.low_bound_value(column)).size());
}
if (!m_core_solver.low_bounds_are_set()) return;
w = std::max(w, (unsigned)T_to_string(m_core_solver.low_bound_value(column)).size());
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::adjust_width_with_upper_bound(unsigned column, unsigned & w) {
w = std::max(w, (unsigned)T_to_string(m_core_solver.upper_bound_value(column)).size());
}
w = std::max(w, (unsigned)T_to_string(m_core_solver.upper_bound_value(column)).size());
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::adjust_width_with_bounds(unsigned column, unsigned & w) {
switch (m_core_solver.get_column_type(column)) {
case fixed:
case boxed:
adjust_width_with_low_bound(column, w);
adjust_width_with_upper_bound(column, w);
break;
case low_bound:
adjust_width_with_low_bound(column, w);
break;
case upper_bound:
adjust_width_with_upper_bound(column, w);
break;
case free_column:
break;
default:
lean_assert(false);
break;
}
switch (m_core_solver.get_column_type(column)) {
case fixed:
case boxed:
adjust_width_with_low_bound(column, w);
adjust_width_with_upper_bound(column, w);
break;
case low_bound:
adjust_width_with_low_bound(column, w);
break;
case upper_bound:
adjust_width_with_upper_bound(column, w);
break;
case free_column:
break;
default:
lean_assert(false);
break;
}
}
template <typename T, typename X> unsigned core_solver_pretty_printer<T, X>:: get_column_width(unsigned column) {
unsigned w = std::max(m_costs[column].size(), T_to_string(m_core_solver.m_x[column]).size());
adjust_width_with_bounds(column, w);
adjust_width_with_basis_heading(column, w);
for (unsigned i = 0; i < nrows(); i++) {
unsigned cellw = m_A[i][column].size();
if (cellw > w) {
w = cellw;
}
unsigned w = std::max(m_costs[column].size(), T_to_string(m_core_solver.m_x[column]).size());
adjust_width_with_bounds(column, w);
adjust_width_with_basis_heading(column, w);
for (unsigned i = 0; i < nrows(); i++) {
unsigned cellw = m_A[i][column].size();
if (cellw > w) {
w = cellw;
}
w = std::max(w, (unsigned)T_to_string(m_exact_column_norms[column]).size());
w = std::max(w, (unsigned)T_to_string(m_core_solver.m_column_norms[column]).size());
return w;
}
w = std::max(w, (unsigned)T_to_string(m_exact_column_norms[column]).size());
w = std::max(w, (unsigned)T_to_string(m_core_solver.m_column_norms[column]).size());
return w;
}
template <typename T, typename X> std::string core_solver_pretty_printer<T, X>::regular_cell_string(unsigned row, unsigned column, std::string name) {
T t = fabs(m_core_solver.m_ed[row]);
if ( t == 1) return name;
return T_to_string(t) + name;
}
T t = fabs(m_core_solver.m_ed[row]);
if ( t == 1) return name;
return T_to_string(t) + name;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::set_coeff(vector<string>& row, vector<string> & row_signs, unsigned col, const T & t, string name) {
if (numeric_traits<T>::is_zero(t)) {
return;
if (numeric_traits<T>::is_zero(t)) {
return;
}
if (col > 0) {
if (t > 0) {
row_signs[col] = "+";
row[col] = t != 1? T_to_string(t) + name : name;
} else {
row_signs[col] = "-";
row[col] = t != -1? T_to_string(-t) + name: name;
}
if (col > 0) {
if (t > 0) {
row_signs[col] = "+";
row[col] = t != 1? T_to_string(t) + name : name;
} else {
row_signs[col] = "-";
row[col] = t != -1? T_to_string(-t) + name: name;
}
} else { // col == 0
if (t == -1) {
row[col] = "-" + name;
} else if (t == 1) {
row[col] = name;
} else {
row[col] = T_to_string(t) + name;
}
} else { // col == 0
if (t == -1) {
row[col] = "-" + name;
} else if (t == 1) {
row[col] = name;
} else {
row[col] = T_to_string(t) + name;
}
}
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_x() {
if (ncols() == 0) {
return;
}
int blanks = m_title_width + 1 - m_x_title.size();
std::cout << m_x_title;
print_blanks(blanks);
auto bh = m_core_solver.m_x;
for (unsigned i = 0; i < ncols(); i++) {
string s = T_to_string(bh[i]);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
std::cout << s << " "; // the column interval
}
std::cout << std::endl;
if (ncols() == 0) {
return;
}
int blanks = m_title_width + 1 - m_x_title.size();
std::cout << m_x_title;
print_blanks(blanks);
auto bh = m_core_solver.m_x;
for (unsigned i = 0; i < ncols(); i++) {
string s = T_to_string(bh[i]);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
std::cout << s << " "; // the column interval
}
std::cout << std::endl;
}
template <typename T, typename X> std::string core_solver_pretty_printer<T, X>::get_low_bound_string(unsigned j) {
switch (m_core_solver.get_column_type(j)){
case boxed:
case low_bound:
case fixed:
if (m_core_solver.low_bounds_are_set())
return T_to_string(m_core_solver.low_bound_value(j));
else
return std::string("0");
break;
default:
return std::string();
}
switch (m_core_solver.get_column_type(j)){
case boxed:
case low_bound:
case fixed:
if (m_core_solver.low_bounds_are_set())
return T_to_string(m_core_solver.low_bound_value(j));
else
return std::string("0");
break;
default:
return std::string();
}
}
template <typename T, typename X> std::string core_solver_pretty_printer<T, X>::get_upp_bound_string(unsigned j) {
switch (m_core_solver.get_column_type(j)){
case boxed:
case upper_bound:
case fixed:
return T_to_string(m_core_solver.upper_bound_value(j));
break;
default:
return std::string();
}
switch (m_core_solver.get_column_type(j)){
case boxed:
case upper_bound:
case fixed:
return T_to_string(m_core_solver.upper_bound_value(j));
break;
default:
return std::string();
}
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_lows() {
if (ncols() == 0) {
return;
}
int blanks = m_title_width + 1 - m_low_bounds_title.size();
std::cout << m_low_bounds_title;
print_blanks(blanks);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_low_bound_string(i);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
std::cout << s << " "; // the column interval
}
std::cout << std::endl;
if (ncols() == 0) {
return;
}
int blanks = m_title_width + 1 - m_low_bounds_title.size();
std::cout << m_low_bounds_title;
print_blanks(blanks);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_low_bound_string(i);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
std::cout << s << " "; // the column interval
}
std::cout << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_upps() {
if (ncols() == 0) {
return;
}
int blanks = m_title_width + 1 - m_upp_bounds_title.size();
std::cout << m_upp_bounds_title;
print_blanks(blanks);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_upp_bound_string(i);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
std::cout << s << " "; // the column interval
}
std::cout << std::endl;
if (ncols() == 0) {
return;
}
int blanks = m_title_width + 1 - m_upp_bounds_title.size();
std::cout << m_upp_bounds_title;
print_blanks(blanks);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_upp_bound_string(i);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
std::cout << s << " "; // the column interval
}
std::cout << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_exact_norms() {
int blanks = m_title_width + 1 - m_exact_norm_title.size();
std::cout << m_exact_norm_title;
int blanks = m_title_width + 1 - m_exact_norm_title.size();
std::cout << m_exact_norm_title;
print_blanks(blanks);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_exact_column_norm_string(i);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
for (unsigned i = 0; i < ncols(); i++) {
string s = get_exact_column_norm_string(i);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
std::cout << s << " ";
}
std::cout << std::endl;
std::cout << s << " ";
}
std::cout << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_approx_norms() {
int blanks = m_title_width + 1 - m_approx_norm_title.size();
std::cout << m_approx_norm_title;
int blanks = m_title_width + 1 - m_approx_norm_title.size();
std::cout << m_approx_norm_title;
print_blanks(blanks);
for (unsigned i = 0; i < ncols(); i++) {
string s = T_to_string(m_core_solver.m_column_norms[i]);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
for (unsigned i = 0; i < ncols(); i++) {
string s = T_to_string(m_core_solver.m_column_norms[i]);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
std::cout << s << " ";
}
std::cout << std::endl;
std::cout << s << " ";
}
std::cout << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print() {
for (unsigned i = 0; i < nrows(); i++) {
print_row(i);
}
print_bottom_line();
print_cost();
print_x();
print_basis_heading();
print_lows();
print_upps();
print_exact_norms();
print_approx_norms();
std::cout << std::endl;
for (unsigned i = 0; i < nrows(); i++) {
print_row(i);
}
print_bottom_line();
print_cost();
print_x();
print_basis_heading();
print_lows();
print_upps();
print_exact_norms();
print_approx_norms();
std::cout << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_basis_heading() {
int blanks = m_title_width + 1 - m_basis_heading_title.size();
std::cout << m_basis_heading_title;
print_blanks(blanks);
int blanks = m_title_width + 1 - m_basis_heading_title.size();
std::cout << m_basis_heading_title;
print_blanks(blanks);
if (ncols() == 0) {
return;
}
auto bh = m_core_solver.m_basis_heading;
for (unsigned i = 0; i < ncols(); i++) {
string s = T_to_string(bh[i]);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
std::cout << s << " "; // the column interval
}
std::cout << std::endl;
if (ncols() == 0) {
return;
}
auto bh = m_core_solver.m_basis_heading;
for (unsigned i = 0; i < ncols(); i++) {
string s = T_to_string(bh[i]);
int blanks = m_column_widths[i] - s.size();
print_blanks(blanks);
std::cout << s << " "; // the column interval
}
std::cout << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_cost() {
int blanks = m_title_width + 1 - m_cost_title.size();
std::cout << m_cost_title;
print_blanks(blanks);
print_given_rows(m_costs, m_cost_signs, m_core_solver.get_cost());
}
int blanks = m_title_width + 1 - m_cost_title.size();
std::cout << m_cost_title;
print_blanks(blanks);
print_given_rows(m_costs, m_cost_signs, m_core_solver.get_cost());
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_given_rows(vector<string> & row, vector<string> & signs, X rst) {
for (unsigned col = 0; col < row.size(); col++) {
unsigned width = m_column_widths[col];
string s = row[col];
int number_of_blanks = width - s.size();
lean_assert(number_of_blanks >= 0);
print_blanks(number_of_blanks);
std::cout << s << ' ';
if (col < row.size() - 1) {
std::cout << signs[col + 1] << ' ';
}
for (unsigned col = 0; col < row.size(); col++) {
unsigned width = m_column_widths[col];
string s = row[col];
int number_of_blanks = width - s.size();
lean_assert(number_of_blanks >= 0);
print_blanks(number_of_blanks);
std::cout << s << ' ';
if (col < row.size() - 1) {
std::cout << signs[col + 1] << ' ';
}
std::cout << '=';
string rs = T_to_string(rst);
int nb = m_rs_width - rs.size();
lean_assert(nb >= 0);
print_blanks(nb + 1);
std::cout << rs << std::endl;
}
std::cout << '=';
string rs = T_to_string(rst);
int nb = m_rs_width - rs.size();
lean_assert(nb >= 0);
print_blanks(nb + 1);
std::cout << rs << std::endl;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_row(unsigned i){
print_blanks(m_title_width + 1);
auto row = m_A[i];
auto sign_row = m_signs[i];
auto rs = m_rs[i];
print_given_rows(row, sign_row, rs);
}
print_blanks(m_title_width + 1);
auto row = m_A[i];
auto sign_row = m_signs[i];
auto rs = m_rs[i];
print_given_rows(row, sign_row, rs);
}
}

274
src/util/lp/dense_matrix.cpp
View file

@ -8,193 +8,193 @@
#include "util/lp/dense_matrix.h"
namespace lean {
template <typename T, typename X> dense_matrix<T,X>::dense_matrix(unsigned m, unsigned n) : m_m(m), m_n(n) {
m_values = new T[m * n];
for (unsigned i = 0; i < m * n; i ++)
m_values[i] = numeric_traits<T>::zero();
}
m_values = new T[m * n];
for (unsigned i = 0; i < m * n; i ++)
m_values[i] = numeric_traits<T>::zero();
}
template <typename T, typename X> dense_matrix<T,X> dense_matrix<T,X>::operator*=(matrix<T, X> const & a){
lean_assert(column_count() == a.row_count());
dense_matrix c(row_count(), a.column_count());
for (unsigned i = 0; i < row_count(); i++) {
for (unsigned j = 0; j < a.column_count(); j++) {
T v = numeric_traits<T>::zero();
for (unsigned k = 0; k < a.column_count(); k++) {
v += get_elem(i, k) * a(k, j);
}
c.set_elem(i, j, v);
lean_assert(column_count() == a.row_count());
dense_matrix c(row_count(), a.column_count());
for (unsigned i = 0; i < row_count(); i++) {
for (unsigned j = 0; j < a.column_count(); j++) {
T v = numeric_traits<T>::zero();
for (unsigned k = 0; k < a.column_count(); k++) {
v += get_elem(i, k) * a(k, j);
}
c.set_elem(i, j, v);
}
*this = c;
return *this;
}
*this = c;
return *this;
}
template <typename T, typename X> dense_matrix<T,X>&
dense_matrix<T,X>::operator=(matrix<T, X> const & other){
if ( this == & other)
return *this;
m_values = new T[m_m * m_n];
for (unsigned i = 0; i < m_m; i ++)
for (unsigned j = 0; j < m_n; j++)
m_values[i * m_n + j] = other.get_elem(i, j);
dense_matrix<T,X>::operator=(matrix<T, X> const & other){
if ( this == & other)
return *this;
}
m_values = new T[m_m * m_n];
for (unsigned i = 0; i < m_m; i ++)
for (unsigned j = 0; j < m_n; j++)
m_values[i * m_n + j] = other.get_elem(i, j);
return *this;
}
template <typename T, typename X> dense_matrix<T,X>&
dense_matrix<T,X>::operator=(dense_matrix const & other){
if ( this == & other)
return *this;
m_m = other.m_m;
m_n = other.m_n;
delete [] m_values;
m_values = new T[m_m * m_n];
for (unsigned i = 0; i < m_m; i ++)
for (unsigned j = 0; j < m_n; j++)
m_values[i * m_n + j] = other.get_elem(i, j);
dense_matrix<T,X>::operator=(dense_matrix const & other){
if ( this == & other)
return *this;
}
m_m = other.m_m;
m_n = other.m_n;
delete [] m_values;
m_values = new T[m_m * m_n];
for (unsigned i = 0; i < m_m; i ++)
for (unsigned j = 0; j < m_n; j++)
m_values[i * m_n + j] = other.get_elem(i, j);
return *this;
}
template <typename T, typename X> dense_matrix<T,X>::dense_matrix(dense_matrix<T, X> const & other) : m_m(other.row_count()), m_n(other.column_count()) {
m_values = new T[m_m * m_n];
for (unsigned i = 0; i < m_m; i ++)
for (unsigned j = 0; j < m_n; j++)
m_values[i * m_n + j] = other.get_elem(i, j);
}
m_values = new T[m_m * m_n];
for (unsigned i = 0; i < m_m; i ++)
for (unsigned j = 0; j < m_n; j++)
m_values[i * m_n + j] = other.get_elem(i, j);
}
template <typename T, typename X> dense_matrix<T,X>::dense_matrix(matrix<T, X> const & other) :
m_m(other.row_count()),
m_n(other.column_count()) {
m_values = new T[m_m * m_n];
for (unsigned i = 0; i < m_m; i++)
for (unsigned j = 0; j < m_n; j++)
m_values[i * m_n + j] = other.get_elem(i, j);
}
m_m(other.row_count()),
m_n(other.column_count()) {
m_values = new T[m_m * m_n];
for (unsigned i = 0; i < m_m; i++)
for (unsigned j = 0; j < m_n; j++)
m_values[i * m_n + j] = other.get_elem(i, j);
}
template <typename T, typename X> void dense_matrix<T, X>::apply_from_right(T * w) {
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(j, i);
}
t[i] = v;
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(j, i);
}
for (int i = 0; i < m_m; i++) {
w[i] = t[i];
}
delete [] t;
t[i] = v;
}
for (int i = 0; i < m_m; i++) {
w[i] = t[i];
}
delete [] t;
}
template <typename T, typename X> void dense_matrix<T, X>::apply_from_right(std::vector <T> & w) {
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(j, i);
}
t[i] = v;
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(j, i);
}
for (int i = 0; i < m_m; i++) {
w[i] = t[i];
}
delete [] t;
t[i] = v;
}
for (int i = 0; i < m_m; i++) {
w[i] = t[i];
}
delete [] t;
}
template <typename T, typename X> T* dense_matrix<T, X>::
apply_from_left_with_different_dims(std::vector<T> & w) {
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_n; j++) {
v += w[j]* get_elem(i, j);
}
t[i] = v;
apply_from_left_with_different_dims(std::vector<T> & w) {
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_n; j++) {
v += w[j]* get_elem(i, j);
}
return t;
t[i] = v;
}
return t;
}
template <typename T, typename X> void dense_matrix<T, X>::apply_from_left(std::vector<T> & w) {
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(i, j);
}
t[i] = v;
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(i, j);
}
for (int i = 0; i < m_m; i ++) {
w[i] = t[i];
}
delete [] t;
t[i] = v;
}
for (int i = 0; i < m_m; i ++) {
w[i] = t[i];
}
delete [] t;
}
template <typename T, typename X> void dense_matrix<T, X>::apply_from_left(X * w, lp_settings & ) {
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(i, j);
}
t[i] = v;
T * t = new T[m_m];
for (int i = 0; i < m_m; i ++) {
T v = numeric_traits<T>::zero();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(i, j);
}
for (int i = 0; i < m_m; i ++) {
w[i] = t[i];
}
delete [] t;
t[i] = v;
}
for (int i = 0; i < m_m; i ++) {
w[i] = t[i];
}
delete [] t;
}
template <typename T, typename X> void dense_matrix<T, X>::apply_from_left_to_X(std::vector<X> & w, lp_settings & ) {
std::vector<X> t(m_m);
for (int i = 0; i < m_m; i ++) {
X v = zero_of_type<X>();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(i, j);
}
t[i] = v;
}
for (int i = 0; i < m_m; i ++) {
w[i] = t[i];
std::vector<X> t(m_m);
for (int i = 0; i < m_m; i ++) {
X v = zero_of_type<X>();
for (int j = 0; j < m_m; j++) {
v += w[j]* get_elem(i, j);
}
t[i] = v;
}
// This method pivots row i to row i0 by muliplying row i by
// alpha and adding it to row i0.
for (int i = 0; i < m_m; i ++) {
w[i] = t[i];
}
}
// This method pivots row i to row i0 by muliplying row i by
// alpha and adding it to row i0.
template <typename T, typename X> void dense_matrix<T, X>::pivot_row_to_row(unsigned i, T alpha, unsigned i0,
double & pivot_epsilon) {
thread_local T _0 = numeric_traits<T>::zero();
for (unsigned j = 0; j < m_n; j++) {
m_values[i0 * m_n + j] += m_values[i * m_n + j] * alpha;
if (fabs(m_values[i0 + m_n + j]) < pivot_epsilon) {
m_values[i0 + m_n + j] = _0;
}
double & pivot_epsilon) {
thread_local T _0 = numeric_traits<T>::zero();
for (unsigned j = 0; j < m_n; j++) {
m_values[i0 * m_n + j] += m_values[i * m_n + j] * alpha;
if (fabs(m_values[i0 + m_n + j]) < pivot_epsilon) {
m_values[i0 + m_n + j] = _0;
}
}
}
template <typename T, typename X> void dense_matrix<T, X>::swap_columns(unsigned a, unsigned b) {
for (unsigned i = 0; i < m_m; i++) {
T t = get_elem(i, a);
set_elem(i, a, get_elem(i, b));
set_elem(i, b, t);
}
for (unsigned i = 0; i < m_m; i++) {
T t = get_elem(i, a);
set_elem(i, a, get_elem(i, b));
set_elem(i, b, t);
}
}
template <typename T, typename X> void dense_matrix<T, X>::swap_rows(unsigned a, unsigned b) {
for (unsigned i = 0; i < m_n; i++) {
T t = get_elem(a, i);
set_elem(a, i, get_elem(b, i));
set_elem(b, i, t);
}
for (unsigned i = 0; i < m_n; i++) {
T t = get_elem(a, i);
set_elem(a, i, get_elem(b, i));
set_elem(b, i, t);
}
}
template <typename T, typename X> void dense_matrix<T, X>::multiply_row_by_constant(unsigned row, T & t) {
for (unsigned i = 0; i < m_n; i++) {
set_elem(row, i, t * get_elem(row, i));
}
for (unsigned i = 0; i < m_n; i++) {
set_elem(row, i, t * get_elem(row, i));
}
}
template <typename T, typename X>
dense_matrix<T, X> operator* (matrix<T, X> & a, matrix<T, X> & b){

76
src/util/lp/indexed_vector.cpp
View file

@ -38,57 +38,57 @@ void print_vector(const std::vector<mpq> & t) {
}
template <typename T>
void indexed_vector<T>::resize(unsigned data_size) {
m_index.clear();
m_data.resize(data_size, numeric_traits<T>::zero());
}
void indexed_vector<T>::resize(unsigned data_size) {
m_index.clear();
m_data.resize(data_size, numeric_traits<T>::zero());
}
template <typename T>
void indexed_vector<T>::set_value(T value, unsigned index) {
m_data[index] = value;
m_index.push_back(index);
}
void indexed_vector<T>::set_value(T value, unsigned index) {
m_data[index] = value;
m_index.push_back(index);
}
template <typename T>
void indexed_vector<T>::clear() {
for (unsigned i : m_index)
m_data[i] = numeric_traits<T>::zero();
m_index.clear();
}
void indexed_vector<T>::clear() {
for (unsigned i : m_index)
m_data[i] = numeric_traits<T>::zero();
m_index.clear();
}
template <typename T>
void indexed_vector<T>::clear_all() {
unsigned i = m_data.size();
while (i--) m_data[i] = numeric_traits<T>::zero();
m_index.clear();
}
void indexed_vector<T>::clear_all() {
unsigned i = m_data.size();
while (i--) m_data[i] = numeric_traits<T>::zero();
m_index.clear();
}
template <typename T>
void indexed_vector<T>::erase_from_index(unsigned j) {
auto it = std::find(m_index.begin(), m_index.end(), j);
if (it != m_index.end()) m_index.erase(it);
}
void indexed_vector<T>::erase_from_index(unsigned j) {
auto it = std::find(m_index.begin(), m_index.end(), j);
if (it != m_index.end()) m_index.erase(it);
}
#ifdef LEAN_DEBUG
template <typename T>
bool indexed_vector<T>::is_OK() const {
int size = 0;
for (unsigned i = 0; i < m_data.size(); i++) {
if (!is_zero(m_data[i])) {
if (std::find(m_index.begin(), m_index.end(), i) == m_index.end())
return false;
size++;
}
bool indexed_vector<T>::is_OK() const {
int size = 0;
for (unsigned i = 0; i < m_data.size(); i++) {
if (!is_zero(m_data[i])) {
if (std::find(m_index.begin(), m_index.end(), i) == m_index.end())
return false;
size++;
}
return size == m_index.size();
}
return size == m_index.size();
}
template <typename T>
void indexed_vector<T>::print() {
std::cout << "m_index " << std::endl;
for (unsigned i = 0; i < m_index.size(); i++) {
std::cout << m_index[i] << " ";
}
std::cout << std::endl;
print_vector(m_data);
void indexed_vector<T>::print() {
std::cout << "m_index " << std::endl;
for (unsigned i = 0; i < m_index.size(); i++) {
std::cout << m_index[i] << " ";
}
std::cout << std::endl;
print_vector(m_data);
}
#endif
}

12
src/util/lp/lar_solver.cpp
View file

@ -17,7 +17,7 @@ double conversion_helper <double>::get_low_bound(const column_info<mpq> & ci) {
return ci.get_low_bound().get_double() + eps;
}
double conversion_helper <double>::get_upper_bound(const column_info<mpq> & ci) {
double conversion_helper <double>::get_upper_bound(const column_info<mpq> & ci) {
if (!ci.upper_bound_is_strict())
return ci.get_upper_bound().get_double();
double eps = 0.00001;
@ -285,11 +285,11 @@ template <typename V> V lar_solver::get_column_val(std::vector<V> & low_bound, s
}
lar_solver::~lar_solver() {
std::vector<canonic_left_side*> to_delete;
for (auto it : m_canonic_left_sides)
to_delete.push_back(it);
for (auto t : to_delete)
delete t;
std::vector<canonic_left_side*> to_delete;
for (auto it : m_canonic_left_sides)
to_delete.push_back(it);
for (auto t : to_delete)
delete t;
}

4
src/util/lp/lp_dual_simplex.cpp
View file

@ -257,8 +257,8 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fill_costs_and_bou
}
template <typename T, typename X> void lp_dual_simplex<T, X>::fill_first_stage_solver_fields_for_row_slack_and_artificial(unsigned row,
unsigned & slack_var,
unsigned & artificial) {
unsigned & slack_var,
unsigned & artificial) {
lean_assert(row < this->row_count());
auto & constraint = this->m_constraints[this->m_core_solver_rows_to_external_rows[row]];
// we need to bring the program to the form Ax = b

18
src/util/lp/lp_primal_core_solver.cpp
View file

@ -210,15 +210,15 @@ template <typename T, typename X> X lp_primal_core_solver<T, X>::get_max_boun
// stage1 constructor
template <typename T, typename X> lp_primal_core_solver<T, X>:: lp_primal_core_solver(static_matrix<T, X> & A,
std::vector<X> & b, // the right side vector
std::vector<X> & x, // the number of elements in x needs to be at least as large as the number of columns in A
std::vector<unsigned> & basis,
std::vector<T> & costs,
std::vector<column_type> & column_type_array,
std::vector<X> & low_bound_values,
std::vector<X> & upper_bound_values,
lp_settings & settings,
std::unordered_map<unsigned, std::string> const & column_names):
std::vector<X> & b, // the right side vector
std::vector<X> & x, // the number of elements in x needs to be at least as large as the number of columns in A
std::vector<unsigned> & basis,
std::vector<T> & costs,
std::vector<column_type> & column_type_array,
std::vector<X> & low_bound_values,
std::vector<X> & upper_bound_values,
lp_settings & settings,
std::unordered_map<unsigned, std::string> const & column_names):
lp_core_solver_base<T, X>(A, b,
basis,
x,

6
src/util/lp/lp_primal_core_solver_instances.cpp
View file

@ -7,7 +7,7 @@
#include "util/lp/lar_solver.h"
#include "util/lp/lp_primal_core_solver.cpp"
namespace lean {
//template void lar_solver::find_solution_signature_with_doubles(lar_solution_signature&);
//template void lar_solver::find_solution_signature_with_doubles(lar_solution_signature&);
template void lp_primal_core_solver<double, double>::find_feasible_solution();
template lp_primal_core_solver<double, double>::lp_primal_core_solver(static_matrix<double, double>&, std::vector<double, std::allocator<double> >&, std::vector<double, std::allocator<double> >&, std::vector<unsigned int, std::allocator<unsigned int> >&, std::vector<double, std::allocator<double> >&, std::vector<column_type, std::allocator<column_type> >&, std::vector<double, std::allocator<double> >&, lp_settings&, std::unordered_map<unsigned int, std::string, std::hash<unsigned int>, std::equal_to<unsigned int>, std::allocator<std::pair<unsigned int const, std::string> > > const&);
template lp_primal_core_solver<double, double>::lp_primal_core_solver(static_matrix<double, double>&, std::vector<double, std::allocator<double> >&, std::vector<double, std::allocator<double> >&, std::vector<unsigned int, std::allocator<unsigned int> >&, std::vector<double, std::allocator<double> >&, std::vector<column_type, std::allocator<column_type> >&, std::vector<double, std::allocator<double> >&, std::vector<double, std::allocator<double> >&, lp_settings&, std::unordered_map<unsigned int, std::string, std::hash<unsigned int>, std::equal_to<unsigned int>, std::allocator<std::pair<unsigned int const, std::string> > > const&);
@ -15,6 +15,6 @@ template lp_primal_core_solver<double, double>::lp_primal_core_solver(static_mat
template unsigned lp_primal_core_solver<double, double>::solve();
template lp_primal_core_solver<mpq, mpq>::lp_primal_core_solver(static_matrix<mpq, mpq>&, std::vector<mpq, std::allocator<mpq> >&, std::vector<mpq, std::allocator<mpq> >&, std::vector<unsigned int, std::allocator<unsigned int> >&, std::vector<mpq, std::allocator<mpq> >&, std::vector<column_type, std::allocator<column_type> >&, std::vector<mpq, std::allocator<mpq> >&, std::vector<mpq, std::allocator<mpq> >&, lp_settings&, std::unordered_map<unsigned int, std::string, std::hash<unsigned int>, std::equal_to<unsigned int>, std::allocator<std::pair<unsigned int const, std::string> > > const&);
template unsigned lp_primal_core_solver<mpq, mpq>::solve();
//template void lp_primal_simplex<double, double>::solve_with_total_inf();
//template void lp_primal_simplex<mpq, mpq>::solve_with_total_inf();
//template void lp_primal_simplex<double, double>::solve_with_total_inf();
//template void lp_primal_simplex<mpq, mpq>::solve_with_total_inf();
}

6
src/util/lp/lp_primal_simplex.cpp
View file

@ -76,9 +76,9 @@ template <typename T, typename X> void lp_primal_simplex<T, X>::set_zero_bound(b
}
template <typename T, typename X> void lp_primal_simplex<T, X>::fill_costs_and_x_for_first_stage_solver_for_row(
int row,
unsigned & slack_var,
unsigned & artificial) {
int row,
unsigned & slack_var,
unsigned & artificial) {
lean_assert(row >= 0 && row < this->row_count());
auto & constraint = this->m_constraints[this->m_core_solver_rows_to_external_rows[row]];
// we need to bring the program to the form Ax = b

2
src/util/lp/lp_settings.cpp
View file

@ -69,7 +69,7 @@ int get_millisecond_span(int start_time) {
}
void my_random_init(unsigned * seed) {
#ifdef LEAN_WINDOWS
srand(*seed);
srand(*seed);
#else
rand_r(seed);
#endif

34
src/util/lp/matrix.cpp
View file

@ -9,26 +9,26 @@
namespace lean {
template <typename T, typename X>
bool matrix<T, X>::is_equal(const matrix<T, X>& other) {
if (other.row_count() != row_count() || other.column_count() != column_count())
return false;
for (unsigned i = 0; i < row_count(); i++) {
for (unsigned j = 0; j < column_count(); j++) {
auto a = get_elem(i, j);
auto b = other.get_elem(i, j);
if (numeric_traits<T>::precise()) {
if (a != b) return false;
} else if (fabs(numeric_traits<T>::get_double(a - b)) > 0.000001) {
// cout << "returning false from operator== of matrix comparison" << endl;
// cout << "this matrix is " << endl;
// print_matrix(*this);
// cout << "other matrix is " << endl;
// print_matrix(other);
return false;
}
if (other.row_count() != row_count() || other.column_count() != column_count())
return false;
for (unsigned i = 0; i < row_count(); i++) {
for (unsigned j = 0; j < column_count(); j++) {
auto a = get_elem(i, j);
auto b = other.get_elem(i, j);
if (numeric_traits<T>::precise()) {
if (a != b) return false;
} else if (fabs(numeric_traits<T>::get_double(a - b)) > 0.000001) {
// cout << "returning false from operator== of matrix comparison" << endl;
// cout << "this matrix is " << endl;
// print_matrix(*this);
// cout << "other matrix is " << endl;
// print_matrix(other);
return false;
}
}
return true;
}
return true;
}
template <typename T, typename X>
void apply_to_vector(matrix<T, X> & m, T * w) {

248
src/util/lp/row_eta_matrix.cpp
View file

@ -9,146 +9,146 @@
namespace lean {
template <typename T, typename X>
void row_eta_matrix<T, X>::apply_from_left(std::vector<X> & w, lp_settings &) {
// #ifdef LEAN_DEBUG
// dense_matrix<T> deb(*this);
// auto clone_w = clone_vector<T>(w, m_dimension);
// deb.apply_from_left(clone_w, settings);
// #endif
auto w_at_row = w[m_row];
for (auto & it :m_row_vector.m_data) {
w_at_row += w[it.first] * it.second;
// #ifdef LEAN_DEBUG
// dense_matrix<T> deb(*this);
// auto clone_w = clone_vector<T>(w, m_dimension);
// deb.apply_from_left(clone_w, settings);
// #endif
auto w_at_row = w[m_row];
for (auto & it :m_row_vector.m_data) {
w_at_row += w[it.first] * it.second;
}
w[m_row] = w_at_row;
// #ifdef LEAN_DEBUG
// lean_assert(vectors_are_equal<T>(clone_w, w, m_dimension));
// delete [] clone_w;
// #endif
}
template <typename T, typename X>
void row_eta_matrix<T, X>::apply_from_left_local_to_T(indexed_vector<T> & w, lp_settings & settings) {
auto w_at_row = w[m_row];
bool was_zero_at_m_row = is_zero(w_at_row);
for (auto & it : m_row_vector.m_data) {
w_at_row += w[it.first] * it.second;
}
if (!settings.abs_val_is_smaller_than_drop_tolerance(w_at_row)){
if (was_zero_at_m_row) {
w.m_index.push_back(m_row);
}
w[m_row] = w_at_row;
// #ifdef LEAN_DEBUG
// lean_assert(vectors_are_equal<T>(clone_w, w, m_dimension));
// delete [] clone_w;
// #endif
} else if (!was_zero_at_m_row){
w[m_row] = zero_of_type<T>();
auto it = std::find(w.m_index.begin(), w.m_index.end(), m_row);
w.m_index.erase(it);
}
template <typename T, typename X>
void row_eta_matrix<T, X>::apply_from_left_local_to_T(indexed_vector<T> & w, lp_settings & settings) {
auto w_at_row = w[m_row];
bool was_zero_at_m_row = is_zero(w_at_row);
for (auto & it : m_row_vector.m_data) {
w_at_row += w[it.first] * it.second;
}
if (!settings.abs_val_is_smaller_than_drop_tolerance(w_at_row)){
if (was_zero_at_m_row) {
w.m_index.push_back(m_row);
}
w[m_row] = w_at_row;
} else if (!was_zero_at_m_row){
w[m_row] = zero_of_type<T>();
auto it = std::find(w.m_index.begin(), w.m_index.end(), m_row);
w.m_index.erase(it);
}
lean_assert(check_vector_for_small_values(w, settings));
lean_assert(check_vector_for_small_values(w, settings));
}
template <typename T, typename X>
void row_eta_matrix<T, X>::apply_from_left_local_to_X(indexed_vector<X> & w, lp_settings & settings) {
auto w_at_row = w[m_row];
bool was_zero_at_m_row = is_zero(w_at_row);
void row_eta_matrix<T, X>::apply_from_left_local_to_X(indexed_vector<X> & w, lp_settings & settings) {
auto w_at_row = w[m_row];
bool was_zero_at_m_row = is_zero(w_at_row);
for (auto & it : m_row_vector.m_data) {
w_at_row += w[it.first] * it.second;
}
for (auto & it : m_row_vector.m_data) {
w_at_row += w[it.first] * it.second;
}
if (!settings.abs_val_is_smaller_than_drop_tolerance(w_at_row)){
if (was_zero_at_m_row) {
w.m_index.push_back(m_row);
}
w[m_row] = w_at_row;
} else if (!was_zero_at_m_row){
w[m_row] = zero_of_type<X>();
auto it = std::find(w.m_index.begin(), w.m_index.end(), m_row);
w.m_index.erase(it);
if (!settings.abs_val_is_smaller_than_drop_tolerance(w_at_row)){
if (was_zero_at_m_row) {
w.m_index.push_back(m_row);
}
lean_assert(check_vector_for_small_values(w, settings));
w[m_row] = w_at_row;
} else if (!was_zero_at_m_row){
w[m_row] = zero_of_type<X>();
auto it = std::find(w.m_index.begin(), w.m_index.end(), m_row);
w.m_index.erase(it);
}
lean_assert(check_vector_for_small_values(w, settings));
}
template <typename T, typename X>
void row_eta_matrix<T, X>::apply_from_right(std::vector<T> & w) {
T w_row = w[m_row];
if (numeric_traits<T>::is_zero(w_row)) return;
void row_eta_matrix<T, X>::apply_from_right(std::vector<T> & w) {
T w_row = w[m_row];
if (numeric_traits<T>::is_zero(w_row)) return;
#ifdef LEAN_DEBUG
// dense_matrix<T> deb(*this);
// auto clone_w = clone_vector<T>(w, m_dimension);
// deb.apply_from_right(clone_w);
#endif
for (auto & it : m_row_vector.m_data) {
w[it.first] += w_row * it.second;
}
#ifdef LEAN_DEBUG
// lean_assert(vectors_are_equal<T>(clone_w, w, m_dimension));
// delete clone_w;
#endif
}
template <typename T, typename X>
void row_eta_matrix<T, X>::apply_from_right(indexed_vector<T> & w) {
T w_row = w[m_row];
if (numeric_traits<T>::is_zero(w_row)) return;
#ifdef LEAN_DEBUG
// dense_matrix<T> deb(*this);
// auto clone_w = clone_vector<T>(w.m_data, m_dimension);
// deb.apply_from_right(clone_w);
#endif
for (auto & it : m_row_vector.m_data) {
T old_val = w[it.first];
T v = w[it.index()] += w_row * it.second;
if (numeric_traits<T>::is_zero(old_val)) {
w.m_index.push_back(it.index());
} else if (numeric_traits<T>::is_zero(v)) { // it is a very rare case
auto w_it = std::find(w.m_index.begin(), w.m_index.end(), it.index());
lean_assert(w_it != w.m_index.end());
w.m_index.erase(w_it);
}
}
#ifdef LEAN_DEBUG
// lean_assert(vectors_are_equal<T>(clone_w, w.m_data, m_dimension));
// for (unsigned i = 0; i < m_dimension; i++) {
// if (!numeric_traits<T>::is_zero(w.m_data[i])) {
// lean_assert(std::find(w.m_index.begin(), w.m_index.end(), i) != w.m_index.end());
// }
// }
// delete clone_w;
#endif
}
template <typename T, typename X>
void row_eta_matrix<T, X>::conjugate_by_permutation(permutation_matrix<T, X> & p) {
// this = p * this * p(-1)
#ifdef LEAN_DEBUG
// auto rev = p.get_reverse();
// auto deb = ((*this) * rev);
// deb = p * deb;
#endif
m_row = p.apply_reverse(m_row);
// copy aside the column indices
std::vector<unsigned> columns;
for (auto & it : m_row_vector.m_data)
columns.push_back(it.first);
for (unsigned i = columns.size(); i-- > 0;)
m_row_vector.m_data[i].first = p.get_rev(columns[i]);
#ifdef LEAN_DEBUG
// lean_assert(deb == *this);
// dense_matrix<T> deb(*this);
// auto clone_w = clone_vector<T>(w, m_dimension);
// deb.apply_from_right(clone_w);
#endif
for (auto & it : m_row_vector.m_data) {
w[it.first] += w_row * it.second;
}
#ifdef LEAN_DEBUG
template <typename T, typename X>
T row_eta_matrix<T, X>::get_elem(unsigned row, unsigned col) const {
if (row == m_row){
if (col == row) {
return numeric_traits<T>::one();
}
return m_row_vector[col];
}
return col == row ? numeric_traits<T>::one() : numeric_traits<T>::zero();
}
// lean_assert(vectors_are_equal<T>(clone_w, w, m_dimension));
// delete clone_w;
#endif
}
template <typename T, typename X>
void row_eta_matrix<T, X>::apply_from_right(indexed_vector<T> & w) {
T w_row = w[m_row];
if (numeric_traits<T>::is_zero(w_row)) return;
#ifdef LEAN_DEBUG
// dense_matrix<T> deb(*this);
// auto clone_w = clone_vector<T>(w.m_data, m_dimension);
// deb.apply_from_right(clone_w);
#endif
for (auto & it : m_row_vector.m_data) {
T old_val = w[it.first];
T v = w[it.index()] += w_row * it.second;
if (numeric_traits<T>::is_zero(old_val)) {
w.m_index.push_back(it.index());
} else if (numeric_traits<T>::is_zero(v)) { // it is a very rare case
auto w_it = std::find(w.m_index.begin(), w.m_index.end(), it.index());
lean_assert(w_it != w.m_index.end());
w.m_index.erase(w_it);
}
}
#ifdef LEAN_DEBUG
// lean_assert(vectors_are_equal<T>(clone_w, w.m_data, m_dimension));
// for (unsigned i = 0; i < m_dimension; i++) {
// if (!numeric_traits<T>::is_zero(w.m_data[i])) {
// lean_assert(std::find(w.m_index.begin(), w.m_index.end(), i) != w.m_index.end());
// }
// }
// delete clone_w;
#endif
}
template <typename T, typename X>
void row_eta_matrix<T, X>::conjugate_by_permutation(permutation_matrix<T, X> & p) {
// this = p * this * p(-1)
#ifdef LEAN_DEBUG
// auto rev = p.get_reverse();
// auto deb = ((*this) * rev);
// deb = p * deb;
#endif
m_row = p.apply_reverse(m_row);
// copy aside the column indices
std::vector<unsigned> columns;
for (auto & it : m_row_vector.m_data)
columns.push_back(it.first);
for (unsigned i = columns.size(); i-- > 0;)
m_row_vector.m_data[i].first = p.get_rev(columns[i]);
#ifdef LEAN_DEBUG
// lean_assert(deb == *this);
#endif
}
#ifdef LEAN_DEBUG
template <typename T, typename X>
T row_eta_matrix<T, X>::get_elem(unsigned row, unsigned col) const {
if (row == m_row){
if (col == row) {
return numeric_traits<T>::one();
}
return m_row_vector[col];
}
return col == row ? numeric_traits<T>::one() : numeric_traits<T>::zero();
}
#endif
}

34
src/util/lp/sparse_matrix_instances.cpp
View file

@ -72,23 +72,23 @@ template void sparse_matrix<mpq, mpq>::double_solve_U_y<mpq>(std::vector<mpq>&);
template void sparse_matrix<mpq, numeric_pair<mpq>>::double_solve_U_y<mpq>(std::vector<mpq>&);
template void sparse_matrix<mpq, numeric_pair<mpq> >::double_solve_U_y<numeric_pair<mpq> >(std::vector<numeric_pair<mpq>>&);
/////////////////////
/*
template void lu<double, double>::create_initial_factorization();
template void lu<mpq, mpq>::create_initial_factorization();
template void lu<mpq, mpq>::replace_column(unsigned int, mpq, indexed_vector<mpq>&);
template void lu<mpq, numeric_pair<mpq> >::create_initial_factorization();
template void lu<mpq, numeric_pair<mpq> >::replace_column(unsigned int, mpq, indexed_vector<mpq>&);
template void lu<double, double>::init_vector_w(unsigned int, indexed_vector<double>&);
template void lu<mpq, mpq>::find_error_of_yB(std::vector<mpq, std::allocator<mpq> >&, std::vector<mpq, std::allocator<mpq> > const&);
template void lu<mpq, mpq>::init_vector_w(unsigned int, indexed_vector<mpq>&);
template void lu<mpq, numeric_pair<mpq> >::find_error_of_yB(std::vector<mpq, std::allocator<mpq> >&, std::vector<mpq, std::allocator<mpq> > const&);
template void lu<mpq, numeric_pair<mpq> >::init_vector_w(unsigned int, indexed_vector<mpq>&);
template void lu<double, double>::find_error_of_yB(std::vector<double, std::allocator<double> >&, std::vector<double, std::allocator<double> > const&);
#ifdef LEAN_DEBUG
template void print_matrix<double, double>(static_matrix<double, double>&);
#endif
*/
/////////////////////
/*
template void lu<double, double>::create_initial_factorization();
template void lu<mpq, mpq>::create_initial_factorization();
template void lu<mpq, mpq>::replace_column(unsigned int, mpq, indexed_vector<mpq>&);
template void lu<mpq, numeric_pair<mpq> >::create_initial_factorization();
template void lu<mpq, numeric_pair<mpq> >::replace_column(unsigned int, mpq, indexed_vector<mpq>&);
template void lu<double, double>::init_vector_w(unsigned int, indexed_vector<double>&);
template void lu<mpq, mpq>::find_error_of_yB(std::vector<mpq, std::allocator<mpq> >&, std::vector<mpq, std::allocator<mpq> > const&);
template void lu<mpq, mpq>::init_vector_w(unsigned int, indexed_vector<mpq>&);
template void lu<mpq, numeric_pair<mpq> >::find_error_of_yB(std::vector<mpq, std::allocator<mpq> >&, std::vector<mpq, std::allocator<mpq> > const&);
template void lu<mpq, numeric_pair<mpq> >::init_vector_w(unsigned int, indexed_vector<mpq>&);
template void lu<double, double>::find_error_of_yB(std::vector<double, std::allocator<double> >&, std::vector<double, std::allocator<double> > const&);
#ifdef LEAN_DEBUG
template void print_matrix<double, double>(static_matrix<double, double>&);
#endif
*/
#ifdef LEAN_DEBUG
template bool sparse_matrix<double, double>::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const;
template bool sparse_matrix<mpq, mpq>::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const;

18
src/util/lp/static_matrix.cpp
View file

@ -20,7 +20,7 @@ void static_matrix<T, X>::init_row_columns(unsigned m, unsigned n) {
}
}
// constructor that copies columns of the basis from A
// constructor that copies columns of the basis from A
template <typename T, typename X>
static_matrix<T, X>::static_matrix(static_matrix const &A, unsigned * basis) :
m_work_pivot_vector(A.row_count(), numeric_traits<T>::zero()) {
@ -52,7 +52,7 @@ template <typename T, typename X> void static_matrix<T, X>::init_empty_matrix
}
template <typename T, typename X>
template <typename L>
L static_matrix<T, X>::dot_product_with_row(unsigned row, const std::vector<L> & w) {
L static_matrix<T, X>::dot_product_with_row(unsigned row, const std::vector<L> & w) {
L ret = zero_of_type<L>();
lean_assert(row < m_rows.size());
for (auto & it : m_rows[row]) {
@ -107,7 +107,7 @@ template <typename T, typename X> void static_matrix<T, X>::remove_last_colum
break;
}
offset--;
}
}
}
m_columns.pop_back();
}
@ -144,7 +144,7 @@ template <typename T, typename X> void static_matrix<T, X>::regen_domain() {
for (auto & t : m_rows[i]) {
m_domain.insert(std::make_pair(i, t.m_j));
}
}
}
}
#endif
@ -183,12 +183,12 @@ template <typename T, typename X> void static_matrix<T, X>::copy_column_to_ve
}
template <typename T, typename X> void static_matrix<T, X>::copy_column_to_vector (unsigned j, std::vector<T> & v) const {
v.resize(row_count(), numeric_traits<T>::zero());
for (auto & it : m_columns[j]) {
if (!is_zero(it.m_value))
v[it.m_i]=it.m_value;
}
v.resize(row_count(), numeric_traits<T>::zero());
for (auto & it : m_columns[j]) {
if (!is_zero(it.m_value))
v[it.m_i]=it.m_value;
}
}
template <typename T, typename X> void static_matrix<T, X>::add_column_to_vector (const T & a, unsigned j, T * v) const {
for (auto & it : m_columns[j]) {