views.hpp

Go to the documentation of this file.

//
// Copyright (c) 2022 INRIA
//
#ifndef PROXSUITE_PROXQP_DENSE_VIEWS_HPP
#define PROXSUITE_PROXQP_DENSE_VIEWS_HPP
 
#include <proxsuite/linalg/veg/type_traits/core.hpp>
#include <proxsuite/linalg/veg/util/dbg.hpp>
#include <cstring>
#include <type_traits>
#include <Eigen/Core>
 
#define LDLT_CONCEPT(...)                                                      \
  VEG_CONCEPT_MACRO(::proxsuite::proxqp::concepts, __VA_ARGS__)
#define LDLT_CHECK_CONCEPT(...)                                                \
  VEG_CHECK_CONCEPT_MACRO(::proxqp::concepts, __VA_ARGS__)
 
namespace proxsuite {
namespace proxqp {
 
using usize = decltype(sizeof(0));
namespace detail {
template<typename Fn>
struct FnInfo;
template<typename Ret_, typename... Args>
struct FnInfo<auto(Args...)->Ret_>
{
  template<usize I>
  using Arg = proxsuite::linalg::veg::ith<I, Args...>;
  using Ret = Ret_;
};
} // namespace detail
 
#define LDLT_IMPL_GET_PARAM(Fn, Idx)                                           \
  typename ::proxsuite::proxqp::detail::FnInfo<                                \
    decltype Fn /* NOLINT */>::template Arg<(Idx)>,
 
#define LDLT_IMPL_GET_PARAMS_0(NParams, ...)                                   \
  __VEG_PP_TUPLE_FOR_EACH(LDLT_IMPL_GET_PARAM,                                 \
                          (__VA_ARGS__),                                       \
                          __VEG_PP_MAKE_TUPLE(__VEG_IMPL_PP_DEC(NParams)))
 
#define LDLT_IMPL_GET_PARAMS_1(NParams, ...)
 
#define LDLT_IMPL_GET_PARAMS(NParams, ...)                                     \
  __VEG_PP_CAT2(LDLT_IMPL_GET_PARAMS_, __VEG_IMPL_PP_IS_1(NParams))            \
  (NParams, __VA_ARGS__)
 
#define LDLT_EXPLICIT_TPL_DEF(NParams, ...)                                    \
  template auto __VA_ARGS__(                                                   \
    LDLT_IMPL_GET_PARAMS(NParams, __VA_ARGS__)                                 \
      typename ::proxsuite::proxqp::detail::FnInfo<                            \
        decltype(__VA_ARGS__)>::template Arg<(NParams) - 1>) ->                \
    typename ::proxsuite::proxqp::detail::FnInfo<decltype(__VA_ARGS__)>::Ret
#define LDLT_EXPLICIT_TPL_DECL(NParams, ...)                                   \
  extern LDLT_EXPLICIT_TPL_DEF(NParams, __VA_ARGS__)
 
using proxsuite::linalg::veg::i32;
using proxsuite::linalg::veg::i64;
using proxsuite::linalg::veg::isize;
using proxsuite::linalg::veg::u32;
using proxsuite::linalg::veg::u64;
using proxsuite::linalg::veg::usize;
using f32 = float;
using f64 = double;
namespace detail {
 
struct NoCopy
{
  NoCopy() = default;
  ~NoCopy() = default;
 
  NoCopy(NoCopy const&) = delete;
  NoCopy(NoCopy&&) = delete;
  auto operator=(NoCopy const&) -> NoCopy& = delete;
  auto operator=(NoCopy&&) -> NoCopy& = delete;
};
 
template<typename Fn>
struct Defer /* NOLINT */
{
  Fn fn;
  NoCopy _;
 
  VEG_INLINE ~Defer() noexcept(noexcept(VEG_FWD(fn)())) { VEG_FWD(fn)(); }
};
 
namespace nb {
struct defer
{
  template<typename Fn>
  VEG_INLINE constexpr auto operator()(Fn fn) const -> Defer<Fn>
  {
    return { VEG_FWD(fn), {} };
  }
};
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
 
template<typename T>
constexpr auto
min_list_impl(T init, T const* arr, usize n) noexcept -> T
{
  return (n == 0)
           ? init
           : nb::min2{}(init, detail::min_list_impl(*arr, arr + 1, n - 1));
}
template<typename T, usize N>
constexpr auto
cx_min_list(T const (&arr)[N]) noexcept -> T
{
  return detail::min_list_impl( //
    arr[0],
    arr + 1,
    N - 1);
}
 
namespace nb {
struct max_list
{
  template<typename T>
  VEG_INLINE auto operator()(std::initializer_list<T> list) const -> T
  {
    T const* data = list.begin();
    isize len = isize(list.size());
 
    T current_max = data[0];
    for (isize i = 1; i < len; ++i) {
      if (data[i] > current_max) {
        current_max = data[i];
      }
    }
    return current_max;
  }
};
} // namespace nb
VEG_NIEBLOID(defer);
VEG_NIEBLOID(max2);
VEG_NIEBLOID(min2);
VEG_NIEBLOID(max_list);
 
template<typename T, bool = std::is_floating_point<T>::value>
struct SetZeroImpl
{
  static void fn(T* dest, usize n)
  {
    for (usize i = 0; i < n; ++i) {
      *dest = 0;
    }
  }
};
 
template<typename T>
struct SetZeroImpl<T, true>
{
  static void fn(T* dest, usize n)
  {
    // TODO: assert bit representation is zero
    std::memset(dest, 0, n * sizeof(T));
  }
};
 
template<typename T>
void
set_zero(T* dest, usize n)
{
  SetZeroImpl<T>::fn(dest, n);
}
 
constexpr auto
round_up(isize n, isize k) noexcept -> isize
{
  return (n + k - 1) / k * k;
}
constexpr auto
uround_up(usize n, usize k) noexcept -> usize
{
  return (n + k - 1) / k * k;
}
 
inline auto
bytes_to_prev_aligned(void* ptr, usize align) noexcept -> isize
{
  using UPtr = std::uintptr_t;
 
  UPtr mask = align - 1;
  UPtr iptr = UPtr(ptr);
  UPtr aligned_ptr = iptr & ~mask;
  return isize(aligned_ptr - iptr);
}
inline auto
bytes_to_next_aligned(void* ptr, usize align) noexcept -> isize
{
  using UPtr = std::uintptr_t;
 
  UPtr mask = align - 1;
  UPtr iptr = UPtr(ptr);
  UPtr aligned_ptr = (iptr + mask) & ~mask;
  return isize(aligned_ptr - iptr);
}
 
inline auto
next_aligned(void* ptr, usize align) noexcept -> void*
{
  using BytePtr = unsigned char*;
  using VoidPtr = void*;
  return VoidPtr(BytePtr(ptr) + detail::bytes_to_next_aligned(ptr, align));
}
inline auto
prev_aligned(void* ptr, usize align) noexcept -> void*
{
  using BytePtr = unsigned char*;
  using VoidPtr = void*;
  return VoidPtr(BytePtr(ptr) + detail::bytes_to_prev_aligned(ptr, align));
}
 
} // namespace detail
 
enum struct Layout : unsigned char
{
  colmajor = 0,
  rowmajor = 1,
};
 
constexpr Layout colmajor = Layout::colmajor;
constexpr Layout rowmajor = Layout::rowmajor;
 
constexpr auto
flip_layout(Layout l) noexcept -> Layout
{
  return Layout(1 - u32(l));
}
constexpr auto
to_eigen_layout(Layout l) -> int
{
  return l == colmajor ? Eigen::ColMajor : Eigen::RowMajor;
}
constexpr auto
from_eigen_layout(int l) -> Layout
{
  return (unsigned(l) & Eigen::RowMajorBit) == Eigen::RowMajor ? rowmajor
                                                               : colmajor;
}
 
static_assert(to_eigen_layout(from_eigen_layout(Eigen::ColMajor)) ==
                Eigen::ColMajor,
              ".");
static_assert(to_eigen_layout(from_eigen_layout(Eigen::RowMajor)) ==
                Eigen::RowMajor,
              ".");
 
namespace detail {
template<Layout L>
struct ElementAccess;
 
template<>
struct ElementAccess<Layout::colmajor>
{
  template<typename T>
  VEG_INLINE static constexpr auto offset(T* ptr,
                                          isize row,
                                          isize col,
                                          isize outer_stride) noexcept -> T*
  {
    return ptr + (usize(row) + usize(col) * usize(outer_stride));
  }
 
  using NextRowStride = Eigen::Stride<0, 0>;
  using NextColStride = Eigen::InnerStride<Eigen::Dynamic>;
  VEG_INLINE static auto next_row_stride(isize outer_stride) noexcept
    -> NextRowStride
  {
    (void)outer_stride;
    return NextRowStride{};
  }
  VEG_INLINE static auto next_col_stride(isize outer_stride) noexcept
    -> NextColStride
  {
    return NextColStride /* NOLINT(modernize-return-braced-init-list) */ (
      outer_stride);
  }
 
  template<typename T>
  VEG_INLINE static void transpose_if_rowmajor(T* ptr,
                                               isize dim,
                                               isize outer_stride)
  {
    (void)ptr, (void)dim, (void)outer_stride;
  }
};
 
template<>
struct ElementAccess<Layout::rowmajor>
{
  template<typename T>
  VEG_INLINE static constexpr auto offset(T* ptr,
                                          isize row,
                                          isize col,
                                          isize outer_stride) noexcept -> T*
  {
    return ptr + (usize(col) + usize(row) * usize(outer_stride));
  }
 
  using NextColStride = Eigen::Stride<0, 0>;
  using NextRowStride = Eigen::InnerStride<Eigen::Dynamic>;
  VEG_INLINE static auto next_col_stride(isize outer_stride) noexcept
    -> NextColStride
  {
    (void)outer_stride;
    return NextColStride{};
  }
  VEG_INLINE static auto next_row_stride(isize outer_stride) noexcept
    -> NextRowStride
  {
    return NextRowStride /* NOLINT(modernize-return-braced-init-list) */ (
      outer_stride);
  }
 
  template<typename T>
  VEG_INLINE static void transpose_if_rowmajor(T* ptr,
                                               isize dim,
                                               isize outer_stride)
  {
    Eigen::Map<                          //
      Eigen::Matrix<                     //
        T,                               //
        Eigen::Dynamic,                  //
        Eigen::Dynamic                   //
        >,                               //
      Eigen::Unaligned,                  //
      Eigen::OuterStride<Eigen::Dynamic> //
      >{
      ptr,
      dim,
      dim,
      Eigen::OuterStride<Eigen::Dynamic>(outer_stride),
    }
      .transposeInPlace();
  }
};
} // namespace detail
 
namespace detail {
template<typename T>
struct unlref
{
  using Type = T;
};
template<typename T>
struct unlref<T&>
{
  using Type = T;
};
 
template<typename T>
auto
is_eigen_matrix_base_impl(Eigen::MatrixBase<T> const volatile*)
  -> proxsuite::linalg::veg::meta::true_type;
auto
is_eigen_matrix_base_impl(void const volatile*)
  -> proxsuite::linalg::veg::meta::false_type;
 
template<typename T>
auto
is_eigen_owning_matrix_base_impl(Eigen::PlainObjectBase<T> const volatile*)
  -> proxsuite::linalg::veg::meta::true_type;
auto
is_eigen_owning_matrix_base_impl(void const volatile*)
  -> proxsuite::linalg::veg::meta::false_type;
 
template<typename... Ts>
using Void = void;
 
template<typename Mat, typename T>
using DataExpr = decltype(static_cast<T*>(VEG_DECLVAL(Mat&).data()));
 
template<typename Dummy,
         typename Fallback,
         template<typename...> class F,
         typename... Ts>
struct DetectedImpl : proxsuite::linalg::veg::meta::false_type
{
  using Type = Fallback;
};
 
template<typename Fallback, template<typename...> class F, typename... Ts>
struct DetectedImpl<Void<F<Ts...>>, Fallback, F, Ts...>
  : proxsuite::linalg::veg::meta::true_type
{
  using Type = F<Ts...>;
};
 
template<typename Fallback, template<typename...> class F, typename... Ts>
using Detected = typename DetectedImpl<void, Fallback, F, Ts...>::Type;
 
template<typename T>
using CompTimeColsImpl =
  proxsuite::linalg::veg::meta::constant<isize, isize(T::ColsAtCompileTime)>;
template<typename T>
using CompTimeRowsImpl =
  proxsuite::linalg::veg::meta::constant<isize, isize(T::RowsAtCompileTime)>;
template<typename T>
using CompTimeInnerStrideImpl =
  proxsuite::linalg::veg::meta::constant<isize,
                                         isize(T::InnerStrideAtCompileTime)>;
template<typename T>
using LayoutImpl = proxsuite::linalg::veg::meta::
  constant<Layout, (bool(T::IsRowMajor) ? rowmajor : colmajor)>;
 
template<typename T, Layout L>
using EigenMatMap = Eigen::Map<      //
  Eigen::Matrix<                     //
    T,                               //
    Eigen::Dynamic,                  //
    Eigen::Dynamic,                  //
    (L == colmajor)                  //
      ? Eigen::ColMajor              //
      : Eigen::RowMajor              //
    > const,                         //
  Eigen::Unaligned,                  //
  Eigen::OuterStride<Eigen::Dynamic> //
  >;
template<typename T, Layout L>
using EigenMatMapMut = Eigen::Map<   //
  Eigen::Matrix<                     //
    T,                               //
    Eigen::Dynamic,                  //
    Eigen::Dynamic,                  //
    (L == colmajor)                  //
      ? Eigen::ColMajor              //
      : Eigen::RowMajor              //
    >,                               //
  Eigen::Unaligned,                  //
  Eigen::OuterStride<Eigen::Dynamic> //
  >;
 
template<typename T, typename Stride>
using EigenVecMap = Eigen::Map< //
  Eigen::Matrix<                //
    T,                          //
    Eigen::Dynamic,             //
    1                           //
    > const,                    //
  Eigen::Unaligned,             //
  Stride                        //
  >;
template<typename T, typename Stride>
using EigenVecMapMut = Eigen::Map< //
  Eigen::Matrix<                   //
    T,                             //
    Eigen::Dynamic,                //
    1                              //
    >,                             //
  Eigen::Unaligned,                //
  Stride                           //
  >;
 
template<typename T, Layout L>
using ColToVec = EigenVecMap<T, typename ElementAccess<L>::NextRowStride>;
template<typename T, Layout L>
using RowToVec = EigenVecMap<T, typename ElementAccess<L>::NextColStride>;
template<typename T, Layout L>
using ColToVecMut = EigenVecMapMut<T, typename ElementAccess<L>::NextRowStride>;
template<typename T, Layout L>
using RowToVecMut = EigenVecMapMut<T, typename ElementAccess<L>::NextColStride>;
 
template<typename T>
using VecMap = EigenVecMap<T, Eigen::Stride<0, 0>>;
template<typename T>
using VecMapMut = EigenVecMapMut<T, Eigen::Stride<0, 0>>;
 
} // namespace detail
 
template<typename T>
using unref = typename detail::unlref<T&>::Type;
 
namespace eigen {
template<typename T>
using CompTimeCols =
  detail::Detected<proxsuite::linalg::veg::meta::constant<isize, 0>,
                   detail::CompTimeColsImpl,
                   T>;
template<typename T>
using CompTimeRows =
  detail::Detected<proxsuite::linalg::veg::meta::constant<isize, 0>,
                   detail::CompTimeRowsImpl,
                   T>;
template<typename T>
using CompTimeInnerStride =
  detail::Detected<proxsuite::linalg::veg::meta::constant<isize, 0>,
                   detail::CompTimeInnerStrideImpl,
                   T>;
template<typename T>
using GetLayout =
  detail::Detected<proxsuite::linalg::veg::meta::
                     constant<Layout, Layout(static_cast<unsigned char>(-1))>,
                   detail::LayoutImpl,
                   T>;
} // namespace eigen
 
namespace concepts {
VEG_DEF_CONCEPT(typename T, rvalue_ref, std::is_rvalue_reference<T>::value);
VEG_DEF_CONCEPT(typename T, lvalue_ref, std::is_lvalue_reference<T>::value);
VEG_DEF_CONCEPT((template<typename...> class F, typename... Ts),
                detected,
                detail::DetectedImpl<void, void, F, Ts...>::value);
 
namespace aux {
VEG_DEF_CONCEPT((typename Mat, typename T),
                has_data_expr,
                LDLT_CONCEPT(detected<detail::DataExpr, Mat, T>));
 
VEG_DEF_CONCEPT((typename Mat),
                matrix_base,
                decltype(detail::is_eigen_matrix_base_impl(
                  static_cast<Mat*>(nullptr)))::value);
 
VEG_DEF_CONCEPT((typename Mat),
                is_plain_object_base,
                decltype(detail::is_eigen_owning_matrix_base_impl(
                  static_cast<Mat*>(nullptr)))::value);
 
VEG_DEF_CONCEPT((typename Mat),
                tmp_matrix,
                (LDLT_CONCEPT(aux::is_plain_object_base<unref<Mat>>) &&
                 !LDLT_CONCEPT(lvalue_ref<Mat>)));
} // namespace aux
 
VEG_DEF_CONCEPT((typename Mat, typename T),
                eigen_view,
                (LDLT_CONCEPT(aux::matrix_base<unref<Mat>>) &&
                 LDLT_CONCEPT(aux::has_data_expr<Mat, T const>)));
 
VEG_DEF_CONCEPT((typename Mat, typename T),
                eigen_view_mut,
                (LDLT_CONCEPT(aux::matrix_base<unref<Mat>>) &&
                 LDLT_CONCEPT(aux::has_data_expr<Mat, T>) &&
                 !LDLT_CONCEPT(aux::tmp_matrix<Mat>)));
 
VEG_DEF_CONCEPT((typename Mat, typename T),
                eigen_strided_vector_view,
                (LDLT_CONCEPT(eigen_view<Mat, T>) &&
                 (eigen::CompTimeCols<unref<Mat>>::value == 1)));
 
VEG_DEF_CONCEPT((typename Mat, typename T),
                eigen_strided_vector_view_mut,
                (LDLT_CONCEPT(eigen_view_mut<Mat, T>) &&
                 (eigen::CompTimeCols<unref<Mat>>::value == 1)));
 
VEG_DEF_CONCEPT((typename Mat, typename T),
                eigen_vector_view,
                (LDLT_CONCEPT(eigen_strided_vector_view<Mat, T>) &&
                 (eigen::CompTimeInnerStride<unref<Mat>>::value == 1)));
 
VEG_DEF_CONCEPT((typename Mat, typename T),
                eigen_vector_view_mut,
                (LDLT_CONCEPT(eigen_strided_vector_view_mut<Mat, T>) &&
                 (eigen::CompTimeInnerStride<unref<Mat>>::value == 1)));
} // namespace concepts
 
inline namespace tags {
VEG_TAG(from_ptr_size, FromPtrSize);
VEG_TAG(from_ptr_size_stride, FromPtrSizeStride);
VEG_TAG(from_ptr_rows_cols_stride, FromPtrRowsColsStride);
VEG_TAG(from_eigen, FromEigen);
} // namespace tags
 
template<typename T>
struct VectorView
{
  T const* data;
  isize dim;
 
  VEG_INLINE
  VectorView(FromPtrSize /*tag*/, T const* _data, isize _dim) noexcept
    : data(_data)
    , dim(_dim)
  {
  }
 
  VEG_TEMPLATE(typename Vec,
               requires(LDLT_CONCEPT(eigen_vector_view<Vec, T>)),
               VEG_INLINE VectorView,
               (/*tag*/, FromEigen),
               (vec, Vec const&)) noexcept
    : data(vec.data())
    , dim(vec.rows())
  {
  }
 
  VEG_INLINE auto ptr(isize index) const noexcept -> T const*
  {
    return data + index;
  }
  VEG_INLINE auto operator()(isize index) const noexcept -> T const&
  {
    return *ptr(index);
  }
  VEG_INLINE auto segment(isize i, isize size) const noexcept -> VectorView
  {
    return {
      from_ptr_size,
      data + i,
      size,
    };
  }
  VEG_INLINE auto to_eigen() const -> detail::VecMap<T>
  {
    return detail::VecMap<T>(data, Eigen::Index(dim));
  }
};
 
template<typename T>
struct VectorViewMut
{
  T* data;
  isize dim;
 
  VEG_INLINE
  VectorViewMut(FromPtrSize /*tag*/, T* _data, isize _dim) noexcept
    : data(_data)
    , dim(_dim)
  {
  }
 
  VEG_TEMPLATE(typename Vec,
               requires(LDLT_CONCEPT(eigen_vector_view_mut<Vec, T>)),
               VEG_INLINE VectorViewMut,
               (/*tag*/, FromEigen),
               (vec, Vec&&)) noexcept
    : data(vec.data())
    , dim(vec.rows())
  {
  }
 
  VEG_INLINE auto as_const() const noexcept -> VectorView<T>
  {
    return {
      from_ptr_size,
      data,
      dim,
    };
  }
  VEG_INLINE auto ptr(isize index) const noexcept -> T* { return data + index; }
  VEG_INLINE auto operator()(isize index) const noexcept -> T&
  {
    return *ptr(index);
  }
  VEG_INLINE auto segment(isize i, isize size) const noexcept -> VectorViewMut
  {
    return {
      from_ptr_size,
      data + i,
      size,
    };
  }
  VEG_INLINE auto to_eigen() const -> detail::VecMapMut<T>
  {
    return detail::VecMapMut<T>(data, Eigen::Index(dim));
  }
};
 
template<typename T>
struct StridedVectorView
{
  T const* data;
  isize dim;
  isize stride;
 
  VEG_INLINE
  StridedVectorView(FromPtrSizeStride /*tag*/,
                    T const* _data,
                    isize _dim,
                    isize _stride) noexcept
    : data(_data)
    , dim(_dim)
    , stride(_stride)
  {
  }
 
  VEG_TEMPLATE(typename Vec,
               requires(LDLT_CONCEPT(eigen_strided_vector_view<Vec, T>)),
               VEG_INLINE StridedVectorView,
               (/*tag*/, FromEigen),
               (vec, Vec const&)) noexcept
    : data(vec.data())
    , dim(vec.rows())
    , stride(vec.innerStride())
  {
  }
 
  VEG_INLINE auto ptr(isize index) const noexcept -> T const*
  {
    return data + stride * index;
  }
  VEG_INLINE auto operator()(isize index) const noexcept -> T const&
  {
    return *ptr(index);
  }
  VEG_INLINE auto segment(isize i, isize size) const noexcept
    -> StridedVectorView
  {
    return {
      from_ptr_size_stride,
      data + stride * i,
      size,
      stride,
    };
  }
  VEG_INLINE auto to_eigen() const
    -> detail::EigenVecMap<T, Eigen::InnerStride<Eigen::Dynamic>>
  {
    return detail::EigenVecMap<T, Eigen::InnerStride<Eigen::Dynamic>>(
      data,
      Eigen::Index(dim),
      Eigen::Index(1),
      Eigen::InnerStride<Eigen::Dynamic>(Eigen::Index(stride)));
  }
};
 
template<typename T>
struct StridedVectorViewMut
{
  T* data;
  isize dim;
  isize stride;
 
  VEG_INLINE
  StridedVectorViewMut(FromPtrSizeStride /*tag*/,
                       T* _data,
                       isize _dim,
                       isize _stride) noexcept
    : data(_data)
    , dim(_dim)
    , stride(_stride)
  {
  }
 
  VEG_TEMPLATE(typename Vec,
               requires(LDLT_CONCEPT(eigen_strided_vector_view_mut<Vec, T>)),
               VEG_INLINE StridedVectorViewMut,
               (/*tag*/, FromEigen),
               (vec, Vec&&)) noexcept
    : data(vec.data())
    , dim(vec.rows())
    , stride(vec.innerStride())
  {
  }
 
  VEG_INLINE auto as_const() const noexcept -> StridedVectorView<T>
  {
    return {
      from_ptr_size_stride,
      data,
      dim,
      stride,
    };
  }
  VEG_INLINE auto ptr(isize index) const noexcept -> T*
  {
    return data + stride * index;
  }
  VEG_INLINE auto operator()(isize index) const noexcept -> T&
  {
    return *ptr(index);
  }
  VEG_INLINE auto segment(isize i, isize size) const noexcept
    -> StridedVectorViewMut
  {
    return {
      from_ptr_size_stride,
      data + stride * i,
      size,
      stride,
    };
  }
  VEG_INLINE auto to_eigen() const
    -> detail::EigenVecMapMut<T, Eigen::InnerStride<Eigen::Dynamic>>
  {
    return detail::EigenVecMapMut<T, Eigen::InnerStride<Eigen::Dynamic>>(
      data,
      Eigen::Index(dim),
      Eigen::Index(1),
      Eigen::InnerStride<Eigen::Dynamic>(Eigen::Index(stride)));
  }
};
 
template<typename T, Layout L>
struct MatrixView
{
  T const* data;
  isize rows;
  isize cols;
  isize outer_stride;
 
  VEG_INLINE MatrixView(FromPtrRowsColsStride /*tag*/,
                        T const* _data,
                        isize _rows,
                        isize _cols,
                        isize _outer_stride) noexcept
    : data(_data)
    , rows(_rows)
    , cols(_cols)
    , outer_stride(_outer_stride)
  {
  }
 
  VEG_TEMPLATE(typename Mat,
               requires(LDLT_CONCEPT(eigen_view<Mat, T>) &&
                        eigen::GetLayout<unref<Mat>>::value == L),
               VEG_INLINE MatrixView,
               (/*tag*/, FromEigen),
               (mat, Mat const&)) noexcept
    : data(mat.data())
    , rows(mat.rows())
    , cols(mat.cols())
    , outer_stride(mat.outerStride())
  {
  }
 
  VEG_INLINE auto ptr(isize row, isize col) const noexcept -> T const*
  {
    return detail::ElementAccess<L>::offset(data, row, col, outer_stride);
  }
  VEG_INLINE auto operator()(isize row, isize col) const noexcept -> T const&
  {
    return *ptr(row, col);
  }
  VEG_INLINE auto block(isize row,
                        isize col,
                        isize nrows,
                        isize ncols) const noexcept -> MatrixView
  {
    return {
      from_ptr_rows_cols_stride,
      detail::ElementAccess<L>::offset(data, row, col, outer_stride),
      nrows,
      ncols,
      outer_stride,
    };
  }
 
private:
  VEG_INLINE auto col_impl(
    proxsuite::linalg::veg::meta::constant<Layout, colmajor> /*tag*/,
    isize c) const noexcept -> VectorView<T>
  {
    return {
      from_ptr_size,
      data + c * outer_stride,
      rows,
    };
  }
  VEG_INLINE auto col_impl(
    proxsuite::linalg::veg::meta::constant<Layout, rowmajor> /*tag*/,
    isize c) const noexcept -> StridedVectorView<T>
  {
    return {
      from_ptr_size_stride,
      data + c,
      rows,
      outer_stride,
    };
  }
 
public:
  VEG_INLINE auto col(isize c) const noexcept -> proxsuite::linalg::veg::meta::
    if_t<(L == colmajor), VectorView<T>, StridedVectorView<T>>
  {
    return col_impl(proxsuite::linalg::veg::meta::constant<Layout, L>{}, c);
  }
  VEG_INLINE auto row(isize r) const noexcept -> proxsuite::linalg::veg::meta::
    if_t<(L == rowmajor), VectorView<T>, StridedVectorView<T>>
  {
    return trans().col(r);
  }
  VEG_INLINE auto trans() const noexcept
    -> MatrixView<T, proxqp::flip_layout(L)>
  {
    return {
      from_ptr_rows_cols_stride, data, cols, rows, outer_stride,
    };
  }
  VEG_INLINE auto to_eigen() const noexcept -> detail::EigenMatMap<T, L>
  {
    return detail::EigenMatMap<T, L>(
      data,
      Eigen::Index(rows),
      Eigen::Index(cols),
      Eigen::OuterStride<Eigen::Dynamic>(Eigen::Index(outer_stride)));
  }
};
 
template<typename T, Layout L>
struct MatrixViewMut
{
  T* data;
  isize rows;
  isize cols;
  isize outer_stride;
 
  VEG_INLINE MatrixViewMut(FromPtrRowsColsStride /*tag*/,
                           T* _data,
                           isize _rows,
                           isize _cols,
                           isize _outer_stride) noexcept
    : data(_data)
    , rows(_rows)
    , cols(_cols)
    , outer_stride(_outer_stride)
  {
  }
 
  VEG_TEMPLATE(typename Mat,
               requires(LDLT_CONCEPT(eigen_view<Mat, T>) &&
                        eigen::GetLayout<unref<Mat>>::value == L),
               VEG_INLINE MatrixViewMut,
               (/*tag*/, FromEigen),
               (mat, Mat&&)) noexcept
    : data(mat.data())
    , rows(mat.rows())
    , cols(mat.cols())
    , outer_stride(mat.outerStride())
  {
  }
 
  VEG_INLINE auto ptr(isize row, isize col) const noexcept -> T*
  {
    return detail::ElementAccess<L>::offset(data, row, col, outer_stride);
  }
  VEG_INLINE auto operator()(isize row, isize col) const noexcept -> T&
  {
    return *ptr(row, col);
  }
  VEG_INLINE auto block(isize row,
                        isize col,
                        isize nrows,
                        isize ncols) const noexcept -> MatrixViewMut
  {
    return {
      from_ptr_rows_cols_stride,
      detail::ElementAccess<L>::offset(data, row, col, outer_stride),
      nrows,
      ncols,
      outer_stride,
    };
  }
 
private:
  VEG_INLINE auto col_impl(
    proxsuite::linalg::veg::meta::constant<Layout, colmajor> /*tag*/,
    isize c) const noexcept -> VectorViewMut<T>
  {
    return {
      from_ptr_size,
      data + c * outer_stride,
      rows,
    };
  }
  VEG_INLINE auto col_impl(
    proxsuite::linalg::veg::meta::constant<Layout, rowmajor> /*tag*/,
    isize c) const noexcept -> StridedVectorViewMut<T>
  {
    return {
      from_ptr_size_stride,
      data + c,
      rows,
      outer_stride,
    };
  }
 
public:
  VEG_INLINE auto col(isize c) const noexcept -> proxsuite::linalg::veg::meta::
    if_t<(L == colmajor), VectorViewMut<T>, StridedVectorViewMut<T>>
  {
    return col_impl(proxsuite::linalg::veg::meta::constant<Layout, L>{}, c);
  }
  VEG_INLINE auto row(isize r) const noexcept -> proxsuite::linalg::veg::meta::
    if_t<(L == rowmajor), VectorViewMut<T>, StridedVectorViewMut<T>>
  {
    return trans().col(r);
  }
  VEG_INLINE auto trans() const noexcept
    -> MatrixViewMut<T, proxqp::flip_layout(L)>
  {
    return {
      from_ptr_rows_cols_stride, data, cols, rows, outer_stride,
    };
  }
  VEG_INLINE auto to_eigen() const noexcept -> detail::EigenMatMapMut<T, L>
  {
    return detail::EigenMatMapMut<T, L>(
      data,
      Eigen::Index(rows),
      Eigen::Index(cols),
      Eigen::OuterStride<Eigen::Dynamic>(Eigen::Index(outer_stride)));
  }
  VEG_INLINE auto as_const() const noexcept -> MatrixView<T, L>
  {
    return {
      from_ptr_rows_cols_stride, data, rows, cols, outer_stride,
    };
  }
};
 
template<typename T>
struct LdltView
{
private:
  MatrixView<T, colmajor> ld;
 
public:
  explicit LdltView(MatrixView<T, colmajor> ld) noexcept
    : ld(ld)
  {
    VEG_DEBUG_ASSERT(ld.rows == ld.cols);
  }
 
  VEG_INLINE auto l() const noexcept -> MatrixView<T, colmajor> { return ld; }
  VEG_INLINE auto d() const noexcept -> StridedVectorView<T>
  {
    return { from_ptr_size_stride, ld.data, ld.rows, ld.outer_stride + 1 };
  }
 
  VEG_INLINE auto head(isize k) const -> LdltView
  {
    return LdltView{ ld.block(0, 0, k, k) };
  }
  VEG_INLINE auto tail(isize k) const -> LdltView
  {
    isize n = ld.rows;
    return LdltView{ ld.block(n - k, n - k, k, k) };
  }
};
template<typename T>
struct LdltViewMut
{
private:
  MatrixViewMut<T, colmajor> ld;
 
public:
  explicit LdltViewMut(MatrixViewMut<T, colmajor> ld) noexcept
    : ld(ld)
  {
    VEG_DEBUG_ASSERT(ld.rows == ld.cols);
  }
 
  VEG_INLINE auto l() const noexcept -> MatrixView<T, colmajor>
  {
    return ld.as_const();
  }
  VEG_INLINE auto l_mut() const noexcept -> MatrixViewMut<T, colmajor>
  {
    return ld;
  }
  VEG_INLINE auto d() const noexcept -> StridedVectorView<T>
  {
    return { from_ptr_size_stride, ld.data, ld.rows, ld.outer_stride + 1 };
  }
  VEG_INLINE auto d_mut() const noexcept -> StridedVectorViewMut<T>
  {
    return { from_ptr_size_stride, ld.data, ld.rows, ld.outer_stride + 1 };
  }
 
  VEG_INLINE auto as_const() const noexcept -> LdltView<T>
  {
    return LdltView<T>{ ld.as_const() };
  }
 
  VEG_INLINE auto head(isize k) const -> LdltViewMut
  {
    return LdltViewMut{ ld.block(0, 0, k, k) };
  }
  VEG_INLINE auto tail(isize k) const -> LdltViewMut
  {
    isize n = ld.rows;
    return LdltViewMut{ ld.block(n - k, n - k, k, k) };
  }
};
 
namespace detail {
template<typename T>
void
noalias_mul_add(MatrixViewMut<T, colmajor> dst,
                MatrixView<T, colmajor> lhs,
                MatrixView<T, colmajor> rhs,
                T factor)
{
 
  if ((dst.cols == 0) || (dst.rows == 0) || (lhs.cols == 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 (dst.cols == 1 && dst.rows == 1) {
    // dot
    auto rhs_col = rhs.col(0);
    auto lhs_row = lhs.row(0);
    auto lhs_as_col = lhs.col(0);
    lhs_as_col.dim = lhs_row.dim;
    if (lhs_row.stride == 1) {
      dst(0, 0) += factor * lhs_as_col.to_eigen().dot(rhs_col.to_eigen());
    } else {
      dst(0, 0) += factor * lhs_row.to_eigen().dot(rhs_col.to_eigen());
    }
  } else if (dst.cols == 1) {
    // gemv
    auto rhs_col = rhs.col(0);
    auto dst_col = dst.col(0);
    dst_col.to_eigen().noalias().operator+=(
      factor * LAZY_PRODUCT(lhs.to_eigen(), rhs_col.to_eigen()));
  }
 
#if !EIGEN_VERSION_AT_LEAST(3, 3, 8)
  else 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 Stride = Eigen::OuterStride<Eigen::Dynamic>;
    using Mat =
      Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor, 20, 20>;
    using MapMut = Eigen::Map<Mat, Eigen::Unaligned, Stride>;
    using Map = Eigen::Map<Mat const, Eigen::Unaligned, Stride>;
 
    MapMut(dst.data, dst.rows, dst.cols, Stride(dst.outer_stride))
      .noalias()
      .operator+=(
        factor *
        LAZY_PRODUCT(
          Map(lhs.data, lhs.rows, lhs.cols, Stride(lhs.outer_stride)),
          Map(rhs.data, rhs.rows, rhs.cols, Stride(rhs.outer_stride))));
  }
#endif
 
  else {
    // gemm
    dst.to_eigen().noalias().operator+=(
      factor * LAZY_PRODUCT(lhs.to_eigen(), rhs.to_eigen()));
  }
 
#undef LAZY_PRODUCT
}
 
template<typename T>
void
noalias_mul_add_vec(VectorViewMut<T> dst,
                    MatrixView<T, colmajor> lhs,
                    VectorView<T> rhs,
                    T factor)
{
  detail::noalias_mul_add<T>(
    {
      from_ptr_rows_cols_stride,
      dst.data,
      dst.dim,
      1,
      0,
    },
    lhs,
    {
      from_ptr_rows_cols_stride,
      rhs.data,
      rhs.dim,
      1,
      0,
    },
    VEG_FWD(factor));
}
 
template<typename T>
auto
dot(StridedVectorView<T> lhs, VectorView<T> rhs) -> T
{
  auto out = T(0);
  detail::noalias_mul_add<T>(
    {
      from_ptr_rows_cols_stride,
      std::addressof(out),
      1,
      1,
      0,
    },
    {
      from_ptr_rows_cols_stride,
      lhs.data,
      1,
      lhs.dim,
      lhs.stride,
    },
    {
      from_ptr_rows_cols_stride,
      rhs.data,
      rhs.dim,
      1,
      0,
    },
    1);
  return out;
}
template<typename T>
void
assign_cwise_prod(VectorViewMut<T> out,
                  StridedVectorView<T> lhs,
                  StridedVectorView<T> rhs)
{
  out.to_eigen() = lhs.to_eigen().cwiseProduct(rhs.to_eigen());
}
template<typename T>
void
assign_scalar_prod(VectorViewMut<T> out, T factor, VectorView<T> in)
{
  out.to_eigen() = in.to_eigen().operator*(factor);
}
 
template<typename T>
void
trans_tr_unit_up_solve_in_place_on_right(MatrixView<T, colmajor> tr,
                                         MatrixViewMut<T, colmajor> rhs)
{
  if (rhs.cols == 1) {
    tr.to_eigen()
      .transpose()
      .template triangularView<Eigen::UnitUpper>()
      .template solveInPlace<Eigen::OnTheRight>(rhs.col(0).to_eigen());
  } else {
    tr.to_eigen()
      .transpose()
      .template triangularView<Eigen::UnitUpper>()
      .template solveInPlace<Eigen::OnTheRight>(rhs.to_eigen());
  }
}
 
template<typename T>
void
apply_diag_inv_on_right(MatrixViewMut<T, colmajor> out,
                        StridedVectorView<T> d,
                        MatrixView<T, colmajor> in)
{
  if (out.cols == 1) {
    out.col(0).to_eigen() =
      in.col(0).to_eigen().operator*(d.to_eigen().asDiagonal().inverse());
  } else {
    out.to_eigen() =
      in.to_eigen().operator*(d.to_eigen().asDiagonal().inverse());
  }
}
template<typename T>
void
apply_diag_on_right(MatrixViewMut<T, colmajor> out,
                    StridedVectorView<T> d,
                    MatrixView<T, colmajor> in)
{
  if (out.cols == 1) {
    out.col(0).to_eigen() =
      in.col(0).to_eigen().operator*(d.to_eigen().asDiagonal());
  } else {
    out.to_eigen() = in.to_eigen().operator*(d.to_eigen().asDiagonal());
  }
}
 
template<typename T>
void
noalias_mul_sub_tr_lo(MatrixViewMut<T, colmajor> out,
                      MatrixView<T, colmajor> lhs,
                      MatrixView<T, rowmajor> rhs)
{
  if (lhs.cols == 1) {
    out.to_eigen().template triangularView<Eigen::Lower>().operator-=(
      lhs.col(0).to_eigen().operator*(
        Eigen::Map<Eigen::Matrix<T, 1, Eigen::Dynamic> const>(
          rhs.data, 1, rhs.cols)));
  } else {
    out.to_eigen().template triangularView<Eigen::Lower>().operator-=(
      lhs.to_eigen().operator*(rhs.to_eigen()));
  }
}
 
} // namespace detail
} // namespace proxqp
} // namespace proxsuite
 
namespace proxsuite {
namespace proxqp {
 
namespace dense {
 
struct EigenAllowAlloc
{
  bool alloc_was_allowed;
  EigenAllowAlloc(EigenAllowAlloc&&) = delete;
  EigenAllowAlloc(EigenAllowAlloc const&) = delete;
  auto operator=(EigenAllowAlloc&&) -> EigenAllowAlloc& = delete;
  auto operator=(EigenAllowAlloc const&) -> EigenAllowAlloc& = delete;
 
#if defined(EIGEN_RUNTIME_NO_MALLOC)
  EigenAllowAlloc() noexcept
    : alloc_was_allowed(Eigen::internal::is_malloc_allowed())
  {
    Eigen::internal::set_is_malloc_allowed(true);
  }
  ~EigenAllowAlloc() noexcept
  {
    Eigen::internal::set_is_malloc_allowed(alloc_was_allowed);
  }
#else
  EigenAllowAlloc() = default;
#endif
};
 
template<typename T>
struct QpView
{
 
  static constexpr Layout layout = rowmajor;
 
  MatrixView<T, layout> H;
  VectorView<T> g;
 
  MatrixView<T, layout> A;
  VectorView<T> b;
  MatrixView<T, layout> C;
  VectorView<T> d;
};
 
template<typename Scalar>
struct QpViewBox
{
  static constexpr Layout layout = rowmajor;
 
  MatrixView<Scalar, layout> H;
  VectorView<Scalar> g;
 
  MatrixView<Scalar, layout> A;
  VectorView<Scalar> b;
  MatrixView<Scalar, layout> C;
  VectorView<Scalar> u;
  VectorView<Scalar> l;
  VectorView<Scalar> I;
  VectorView<Scalar> u_box;
  VectorView<Scalar> l_box;
};
 
template<typename T>
struct QpViewMut
{
  static constexpr Layout layout = rowmajor;
 
  MatrixViewMut<T, layout> H;
  VectorViewMut<T> g;
 
  MatrixViewMut<T, layout> A;
  VectorViewMut<T> b;
  MatrixViewMut<T, layout> C;
  VectorViewMut<T> d;
 
  VEG_INLINE constexpr auto as_const() const noexcept -> QpView<T>
  {
    return {
      H.as_const(), g.as_const(), A.as_const(),
      b.as_const(), C.as_const(), d.as_const(),
    };
  }
};
 
template<typename Scalar>
struct QpViewBoxMut
{
  static constexpr Layout layout = rowmajor;
 
  MatrixViewMut<Scalar, layout> H;
  VectorViewMut<Scalar> g;
 
  MatrixViewMut<Scalar, layout> A;
  VectorViewMut<Scalar> b;
  MatrixViewMut<Scalar, layout> C;
  VectorViewMut<Scalar> u;
  VectorViewMut<Scalar> l;
  VectorViewMut<Scalar> I;
  VectorViewMut<Scalar> l_box;
  VectorViewMut<Scalar> u_box;
 
  VEG_INLINE constexpr auto as_const() const noexcept -> QpViewBox<Scalar>
  {
    return { H.as_const(),     g.as_const(),    A.as_const(), b.as_const(),
             C.as_const(),     u.as_const(),    l.as_const(), I.as_const(),
             u_box.as_const(), l_box.as_const() };
  }
};
 
namespace nb {
struct pow
{
  template<typename T>
  auto operator()(T x, T y) const -> T
  {
    using std::pow;
    return pow(x, y);
  }
};
struct infty_norm
{
  template<typename D>
  auto operator()(Eigen::MatrixBase<D> const& mat) const -> typename D::Scalar
  {
    if (mat.rows() == 0 || mat.cols() == 0) {
      return typename D::Scalar(0);
    } else {
      return mat.template lpNorm<Eigen::Infinity>();
    }
  }
};
struct sqrt
{
  template<typename T>
  auto operator()(T x) const -> T
  {
    using std::sqrt;
    return sqrt(x);
  }
};
struct fabs
{
  template<typename T>
  auto operator()(T x) const -> T
  {
    using std::fabs;
    return fabs(x);
  }
};
} // namespace nb
VEG_NIEBLOID(fabs);
VEG_NIEBLOID(sqrt);
VEG_NIEBLOID(pow);
VEG_NIEBLOID(infty_norm);
} // namespace dense
} // namespace proxqp
} // namespace proxsuite
 
#endif /* end of include guard PROXSUITE_PROXQP_DENSE_VIEWS_HPP */