10#ifndef TLAPACK_LEGACYARRAY_HH
11#define TLAPACK_LEGACYARRAY_HH
28 template <
typename T,
class idx_t, Layout L>
31 template <
typename T,
class idx_t,
class int_t, Direction D>
32 std::true_type is_legacy_type_f(
35 std::false_type is_legacy_type_f(
const void*);
42 decltype(internal::is_legacy_type_f(std::declval<T*>()))::value;
50 template <
typename T,
class idx_t, Layout L>
52 static constexpr Layout value = L;
56 template <
typename T,
class idx_t,
typename int_t, Direction D>
58 static constexpr Layout value = Layout::Strided;
61 template <
class T,
class idx_t,
typename int_t, Direction D>
66 template <
typename T,
class idx_t, Layout layout>
71 template <
class T,
class idx_t,
typename int_t, Direction D>
76 template <
typename T,
class idx_t, Layout layout>
82 template <
class U,
class idx_t, Layout layout>
85 constexpr auto operator()(std::vector<T>&
v, idx_t m, idx_t n)
const
94 template <
class U,
class idx_t, Layout layout,
int m,
int n>
96 static_assert(m >= 0 && n >= 0);
106 template <
class U,
class idx_t,
typename int_t, Direction D>
109 constexpr auto operator()(std::vector<T>&
v, idx_t n)
const
118 template <
class U,
class idx_t,
typename int_t, Direction D,
int n>
120 static_assert(n >= 0);
122 template <
typename T>
134template <
typename T,
class idx_t, Layout layout>
141template <
typename T,
class idx_t, Layout layout>
148template <
typename T,
class idx_t, Layout layout>
155template <
typename T,
class idx_t,
typename int_t, Direction direction>
162template <
typename T,
class idx_t>
163constexpr auto nrows(
const LegacyBandedMatrix<T, idx_t>& A)
noexcept
169template <
typename T,
class idx_t>
170constexpr auto ncols(
const LegacyBandedMatrix<T, idx_t>& A)
noexcept
176template <
typename T,
class idx_t>
177constexpr auto size(
const LegacyBandedMatrix<T, idx_t>& A)
noexcept
183template <
typename T,
class idx_t>
184constexpr auto lowerband(
const LegacyBandedMatrix<T, idx_t>& A)
noexcept
190template <
typename T,
class idx_t>
191constexpr auto upperband(
const LegacyBandedMatrix<T, idx_t>& A)
noexcept
199#define isSlice(SliceSpec) !std::is_convertible<SliceSpec, idx_t>::value
208 typename std::enable_if<isSlice(SliceSpecRow) && isSlice(SliceSpecCol),
212 SliceSpecCol&& cols)
noexcept
214 assert((rows.first >= 0 and (idx_t) rows.first < nrows(A)) ||
215 rows.first == rows.second);
216 assert(rows.second >= 0 and (idx_t) rows.second <= nrows(A));
217 assert(rows.first <= rows.second);
218 assert((cols.first >= 0 and (idx_t) cols.first < ncols(A)) ||
219 cols.first == cols.second);
220 assert(cols.second >= 0 and (idx_t) cols.second <= ncols(A));
221 assert(cols.first <= cols.second);
223 rows.second - rows.first, cols.second - cols.first,
224 (layout == Layout::ColMajor) ? &A.ptr[rows.first + cols.first * A.ldim]
225 : &A.ptr[rows.first * A.ldim + cols.first],
232template <
typename T,
class idx_t, Layout layout,
class SliceSpecCol>
235 SliceSpecCol&& cols)
noexcept
237 assert((cols.first >= 0 and (idx_t) cols.first < ncols(A)) ||
238 cols.first == cols.second);
239 assert(cols.second >= 0 and (idx_t) cols.second <= ncols(A));
240 assert(cols.first <= cols.second);
241 assert(rowIdx >= 0 and rowIdx < nrows(A));
243 cols.second - cols.first,
244 (layout == Layout::ColMajor) ? &A.ptr[rowIdx + cols.first * A.ldim]
245 : &A.ptr[rowIdx * A.ldim + cols.first],
250template <
typename T,
class idx_t, Layout layout,
class SliceSpecRow>
255 assert((rows.first >= 0 and (idx_t) rows.first < nrows(A)) ||
256 rows.first == rows.second);
257 assert(rows.second >= 0 and (idx_t) rows.second <= nrows(A));
258 assert(rows.first <= rows.second);
259 assert(colIdx >= 0 and colIdx < ncols(A));
261 rows.second - rows.first,
262 (layout == Layout::ColMajor) ? &A.ptr[rows.first + colIdx * A.ldim]
263 : &A.ptr[rows.first * A.ldim + colIdx],
268template <
typename T,
class idx_t, Layout layout,
class SliceSpec>
270 SliceSpec&& rows)
noexcept
272 assert((rows.first >= 0 and (idx_t) rows.first < nrows(A)) ||
273 rows.first == rows.second);
274 assert(rows.second >= 0 and (idx_t) rows.second <= nrows(A));
275 assert(rows.first <= rows.second);
277 rows.second - rows.first, A.n,
278 (layout == Layout::ColMajor) ? &A.ptr[rows.first]
279 : &A.ptr[rows.first * A.ldim],
284template <
typename T,
class idx_t>
288 assert(rowIdx >= 0 and rowIdx < nrows(A));
293template <
typename T,
class idx_t>
298 assert(rowIdx >= 0 and rowIdx < nrows(A));
303template <
typename T,
class idx_t, Layout layout,
class SliceSpec>
305 SliceSpec&& cols)
noexcept
307 assert((cols.first >= 0 and (idx_t) cols.first < ncols(A)) ||
308 cols.first == cols.second);
309 assert(cols.second >= 0 and (idx_t) cols.second <= ncols(A));
310 assert(cols.first <= cols.second);
312 A.m, cols.second - cols.first,
313 (layout == Layout::ColMajor) ? &A.ptr[cols.first * A.ldim]
314 : &A.ptr[cols.first],
319template <
typename T,
class idx_t>
323 assert(colIdx >= 0 and colIdx < ncols(A));
328template <
typename T,
class idx_t>
333 assert(colIdx >= 0 and colIdx < ncols(A));
338template <
typename T,
class idx_t, Layout layout>
340 int diagIdx = 0) noexcept
342 assert(diagIdx >= 0 || (idx_t)(-diagIdx) < nrows(A));
343 assert(diagIdx <= 0 || (idx_t)diagIdx < ncols(A));
345 (
layout == Layout::ColMajor)
346 ? ((diagIdx >= 0) ? &A.ptr[diagIdx * A.ldim] : &A.ptr[-diagIdx])
347 : ((diagIdx >= 0) ? &A.ptr[diagIdx] : &A.ptr[-diagIdx * A.ldim]);
348 idx_t n = (diagIdx >= 0) ? std::min(A.m + diagIdx, A.n) - (idx_t)diagIdx
349 : std::min(A.m, A.n - diagIdx) + (idx_t)diagIdx;
355template <
typename T,
class idx_t>
356constexpr auto transpose_view(
362template <
typename T,
class idx_t>
363constexpr auto transpose_view(
377 SliceSpec&& rows)
noexcept
379 assert((rows.first >= 0 and (idx_t) rows.first < size(v)) ||
380 rows.first == rows.second);
381 assert(rows.second >= 0 and (idx_t) rows.second <= size(v));
382 assert(rows.first <= rows.second);
384 rows.second - rows.first, &v.ptr[rows.first * v.inc], v.inc);
390#define isSlice(SliceSpec) !std::is_convertible<SliceSpec, idx_t>::value
399 typename std::enable_if<isSlice(SliceSpecRow) && isSlice(SliceSpecCol),
403 SliceSpecCol&& cols)
noexcept
405 assert((rows.first >= 0 and (idx_t) rows.first < nrows(A)) ||
406 rows.first == rows.second);
407 assert(rows.second >= 0 and (idx_t) rows.second <= nrows(A));
408 assert(rows.first <= rows.second);
409 assert((cols.first >= 0 and (idx_t) cols.first < ncols(A)) ||
410 cols.first == cols.second);
411 assert(cols.second >= 0 and (idx_t) cols.second <= ncols(A));
412 assert(cols.first <= cols.second);
414 rows.second - rows.first, cols.second - cols.first,
415 (layout == Layout::ColMajor) ? &A.ptr[rows.first + cols.first * A.ldim]
416 : &A.ptr[rows.first * A.ldim + cols.first],
423template <
typename T,
class idx_t, Layout layout,
class SliceSpecCol>
426 SliceSpecCol&& cols)
noexcept
428 assert((cols.first >= 0 and (idx_t) cols.first < ncols(A)) ||
429 cols.first == cols.second);
430 assert(cols.second >= 0 and (idx_t) cols.second <= ncols(A));
431 assert(cols.first <= cols.second);
432 assert(rowIdx >= 0 and rowIdx < nrows(A));
434 cols.second - cols.first,
435 (layout == Layout::ColMajor) ? &A.ptr[rowIdx + cols.first * A.ldim]
436 : &A.ptr[rowIdx * A.ldim + cols.first],
441template <
typename T,
class idx_t, Layout layout,
class SliceSpecRow>
446 assert((rows.first >= 0 and (idx_t) rows.first < nrows(A)) ||
447 rows.first == rows.second);
448 assert(rows.second >= 0 and (idx_t) rows.second <= nrows(A));
449 assert(rows.first <= rows.second);
450 assert(colIdx >= 0 and colIdx < ncols(A));
452 rows.second - rows.first,
453 (layout == Layout::ColMajor) ? &A.ptr[rows.first + colIdx * A.ldim]
454 : &A.ptr[rows.first * A.ldim + colIdx],
459template <
typename T,
class idx_t, Layout layout,
class SliceSpec>
461 SliceSpec&& rows)
noexcept
463 assert((rows.first >= 0 and (idx_t) rows.first < nrows(A)) ||
464 rows.first == rows.second);
465 assert(rows.second >= 0 and (idx_t) rows.second <= nrows(A));
466 assert(rows.first <= rows.second);
468 (layout == Layout::ColMajor)
470 : &A.ptr[rows.first * A.ldim],
475template <
typename T,
class idx_t>
479 assert(rowIdx >= 0 and rowIdx < nrows(A));
484template <
typename T,
class idx_t>
489 assert(rowIdx >= 0 and rowIdx < nrows(A));
494template <
typename T,
class idx_t, Layout layout,
class SliceSpec>
496 SliceSpec&& cols)
noexcept
498 assert((cols.first >= 0 and (idx_t) cols.first < ncols(A)) ||
499 cols.first == cols.second);
500 assert(cols.second >= 0 and (idx_t) cols.second <= ncols(A));
501 assert(cols.first <= cols.second);
503 (layout == Layout::ColMajor)
504 ? &A.ptr[cols.first * A.ldim]
505 : &A.ptr[cols.first],
510template <
typename T,
class idx_t>
514 assert(colIdx >= 0 and colIdx < ncols(A));
519template <
typename T,
class idx_t>
524 assert(colIdx >= 0 and colIdx < ncols(A));
529template <
typename T,
class idx_t, Layout layout>
532 assert(diagIdx >= 0 || (idx_t)(-diagIdx) < nrows(A));
533 assert(diagIdx <= 0 || (idx_t)diagIdx < ncols(A));
535 (
layout == Layout::ColMajor)
536 ? ((diagIdx >= 0) ? &A.ptr[diagIdx * A.ldim] : &A.ptr[-diagIdx])
537 : ((diagIdx >= 0) ? &A.ptr[diagIdx] : &A.ptr[-diagIdx * A.ldim]);
538 idx_t n = (diagIdx >= 0) ? std::min(A.m + diagIdx, A.n) - (idx_t)diagIdx
539 : std::min(A.m, A.n - diagIdx) + (idx_t)diagIdx;
545template <
typename T,
class idx_t>
546constexpr auto transpose_view(
551template <
typename T,
class idx_t>
552constexpr auto transpose_view(
565 SliceSpec&& rows)
noexcept
567 assert((rows.first >= 0 and (idx_t) rows.first < size(v)) ||
568 rows.first == rows.second);
569 assert(rows.second >= 0 and (idx_t) rows.second <= size(v));
570 assert(rows.first <= rows.second);
572 rows.second - rows.first, &v.ptr[rows.first * v.inc], v.inc);
576template <
typename T,
class idx_t, Layout layout>
584 const idx_t size = A.m * A.n;
585 const idx_t new_size = m * n;
586 const bool is_contiguous =
588 (layout == Layout::ColMajor && (A.ldim == A.m || A.n <= 1)) ||
589 (layout == Layout::RowMajor && (A.ldim == A.n || A.m <= 1));
593 throw std::domain_error(
"New size is larger than current size");
596 return std::make_pair(
597 matrix_t(m, n, &A.ptr[0]),
598 (layout == Layout::ColMajor)
599 ? matrix_t(size - new_size, 1, &A.ptr[0] + new_size)
600 : matrix_t(1, size - new_size, &A.ptr[0] + new_size));
603 if (m == A.m || n == 0)
604 return std::make_pair(cols(A, std::pair{(idx_t)0, n}),
605 cols(A, std::pair{n, A.n}));
606 else if (n == A.n || m == 0)
607 return std::make_pair(rows(A, std::pair{(idx_t)0, m}),
608 rows(A, std::pair{m, A.m}));
610 throw std::domain_error(
611 "Cannot reshape to non-contiguous matrix into a matrix if both "
612 "the number of rows and columns are different from the new "
616template <
typename T,
class idx_t, Layout layout>
624 const idx_t size = A.m * A.n;
625 const bool is_contiguous =
627 (layout == Layout::ColMajor && (A.ldim == A.m || A.n <= 1)) ||
628 (layout == Layout::RowMajor && (A.ldim == A.n || A.m <= 1));
632 throw std::domain_error(
"New size is larger than current size");
635 return std::make_pair(vector_t(n, &A.ptr[0]),
636 (layout == Layout::ColMajor)
637 ? matrix_t(size - n, 1, &A.ptr[0] + n)
638 : matrix_t(1, size - n, &A.ptr[0] + n));
642 return std::make_pair(vector_t(0, &A.ptr[0]), A);
645 if constexpr (
layout == Layout::ColMajor)
646 return std::make_pair(
647 vector_t(n, &A.ptr[0], 1),
648 matrix_t(A.m, A.n - 1, &A.ptr[0] + A.ldim, A.ldim));
650 return std::make_pair(
651 vector_t(n, &A.ptr[0], A.ldim),
652 matrix_t(A.m, A.n - 1, &A.ptr[0] + 1, A.ldim));
655 if constexpr (
layout == Layout::ColMajor)
656 return std::make_pair(
657 vector_t(n, &A.ptr[0], A.ldim),
658 matrix_t(A.m - 1, A.n, &A.ptr[0] + 1, A.ldim));
660 return std::make_pair(
661 vector_t(n, &A.ptr[0], 1),
662 matrix_t(A.m - 1, A.n, &A.ptr[0] + A.ldim, A.ldim));
665 throw std::domain_error(
666 "Cannot reshape to non-contiguous matrix into a vector if the "
667 "new size is different from the number of rows and columns.");
673template <
typename T,
class idx_t,
typename int_t, Direction direction>
682 const idx_t new_size = m * n;
683 const idx_t s = v.n - new_size;
684 const bool is_contiguous = (v.inc == 1) || (v.n <= 1);
688 throw std::domain_error(
"New size is larger than current size");
689 if (!is_contiguous && m > 1)
690 throw std::domain_error(
691 "New sizes are not compatible with the current vector.");
694 return std::make_pair(matrix_t(m, n, &v.ptr[0]),
695 vector_t(s, &v.ptr[0] + new_size));
698 return std::make_pair(matrix_t(m, n, &v.ptr[0], v.inc),
699 vector_t(s, &v.ptr[0] + new_size * v.inc, v.inc));
702template <
typename T,
class idx_t,
typename int_t, Direction direction>
708 throw std::domain_error(
"New size is larger than current size");
710 return std::make_pair(slice(v, std::pair{(idx_t)0, n}),
711 slice(v, std::pair{n, v.n}));
719 template <
typename T>
720 constexpr bool cast_to_legacy_type =
721 is_legacy_type<T> || is_stdvector_type<T>
722#ifdef TLAPACK_EIGEN_HH
725#ifdef TLAPACK_MDSPAN_HH
732 template <
class matrixA_t,
typename matrixB_t>
737 ((is_legacy_type<matrixA_t> ||
738 is_legacy_type<matrixB_t>)&&cast_to_legacy_type<matrixA_t> &&
739 cast_to_legacy_type<matrixB_t>),
746 static constexpr Layout L =
747 ((LA == Layout::RowMajor) && (LB == Layout::RowMajor))
756 template <
typename vecA_t,
typename vecB_t>
761 ((is_legacy_type<vecA_t> ||
762 is_legacy_type<vecB_t>)&&cast_to_legacy_type<vecA_t> &&
763 cast_to_legacy_type<vecB_t>),
771#if !defined(TLAPACK_EIGEN_HH) && !defined(TLAPACK_MDSPAN_HH)
772 template <
class vecA_t,
class vecB_t>
776 std::enable_if_t<traits::is_stdvector_type<vecA_t> &&
777 traits::is_stdvector_type<vecB_t>,
783 template <
class vecA_t,
class vecB_t>
784 struct vector_type_traits<
787 std::enable_if_t<traits::is_stdvector_type<vecA_t> &&
788 traits::is_stdvector_type<vecB_t>,
800template <
typename T,
class idx_t, Layout layout>
803 return legacy::Matrix<T, idx_t>{
layout, A.m, A.n, A.ptr, A.ldim};
806template <
class T,
class idx_t,
class int_t, Direction direction>
807constexpr auto legacy_matrix(
810 return legacy::Matrix<T, idx_t>{Layout::ColMajor, 1, v.n, v.ptr,
814template <
typename T,
class idx_t,
typename int_t, Direction direction>
815constexpr auto legacy_vector(
818 assert(direction == Direction::Forward || std::is_signed<idx_t>::value ||
821 return legacy::Vector<T, idx_t>{
823 (direction == Direction::Forward) ? idx_t(v.inc) : idx_t(-v.inc)};
constexpr Layout layout
Layout of a matrix or vector.
Definition arrayTraits.hpp:232
Layout
Definition types.hpp:29
@ ColMajor
Column-major layout.
@ RowMajor
Row-major layout.
constexpr auto diag(T &A, int diagIdx=0) noexcept
Get the Diagonal of an Eigen Matrix.
Definition eigen.hpp:576
Concept for types that represent tlapack::Direction.
Concept for matrices that can be converted to a legacy matrix.
Concept for vectors that can be converted to a legacy vector.
constexpr bool is_legacy_type
True if T is a legacy array.
Definition legacyArray.hpp:41
typename traits::real_type_traits< Types..., int >::type real_type
The common real type of the list of types.
Definition scalar_type_traits.hpp:113
Legacy matrix.
Definition LegacyMatrix.hpp:34
Legacy vector.
Definition LegacyVector.hpp:42
Functor for data creation.
Definition arrayTraits.hpp:89
constexpr auto operator()(std::vector< T > &v, idx_t m, idx_t n=1) const
Creates a m-by-n matrix with entries of type T.
Definition arrayTraits.hpp:105
Functor for data creation with static size.
Definition arrayTraits.hpp:141
constexpr auto operator()(T *v) const
Creates a m-by-n matrix or, if n == -1, a vector of size m.
Definition arrayTraits.hpp:157
Complex type traits for the list of types Types.
Definition scalar_type_traits.hpp:145
Trait to determine the layout of a given data structure.
Definition arrayTraits.hpp:75
Matrix type deduction.
Definition arrayTraits.hpp:176
Real type traits for the list of types Types.
Definition scalar_type_traits.hpp:71
Vector type deduction.
Definition arrayTraits.hpp:203