diff --git a/src/util/lp/binary_heap_priority_queue.cpp b/src/util/lp/binary_heap_priority_queue.cpp index 11abf41a4..f97ab50c0 100644 --- a/src/util/lp/binary_heap_priority_queue.cpp +++ b/src/util/lp/binary_heap_priority_queue.cpp @@ -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 void binary_heap_priority_queue::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 void binary_heap_priority_queue::put_at(unsigned i, unsigned h) { @@ -19,183 +19,183 @@ template void binary_heap_priority_queue::put_at(unsigned i, uns } template void binary_heap_priority_queue::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 bool binary_heap_priority_queue::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(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(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 void binary_heap_priority_queue::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 binary_heap_priority_queue::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 void binary_heap_priority_queue::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 void binary_heap_priority_queue::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 void binary_heap_priority_queue::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 void binary_heap_priority_queue::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 void binary_heap_priority_queue::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 unsigned binary_heap_priority_queue::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 void binary_heap_priority_queue::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 void binary_heap_priority_queue::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 unsigned binary_heap_priority_queue::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 void binary_heap_priority_queue::print() { - std::vector index; - std::vector 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 index; + std::vector 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 } diff --git a/src/util/lp/binary_heap_upair_queue.cpp b/src/util/lp/binary_heap_upair_queue.cpp index 2c0e34c47..8cf51c79a 100644 --- a/src/util/lp/binary_heap_upair_queue.cpp +++ b/src/util/lp/binary_heap_upair_queue.cpp @@ -9,242 +9,108 @@ #include "util/lp/binary_heap_upair_queue.h" namespace lean { template binary_heap_upair_queue::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 unsigned - binary_heap_upair_queue::dequeue_available_spot() { +binary_heap_upair_queue::dequeue_available_spot() { lean_assert(m_available_spots.empty() == false); unsigned ret = m_available_spots.back(); m_available_spots.pop_back(); return ret; } + template void binary_heap_upair_queue::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 bool binary_heap_upair_queue::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 void binary_heap_upair_queue::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_indexsecond; + 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_indexsecond; } + m_q.enqueue(ij_index, priority); +} template void binary_heap_upair_queue::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 T binary_heap_upair_queue::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 bool binary_heap_upair_queue::pair_to_index_is_a_bijection() const { - std::set 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 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 bool binary_heap_upair_queue::available_spots_are_correct() const { - std::set 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 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 -#include -#include -#include -#include -typedef std::pair upair; - -namespace lean { -template -class binary_heap_upair_queue { - binary_heap_priority_queue m_q; - std::unordered_map m_pairs_to_index; - std::vector m_pairs; // inverse to index - std::vector 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_indexsecond; - } - 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 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 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(); - } -}; -} -*/ diff --git a/src/util/lp/core_solver_pretty_printer.cpp b/src/util/lp/core_solver_pretty_printer.cpp index 6baa16400..865979ac1 100644 --- a/src/util/lp/core_solver_pretty_printer.cpp +++ b/src/util/lp/core_solver_pretty_printer.cpp @@ -11,329 +11,329 @@ namespace lean { template core_solver_pretty_printer::core_solver_pretty_printer(lp_core_solver_base & 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(core_solver.m_A.column_count(), "")), - m_signs(core_solver.m_A.row_count(), vector(core_solver.m_A.column_count(), " ")), - m_costs(ncols(), ""), - m_cost_signs(ncols(), " "), - m_rs(ncols(), zero_of_type()) { - 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(core_solver.m_A.column_count(), "")), + m_signs(core_solver.m_A.row_count(), vector(core_solver.m_A.column_count(), " ")), + m_costs(ncols(), ""), + m_cost_signs(ncols(), " "), + m_rs(ncols(), zero_of_type()) { + 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 void core_solver_pretty_printer::init_costs() { - vector 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 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 core_solver_pretty_printer::~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 void core_solver_pretty_printer::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 T core_solver_pretty_printer::current_column_norm() { - T ret = zero_of_type(); - for (T & ed : m_core_solver.m_ed) - ret += ed * ed; - return ret; - } + T ret = zero_of_type(); + for (T & ed : m_core_solver.m_ed) + ret += ed * ed; + return ret; +} template void core_solver_pretty_printer::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 void core_solver_pretty_printer::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 void core_solver_pretty_printer::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 void core_solver_pretty_printer::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 void core_solver_pretty_printer::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 unsigned core_solver_pretty_printer:: 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 std::string core_solver_pretty_printer::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 void core_solver_pretty_printer::set_coeff(vector& row, vector & row_signs, unsigned col, const T & t, string name) { - if (numeric_traits::is_zero(t)) { - return; + if (numeric_traits::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 void core_solver_pretty_printer::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 std::string core_solver_pretty_printer::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 std::string core_solver_pretty_printer::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 void core_solver_pretty_printer::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 void core_solver_pretty_printer::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 void core_solver_pretty_printer::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 void core_solver_pretty_printer::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 void core_solver_pretty_printer::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 void core_solver_pretty_printer::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 void core_solver_pretty_printer::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 void core_solver_pretty_printer::print_given_rows(vector & row, vector & 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 void core_solver_pretty_printer::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); +} } diff --git a/src/util/lp/dense_matrix.cpp b/src/util/lp/dense_matrix.cpp index 8d1b0eb3f..bcb67edc1 100644 --- a/src/util/lp/dense_matrix.cpp +++ b/src/util/lp/dense_matrix.cpp @@ -8,193 +8,193 @@ #include "util/lp/dense_matrix.h" namespace lean { template dense_matrix::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::zero(); - } + m_values = new T[m * n]; + for (unsigned i = 0; i < m * n; i ++) + m_values[i] = numeric_traits::zero(); +} template dense_matrix dense_matrix::operator*=(matrix 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::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::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 dense_matrix& - dense_matrix::operator=(matrix 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::operator=(matrix 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 dense_matrix& - dense_matrix::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::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 dense_matrix::dense_matrix(dense_matrix 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 dense_matrix::dense_matrix(matrix 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 void dense_matrix::apply_from_right(T * w) { - T * t = new T[m_m]; - for (int i = 0; i < m_m; i ++) { - T v = numeric_traits::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::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 void dense_matrix::apply_from_right(std::vector & w) { - T * t = new T[m_m]; - for (int i = 0; i < m_m; i ++) { - T v = numeric_traits::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::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 T* dense_matrix:: - apply_from_left_with_different_dims(std::vector & w) { - T * t = new T[m_m]; - for (int i = 0; i < m_m; i ++) { - T v = numeric_traits::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 & w) { + T * t = new T[m_m]; + for (int i = 0; i < m_m; i ++) { + T v = numeric_traits::zero(); + for (int j = 0; j < m_n; j++) { + v += w[j]* get_elem(i, j); } - - return t; + t[i] = v; } + return t; +} + template void dense_matrix::apply_from_left(std::vector & w) { - T * t = new T[m_m]; - for (int i = 0; i < m_m; i ++) { - T v = numeric_traits::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::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 void dense_matrix::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::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::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 void dense_matrix::apply_from_left_to_X(std::vector & w, lp_settings & ) { - std::vector t(m_m); - for (int i = 0; i < m_m; i ++) { - X v = zero_of_type(); - 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 t(m_m); + for (int i = 0; i < m_m; i ++) { + X v = zero_of_type(); + 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 void dense_matrix::pivot_row_to_row(unsigned i, T alpha, unsigned i0, - double & pivot_epsilon) { - thread_local T _0 = numeric_traits::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::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 void dense_matrix::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 void dense_matrix::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 void dense_matrix::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 dense_matrix operator* (matrix & a, matrix & b){ diff --git a/src/util/lp/indexed_vector.cpp b/src/util/lp/indexed_vector.cpp index fb5666939..eccd782a0 100644 --- a/src/util/lp/indexed_vector.cpp +++ b/src/util/lp/indexed_vector.cpp @@ -38,57 +38,57 @@ void print_vector(const std::vector & t) { } template - void indexed_vector::resize(unsigned data_size) { - m_index.clear(); - m_data.resize(data_size, numeric_traits::zero()); - } +void indexed_vector::resize(unsigned data_size) { + m_index.clear(); + m_data.resize(data_size, numeric_traits::zero()); +} template - void indexed_vector::set_value(T value, unsigned index) { - m_data[index] = value; - m_index.push_back(index); - } +void indexed_vector::set_value(T value, unsigned index) { + m_data[index] = value; + m_index.push_back(index); +} template - void indexed_vector::clear() { - for (unsigned i : m_index) - m_data[i] = numeric_traits::zero(); - m_index.clear(); - } +void indexed_vector::clear() { + for (unsigned i : m_index) + m_data[i] = numeric_traits::zero(); + m_index.clear(); +} template - void indexed_vector::clear_all() { - unsigned i = m_data.size(); - while (i--) m_data[i] = numeric_traits::zero(); - m_index.clear(); - } +void indexed_vector::clear_all() { + unsigned i = m_data.size(); + while (i--) m_data[i] = numeric_traits::zero(); + m_index.clear(); +} template - void indexed_vector::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::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 - bool indexed_vector::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::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 - void indexed_vector::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::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 } diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 30b60425a..c234e2430 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -17,7 +17,7 @@ double conversion_helper ::get_low_bound(const column_info & ci) { return ci.get_low_bound().get_double() + eps; } - double conversion_helper ::get_upper_bound(const column_info & ci) { +double conversion_helper ::get_upper_bound(const column_info & ci) { if (!ci.upper_bound_is_strict()) return ci.get_upper_bound().get_double(); double eps = 0.00001; @@ -285,11 +285,11 @@ template V lar_solver::get_column_val(std::vector & low_bound, s } lar_solver::~lar_solver() { - std::vector to_delete; - for (auto it : m_canonic_left_sides) - to_delete.push_back(it); - for (auto t : to_delete) - delete t; + std::vector to_delete; + for (auto it : m_canonic_left_sides) + to_delete.push_back(it); + for (auto t : to_delete) + delete t; } diff --git a/src/util/lp/lp_dual_simplex.cpp b/src/util/lp/lp_dual_simplex.cpp index c90d325e3..e6b572e4a 100644 --- a/src/util/lp/lp_dual_simplex.cpp +++ b/src/util/lp/lp_dual_simplex.cpp @@ -257,8 +257,8 @@ template void lp_dual_simplex::fill_costs_and_bou } template void lp_dual_simplex::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 diff --git a/src/util/lp/lp_primal_core_solver.cpp b/src/util/lp/lp_primal_core_solver.cpp index 44ccf0c6d..7ae9540ea 100644 --- a/src/util/lp/lp_primal_core_solver.cpp +++ b/src/util/lp/lp_primal_core_solver.cpp @@ -210,15 +210,15 @@ template X lp_primal_core_solver::get_max_boun // stage1 constructor template lp_primal_core_solver:: lp_primal_core_solver(static_matrix & A, - std::vector & b, // the right side vector - std::vector & x, // the number of elements in x needs to be at least as large as the number of columns in A - std::vector & basis, - std::vector & costs, - std::vector & column_type_array, - std::vector & low_bound_values, - std::vector & upper_bound_values, - lp_settings & settings, - std::unordered_map const & column_names): + std::vector & b, // the right side vector + std::vector & x, // the number of elements in x needs to be at least as large as the number of columns in A + std::vector & basis, + std::vector & costs, + std::vector & column_type_array, + std::vector & low_bound_values, + std::vector & upper_bound_values, + lp_settings & settings, + std::unordered_map const & column_names): lp_core_solver_base(A, b, basis, x, diff --git a/src/util/lp/lp_primal_core_solver_instances.cpp b/src/util/lp/lp_primal_core_solver_instances.cpp index 90eef637d..ef5afc832 100644 --- a/src/util/lp/lp_primal_core_solver_instances.cpp +++ b/src/util/lp/lp_primal_core_solver_instances.cpp @@ -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::find_feasible_solution(); template lp_primal_core_solver::lp_primal_core_solver(static_matrix&, std::vector >&, std::vector >&, std::vector >&, std::vector >&, std::vector >&, std::vector >&, lp_settings&, std::unordered_map, std::equal_to, std::allocator > > const&); template lp_primal_core_solver::lp_primal_core_solver(static_matrix&, std::vector >&, std::vector >&, std::vector >&, std::vector >&, std::vector >&, std::vector >&, std::vector >&, lp_settings&, std::unordered_map, std::equal_to, std::allocator > > const&); @@ -15,6 +15,6 @@ template lp_primal_core_solver::lp_primal_core_solver(static_mat template unsigned lp_primal_core_solver::solve(); template lp_primal_core_solver::lp_primal_core_solver(static_matrix&, std::vector >&, std::vector >&, std::vector >&, std::vector >&, std::vector >&, std::vector >&, std::vector >&, lp_settings&, std::unordered_map, std::equal_to, std::allocator > > const&); template unsigned lp_primal_core_solver::solve(); - //template void lp_primal_simplex::solve_with_total_inf(); - //template void lp_primal_simplex::solve_with_total_inf(); +//template void lp_primal_simplex::solve_with_total_inf(); +//template void lp_primal_simplex::solve_with_total_inf(); } diff --git a/src/util/lp/lp_primal_simplex.cpp b/src/util/lp/lp_primal_simplex.cpp index 0ec7f0930..440c46808 100644 --- a/src/util/lp/lp_primal_simplex.cpp +++ b/src/util/lp/lp_primal_simplex.cpp @@ -76,9 +76,9 @@ template void lp_primal_simplex::set_zero_bound(b } template void lp_primal_simplex::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 diff --git a/src/util/lp/lp_settings.cpp b/src/util/lp/lp_settings.cpp index ef1c13421..342ec65e6 100644 --- a/src/util/lp/lp_settings.cpp +++ b/src/util/lp/lp_settings.cpp @@ -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 diff --git a/src/util/lp/matrix.cpp b/src/util/lp/matrix.cpp index a27ed0eae..3656b7103 100644 --- a/src/util/lp/matrix.cpp +++ b/src/util/lp/matrix.cpp @@ -9,26 +9,26 @@ namespace lean { template bool matrix::is_equal(const matrix& 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::precise()) { - if (a != b) return false; - } else if (fabs(numeric_traits::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::precise()) { + if (a != b) return false; + } else if (fabs(numeric_traits::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 void apply_to_vector(matrix & m, T * w) { diff --git a/src/util/lp/row_eta_matrix.cpp b/src/util/lp/row_eta_matrix.cpp index 1968a4108..e84ce96e0 100644 --- a/src/util/lp/row_eta_matrix.cpp +++ b/src/util/lp/row_eta_matrix.cpp @@ -9,146 +9,146 @@ namespace lean { template void row_eta_matrix::apply_from_left(std::vector & w, lp_settings &) { -// #ifdef LEAN_DEBUG -// dense_matrix deb(*this); -// auto clone_w = clone_vector(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 deb(*this); + // auto clone_w = clone_vector(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(clone_w, w, m_dimension)); + // delete [] clone_w; + // #endif +} + +template +void row_eta_matrix::apply_from_left_local_to_T(indexed_vector & 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(clone_w, w, m_dimension)); -// delete [] clone_w; -// #endif + } else if (!was_zero_at_m_row){ + w[m_row] = zero_of_type(); + auto it = std::find(w.m_index.begin(), w.m_index.end(), m_row); + w.m_index.erase(it); } - -template - void row_eta_matrix::apply_from_left_local_to_T(indexed_vector & 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(); - 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 - void row_eta_matrix::apply_from_left_local_to_X(indexed_vector & 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::apply_from_left_local_to_X(indexed_vector & 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(); - 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(); + 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 - void row_eta_matrix::apply_from_right(std::vector & w) { - T w_row = w[m_row]; - if (numeric_traits::is_zero(w_row)) return; +void row_eta_matrix::apply_from_right(std::vector & w) { + T w_row = w[m_row]; + if (numeric_traits::is_zero(w_row)) return; #ifdef LEAN_DEBUG - // dense_matrix deb(*this); - // auto clone_w = clone_vector(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(clone_w, w, m_dimension)); - // delete clone_w; -#endif - } - -template - void row_eta_matrix::apply_from_right(indexed_vector & w) { - T w_row = w[m_row]; - if (numeric_traits::is_zero(w_row)) return; -#ifdef LEAN_DEBUG - // dense_matrix deb(*this); - // auto clone_w = clone_vector(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::is_zero(old_val)) { - w.m_index.push_back(it.index()); - } else if (numeric_traits::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(clone_w, w.m_data, m_dimension)); - // for (unsigned i = 0; i < m_dimension; i++) { - // if (!numeric_traits::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 - void row_eta_matrix::conjugate_by_permutation(permutation_matrix & 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 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 deb(*this); + // auto clone_w = clone_vector(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 - T row_eta_matrix::get_elem(unsigned row, unsigned col) const { - if (row == m_row){ - if (col == row) { - return numeric_traits::one(); - } - return m_row_vector[col]; - } - - return col == row ? numeric_traits::one() : numeric_traits::zero(); - } + // lean_assert(vectors_are_equal(clone_w, w, m_dimension)); + // delete clone_w; +#endif +} + +template +void row_eta_matrix::apply_from_right(indexed_vector & w) { + T w_row = w[m_row]; + if (numeric_traits::is_zero(w_row)) return; +#ifdef LEAN_DEBUG + // dense_matrix deb(*this); + // auto clone_w = clone_vector(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::is_zero(old_val)) { + w.m_index.push_back(it.index()); + } else if (numeric_traits::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(clone_w, w.m_data, m_dimension)); + // for (unsigned i = 0; i < m_dimension; i++) { + // if (!numeric_traits::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 +void row_eta_matrix::conjugate_by_permutation(permutation_matrix & 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 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 +T row_eta_matrix::get_elem(unsigned row, unsigned col) const { + if (row == m_row){ + if (col == row) { + return numeric_traits::one(); + } + return m_row_vector[col]; + } + + return col == row ? numeric_traits::one() : numeric_traits::zero(); +} #endif } diff --git a/src/util/lp/sparse_matrix_instances.cpp b/src/util/lp/sparse_matrix_instances.cpp index 67dc55deb..49ec3db7f 100644 --- a/src/util/lp/sparse_matrix_instances.cpp +++ b/src/util/lp/sparse_matrix_instances.cpp @@ -72,23 +72,23 @@ template void sparse_matrix::double_solve_U_y(std::vector&); template void sparse_matrix>::double_solve_U_y(std::vector&); template void sparse_matrix >::double_solve_U_y >(std::vector>&); - ///////////////////// - /* -template void lu::create_initial_factorization(); -template void lu::create_initial_factorization(); -template void lu::replace_column(unsigned int, mpq, indexed_vector&); -template void lu >::create_initial_factorization(); -template void lu >::replace_column(unsigned int, mpq, indexed_vector&); -template void lu::init_vector_w(unsigned int, indexed_vector&); -template void lu::find_error_of_yB(std::vector >&, std::vector > const&); -template void lu::init_vector_w(unsigned int, indexed_vector&); -template void lu >::find_error_of_yB(std::vector >&, std::vector > const&); -template void lu >::init_vector_w(unsigned int, indexed_vector&); -template void lu::find_error_of_yB(std::vector >&, std::vector > const&); -#ifdef LEAN_DEBUG -template void print_matrix(static_matrix&); -#endif - */ +///////////////////// +/* + template void lu::create_initial_factorization(); + template void lu::create_initial_factorization(); + template void lu::replace_column(unsigned int, mpq, indexed_vector&); + template void lu >::create_initial_factorization(); + template void lu >::replace_column(unsigned int, mpq, indexed_vector&); + template void lu::init_vector_w(unsigned int, indexed_vector&); + template void lu::find_error_of_yB(std::vector >&, std::vector > const&); + template void lu::init_vector_w(unsigned int, indexed_vector&); + template void lu >::find_error_of_yB(std::vector >&, std::vector > const&); + template void lu >::init_vector_w(unsigned int, indexed_vector&); + template void lu::find_error_of_yB(std::vector >&, std::vector > const&); + #ifdef LEAN_DEBUG + template void print_matrix(static_matrix&); + #endif +*/ #ifdef LEAN_DEBUG template bool sparse_matrix::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const; template bool sparse_matrix::is_upper_triangular_and_maximums_are_set_correctly_in_rows(lp_settings&) const; diff --git a/src/util/lp/static_matrix.cpp b/src/util/lp/static_matrix.cpp index d4bdd4c34..2ca7578eb 100644 --- a/src/util/lp/static_matrix.cpp +++ b/src/util/lp/static_matrix.cpp @@ -20,7 +20,7 @@ void static_matrix::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 static_matrix::static_matrix(static_matrix const &A, unsigned * basis) : m_work_pivot_vector(A.row_count(), numeric_traits::zero()) { @@ -52,7 +52,7 @@ template void static_matrix::init_empty_matrix } template template - L static_matrix::dot_product_with_row(unsigned row, const std::vector & w) { +L static_matrix::dot_product_with_row(unsigned row, const std::vector & w) { L ret = zero_of_type(); lean_assert(row < m_rows.size()); for (auto & it : m_rows[row]) { @@ -107,7 +107,7 @@ template void static_matrix::remove_last_colum break; } offset--; - } + } } m_columns.pop_back(); } @@ -144,7 +144,7 @@ template void static_matrix::regen_domain() { for (auto & t : m_rows[i]) { m_domain.insert(std::make_pair(i, t.m_j)); } - } + } } #endif @@ -183,12 +183,12 @@ template void static_matrix::copy_column_to_ve } template void static_matrix::copy_column_to_vector (unsigned j, std::vector & v) const { - v.resize(row_count(), numeric_traits::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::zero()); + for (auto & it : m_columns[j]) { + if (!is_zero(it.m_value)) + v[it.m_i]=it.m_value; } +} template void static_matrix::add_column_to_vector (const T & a, unsigned j, T * v) const { for (auto & it : m_columns[j]) {