Program Listing for File rowmod.hpp
↰ Return to documentation for file (/tmp/ws/src/proxsuite/include/proxsuite/linalg/sparse/rowmod.hpp
)
//
// Copyright (c) 2022 INRIA
//
#ifndef PROXSUITE_LINALG_SPARSE_LDLT_ROWMOD_HPP
#define PROXSUITE_LINALG_SPARSE_LDLT_ROWMOD_HPP
#include "proxsuite/linalg/sparse/update.hpp"
#include <algorithm>
namespace proxsuite {
namespace linalg {
namespace sparse {
template<typename T, typename I>
auto
delete_row_req( //
proxsuite::linalg::veg::Tag<T> /*tag*/,
proxsuite::linalg::veg::Tag<I> /*tag*/,
isize n,
isize max_nnz) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
{
return sparse::rank1_update_req(proxsuite::linalg::veg::Tag<T>{},
proxsuite::linalg::veg::Tag<I>{},
n,
true,
max_nnz);
}
template<typename T, typename I>
auto
delete_row(MatMut<T, I> ld,
I* etree,
I const* perm_inv,
isize pos,
DynStackMut stack) noexcept(false) -> MatMut<T, I>
{
// step 1: delete row k from each column
VEG_ASSERT(!ld.is_compressed());
// we're actually deleting perm_inv[k], so that k is deleted in the permuted
// matrix
usize permuted_pos =
perm_inv == nullptr ? usize(pos) : util::zero_extend(perm_inv[pos]);
auto petree = etree;
I* pldi = ld.row_indices_mut();
T* pldx = ld.values_mut();
I* pldnz = ld.nnz_per_col_mut();
for (usize j = 0; j < permuted_pos; ++j) {
auto col_start = ld.col_start(j) + 1;
auto col_end = ld.col_end(j);
// search for the first row in column j greater than or equal to k
auto it =
std::lower_bound(pldi + col_start, pldi + col_end, I(permuted_pos));
// if an element was found, and it is equal to k
if ((it != (pldi + col_end)) && *it == I(permuted_pos)) {
usize it_pos = usize(it - (pldi + col_start));
usize count = (col_end - col_start - it_pos);
// shift all the row indices back by one position
// to delete row k
std::memmove(it, it + 1, count * sizeof(I));
T* itx = pldx + col_start + it_pos;
VEG_CHECK_CONCEPT(trivially_copyable<T>);
// shift all the values back by one position
std::memmove(itx, itx + 1, count * sizeof(T));
// decrement the non zero count
--pldnz[j];
ld._set_nnz(ld.nnz() - 1);
// adjust the parent of j in the elimination tree if necessary
if (petree[j] == I(permuted_pos)) {
VEG_ASSERT(it_pos == 0);
if (pldnz[j] > 1) {
petree[j] = *it;
} else {
petree[j] = I(-1);
}
}
}
}
// step 2: set d_kk = 1
T d_old = ld.values()[ld.col_start(permuted_pos)];
ld.values_mut()[ld.col_start(permuted_pos)] = 1;
// step 3: perform rank update
isize len = isize(util::zero_extend(ld.nnz_per_col()[permuted_pos])) - 1;
ld = sparse::rank1_update<T, I>( //
ld,
etree,
static_cast<I const*>(nullptr),
VecRef<T, I>{
from_raw_parts,
ld.nrows(),
len,
pldi + ld.col_start(permuted_pos) + 1,
pldx + ld.col_start(permuted_pos) + 1,
},
d_old,
stack);
// step 4: delete col k_
ld.nnz_per_col_mut()[permuted_pos] = 1;
petree[permuted_pos] = I(-1);
return ld;
}
template<typename T, typename I>
auto
add_row_req( //
proxsuite::linalg::veg::Tag<T> /*tag*/,
proxsuite::linalg::veg::Tag<I> /*tag*/,
isize n,
bool id_perm,
isize nnz,
isize max_nnz) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
{
using proxsuite::linalg::veg::dynstack::StackReq;
auto numerical_work = StackReq{ n * isize{ sizeof(T) }, isize{ alignof(T) } };
auto permuted_indices =
StackReq{ (id_perm ? 0 : nnz) * isize{ sizeof(I) }, isize{ alignof(I) } };
auto pattern_diff = StackReq{ n * isize{ sizeof(I) }, isize{ alignof(I) } };
auto merge =
merge_second_col_into_first_req(proxsuite::linalg::veg::Tag<I>{}, n);
auto update = sparse::rank1_update_req(proxsuite::linalg::veg::Tag<T>{},
proxsuite::linalg::veg::Tag<I>{},
n,
true,
max_nnz);
auto req = numerical_work;
req = req & permuted_indices;
req = req & pattern_diff;
req = req & merge;
req = req | update;
return req;
}
template<typename T, typename I>
auto
add_row(MatMut<T, I> ld,
I* etree,
I const* perm_inv,
isize pos,
VecRef<T, I> new_col,
proxsuite::linalg::veg::DoNotDeduce<T> diag_element,
DynStackMut stack) noexcept(false) -> MatMut<T, I>
{
VEG_ASSERT(!ld.is_compressed());
bool id_perm = perm_inv == nullptr;
auto zx = util::zero_extend;
I* pldp = ld.col_ptrs_mut();
I* pldnz = ld.nnz_per_col_mut();
I* pldi = ld.row_indices_mut();
T* pldx = ld.values_mut();
// actually inserting in the position perm_inv[k] so that row k is added in
// the permuted matrix
usize permuted_pos = id_perm ? usize(pos) : zx(perm_inv[pos]);
VEG_ASSERT(pldnz[permuted_pos] == 1);
{
// allocate workspace for numerical step, storage for the k-th row and k-th
// column of the new matrix
auto _lx2_storage = stack.make_new_for_overwrite(
proxsuite::linalg::veg::Tag<T>{}, ld.nrows());
auto plx2_storage = _lx2_storage.ptr_mut();
// allocate workspace for permuted row indices of the new column if
// necessary
auto _new_col_permuted_indices = stack.make_new_for_overwrite(
proxsuite::linalg::veg::Tag<I>{}, id_perm ? isize(0) : new_col.nnz());
auto new_col_permuted_indices =
id_perm ? new_col.row_indices() : _new_col_permuted_indices.ptr();
// copy and sort permuted row indices
if (!id_perm) {
I* pnew_col_permuted_indices = _new_col_permuted_indices.ptr_mut();
for (usize k = 0; k < usize(new_col.nnz()); ++k) {
usize i = zx(new_col.row_indices()[k]);
pnew_col_permuted_indices[k] = perm_inv[i];
}
std::sort(pnew_col_permuted_indices,
pnew_col_permuted_indices + new_col.nnz());
}
// allocate workspace for non-zero pattern of k-th row
auto _l12_nnz_pattern = stack.make_new_for_overwrite(
proxsuite::linalg::veg::Tag<I>{}, isize(permuted_pos));
auto _difference = stack.make_new_for_overwrite(
proxsuite::linalg::veg::Tag<I>{}, ld.nrows() - isize(permuted_pos));
auto pdifference = _difference.ptr_mut();
auto pl12_nnz_pattern = _l12_nnz_pattern.ptr_mut();
usize l12_nnz_pattern_count = 0;
// the non-zero pattern is the set of columns reachable from the non-zero
// pattern of the added column through graph of L_{1..k,1..k}
// instead of graph traversal, we can use the k-th elimination subtree as we
// did in the initial factorization step
// for each row in the added column
{
auto _visited = stack.make_new(proxsuite::linalg::veg::Tag<bool>{},
isize(permuted_pos));
bool* visited = _visited.ptr_mut();
for (usize p = 0; p < usize(new_col.nnz()); ++p) {
auto j = zx(new_col_permuted_indices[p]);
if (j >= permuted_pos) {
break;
}
// add the ancestors of the corresponding column
// ancestors are not sorted, but they are added in topological order,
// which suffices for the triangular solve
while (true) {
if (visited[j]) {
break;
}
visited[j] = true;
pl12_nnz_pattern[l12_nnz_pattern_count] = I(j);
++l12_nnz_pattern_count;
j = util::sign_extend(etree[j]);
if (j == usize(-1) || j >= permuted_pos || visited[j]) {
break;
}
}
}
}
std::sort(pl12_nnz_pattern, pl12_nnz_pattern + l12_nnz_pattern_count);
// zero the elements in the non-zero pattern of the solution (new k-th row)
for (usize p = 0; p < l12_nnz_pattern_count; ++p) {
plx2_storage[zx(pl12_nnz_pattern[p])] = 0;
}
// insert the rhs of the k-th row triangular system in the top part of the
// storage, and the bottom part of the added column in the bottom part of
// the storage
for (usize p = 0; p < usize(new_col.nnz()); ++p) {
auto j = zx(new_col.row_indices()[p]);
auto permuted_j = id_perm ? j : zx(perm_inv[j]);
plx2_storage[permuted_j] = new_col.values()[p];
// add the row indices of the bottom part of the added column, to the
// k-th column of L
if (permuted_j > permuted_pos) {
usize nz = zx(pldnz[permuted_pos]);
VEG_ASSERT(nz < (zx(pldp[permuted_pos + 1]) - zx(pldp[permuted_pos])));
pldi[zx(pldp[permuted_pos]) + nz] = I(permuted_j);
++pldnz[permuted_pos];
ld._set_nnz(ld.nnz() + 1);
}
}
// sort the added row indices
std::sort(pldi + zx(pldp[permuted_pos]) + 1,
pldi + zx(pldp[permuted_pos]) + zx(pldnz[permuted_pos]));
// TODO: fuse loops?
for (usize p = 0; p < l12_nnz_pattern_count; ++p) {
usize j = zx(pl12_nnz_pattern[p]);
auto col_start = ld.col_start(j);
auto col_end = ld.col_end(j);
// update the pattern of the k-th column of L, with that of the bottom
// part of the j-th column of L, ignoring the elements less than or equal
// to k
VEG_BIND(auto,
(_, new_current_col, computed_difference),
sparse::merge_second_col_into_first(
pdifference,
static_cast<T*>(nullptr),
pldi + (zx(pldp[permuted_pos]) + 1),
isize(zx(pldp[permuted_pos + 1]) - zx(pldp[permuted_pos])) - 1,
pldnz[permuted_pos] - 1,
{
unsafe,
from_raw_parts,
pldi + (zx(pldp[j]) + 1),
isize(zx(pldnz[j])) - 1,
},
I(permuted_pos),
false,
stack));
(void)_;
(void)new_current_col;
// update column and global non-zero count
pldnz[permuted_pos] += I(computed_difference.len());
ld._set_nnz(ld.nnz() + computed_difference.len());
for (usize q = 0; q < usize(computed_difference.len()); ++q) {
plx2_storage[zx(computed_difference.ptr()[q])] = 0;
}
// perform triangular solve and matrix vector product simultaneously
auto const xj = plx2_storage[j];
for (usize q = col_start + 1; q < col_end; ++q) {
auto i = zx(pldi[q]);
plx2_storage[i] -= pldx[q] * xj;
}
}
// insert the k-th row into L
for (usize p = 0; p < l12_nnz_pattern_count; ++p) {
// for each column in the non-zero pattern of the k-th row
usize j = zx(pl12_nnz_pattern[p]);
auto col_start = ld.col_start(j);
auto col_end = ld.col_end(j);
T d = pldx[col_start];
T l12_elem = plx2_storage[j];
diag_element -= l12_elem * l12_elem / d;
// check that we have enough space to insert one element
VEG_ASSERT(zx(pldnz[j]) < (zx(pldp[j + 1]) - zx(pldp[j])));
// find the first element greater than k
auto it =
std::lower_bound(pldi + col_start, pldi + col_end, I(permuted_pos));
// if it is the first element, update the elimination tree so that k is
// the new parent of column j
if (it == (pldi + col_start + 1)) {
etree[j] = I(permuted_pos);
}
// shift the row indices up by one position to provide enough space for
// the new element
std::memmove( //
it + 1,
it,
usize((pldi + col_end) - it) * sizeof(I));
VEG_CHECK_CONCEPT(trivially_copyable<T>);
// shift the values up by one position to provide enough space for the
// new element
std::memmove( //
pldx + (it - pldi) + 1,
pldx + (it - pldi),
usize((pldi + col_end) - it) * sizeof(T));
// insert the new row index k
*it = I(permuted_pos);
// insert the new corresponding value
*(pldx + (it - pldi)) = l12_elem / d;
// update the non-zero count
++pldnz[j];
ld._set_nnz(ld.nnz() + 1);
}
// insert the k-th column of L
{
usize col_start = ld.col_start(permuted_pos);
usize col_end = ld.col_end(permuted_pos);
pldx[col_start] = diag_element;
for (usize p = col_start + 1; p < col_end; ++p) {
pldx[p] = plx2_storage[zx(pldi[p])] / diag_element;
}
}
}
// set the parent of the k-th column of L
if (pldnz[permuted_pos] > 1) {
etree[permuted_pos] = pldi[ld.col_start(permuted_pos) + 1];
}
isize len = isize(util::zero_extend(ld.nnz_per_col()[permuted_pos])) - 1;
// perform the rank update with the newly added column
ld = sparse::rank1_update<T, I>(ld,
etree,
static_cast<I const*>(nullptr),
VecRef<T, I>{
from_raw_parts,
ld.nrows(),
len,
pldi + ld.col_start(permuted_pos) + 1,
pldx + ld.col_start(permuted_pos) + 1,
},
-diag_element,
stack);
return ld;
}
} // namespace sparse
} // namespace linalg
} // namespace proxsuite
#endif /* end of include guard PROXSUITE_LINALG_SPARSE_LDLT_ROWMOD_HPP */