.. _program_listing_file__tmp_ws_src_proxsuite_include_proxsuite_linalg_dense_ldlt.hpp:

Program Listing for File ldlt.hpp
=================================

|exhale_lsh| :ref:`Return to documentation for file <file__tmp_ws_src_proxsuite_include_proxsuite_linalg_dense_ldlt.hpp>` (``/tmp/ws/src/proxsuite/include/proxsuite/linalg/dense/ldlt.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   //
   // Copyright (c) 2022 INRIA
   //
   #ifndef PROXSUITE_LINALG_DENSE_LDLT_LDLT_HPP
   #define PROXSUITE_LINALG_DENSE_LDLT_LDLT_HPP
   
   #include "proxsuite/linalg/dense/factorize.hpp"
   #include "proxsuite/linalg/dense/update.hpp"
   #include "proxsuite/linalg/dense/modify.hpp"
   #include "proxsuite/linalg/dense/solve.hpp"
   #include <proxsuite/linalg/veg/vec.hpp>
   
   namespace proxsuite {
   namespace linalg {
   namespace dense {
   namespace _detail {
   struct SimdAlignedSystemAlloc
   {
     friend auto operator==(SimdAlignedSystemAlloc /*unused*/,
                            SimdAlignedSystemAlloc /*unused*/) noexcept -> bool
     {
       return true;
     }
   };
   } // namespace _detail
   } // namespace dense
   } // namespace linalg
   } // namespace proxsuite
   
   template<>
   struct proxsuite::linalg::veg::mem::Alloc<
     proxsuite::linalg::dense::_detail::SimdAlignedSystemAlloc>
   {
   #ifdef PROXSUITE_VECTORIZE
     static constexpr usize min_align = alignof(std::max_align_t) >
                                            SIMDE_NATURAL_VECTOR_SIZE / 8
                                          ? SIMDE_NATURAL_VECTOR_SIZE / 8
                                          : alignof(std::max_align_t);
   #else
     static constexpr usize min_align = 0;
   #endif
   
     using RefMut = proxsuite::linalg::veg::RefMut<
       proxsuite::linalg::dense::_detail::SimdAlignedSystemAlloc>;
   
     VEG_INLINE static auto adjusted_layout(Layout l) noexcept -> Layout
     {
       if (l.align < min_align) {
         l.align = min_align;
       }
       return l;
     }
   
     VEG_INLINE static void dealloc(RefMut /*alloc*/, void* ptr, Layout l) noexcept
     {
       return Alloc<SystemAlloc>::dealloc(
         mut(SystemAlloc{}), ptr, adjusted_layout(l));
     }
   
     VEG_NODISCARD VEG_INLINE static auto alloc(RefMut /*alloc*/,
                                                Layout l) noexcept
       -> mem::AllocBlock
     {
       return Alloc<SystemAlloc>::alloc(mut(SystemAlloc{}), adjusted_layout(l));
     }
   
     VEG_NODISCARD VEG_INLINE static auto grow(RefMut /*alloc*/,
                                               void* ptr,
                                               Layout l,
                                               usize new_size,
                                               RelocFn reloc) noexcept
       -> mem::AllocBlock
     {
       return Alloc<SystemAlloc>::grow(
         mut(SystemAlloc{}), ptr, adjusted_layout(l), new_size, reloc);
     }
     VEG_NODISCARD VEG_INLINE static auto shrink(RefMut /*alloc*/,
                                                 void* ptr,
                                                 Layout l,
                                                 usize new_size,
                                                 RelocFn reloc) noexcept
       -> mem::AllocBlock
     {
       return Alloc<SystemAlloc>::shrink(
         mut(SystemAlloc{}), ptr, adjusted_layout(l), new_size, reloc);
     }
   };
   
   namespace proxsuite {
   namespace linalg {
   namespace dense {
   template<typename T>
   struct Ldlt
   {
   private:
     static constexpr auto DYN = Eigen::Dynamic;
     using ColMat = Eigen::Matrix<T, DYN, DYN, Eigen::ColMajor>;
     using RowMat = Eigen::Matrix<T, DYN, DYN, Eigen::RowMajor>;
     using Vec = Eigen::Matrix<T, DYN, 1>;
   
     using LView = Eigen::TriangularView<Eigen::Map< //
                                           ColMat const,
                                           Eigen::Unaligned,
                                           Eigen::OuterStride<DYN>>,
                                         Eigen::UnitLower>;
     using LViewMut = Eigen::TriangularView<Eigen::Map< //
                                              ColMat,
                                              Eigen::Unaligned,
                                              Eigen::OuterStride<DYN>>,
                                            Eigen::UnitLower>;
   
     using LTView = Eigen::TriangularView<Eigen::Map< //
                                            RowMat const,
                                            Eigen::Unaligned,
                                            Eigen::OuterStride<DYN>>,
                                          Eigen::UnitUpper>;
     using LTViewMut = Eigen::TriangularView<Eigen::Map< //
                                               RowMat,
                                               Eigen::Unaligned,
                                               Eigen::OuterStride<DYN>>,
                                             Eigen::UnitUpper>;
   
     using DView =
       Eigen::Map<Vec const, Eigen::Unaligned, Eigen::InnerStride<DYN>>;
     using DViewMut = Eigen::Map<Vec, Eigen::Unaligned, Eigen::InnerStride<DYN>>;
   
     using VecMapISize = Eigen::Map<Eigen::Matrix<isize, DYN, 1> const>;
     using Perm = Eigen::PermutationWrapper<VecMapISize>;
   
     using StorageSimdVec =
       proxsuite::linalg::veg::Vec<T,
                                   proxsuite::linalg::veg::meta::if_t<
                                     _detail::should_vectorize<T>::value,
                                     _detail::SimdAlignedSystemAlloc,
                                     proxsuite::linalg::veg::mem::SystemAlloc>>;
   
     StorageSimdVec ld_storage;
     isize stride{};
     proxsuite::linalg::veg::Vec<isize> perm;
     proxsuite::linalg::veg::Vec<isize> perm_inv;
   
     // sorted on a best effort basis
     proxsuite::linalg::veg::Vec<T> maybe_sorted_diag;
   
     VEG_REFLECT(Ldlt, ld_storage, stride, perm, perm_inv, maybe_sorted_diag);
   
     static auto adjusted_stride(isize n) noexcept -> isize
     {
       return _detail::adjusted_stride<T>(n);
     }
   
     // soft invariants:
     // - perm.len() == perm_inv.len() == dim
     // - dim < stride
     // - ld_storage.len() >= dim * stride
   public:
     Ldlt() = default;
   
     void reserve_uninit(isize cap) noexcept
     {
       static_assert(VEG_CONCEPT(nothrow_constructible<T>), ".");
   
       auto new_stride = adjusted_stride(cap);
       if (cap <= stride && cap * new_stride <= ld_storage.len()) {
         return;
       }
   
       ld_storage.reserve_exact(cap * new_stride);
       perm.reserve_exact(cap);
       perm_inv.reserve_exact(cap);
       maybe_sorted_diag.reserve_exact(cap);
   
       ld_storage.resize_for_overwrite(cap * new_stride);
       stride = new_stride;
     }
   
     void reserve(isize cap) noexcept
     {
       auto new_stride = adjusted_stride(cap);
       if (cap <= stride && cap * new_stride <= ld_storage.len()) {
         return;
       }
       auto n = dim();
   
       ld_storage.reserve_exact(cap * new_stride);
       perm.reserve_exact(cap);
       perm_inv.reserve_exact(cap);
       maybe_sorted_diag.reserve_exact(cap);
   
       ld_storage.resize_for_overwrite(cap * new_stride);
   
       for (isize i = 0; i < n; ++i) {
         auto col = n - i - 1;
         T* ptr = ld_col_mut().data();
         std::move_backward( //
           ptr + col * stride,
           ptr + col * stride + n,
           ptr + col * new_stride + n);
       }
       stride = new_stride;
     }
   
     static auto rank_r_update_req(isize n, isize r) noexcept
       -> proxsuite::linalg::veg::dynstack::StackReq
     {
       auto w_req = proxsuite::linalg::veg::dynstack::StackReq{
         _detail::adjusted_stride<T>(n) * r * isize{ sizeof(T) },
         _detail::align<T>(),
       };
       auto alpha_req = proxsuite::linalg::veg::dynstack::StackReq{
         r * isize{ sizeof(T) },
         alignof(T),
       };
       return w_req & alpha_req;
     }
   
     static auto delete_at_req(isize n, isize r) noexcept
       -> proxsuite::linalg::veg::dynstack::StackReq
     {
       return proxsuite::linalg::veg::dynstack::StackReq{
         r * isize{ sizeof(isize) },
         alignof(isize),
       } &
              proxsuite::linalg::dense::ldlt_delete_rows_and_cols_req(
                proxsuite::linalg::veg::Tag<T>{}, n, r);
     }
   
     void delete_at(isize const* indices,
                    isize r,
                    proxsuite::linalg::veg::dynstack::DynStackMut stack)
     {
       if (r == 0) {
         return;
       }
   
       VEG_ASSERT(std::is_sorted(indices, indices + r));
   
       isize n = dim();
   
       auto _indices_actual =
         stack.make_new_for_overwrite(proxsuite::linalg::veg::Tag<isize>{}, r);
       auto* indices_actual = _indices_actual.ptr_mut();
   
       for (isize k = 0; k < r; ++k) {
         indices_actual[k] = perm_inv[indices[k]];
       }
   
       proxsuite::linalg::dense::ldlt_delete_rows_and_cols_sort_indices( //
         ld_col_mut(),
         indices_actual,
         r,
         stack);
   
       // PERF: do this in one pass
       for (isize k = 0; k < r; ++k) {
         auto i_actual = indices_actual[r - 1 - k];
         auto i = indices[r - 1 - k];
   
         perm.pop_mid(i_actual);
         perm_inv.pop_mid(i);
         maybe_sorted_diag.pop_mid(i_actual);
   
         for (isize j = 0; j < n - 1 - k; ++j) {
           auto& p_j = perm[j];
           auto& pinv_j = perm_inv[j];
   
           if (p_j > i) {
             --p_j;
           }
           if (pinv_j > i_actual) {
             --pinv_j;
           }
         }
       }
     }
   
     auto choose_insertion_position(isize i, Eigen::Ref<Vec const> a) -> isize
     {
       isize n = dim();
       auto diag_elem = a[i];
   
       isize pos = 0;
       for (; pos < n; ++pos) {
         if (diag_elem >= maybe_sorted_diag[pos]) {
           break;
         }
       }
       return pos;
     }
   
     static auto insert_block_at_req(isize n, isize r) noexcept
       -> proxsuite::linalg::veg::dynstack::StackReq
     {
       using proxsuite::linalg::veg::dynstack::StackReq;
       return StackReq{
         isize{ sizeof(T) } * (adjusted_stride(n + r) * r),
         _detail::align<T>(),
       } &
              proxsuite::linalg::dense::ldlt_insert_rows_and_cols_req(
                proxsuite::linalg::veg::Tag<T>{}, n, r);
     }
   
     void insert_block_at(isize i,
                          Eigen::Ref<ColMat const> a,
                          proxsuite::linalg::veg::dynstack::DynStackMut stack)
     {
   
       isize n = dim();
       isize r = a.cols();
   
       if (r == 0) {
         return;
       }
   
       reserve(n + r);
   
       isize i_actual = choose_insertion_position(i, a.col(0));
   
       for (isize j = 0; j < n; ++j) {
         auto& p_j = perm[j];
         auto& pinv_j = perm_inv[j];
   
         if (p_j >= i) {
           p_j += r;
         }
         if (pinv_j >= i_actual) {
           pinv_j += r;
         }
       }
   
       for (isize k = 0; k < r; ++k) {
         perm.push_mid(i + k, i_actual + k);
         perm_inv.push_mid(i_actual + k, i + k);
         maybe_sorted_diag.push_mid(a(i + k, k), i_actual + k);
       }
   
       LDLT_TEMP_MAT_UNINIT(T, permuted_a, n + r, r, stack);
   
       for (isize k = 0; k < r; ++k) {
         for (isize j = 0; j < n + r; ++j) {
           permuted_a(j, k) = a(perm[j], k);
         }
       }
   
       proxsuite::linalg::dense::ldlt_insert_rows_and_cols(
         ld_col_mut(), i_actual, permuted_a, stack);
     }
   
     static auto diagonal_update_req(isize n, isize r) noexcept
       -> proxsuite::linalg::veg::dynstack::StackReq
     {
       using proxsuite::linalg::veg::dynstack::StackReq;
       auto algo_req = StackReq{
         2 * r * isize{ sizeof(isize) },
         alignof(isize),
       };
       auto w_req = StackReq{
         _detail::adjusted_stride<T>(n) * r * isize{ sizeof(T) },
         _detail::align<T>(),
       };
       auto alpha_req = StackReq{
         r * isize{ sizeof(T) },
         alignof(T),
       };
       return algo_req & w_req & alpha_req;
     }
   
     void diagonal_update_clobber_indices( //
       isize* indices,
       isize r,
       Eigen::Ref<Vec const> alpha,
       proxsuite::linalg::veg::dynstack::DynStackMut stack)
     {
   
       if (r == 0) {
         return;
       }
   
       auto _positions =
         stack.make_new_for_overwrite(proxsuite::linalg::veg::Tag<isize>{}, r);
       auto _sorted_indices =
         stack.make_new_for_overwrite(proxsuite::linalg::veg::Tag<isize>{}, r);
       auto* positions = _positions.ptr_mut();
       auto* sorted_indices = _sorted_indices.ptr_mut();
   
       for (isize k = 0; k < r; ++k) {
         indices[k] = perm_inv[indices[k]];
         positions[k] = k;
       }
   
       std::sort(
         positions, positions + r, [indices](isize i, isize j) noexcept -> bool {
           return indices[i] < indices[j];
         });
   
       for (isize k = 0; k < r; ++k) {
         sorted_indices[k] = indices[positions[k]];
       }
   
       auto first = sorted_indices[0];
       auto n = dim() - first;
   
       LDLT_TEMP_MAT(T, _w, n, r, stack);
       LDLT_TEMP_VEC_UNINIT(T, _alpha, r, stack);
   
       for (isize k = 0; k < r; ++k) {
         _alpha(k) = alpha(positions[k]);
         _w(sorted_indices[k] - first, k) = 1;
       }
   
       proxsuite::linalg::dense::_detail::rank_r_update_clobber_w_impl(
         util::submatrix(ld_col_mut(), first, first, n, n),
         _w.data(),
         _w.outerStride(),
         _alpha.data(),
         _detail::IndicesR{
           first,
           0,
           r,
           sorted_indices,
         });
     }
   
     void rank_r_update( //
       Eigen::Ref<ColMat const> w,
       Eigen::Ref<Vec const> alpha,
       proxsuite::linalg::veg::dynstack::DynStackMut stack)
     {
   
       auto n = dim();
       auto r = w.cols();
       if (r == 0) {
         return;
       }
   
       VEG_ASSERT(w.rows() == n);
   
       LDLT_TEMP_MAT_UNINIT(T, _w, n, r, stack);
       LDLT_TEMP_VEC_UNINIT(T, _alpha, r, stack);
   
       for (isize k = 0; k < r; ++k) {
         auto alpha_tmp = alpha(k);
         _alpha(k) = alpha_tmp;
         for (isize i = 0; i < n; ++i) {
           auto w_tmp = w(perm[i], k);
           _w(i, k) = w_tmp;
           maybe_sorted_diag[i] += alpha_tmp * (w_tmp * w_tmp);
         }
       }
   
       proxsuite::linalg::dense::rank_r_update_clobber_inputs(
         ld_col_mut(), _w, _alpha);
     }
   
     auto dim() const noexcept -> isize { return perm.len(); }
   
     auto ld_col() const noexcept -> Eigen::Map< //
       ColMat const,
       Eigen::Unaligned,
       Eigen::OuterStride<DYN>>
     {
       return { ld_storage.ptr(), dim(), dim(), stride };
     }
     auto ld_col_mut() noexcept -> Eigen::Map< //
       ColMat,
       Eigen::Unaligned,
       Eigen::OuterStride<DYN>>
     {
       return { ld_storage.ptr_mut(), dim(), dim(), stride };
     }
     auto ld_row() const noexcept -> Eigen::Map< //
       RowMat const,
       Eigen::Unaligned,
       Eigen::OuterStride<DYN>>
     {
       return {
         ld_storage.ptr(),
         dim(),
         dim(),
         Eigen::OuterStride<DYN>{ stride },
       };
     }
     auto ld_row_mut() noexcept -> Eigen::Map< //
       RowMat,
       Eigen::Unaligned,
       Eigen::OuterStride<DYN>>
     {
       return {
         ld_storage.ptr_mut(),
         dim(),
         dim(),
         Eigen::OuterStride<DYN>{ stride },
       };
     }
   
     auto l() const noexcept -> LView
     {
       return ld_col().template triangularView<Eigen::UnitLower>();
     }
     auto l_mut() noexcept -> LViewMut
     {
       return ld_col_mut().template triangularView<Eigen::UnitLower>();
     }
     auto lt() const noexcept -> LTView
     {
       return ld_row().template triangularView<Eigen::UnitUpper>();
     }
     auto lt_mut() noexcept -> LTViewMut
     {
       return ld_row_mut().template triangularView<Eigen::UnitUpper>();
     }
   
     auto d() const noexcept -> DView
     {
       return {
         ld_storage.ptr(),
         dim(),
         1,
         Eigen::InnerStride<DYN>{ stride + 1 },
       };
     }
     auto d_mut() noexcept -> DViewMut
     {
       return {
         ld_storage.ptr_mut(),
         dim(),
         1,
         Eigen::InnerStride<DYN>{ stride + 1 },
       };
     }
     auto p() const -> Perm { return { VecMapISize(perm.ptr(), dim()) }; }
     auto pt() const -> Perm { return { VecMapISize(perm_inv.ptr(), dim()) }; }
   
     static auto factorize_req(isize n)
       -> proxsuite::linalg::veg::dynstack::StackReq
     {
       return proxsuite::linalg::veg::dynstack::StackReq{
         n * adjusted_stride(n) * isize{ sizeof(T) },
         _detail::align<T>(),
       } |
              proxsuite::linalg::dense::factorize_req(
                proxsuite::linalg::veg::Tag<T>{}, n);
     }
   
     void factorize(Eigen::Ref<ColMat const> mat /* NOLINT */,
                    proxsuite::linalg::veg::dynstack::DynStackMut stack)
     {
       VEG_ASSERT(mat.rows() == mat.cols());
       isize n = mat.rows();
       reserve_uninit(n);
   
       perm.resize_for_overwrite(n);
       perm_inv.resize_for_overwrite(n);
       maybe_sorted_diag.resize_for_overwrite(n);
   
       proxsuite::linalg::dense::_detail::compute_permutation( //
         perm.ptr_mut(),
         perm_inv.ptr_mut(),
         util::diagonal(mat));
   
       {
         LDLT_TEMP_MAT_UNINIT(T, work, n, n, stack);
         ld_col_mut() = mat;
         proxsuite::linalg::dense::_detail::apply_permutation_tri_lower(
           ld_col_mut(), work, perm.ptr());
       }
   
       for (isize i = 0; i < n; ++i) {
         maybe_sorted_diag[i] = ld_col()(i, i);
       }
   
       proxsuite::linalg::dense::factorize(ld_col_mut(), stack);
     }
   
     static auto solve_in_place_req(isize n)
       -> proxsuite::linalg::veg::dynstack::StackReq
     {
       return {
         n * isize{ sizeof(T) },
         _detail::align<T>(),
       };
     }
   
     void solve_in_place(Eigen::Ref<Vec> rhs,
                         proxsuite::linalg::veg::dynstack::DynStackMut stack) const
     {
       isize n = rhs.rows();
       LDLT_TEMP_VEC_UNINIT(T, work, n, stack);
   
       for (isize i = 0; i < n; ++i) {
         work[i] = rhs[perm[i]];
       }
   
       proxsuite::linalg::dense::solve(ld_col(), work);
   
       for (isize i = 0; i < n; ++i) {
         rhs[i] = work[perm_inv[i]];
       }
     }
   
     auto dbg_reconstructed_matrix_internal() const -> ColMat
     {
       isize n = dim();
       auto tmp = ColMat(n, n);
       tmp = l();
       tmp = tmp * d().asDiagonal();
       auto A = ColMat(tmp * lt());
       return A;
     }
   
     auto dbg_reconstructed_matrix() const -> ColMat
     {
       isize n = dim();
       auto tmp = ColMat(n, n);
       tmp = l();
       tmp = tmp * d().asDiagonal();
       auto A = ColMat(tmp * lt());
   
       for (isize i = 0; i < n; i++) {
         tmp.row(i) = A.row(perm_inv[i]);
       }
       for (isize i = 0; i < n; i++) {
         A.col(i) = tmp.col(perm_inv[i]);
       }
       return A;
     }
   };
   } // namespace dense
   } // namespace linalg
   } // namespace proxsuite
   
   #endif /* end of include guard PROXSUITE_LINALG_DENSE_LDLT_LDLT_HPP */