Function ldlt_with_perm

test/src/sparse_factorization.cpp:73–104

template<typename T, typename I>
auto
ldlt_with_perm(proxsuite::linalg::veg::Slice<I> perm_inv,
               proxsuite::linalg::sparse::MatRef<T, I> a)
  -> Eigen::Matrix<T, -1, -1, Eigen::ColMajor>
{
  using Mat = Eigen::Matrix<T, -1, -1, Eigen::ColMajor>;

  VEG_ASSERT(a.nrows() == a.ncols());

  proxsuite::linalg::veg::isize n = a.nrows();

  Eigen::PermutationMatrix<-1, -1, I> perm_inv_eigen = to_eigen_perm(perm_inv);

  auto ld_perm_eigen =
    Mat((perm_inv_eigen *
         Mat(::to_eigen(a).template selfadjointView<Eigen::Upper>()) *
         perm_inv_eigen.inverse())
          .template triangularView<Eigen::Lower>());
  {
    for (proxsuite::linalg::veg::isize i = 0; i < n; ++i) {
      auto a12 = ld_perm_eigen.row(i).head(i).transpose();
      auto l11 = ld_perm_eigen.topLeftCorner(i, i)
                   .template triangularView<Eigen::UnitLower>();
      auto l12 = ld_perm_eigen.row(i).head(i).transpose();
      auto d1 = ld_perm_eigen.diagonal().head(i);
      l12 = l11.solve(a12);
      l12 = d1.asDiagonal().inverse() * l12;
      ld_perm_eigen(i, i) -= l12.dot(d1.asDiagonal() * l12);
    }
  }
  return ld_perm_eigen;
}

using namespace proxsuite::linalg::sparse;
using namespace proxsuite::linalg::veg;