core.hpp

Go to the documentation of this file.

 
//
// Copyright (c) 2022 INRIA
//
#ifndef PROXSUITE_LINALG_DENSE_LDLT_CORE_HPP
#define PROXSUITE_LINALG_DENSE_LDLT_CORE_HPP
 
#include <proxsuite/linalg/veg/util/dbg.hpp>
#include <proxsuite/linalg/veg/util/assert.hpp>
#include <proxsuite/linalg/veg/memory/dynamic_stack.hpp>
 
#if !(defined(__aarch64__) || defined(__PPC64__) || defined(__ppc64__) ||      \
      defined(_ARCH_PPC64))
#include <immintrin.h>
#endif
 
#ifdef PROXSUITE_VECTORIZE
#include <cmath> // to avoid error of the type no member named 'isnan' in namespace 'std';
#include <simde/x86/avx2.h>
#include <simde/x86/fma.h>
#endif
 
#include <Eigen/Core>
 
#define LDLT_ID(id) __VEG_PP_CAT(id, __LINE__)
 
#define __LDLT_TEMP_VEC_IMPL(Type, Name, Rows, Stack, Make)                    \
  auto LDLT_ID(vec_storage) = (Stack).Make(                                    \
    ::proxsuite::linalg::veg::Tag<__VEG_PP_REMOVE_PAREN(Type)>{},              \
    (Rows),                                                                    \
    ::proxsuite::linalg::dense::_detail::align<__VEG_PP_REMOVE_PAREN(          \
      Type)>());                                                               \
  auto Name /* NOLINT */ = ::Eigen::Map<                                       \
    ::Eigen::Matrix<__VEG_PP_REMOVE_PAREN(Type), ::Eigen::Dynamic, 1>,         \
    ::Eigen::Unaligned,                                                        \
    ::Eigen::Stride<::Eigen::Dynamic, 1>>{                                     \
    LDLT_ID(vec_storage).ptr_mut(),                                            \
    LDLT_ID(vec_storage).len(),                                                \
    ::Eigen::Stride<::Eigen::Dynamic, 1>{                                      \
      LDLT_ID(vec_storage).len(),                                              \
      1,                                                                       \
    },                                                                         \
  };                                                                           \
  static_assert(true, ".")
 
#define __LDLT_TEMP_MAT_IMPL(Type, Name, Rows, Cols, Stack, Make)              \
  ::proxsuite::linalg::veg::isize LDLT_ID(rows) = (Rows);                      \
  ::proxsuite::linalg::veg::isize LDLT_ID(cols) = (Cols);                      \
  ::proxsuite::linalg::veg::isize LDLT_ID(stride) =                            \
    ::proxsuite::linalg::dense::_detail::adjusted_stride<                      \
      __VEG_PP_REMOVE_PAREN(Type)>(LDLT_ID(rows));                             \
  auto LDLT_ID(vec_storage) = (Stack).Make(                                    \
    ::proxsuite::linalg::veg::Tag<__VEG_PP_REMOVE_PAREN(Type)>{},              \
    LDLT_ID(stride) * LDLT_ID(cols),                                           \
    ::proxsuite::linalg::dense::_detail::align<__VEG_PP_REMOVE_PAREN(          \
      Type)>());                                                               \
  auto Name /* NOLINT */ =                                                     \
    ::Eigen::Map<::Eigen::Matrix<__VEG_PP_REMOVE_PAREN(Type),                  \
                                 ::Eigen::Dynamic,                             \
                                 ::Eigen::Dynamic,                             \
                                 ::Eigen::ColMajor>,                           \
                 ::Eigen::Unaligned,                                           \
                 ::Eigen::Stride<::Eigen::Dynamic, 1>>{                        \
      LDLT_ID(vec_storage).ptr_mut(),                                          \
      LDLT_ID(rows),                                                           \
      LDLT_ID(cols),                                                           \
      ::Eigen::Stride<::Eigen::Dynamic, 1>{                                    \
        LDLT_ID(stride),                                                       \
        1,                                                                     \
      },                                                                       \
    };                                                                         \
  static_assert(true, ".")
 
#define LDLT_TEMP_VEC(Type, Name, Rows, Stack)                                 \
  __LDLT_TEMP_VEC_IMPL(Type, Name, Rows, Stack, make_new)
#define LDLT_TEMP_VEC_UNINIT(Type, Name, Rows, Stack)                          \
  __LDLT_TEMP_VEC_IMPL(Type, Name, Rows, Stack, make_new_for_overwrite)
 
#define LDLT_TEMP_MAT(Type, Name, Rows, Cols, Stack)                           \
  __LDLT_TEMP_MAT_IMPL(Type, Name, Rows, Cols, Stack, make_new)
#define LDLT_TEMP_MAT_UNINIT(Type, Name, Rows, Cols, Stack)                    \
  __LDLT_TEMP_MAT_IMPL(Type, Name, Rows, Cols, Stack, make_new_for_overwrite)
 
namespace proxsuite {
namespace linalg {
namespace dense {
using proxsuite::linalg::veg::i32;
using proxsuite::linalg::veg::isize;
using proxsuite::linalg::veg::u32;
using proxsuite::linalg::veg::usize;
using f32 = float;
using f64 = double;
 
namespace _detail {
namespace _simd {
 
#ifdef __clang__
#define DENSE_LDLT_FP_PRAGMA _Pragma("STDC FP_CONTRACT ON")
#else
#define DENSE_LDLT_FP_PRAGMA
#endif
 
static_assert(static_cast<unsigned char>(-1) == 255, "char should have 8 bits");
static_assert(sizeof(f32) == 4, "f32 should be 32 bits");
static_assert(sizeof(f64) == 8, "f64 should be 64 bits");
 
#define LDLT_FN_IMPL3(Fn, Prefix, Suffix)                                      \
  VEG_INLINE static auto Fn(Pack a, Pack b, Pack c) noexcept -> Pack           \
  {                                                                            \
    return Pack{ simde_mm##Prefix##_##Fn##_##Suffix(                           \
      a.inner, b.inner, c.inner) };                                            \
  }                                                                            \
  VEG_NOM_SEMICOLON
 
#define LDLT_ARITHMETIC_IMPL(Prefix, Suffix)                                   \
  LDLT_FN_IMPL3(fmadd, Prefix, Suffix);  /* (a * b + c) */                     \
  LDLT_FN_IMPL3(fnmadd, Prefix, Suffix); /* (-a * b + c) */
 
#define LDLT_LOAD_STORE(Prefix, Suffix)                                        \
  VEG_INLINE static auto load_unaligned(ScalarType const* ptr) noexcept        \
    -> Pack                                                                    \
  {                                                                            \
    return Pack{ simde_mm##Prefix##_loadu_##Suffix(ptr) };                     \
  }                                                                            \
  VEG_INLINE static auto broadcast(ScalarType value) noexcept -> Pack          \
  {                                                                            \
    return Pack{ simde_mm##Prefix##_set1_##Suffix(value) };                    \
  }                                                                            \
  VEG_INLINE void store_unaligned(ScalarType* ptr) const noexcept              \
  {                                                                            \
    simde_mm##Prefix##_storeu_##Suffix(ptr, inner);                            \
  }                                                                            \
  VEG_NOM_SEMICOLON
 
template<typename T, usize N>
struct Pack;
 
template<typename T>
struct Pack<T, 1>
{
  using ScalarType = T;
 
  T inner;
 
  VEG_INLINE static auto fmadd(Pack a, Pack b, Pack c) noexcept -> Pack
  {
    DENSE_LDLT_FP_PRAGMA
    return { a.inner * b.inner + c.inner };
  }
  VEG_INLINE static auto fnmadd(Pack a, Pack b, Pack c) noexcept -> Pack
  {
    return fmadd({ -a.inner }, b, c);
  }
  VEG_INLINE static auto load_unaligned(ScalarType const* ptr) noexcept -> Pack
  {
    return { *ptr };
  }
  VEG_INLINE static auto broadcast(ScalarType value) noexcept -> Pack
  {
    return { value };
  }
  VEG_INLINE void store_unaligned(ScalarType* ptr) const noexcept
  {
    *ptr = inner;
  }
};
 
#ifdef PROXSUITE_VECTORIZE
template<>
struct Pack<f32, 4>
{
  using ScalarType = f32;
 
  simde__m128 inner;
  LDLT_ARITHMETIC_IMPL(, ps)
  LDLT_LOAD_STORE(, ps);
};
 
template<>
struct Pack<f32, 8>
{
  using ScalarType = f32;
 
  simde__m256 inner;
  LDLT_ARITHMETIC_IMPL(256, ps)
  LDLT_LOAD_STORE(256, ps);
};
 
#ifdef __AVX512F__
template<>
struct Pack<f32, 16>
{
  using ScalarType = f32;
 
  __m512 inner;
  VEG_INLINE static auto fmadd(Pack a, Pack b, Pack c) noexcept -> Pack
  {
    DENSE_LDLT_FP_PRAGMA
    return { _mm512_fmadd_ps(a.inner, b.inner, c.inner) };
  }
  VEG_INLINE static auto fnmadd(Pack a, Pack b, Pack c) noexcept -> Pack
  {
    return { _mm512_fnmadd_ps(a.inner, b.inner, c.inner) };
  }
  VEG_INLINE static auto load_unaligned(ScalarType const* ptr) noexcept -> Pack
  {
    return { _mm512_loadu_ps(ptr) };
  }
  VEG_INLINE static auto broadcast(ScalarType value) noexcept -> Pack
  {
    return { _mm512_set1_ps(value) };
  }
  VEG_INLINE void store_unaligned(ScalarType* ptr) const noexcept
  {
    _mm512_storeu_ps(ptr, inner);
  }
};
#endif
 
template<>
struct Pack<f64, 2>
{
  using ScalarType = f64;
 
  simde__m128d inner;
  LDLT_ARITHMETIC_IMPL(, pd)
  LDLT_LOAD_STORE(, pd);
};
template<>
struct Pack<f64, 4>
{
  using ScalarType = f64;
 
  simde__m256d inner;
  LDLT_ARITHMETIC_IMPL(256, pd)
  LDLT_LOAD_STORE(256, pd);
};
 
#ifdef __AVX512F__
template<>
struct Pack<f64, 8>
{
  using ScalarType = f64;
 
  __m512d inner;
  VEG_INLINE static auto fmadd(Pack a, Pack b, Pack c) noexcept -> Pack
  {
    DENSE_LDLT_FP_PRAGMA
    return { _mm512_fmadd_pd(a.inner, b.inner, c.inner) };
  }
  VEG_INLINE static auto fnmadd(Pack a, Pack b, Pack c) noexcept -> Pack
  {
    return { _mm512_fnmadd_pd(a.inner, b.inner, c.inner) };
  }
  VEG_INLINE static auto load_unaligned(ScalarType const* ptr) noexcept -> Pack
  {
    return { _mm512_loadu_pd(ptr) };
  }
  VEG_INLINE static auto broadcast(ScalarType value) noexcept -> Pack
  {
    return { _mm512_set1_pd(value) };
  }
  VEG_INLINE void store_unaligned(ScalarType* ptr) const noexcept
  {
    _mm512_storeu_pd(ptr, inner);
  }
};
#endif
 
#endif
 
template<typename T>
struct NativePackInfo
{
  static constexpr usize N = 1;
  using Type = Pack<f32, N>;
};
 
#ifdef PROXSUITE_VECTORIZE
template<>
struct NativePackInfo<f32>
{
  static constexpr usize N = SIMDE_NATURAL_VECTOR_SIZE / 32;
  using Type = Pack<f32, N>;
};
template<>
struct NativePackInfo<f64>
{
  static constexpr usize N = SIMDE_NATURAL_VECTOR_SIZE / 64;
  using Type = Pack<f64, N>;
};
#endif
 
template<typename T>
using NativePack = typename NativePackInfo<T>::Type;
} // namespace _simd
} // namespace _detail
 
namespace _detail {
using proxsuite::linalg::veg::uncvref_t;
template<bool COND, typename T>
using const_if = proxsuite::linalg::veg::meta::if_t<COND, T const, T>;
 
template<typename T>
using ptr_is_const = proxsuite::linalg::veg::meta::bool_constant<VEG_CONCEPT(
  const_type<proxsuite::linalg::veg::meta::unptr_t<T>>)>;
 
template<typename T>
constexpr auto
round_up(T a, T b) noexcept -> T
{
  return a + (b - 1) / b * b;
}
 
#ifdef PROXSUITE_VECTORIZE
template<typename T>
using should_vectorize =
  proxsuite::linalg::veg::meta::bool_constant<VEG_CONCEPT(same<T, f32>) ||
                                              VEG_CONCEPT(same<T, f64>)>;
#else
template<typename T>
using should_vectorize = proxsuite::linalg::veg::meta::bool_constant<false>;
#endif
 
template<typename T>
auto
adjusted_stride(isize n) noexcept -> isize
{
#ifndef SIMDE_NATURAL_FLOAT_VECTOR_SIZE
  isize simd_stride = 1;
#else
  isize simd_stride =
    (SIMDE_NATURAL_VECTOR_SIZE / CHAR_BIT) / isize{ sizeof(T) };
#endif
  return _detail::should_vectorize<T>::value ? _detail::round_up(n, simd_stride)
                                             : n;
}
template<typename T>
auto
align() noexcept -> isize
{
  return isize{ alignof(T) } * _detail::adjusted_stride<T>(1);
}
 
struct NoCopy
{
  NoCopy() = default;
  ~NoCopy() = default;
 
  NoCopy(NoCopy const&) = delete;
  NoCopy(NoCopy&&) = delete;
  auto operator=(NoCopy const&) -> NoCopy& = delete;
  auto operator=(NoCopy&&) -> NoCopy& = delete;
};
 
namespace nb {
struct max2
{
  template<typename T>
  VEG_INLINE constexpr auto operator()(T const& a, T const& b) const -> T const&
  {
    return a > b ? a : b;
  }
};
struct min2
{
  template<typename T>
  VEG_INLINE constexpr auto operator()(T a, T b) const -> T
  {
    return (a < b) ? a : b;
  }
};
} // namespace nb
VEG_NIEBLOID(min2);
VEG_NIEBLOID(max2);
 
template<typename T>
void
set_zero(T* dest, usize n)
{
  for (usize i = 0; i < n; ++i) {
    *dest = 0;
  }
}
 
template<typename T>
using OwnedMatrix =
  Eigen::Matrix<typename T::Scalar,
                Eigen::Dynamic,
                Eigen::Dynamic,
                bool(T::IsRowMajor) ? Eigen::RowMajor : Eigen::ColMajor>;
 
template<typename T>
using OwnedAll =
  Eigen::Matrix<typename T::Scalar,
                T::RowsAtCompileTime,
                T::ColsAtCompileTime,
                bool(T::IsRowMajor) ? Eigen::RowMajor : Eigen::ColMajor>;
 
template<typename T>
using OwnedRows =
  Eigen::Matrix<typename T::Scalar,
                Eigen::Dynamic,
                T::ColsAtCompileTime,
                bool(T::IsRowMajor) ? Eigen::RowMajor : Eigen::ColMajor>;
 
template<typename T>
using OwnedCols =
  Eigen::Matrix<typename T::Scalar,
                T::RowsAtCompileTime,
                Eigen::Dynamic,
                bool(T::IsRowMajor) ? Eigen::RowMajor : Eigen::ColMajor>;
 
template<typename T>
using OwnedColVector = Eigen::Matrix< //
  typename T::Scalar,
  Eigen::Dynamic,
  1,
  Eigen::ColMajor>;
template<typename T>
using OwnedRowVector = Eigen::Matrix< //
  typename T::Scalar,
  1,
  Eigen::Dynamic,
  Eigen::RowMajor>;
 
template<bool ROWMAJOR>
struct ElemAddrImpl;
 
template<>
struct ElemAddrImpl<false>
{
  template<typename T>
  static auto fn(T* ptr,
                 isize row,
                 isize col,
                 isize outer_stride,
                 isize inner_stride) noexcept -> T*
  {
    return ptr + (inner_stride * row + outer_stride * col);
  }
};
template<>
struct ElemAddrImpl<true>
{
  template<typename T>
  static auto fn(T* ptr,
                 isize row,
                 isize col,
                 isize outer_stride,
                 isize inner_stride) noexcept -> T*
  {
    return ptr + (outer_stride * row + inner_stride * col);
  }
};
 
template<bool COLMAJOR>
struct RowColAccessImpl;
 
template<>
struct RowColAccessImpl<true>
{
  template<typename T>
  using Col =
    Eigen::Map<const_if<ptr_is_const<decltype(VEG_DECLVAL(T&&).data())>::value,
                        OwnedColVector<uncvref_t<T>>>,
               Eigen::Unaligned,
               Eigen::InnerStride<uncvref_t<T>::InnerStrideAtCompileTime>>;
  template<typename T>
  using Row =
    Eigen::Map<const_if<ptr_is_const<decltype(VEG_DECLVAL(T&&).data())>::value,
                        OwnedRowVector<uncvref_t<T>>>,
               Eigen::Unaligned,
               Eigen::InnerStride<uncvref_t<T>::OuterStrideAtCompileTime>>;
 
  template<typename T>
  static auto col(T&& mat, isize col_idx) noexcept -> Col<T>
  {
    return {
      mat.data() + col_idx * mat.outerStride(),
      mat.rows(),
      1,
      Eigen::InnerStride<uncvref_t<T>::InnerStrideAtCompileTime>{
        mat.innerStride(),
      },
    };
  }
  template<typename T>
  static auto row(T&& mat, isize row_idx) noexcept -> Row<T>
  {
    return {
      mat.data() + row_idx * mat.innerStride(),
      1,
      mat.cols(),
      Eigen::InnerStride<uncvref_t<T>::OuterStrideAtCompileTime>{
        mat.outerStride(),
      },
    };
  }
};
template<>
struct RowColAccessImpl<false>
{
  template<typename T>
  using Col =
    Eigen::Map<const_if<ptr_is_const<decltype(VEG_DECLVAL(T&&).data())>::value,
                        OwnedColVector<uncvref_t<T>>>,
               Eigen::Unaligned,
               Eigen::InnerStride<uncvref_t<T>::OuterStrideAtCompileTime>>;
  template<typename T>
  using Row =
    Eigen::Map<const_if<ptr_is_const<decltype(VEG_DECLVAL(T&&).data())>::value,
                        OwnedRowVector<uncvref_t<T>>>,
               Eigen::Unaligned,
               Eigen::InnerStride<uncvref_t<T>::InnerStrideAtCompileTime>>;
 
  template<typename T>
  static auto col(T&& mat, isize col_idx) noexcept -> Col<T>
  {
    return {
      mat.data() + col_idx * mat.innerStride(),
      mat.rows(),
      1,
      Eigen::InnerStride<uncvref_t<T>::OuterStrideAtCompileTime>{
        mat.outerStride(),
      },
    };
  }
  template<typename T>
  static auto row(T&& mat, isize row_idx) noexcept -> Row<T>
  {
    return {
      mat.data() + row_idx * mat.outerStride(),
      1,
      mat.cols(),
      Eigen::InnerStride<uncvref_t<T>::InnerStrideAtCompileTime>{
        mat.innerStride(),
      },
    };
  }
};
template<typename T>
using StrideOf = Eigen::Stride< //
  T::OuterStrideAtCompileTime,
  T::InnerStrideAtCompileTime>;
} // namespace _detail
 
namespace util {
template<bool COLMAJOR, typename T>
auto
elem_addr(T* ptr,
          isize row,
          isize col,
          isize outer_stride,
          isize inner_stride) noexcept -> T*
{
  return _detail::ElemAddrImpl<!COLMAJOR>::fn(
    ptr, row, col, outer_stride, inner_stride);
}
 
template<typename Mat>
auto
matrix_elem_addr(Mat&& mat, isize row, isize col) noexcept
  -> decltype(mat.data())
{
  return util::elem_addr<!bool(
    proxsuite::linalg::veg::uncvref_t<Mat>::IsRowMajor)>( //
    mat.data(),
    row,
    col,
    mat.outerStride(),
    mat.innerStride());
}
 
template<typename T>
auto
col(T&& mat, isize col_idx) noexcept -> typename _detail::RowColAccessImpl<
  !bool(proxsuite::linalg::veg::uncvref_t<T>::IsRowMajor)>::template Col<T>
{
  return _detail::RowColAccessImpl<!bool(
    proxsuite::linalg::veg::uncvref_t<T>::IsRowMajor)>::col(mat, col_idx);
}
template<typename T>
auto
row(T&& mat, isize row_idx) noexcept -> typename _detail::RowColAccessImpl<
  !bool(proxsuite::linalg::veg::uncvref_t<T>::IsRowMajor)>::template Row<T>
{
  return _detail::RowColAccessImpl<!bool(
    proxsuite::linalg::veg::uncvref_t<T>::IsRowMajor)>::row(mat, row_idx);
}
 
template<typename Mat>
auto
trans(Mat&& mat) noexcept -> Eigen::Map< //
  _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
                    Eigen::Matrix< //
                      typename proxsuite::linalg::veg::uncvref_t<Mat>::Scalar,
                      proxsuite::linalg::veg::uncvref_t<Mat>::ColsAtCompileTime,
                      proxsuite::linalg::veg::uncvref_t<Mat>::RowsAtCompileTime,
                      bool(proxsuite::linalg::veg::uncvref_t<Mat>::IsRowMajor)
                        ? Eigen::ColMajor
                        : Eigen::RowMajor>>,
  Eigen::Unaligned,
  _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>>
{
  return {
    mat.data(),
    mat.cols(),
    mat.rows(),
    _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>{
      mat.outerStride(),
      mat.innerStride(),
    },
  };
}
 
template<typename Mat>
auto
diagonal(Mat&& mat) noexcept -> Eigen::Map< //
  _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
                    Eigen::Matrix< //
                      typename proxsuite::linalg::veg::uncvref_t<Mat>::Scalar,
                      Eigen::Dynamic,
                      1,
                      Eigen::ColMajor>>,
  Eigen::Unaligned,
  Eigen::InnerStride<Eigen::Dynamic>>
{
  VEG_DEBUG_ASSERT( //
    mat.rows() == mat.cols());
  return { mat.data(),
           mat.rows(),
           1,
           Eigen::InnerStride<Eigen::Dynamic>{ mat.outerStride() + 1 } };
}
 
template<typename Mat>
auto
submatrix(Mat&& mat,
          isize row_start,
          isize col_start,
          isize nrows,
          isize ncols) noexcept
  -> Eigen::Map<_detail::const_if<
                  _detail::ptr_is_const<decltype(mat.data())>::value,
                  _detail::OwnedMatrix<proxsuite::linalg::veg::uncvref_t<Mat>>>,
                Eigen::Unaligned,
                _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>>
{
  return {
    util::elem_addr<!bool(proxsuite::linalg::veg::uncvref_t<Mat>::IsRowMajor)>(
      mat.data(), row_start, col_start, mat.outerStride(), mat.innerStride()),
    nrows,
    ncols,
    _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>{
      mat.outerStride(),
      mat.innerStride(),
    },
  };
}
 
template<typename Mat>
auto
to_view(Mat&& mat) noexcept -> Eigen::Map<
  _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
                    _detail::OwnedAll<proxsuite::linalg::veg::uncvref_t<Mat>>>,
  Eigen::Unaligned,
  _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>>
{
  return {
    mat.data(),
    mat.rows(),
    mat.cols(),
    _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>{
      mat.outerStride(),
      mat.innerStride(),
    },
  };
}
 
template<typename Mat>
auto
to_view_dyn_rows(Mat&& mat) noexcept -> Eigen::Map<
  _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
                    _detail::OwnedRows<proxsuite::linalg::veg::uncvref_t<Mat>>>,
  Eigen::Unaligned,
  _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>>
{
  return {
    mat.data(),
    mat.rows(),
    mat.cols(),
    _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>{
      mat.outerStride(),
      mat.innerStride(),
    },
  };
}
 
template<typename Mat>
auto
to_view_dyn_cols(Mat&& mat) noexcept -> Eigen::Map<
  _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
                    _detail::OwnedCols<proxsuite::linalg::veg::uncvref_t<Mat>>>,
  Eigen::Unaligned,
  _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>>
{
  return {
    mat.data(),
    mat.rows(),
    mat.cols(),
    _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>{
      mat.outerStride(),
      mat.innerStride(),
    },
  };
}
 
template<typename Mat>
auto
to_view_dyn(Mat&& mat) noexcept
  -> Eigen::Map<_detail::const_if<
                  _detail::ptr_is_const<decltype(mat.data())>::value,
                  _detail::OwnedMatrix<proxsuite::linalg::veg::uncvref_t<Mat>>>,
                Eigen::Unaligned,
                _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>>
{
  return {
    mat.data(),
    mat.rows(),
    mat.cols(),
    _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>{
      mat.outerStride(),
      mat.innerStride(),
    },
  };
}
 
template<typename Mat>
auto
subrows(Mat&& mat, isize row_start, isize nrows) noexcept -> Eigen::Map<
  _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
                    _detail::OwnedRows<proxsuite::linalg::veg::uncvref_t<Mat>>>,
  Eigen::Unaligned,
  _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>>
{
  return {
    util::elem_addr<!bool(proxsuite::linalg::veg::uncvref_t<Mat>::IsRowMajor)>(
      mat.data(), row_start, 0, mat.outerStride(), mat.innerStride()),
    nrows,
    mat.cols(),
    _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>{
      mat.outerStride(),
      mat.innerStride(),
    },
  };
}
 
template<typename Mat>
auto
subcols(Mat&& mat, isize col_start, isize ncols) noexcept -> Eigen::Map<
  _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
                    _detail::OwnedCols<proxsuite::linalg::veg::uncvref_t<Mat>>>,
  Eigen::Unaligned,
  _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>>
{
  return {
    util::elem_addr<!bool(proxsuite::linalg::veg::uncvref_t<Mat>::IsRowMajor)>(
      mat.data(), 0, col_start, mat.outerStride(), mat.innerStride()),
    mat.rows(),
    ncols,
    _detail::StrideOf<proxsuite::linalg::veg::uncvref_t<Mat>>{
      mat.outerStride(),
      mat.innerStride(),
    },
  };
}
 
} // namespace util
namespace _detail {
template<typename Dst, typename Lhs, typename Rhs, typename T>
void
noalias_mul_add_impl(Dst dst, Lhs lhs, Rhs rhs, T factor)
{
  VEG_ASSERT_ALL_OF(dst.rows() == lhs.rows(),
                    dst.cols() == rhs.cols(),
                    lhs.cols() == rhs.rows());
 
  isize nrows = dst.rows();
  isize ncols = dst.cols();
  isize depth = lhs.cols();
 
  if (nrows == 0 || ncols == 0 || depth == 0) {
    return;
  }
 
#if !EIGEN_VERSION_AT_LEAST(3, 3, 8)
#define LAZY_PRODUCT(a, b) a.lazyProduct(b)
#else
#define LAZY_PRODUCT(a, b) a.operator*(b)
#endif
 
#if !EIGEN_VERSION_AT_LEAST(3, 3, 8)
  if ((dst.rows() < 20) && (dst.cols() < 20) && (rhs.rows() < 20)) {
    // gemm
    // workaround for eigen 3.3.7 bug:
    // https://gitlab.com/libeigen/eigen/-/issues/1562
    using Mat =
      Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor, 20, 20>;
    using MapMut =
      Eigen::Map<Mat, Eigen::Unaligned, Eigen::OuterStride<Eigen::Dynamic>>;
 
    auto dst_ =
      MapMut(dst.data(), dst.rows(), dst.cols(), { dst.outerStride() });
    dst_.noalias().operator+=(factor * LAZY_PRODUCT(lhs, rhs));
  } else
#endif
  {
    dst.noalias().operator+=(factor * LAZY_PRODUCT(lhs, rhs));
  }
 
#undef LAZY_PRODUCT
}
} // namespace _detail
namespace util {
template<typename Dst, typename Lhs, typename Rhs, typename T>
void
noalias_mul_add(Dst&& dst, Lhs const& lhs, Rhs const& rhs, T factor)
{
  _detail::noalias_mul_add_impl(
    util::submatrix(dst, 0, 0, dst.rows(), dst.cols()),
    util::submatrix(lhs, 0, 0, lhs.rows(), lhs.cols()),
    util::submatrix(rhs, 0, 0, rhs.rows(), rhs.cols()),
    factor);
}
} // namespace util
template<typename T>
auto
temp_mat_req(proxsuite::linalg::veg::Tag<T> /*tag*/,
             isize rows,
             isize cols) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
{
  return {
    _detail::adjusted_stride<T>(rows) * cols * isize{ sizeof(T) },
    _detail::align<T>(),
  };
}
 
template<typename T>
auto
temp_vec_req(proxsuite::linalg::veg::Tag<T> /*tag*/, isize rows) noexcept
  -> proxsuite::linalg::veg::dynstack::StackReq
{
  return {
    rows * isize{ sizeof(T) },
    _detail::align<T>(),
  };
}
} // namespace dense
} // namespace linalg
} // namespace proxsuite
 
#endif /* end of include guard PROXSUITE_LINALG_DENSE_LDLT_CORE_HPP */