proxsuite 0.7.1
The Advanced Proximal Optimization Toolbox
Loading...
Searching...
No Matches
core.hpp
Go to the documentation of this file.
1
2//
3// Copyright (c) 2022 INRIA
4//
5#ifndef PROXSUITE_LINALG_DENSE_LDLT_CORE_HPP
6#define PROXSUITE_LINALG_DENSE_LDLT_CORE_HPP
7
11
12#if !(defined(__aarch64__) || defined(__PPC64__) || defined(__ppc64__) || \
13 defined(_ARCH_PPC64))
14#include <immintrin.h>
15#endif
16
17#ifdef PROXSUITE_VECTORIZE
18#include <cmath> // to avoid error of the type no member named 'isnan' in namespace 'std';
19#include <simde/x86/avx2.h>
20#include <simde/x86/fma.h>
21#endif
22
23#include <Eigen/Core>
24
25#define LDLT_ID(id) __VEG_PP_CAT(id, __LINE__)
26
27#define __LDLT_TEMP_VEC_IMPL(Type, Name, Rows, Stack, Make) \
28 auto LDLT_ID(vec_storage) = (Stack).Make( \
29 ::proxsuite::linalg::veg::Tag<__VEG_PP_REMOVE_PAREN(Type)>{}, \
30 (Rows), \
31 ::proxsuite::linalg::dense::_detail::align<__VEG_PP_REMOVE_PAREN( \
32 Type)>()); \
33 auto Name /* NOLINT */ = ::Eigen::Map< \
34 ::Eigen::Matrix<__VEG_PP_REMOVE_PAREN(Type), ::Eigen::Dynamic, 1>, \
35 ::Eigen::Unaligned, \
36 ::Eigen::Stride<::Eigen::Dynamic, 1>>{ \
37 LDLT_ID(vec_storage).ptr_mut(), \
38 LDLT_ID(vec_storage).len(), \
39 ::Eigen::Stride<::Eigen::Dynamic, 1>{ \
40 LDLT_ID(vec_storage).len(), \
41 1, \
42 }, \
43 }; \
44 static_assert(true, ".")
45
46#define __LDLT_TEMP_MAT_IMPL(Type, Name, Rows, Cols, Stack, Make) \
47 ::proxsuite::linalg::veg::isize LDLT_ID(rows) = (Rows); \
48 ::proxsuite::linalg::veg::isize LDLT_ID(cols) = (Cols); \
49 ::proxsuite::linalg::veg::isize LDLT_ID(stride) = \
50 ::proxsuite::linalg::dense::_detail::adjusted_stride< \
51 __VEG_PP_REMOVE_PAREN(Type)>(LDLT_ID(rows)); \
52 auto LDLT_ID(vec_storage) = (Stack).Make( \
53 ::proxsuite::linalg::veg::Tag<__VEG_PP_REMOVE_PAREN(Type)>{}, \
54 LDLT_ID(stride) * LDLT_ID(cols), \
55 ::proxsuite::linalg::dense::_detail::align<__VEG_PP_REMOVE_PAREN( \
56 Type)>()); \
57 auto Name /* NOLINT */ = \
58 ::Eigen::Map<::Eigen::Matrix<__VEG_PP_REMOVE_PAREN(Type), \
59 ::Eigen::Dynamic, \
60 ::Eigen::Dynamic, \
61 ::Eigen::ColMajor>, \
62 ::Eigen::Unaligned, \
63 ::Eigen::Stride<::Eigen::Dynamic, 1>>{ \
64 LDLT_ID(vec_storage).ptr_mut(), \
65 LDLT_ID(rows), \
66 LDLT_ID(cols), \
67 ::Eigen::Stride<::Eigen::Dynamic, 1>{ \
68 LDLT_ID(stride), \
69 1, \
70 }, \
71 }; \
72 static_assert(true, ".")
73
74#define LDLT_TEMP_VEC(Type, Name, Rows, Stack) \
75 __LDLT_TEMP_VEC_IMPL(Type, Name, Rows, Stack, make_new)
76#define LDLT_TEMP_VEC_UNINIT(Type, Name, Rows, Stack) \
77 __LDLT_TEMP_VEC_IMPL(Type, Name, Rows, Stack, make_new_for_overwrite)
78
79#define LDLT_TEMP_MAT(Type, Name, Rows, Cols, Stack) \
80 __LDLT_TEMP_MAT_IMPL(Type, Name, Rows, Cols, Stack, make_new)
81#define LDLT_TEMP_MAT_UNINIT(Type, Name, Rows, Cols, Stack) \
82 __LDLT_TEMP_MAT_IMPL(Type, Name, Rows, Cols, Stack, make_new_for_overwrite)
83
84namespace proxsuite {
85namespace linalg {
86namespace dense {
91using f32 = float;
92using f64 = double;
93
94namespace _detail {
95namespace _simd {
96
97#ifdef __clang__
98#define DENSE_LDLT_FP_PRAGMA _Pragma("STDC FP_CONTRACT ON")
99#else
100#define DENSE_LDLT_FP_PRAGMA
101#endif
102
103static_assert(static_cast<unsigned char>(-1) == 255, "char should have 8 bits");
104static_assert(sizeof(f32) == 4, "f32 should be 32 bits");
105static_assert(sizeof(f64) == 8, "f64 should be 64 bits");
106
107#define LDLT_FN_IMPL3(Fn, Prefix, Suffix) \
108 VEG_INLINE static auto Fn(Pack a, Pack b, Pack c) noexcept -> Pack \
109 { \
110 return Pack{ simde_mm##Prefix##_##Fn##_##Suffix( \
111 a.inner, b.inner, c.inner) }; \
112 } \
113 VEG_NOM_SEMICOLON
114
115#define LDLT_ARITHMETIC_IMPL(Prefix, Suffix) \
116 LDLT_FN_IMPL3(fmadd, Prefix, Suffix); /* (a * b + c) */ \
117 LDLT_FN_IMPL3(fnmadd, Prefix, Suffix); /* (-a * b + c) */
118
119#define LDLT_LOAD_STORE(Prefix, Suffix) \
120 VEG_INLINE static auto load_unaligned(ScalarType const* ptr) noexcept \
121 -> Pack \
122 { \
123 return Pack{ simde_mm##Prefix##_loadu_##Suffix(ptr) }; \
124 } \
125 VEG_INLINE static auto broadcast(ScalarType value) noexcept -> Pack \
126 { \
127 return Pack{ simde_mm##Prefix##_set1_##Suffix(value) }; \
128 } \
129 VEG_INLINE void store_unaligned(ScalarType* ptr) const noexcept \
130 { \
131 simde_mm##Prefix##_storeu_##Suffix(ptr, inner); \
132 } \
133 VEG_NOM_SEMICOLON
134
135template<typename T, usize N>
136struct Pack;
137
138template<typename T>
139struct Pack<T, 1>
140{
141 using ScalarType = T;
142
144
145 VEG_INLINE static auto fmadd(Pack a, Pack b, Pack c) noexcept -> Pack
146 {
148 return { a.inner * b.inner + c.inner };
149 }
150 VEG_INLINE static auto fnmadd(Pack a, Pack b, Pack c) noexcept -> Pack
151 {
152 return fmadd({ -a.inner }, b, c);
153 }
154 VEG_INLINE static auto load_unaligned(ScalarType const* ptr) noexcept -> Pack
155 {
156 return { *ptr };
157 }
158 VEG_INLINE static auto broadcast(ScalarType value) noexcept -> Pack
159 {
160 return { value };
161 }
162 VEG_INLINE void store_unaligned(ScalarType* ptr) const noexcept
163 {
164 *ptr = inner;
165 }
166};
167
168#ifdef PROXSUITE_VECTORIZE
169template<>
170struct Pack<f32, 4>
171{
172 using ScalarType = f32;
173
174 simde__m128 inner;
176 LDLT_LOAD_STORE(, ps);
177};
178
179template<>
180struct Pack<f32, 8>
181{
182 using ScalarType = f32;
183
184 simde__m256 inner;
185 LDLT_ARITHMETIC_IMPL(256, ps)
186 LDLT_LOAD_STORE(256, ps);
187};
188
189#ifdef __AVX512F__
190template<>
191struct Pack<f32, 16>
192{
193 using ScalarType = f32;
194
195 __m512 inner;
196 VEG_INLINE static auto fmadd(Pack a, Pack b, Pack c) noexcept -> Pack
197 {
199 return { _mm512_fmadd_ps(a.inner, b.inner, c.inner) };
200 }
201 VEG_INLINE static auto fnmadd(Pack a, Pack b, Pack c) noexcept -> Pack
202 {
203 return { _mm512_fnmadd_ps(a.inner, b.inner, c.inner) };
204 }
205 VEG_INLINE static auto load_unaligned(ScalarType const* ptr) noexcept -> Pack
206 {
207 return { _mm512_loadu_ps(ptr) };
208 }
209 VEG_INLINE static auto broadcast(ScalarType value) noexcept -> Pack
210 {
211 return { _mm512_set1_ps(value) };
212 }
213 VEG_INLINE void store_unaligned(ScalarType* ptr) const noexcept
214 {
215 _mm512_storeu_ps(ptr, inner);
216 }
217};
218#endif
219
220template<>
221struct Pack<f64, 2>
222{
223 using ScalarType = f64;
224
225 simde__m128d inner;
227 LDLT_LOAD_STORE(, pd);
228};
229template<>
230struct Pack<f64, 4>
231{
232 using ScalarType = f64;
233
234 simde__m256d inner;
235 LDLT_ARITHMETIC_IMPL(256, pd)
236 LDLT_LOAD_STORE(256, pd);
237};
238
239#ifdef __AVX512F__
240template<>
241struct Pack<f64, 8>
242{
243 using ScalarType = f64;
244
245 __m512d inner;
246 VEG_INLINE static auto fmadd(Pack a, Pack b, Pack c) noexcept -> Pack
247 {
249 return { _mm512_fmadd_pd(a.inner, b.inner, c.inner) };
250 }
251 VEG_INLINE static auto fnmadd(Pack a, Pack b, Pack c) noexcept -> Pack
252 {
253 return { _mm512_fnmadd_pd(a.inner, b.inner, c.inner) };
254 }
255 VEG_INLINE static auto load_unaligned(ScalarType const* ptr) noexcept -> Pack
256 {
257 return { _mm512_loadu_pd(ptr) };
258 }
259 VEG_INLINE static auto broadcast(ScalarType value) noexcept -> Pack
260 {
261 return { _mm512_set1_pd(value) };
262 }
263 VEG_INLINE void store_unaligned(ScalarType* ptr) const noexcept
264 {
265 _mm512_storeu_pd(ptr, inner);
266 }
267};
268#endif
269
270#endif
271
272template<typename T>
274{
275 static constexpr usize N = 1;
277};
278
279#ifdef PROXSUITE_VECTORIZE
280template<>
281struct NativePackInfo<f32>
282{
283 static constexpr usize N = SIMDE_NATURAL_VECTOR_SIZE / 32;
284 using Type = Pack<f32, N>;
285};
286template<>
287struct NativePackInfo<f64>
288{
289 static constexpr usize N = SIMDE_NATURAL_VECTOR_SIZE / 64;
290 using Type = Pack<f64, N>;
291};
292#endif
293
294template<typename T>
296} // namespace _simd
297} // namespace _detail
298
299namespace _detail {
301template<bool COND, typename T>
303
304template<typename T>
307
308template<typename T>
309constexpr auto
310round_up(T a, T b) noexcept -> T
311{
312 return a + (b - 1) / b * b;
313}
314
315#ifdef PROXSUITE_VECTORIZE
316template<typename T>
317using should_vectorize =
319 VEG_CONCEPT(same<T, f64>)>;
320#else
321template<typename T>
323#endif
324
325template<typename T>
326auto
328{
329#ifndef SIMDE_NATURAL_FLOAT_VECTOR_SIZE
330 isize simd_stride = 1;
331#else
332 isize simd_stride =
333 (SIMDE_NATURAL_VECTOR_SIZE / CHAR_BIT) / isize{ sizeof(T) };
334#endif
336 : n;
337}
338template<typename T>
339auto
340align() noexcept -> isize
341{
342 return isize{ alignof(T) } * _detail::adjusted_stride<T>(1);
343}
344
345struct NoCopy
346{
347 NoCopy() = default;
348 ~NoCopy() = default;
349
350 NoCopy(NoCopy const&) = delete;
351 NoCopy(NoCopy&&) = delete;
352 auto operator=(NoCopy const&) -> NoCopy& = delete;
353 auto operator=(NoCopy&&) -> NoCopy& = delete;
354};
355
356namespace nb {
357struct max2
358{
359 template<typename T>
360 VEG_INLINE constexpr auto operator()(T const& a, T const& b) const -> T const&
361 {
362 return a > b ? a : b;
363 }
364};
365struct min2
366{
367 template<typename T>
368 VEG_INLINE constexpr auto operator()(T a, T b) const -> T
369 {
370 return (a < b) ? a : b;
371 }
372};
373} // namespace nb
376
377template<typename T>
378void
379set_zero(T* dest, usize n)
380{
381 for (usize i = 0; i < n; ++i) {
382 *dest = 0;
383 }
384}
385
386template<typename T>
388 Eigen::Matrix<typename T::Scalar,
389 Eigen::Dynamic,
390 Eigen::Dynamic,
391 bool(T::IsRowMajor) ? Eigen::RowMajor : Eigen::ColMajor>;
392
393template<typename T>
394using OwnedAll =
395 Eigen::Matrix<typename T::Scalar,
396 T::RowsAtCompileTime,
397 T::ColsAtCompileTime,
398 bool(T::IsRowMajor) ? Eigen::RowMajor : Eigen::ColMajor>;
399
400template<typename T>
402 Eigen::Matrix<typename T::Scalar,
403 Eigen::Dynamic,
404 T::ColsAtCompileTime,
405 bool(T::IsRowMajor) ? Eigen::RowMajor : Eigen::ColMajor>;
406
407template<typename T>
409 Eigen::Matrix<typename T::Scalar,
410 T::RowsAtCompileTime,
411 Eigen::Dynamic,
412 bool(T::IsRowMajor) ? Eigen::RowMajor : Eigen::ColMajor>;
413
414template<typename T>
415using OwnedColVector = Eigen::Matrix< //
416 typename T::Scalar,
417 Eigen::Dynamic,
418 1,
419 Eigen::ColMajor>;
420template<typename T>
421using OwnedRowVector = Eigen::Matrix< //
422 typename T::Scalar,
423 1,
424 Eigen::Dynamic,
425 Eigen::RowMajor>;
426
427template<bool ROWMAJOR>
429
430template<>
431struct ElemAddrImpl<false>
432{
433 template<typename T>
434 static auto fn(T* ptr,
435 isize row,
436 isize col,
437 isize outer_stride,
438 isize inner_stride) noexcept -> T*
439 {
440 return ptr + (inner_stride * row + outer_stride * col);
441 }
442};
443template<>
444struct ElemAddrImpl<true>
445{
446 template<typename T>
447 static auto fn(T* ptr,
448 isize row,
449 isize col,
450 isize outer_stride,
451 isize inner_stride) noexcept -> T*
452 {
453 return ptr + (outer_stride * row + inner_stride * col);
454 }
455};
456
457template<bool COLMAJOR>
459
460template<>
462{
463 template<typename T>
464 using Col =
465 Eigen::Map<const_if<ptr_is_const<decltype(VEG_DECLVAL(T&&).data())>::value,
466 OwnedColVector<uncvref_t<T>>>,
467 Eigen::Unaligned,
468 Eigen::InnerStride<uncvref_t<T>::InnerStrideAtCompileTime>>;
469 template<typename T>
470 using Row =
471 Eigen::Map<const_if<ptr_is_const<decltype(VEG_DECLVAL(T&&).data())>::value,
472 OwnedRowVector<uncvref_t<T>>>,
473 Eigen::Unaligned,
474 Eigen::InnerStride<uncvref_t<T>::OuterStrideAtCompileTime>>;
475
476 template<typename T>
477 static auto col(T&& mat, isize col_idx) noexcept -> Col<T>
478 {
479 return {
480 mat.data() + col_idx * mat.outerStride(),
481 mat.rows(),
482 1,
483 Eigen::InnerStride<uncvref_t<T>::InnerStrideAtCompileTime>{
484 mat.innerStride(),
485 },
486 };
487 }
488 template<typename T>
489 static auto row(T&& mat, isize row_idx) noexcept -> Row<T>
490 {
491 return {
492 mat.data() + row_idx * mat.innerStride(),
493 1,
494 mat.cols(),
495 Eigen::InnerStride<uncvref_t<T>::OuterStrideAtCompileTime>{
496 mat.outerStride(),
497 },
498 };
499 }
500};
501template<>
502struct RowColAccessImpl<false>
503{
504 template<typename T>
505 using Col =
506 Eigen::Map<const_if<ptr_is_const<decltype(VEG_DECLVAL(T&&).data())>::value,
507 OwnedColVector<uncvref_t<T>>>,
508 Eigen::Unaligned,
509 Eigen::InnerStride<uncvref_t<T>::OuterStrideAtCompileTime>>;
510 template<typename T>
511 using Row =
512 Eigen::Map<const_if<ptr_is_const<decltype(VEG_DECLVAL(T&&).data())>::value,
513 OwnedRowVector<uncvref_t<T>>>,
514 Eigen::Unaligned,
515 Eigen::InnerStride<uncvref_t<T>::InnerStrideAtCompileTime>>;
516
517 template<typename T>
518 static auto col(T&& mat, isize col_idx) noexcept -> Col<T>
519 {
520 return {
521 mat.data() + col_idx * mat.innerStride(),
522 mat.rows(),
523 1,
524 Eigen::InnerStride<uncvref_t<T>::OuterStrideAtCompileTime>{
525 mat.outerStride(),
526 },
527 };
528 }
529 template<typename T>
530 static auto row(T&& mat, isize row_idx) noexcept -> Row<T>
531 {
532 return {
533 mat.data() + row_idx * mat.outerStride(),
534 1,
535 mat.cols(),
536 Eigen::InnerStride<uncvref_t<T>::InnerStrideAtCompileTime>{
537 mat.innerStride(),
538 },
539 };
540 }
541};
542template<typename T>
543using StrideOf = Eigen::Stride< //
544 T::OuterStrideAtCompileTime,
545 T::InnerStrideAtCompileTime>;
546} // namespace _detail
547
548namespace util {
549template<bool COLMAJOR, typename T>
550auto
552 isize row,
553 isize col,
554 isize outer_stride,
555 isize inner_stride) noexcept -> T*
556{
558 ptr, row, col, outer_stride, inner_stride);
559}
560
561template<typename Mat>
562auto
563matrix_elem_addr(Mat&& mat, isize row, isize col) noexcept
564 -> decltype(mat.data())
565{
566 return util::elem_addr<!bool(
568 mat.data(),
569 row,
570 col,
571 mat.outerStride(),
572 mat.innerStride());
573}
574
575template<typename T>
576auto
577col(T&& mat, isize col_idx) noexcept -> typename _detail::RowColAccessImpl<
579{
580 return _detail::RowColAccessImpl<!bool(
582}
583template<typename T>
584auto
585row(T&& mat, isize row_idx) noexcept -> typename _detail::RowColAccessImpl<
587{
588 return _detail::RowColAccessImpl<!bool(
590}
591
592template<typename Mat>
593auto
594trans(Mat&& mat) noexcept -> Eigen::Map< //
595 _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
596 Eigen::Matrix< //
601 ? Eigen::ColMajor
602 : Eigen::RowMajor>>,
603 Eigen::Unaligned,
605{
606 return {
607 mat.data(),
608 mat.cols(),
609 mat.rows(),
611 mat.outerStride(),
612 mat.innerStride(),
613 },
614 };
615}
616
617template<typename Mat>
618auto
619diagonal(Mat&& mat) noexcept -> Eigen::Map< //
620 _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
621 Eigen::Matrix< //
623 Eigen::Dynamic,
624 1,
625 Eigen::ColMajor>>,
626 Eigen::Unaligned,
627 Eigen::InnerStride<Eigen::Dynamic>>
628{
630 mat.rows() == mat.cols());
631 return { mat.data(),
632 mat.rows(),
633 1,
634 Eigen::InnerStride<Eigen::Dynamic>{ mat.outerStride() + 1 } };
635}
636
637template<typename Mat>
638auto
639submatrix(Mat&& mat,
640 isize row_start,
641 isize col_start,
642 isize nrows,
643 isize ncols) noexcept
644 -> Eigen::Map<_detail::const_if<
645 _detail::ptr_is_const<decltype(mat.data())>::value,
647 Eigen::Unaligned,
649{
650 return {
652 mat.data(), row_start, col_start, mat.outerStride(), mat.innerStride()),
653 nrows,
654 ncols,
656 mat.outerStride(),
657 mat.innerStride(),
658 },
659 };
660}
661
662template<typename Mat>
663auto
664to_view(Mat&& mat) noexcept -> Eigen::Map<
665 _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
667 Eigen::Unaligned,
669{
670 return {
671 mat.data(),
672 mat.rows(),
673 mat.cols(),
675 mat.outerStride(),
676 mat.innerStride(),
677 },
678 };
679}
680
681template<typename Mat>
682auto
683to_view_dyn_rows(Mat&& mat) noexcept -> Eigen::Map<
684 _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
686 Eigen::Unaligned,
688{
689 return {
690 mat.data(),
691 mat.rows(),
692 mat.cols(),
694 mat.outerStride(),
695 mat.innerStride(),
696 },
697 };
698}
699
700template<typename Mat>
701auto
702to_view_dyn_cols(Mat&& mat) noexcept -> Eigen::Map<
703 _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
705 Eigen::Unaligned,
707{
708 return {
709 mat.data(),
710 mat.rows(),
711 mat.cols(),
713 mat.outerStride(),
714 mat.innerStride(),
715 },
716 };
717}
718
719template<typename Mat>
720auto
721to_view_dyn(Mat&& mat) noexcept
722 -> Eigen::Map<_detail::const_if<
723 _detail::ptr_is_const<decltype(mat.data())>::value,
725 Eigen::Unaligned,
727{
728 return {
729 mat.data(),
730 mat.rows(),
731 mat.cols(),
733 mat.outerStride(),
734 mat.innerStride(),
735 },
736 };
737}
738
739template<typename Mat>
740auto
741subrows(Mat&& mat, isize row_start, isize nrows) noexcept -> Eigen::Map<
742 _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
744 Eigen::Unaligned,
746{
747 return {
749 mat.data(), row_start, 0, mat.outerStride(), mat.innerStride()),
750 nrows,
751 mat.cols(),
753 mat.outerStride(),
754 mat.innerStride(),
755 },
756 };
757}
758
759template<typename Mat>
760auto
761subcols(Mat&& mat, isize col_start, isize ncols) noexcept -> Eigen::Map<
762 _detail::const_if<_detail::ptr_is_const<decltype(mat.data())>::value,
764 Eigen::Unaligned,
766{
767 return {
769 mat.data(), 0, col_start, mat.outerStride(), mat.innerStride()),
770 mat.rows(),
771 ncols,
773 mat.outerStride(),
774 mat.innerStride(),
775 },
776 };
777}
778
779} // namespace util
780namespace _detail {
781template<typename Dst, typename Lhs, typename Rhs, typename T>
782void
783noalias_mul_add_impl(Dst dst, Lhs lhs, Rhs rhs, T factor)
784{
785 VEG_ASSERT_ALL_OF(dst.rows() == lhs.rows(),
786 dst.cols() == rhs.cols(),
787 lhs.cols() == rhs.rows());
788
789 isize nrows = dst.rows();
790 isize ncols = dst.cols();
791 isize depth = lhs.cols();
792
793 if (nrows == 0 || ncols == 0 || depth == 0) {
794 return;
795 }
796
797#if !EIGEN_VERSION_AT_LEAST(3, 3, 8)
798#define LAZY_PRODUCT(a, b) a.lazyProduct(b)
799#else
800#define LAZY_PRODUCT(a, b) a.operator*(b)
801#endif
802
803#if !EIGEN_VERSION_AT_LEAST(3, 3, 8)
804 if ((dst.rows() < 20) && (dst.cols() < 20) && (rhs.rows() < 20)) {
805 // gemm
806 // workaround for eigen 3.3.7 bug:
807 // https://gitlab.com/libeigen/eigen/-/issues/1562
808 using Mat =
809 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor, 20, 20>;
810 using MapMut =
811 Eigen::Map<Mat, Eigen::Unaligned, Eigen::OuterStride<Eigen::Dynamic>>;
812
813 auto dst_ =
814 MapMut(dst.data(), dst.rows(), dst.cols(), { dst.outerStride() });
815 dst_.noalias().operator+=(factor * LAZY_PRODUCT(lhs, rhs));
816 } else
817#endif
818 {
819 dst.noalias().operator+=(factor * LAZY_PRODUCT(lhs, rhs));
820 }
821
822#undef LAZY_PRODUCT
823}
824} // namespace _detail
825namespace util {
826template<typename Dst, typename Lhs, typename Rhs, typename T>
827void
828noalias_mul_add(Dst&& dst, Lhs const& lhs, Rhs const& rhs, T factor)
829{
831 util::submatrix(dst, 0, 0, dst.rows(), dst.cols()),
832 util::submatrix(lhs, 0, 0, lhs.rows(), lhs.cols()),
833 util::submatrix(rhs, 0, 0, rhs.rows(), rhs.cols()),
834 factor);
835}
836} // namespace util
837template<typename T>
838auto
839temp_mat_req(proxsuite::linalg::veg::Tag<T> /*tag*/,
840 isize rows,
842{
843 return {
844 _detail::adjusted_stride<T>(rows) * cols * isize{ sizeof(T) },
846 };
847}
848
849template<typename T>
850auto
851temp_vec_req(proxsuite::linalg::veg::Tag<T> /*tag*/, isize rows) noexcept
853{
854 return {
855 rows * isize{ sizeof(T) },
857 };
858}
859} // namespace dense
860} // namespace linalg
861} // namespace proxsuite
862
863#endif /* end of include guard PROXSUITE_LINALG_DENSE_LDLT_CORE_HPP */
#define VEG_DEBUG_ASSERT(...)
Definition assert.hpp:38
#define VEG_ASSERT_ALL_OF(...)
#define DENSE_LDLT_FP_PRAGMA
Definition core.hpp:100
#define LAZY_PRODUCT(a, b)
#define LDLT_LOAD_STORE(Prefix, Suffix)
Definition core.hpp:119
#define LDLT_ARITHMETIC_IMPL(Prefix, Suffix)
Definition core.hpp:115
#define VEG_CONCEPT(...)
Definition macros.hpp:1241
#define VEG_NIEBLOID(Name)
Definition macros.hpp:545
#define VEG_INLINE
Definition macros.hpp:118
#define VEG_DECLVAL(...)
Definition macros.hpp:131
typename NativePackInfo< T >::Type NativePack
Definition core.hpp:295
Eigen::Matrix< typename T::Scalar, Eigen::Dynamic, Eigen::Dynamic, bool(T::IsRowMajor) ? Eigen::RowMajor :Eigen::ColMajor > OwnedMatrix
Definition core.hpp:387
constexpr auto round_up(T a, T b) noexcept -> T
Definition core.hpp:310
void set_zero(T *dest, usize n)
Definition core.hpp:379
void noalias_mul_add_impl(Dst dst, Lhs lhs, Rhs rhs, T factor)
Definition core.hpp:783
auto adjusted_stride(isize n) noexcept -> isize
Definition core.hpp:327
Eigen::Matrix< typename T::Scalar, Eigen::Dynamic, 1, Eigen::ColMajor > OwnedColVector
Definition core.hpp:415
Eigen::Matrix< typename T::Scalar, T::RowsAtCompileTime, Eigen::Dynamic, bool(T::IsRowMajor) ? Eigen::RowMajor :Eigen::ColMajor > OwnedCols
Definition core.hpp:408
Eigen::Matrix< typename T::Scalar, Eigen::Dynamic, T::ColsAtCompileTime, bool(T::IsRowMajor) ? Eigen::RowMajor :Eigen::ColMajor > OwnedRows
Definition core.hpp:401
Eigen::Stride< T::OuterStrideAtCompileTime, T::InnerStrideAtCompileTime > StrideOf
Definition core.hpp:543
proxsuite::linalg::veg::meta::if_t< COND, T const, T > const_if
Definition core.hpp:302
Eigen::Matrix< typename T::Scalar, 1, Eigen::Dynamic, Eigen::RowMajor > OwnedRowVector
Definition core.hpp:421
Eigen::Matrix< typename T::Scalar, T::RowsAtCompileTime, T::ColsAtCompileTime, bool(T::IsRowMajor) ? Eigen::RowMajor :Eigen::ColMajor > OwnedAll
Definition core.hpp:394
proxsuite::linalg::veg::meta::bool_constant< false > should_vectorize
Definition core.hpp:322
auto align() noexcept -> isize
Definition core.hpp:340
auto matrix_elem_addr(Mat &&mat, isize row, isize col) noexcept -> decltype(mat.data())
Definition core.hpp:563
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 > > >
Definition core.hpp:594
void noalias_mul_add(Dst &&dst, Lhs const &lhs, Rhs const &rhs, T factor)
Definition core.hpp:828
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 > > >
Definition core.hpp:683
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 > > >
Definition core.hpp:761
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 > >
Definition core.hpp:619
auto row(T &&mat, isize row_idx) noexcept -> typename _detail::RowColAccessImpl< !bool(proxsuite::linalg::veg::uncvref_t< T >::IsRowMajor)>::template Row< T >
Definition core.hpp:585
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 > > >
Definition core.hpp:741
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 > > >
Definition core.hpp:639
auto elem_addr(T *ptr, isize row, isize col, isize outer_stride, isize inner_stride) noexcept -> T *
Definition core.hpp:551
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 > > >
Definition core.hpp:721
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 > > >
Definition core.hpp:664
auto col(T &&mat, isize col_idx) noexcept -> typename _detail::RowColAccessImpl< !bool(proxsuite::linalg::veg::uncvref_t< T >::IsRowMajor)>::template Col< T >
Definition core.hpp:577
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 > > >
Definition core.hpp:702
auto temp_vec_req(proxsuite::linalg::veg::Tag< T >, isize rows) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
Definition core.hpp:851
auto temp_mat_req(proxsuite::linalg::veg::Tag< T >, isize rows, isize cols) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
Definition core.hpp:839
typename _detail::_meta::uncvlref< T & >::type uncvref_t
Definition macros.hpp:887
typename _detail::_meta::is_pointer< T >::type unptr_t
Definition core.hpp:270
typename _detail::_meta::conditional_< B >::template type< T, F > if_t
Definition core.hpp:83
std::uint32_t u32
Definition typedefs.hpp:48
_detail::_meta::make_signed< usize >::Type isize
Definition typedefs.hpp:43
decltype(sizeof(0)) usize
Definition macros.hpp:702
static auto fn(T *ptr, isize row, isize col, isize outer_stride, isize inner_stride) noexcept -> T *
Definition core.hpp:434
static auto fn(T *ptr, isize row, isize col, isize outer_stride, isize inner_stride) noexcept -> T *
Definition core.hpp:447
auto operator=(NoCopy &&) -> NoCopy &=delete
auto operator=(NoCopy const &) -> NoCopy &=delete
static auto col(T &&mat, isize col_idx) noexcept -> Col< T >
Definition core.hpp:518
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 > > Row
Definition core.hpp:511
static auto row(T &&mat, isize row_idx) noexcept -> Row< T >
Definition core.hpp:530
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 > > Col
Definition core.hpp:505
static auto col(T &&mat, isize col_idx) noexcept -> Col< T >
Definition core.hpp:477
static auto row(T &&mat, isize row_idx) noexcept -> Row< T >
Definition core.hpp:489
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 > > Col
Definition core.hpp:464
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 > > Row
Definition core.hpp:470
VEG_INLINE void store_unaligned(ScalarType *ptr) const noexcept
Definition core.hpp:162
static VEG_INLINE auto load_unaligned(ScalarType const *ptr) noexcept -> Pack
Definition core.hpp:154
static VEG_INLINE auto fmadd(Pack a, Pack b, Pack c) noexcept -> Pack
Definition core.hpp:145
static VEG_INLINE auto broadcast(ScalarType value) noexcept -> Pack
Definition core.hpp:158
static VEG_INLINE auto fnmadd(Pack a, Pack b, Pack c) noexcept -> Pack
Definition core.hpp:150
VEG_INLINE constexpr auto operator()(T const &a, T const &b) const -> T const &
Definition core.hpp:360
VEG_INLINE constexpr auto operator()(T a, T b) const -> T
Definition core.hpp:368