<T>LAPACK 0.1.1
C++ Template Linear Algebra PACKage
Loading...
Searching...
No Matches
eigen.hpp
Go to the documentation of this file.
1
3//
4// Copyright (c) 2021-2023, University of Colorado Denver. All rights reserved.
5//
6// This file is part of <T>LAPACK.
7// <T>LAPACK is free software: you can redistribute it and/or modify it under
8// the terms of the BSD 3-Clause license. See the accompanying LICENSE file.
9
10#ifndef TLAPACK_EIGEN_HH
11#define TLAPACK_EIGEN_HH
12
13#include <Eigen/Core>
14#include <cassert>
15
16#include "tlapack/base/arrayTraits.hpp"
17#include "tlapack/plugins/stdvector.hpp"
18
19namespace tlapack {
20
21// -----------------------------------------------------------------------------
22// Helpers
23
24namespace traits {
25 namespace internal {
26 // Auxiliary constexpr routines
27
28 template <class Derived>
29 std::true_type is_eigen_dense_f(const Eigen::DenseBase<Derived>*);
30 std::false_type is_eigen_dense_f(const void*);
31
32 template <class Derived>
33 std::true_type is_eigen_matrix_f(const Eigen::MatrixBase<Derived>*);
34 std::false_type is_eigen_matrix_f(const void*);
35
36 template <class XprType, int BlockRows, int BlockCols, bool InnerPanel>
37 std::true_type is_eigen_block_f(
38 const Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>*);
39 std::false_type is_eigen_block_f(const void*);
40
43 template <class T>
44 constexpr bool is_eigen_dense =
45 decltype(is_eigen_dense_f(std::declval<T*>()))::value;
46
49 template <class T>
50 constexpr bool is_eigen_matrix =
51 decltype(is_eigen_matrix_f(std::declval<T*>()))::value;
52
55 template <class T>
56 constexpr bool is_eigen_block =
57 decltype(is_eigen_block_f(std::declval<T*>()))::value;
58
59 template <class>
60 struct isStdComplex : public std::false_type {};
61 template <class T>
62 struct isStdComplex<std::complex<T>> : public std::true_type {};
63
64 } // namespace internal
65
68 template <class T>
69 constexpr bool is_eigen_type = internal::is_eigen_dense<T>;
70
71} // namespace traits
72
73// -----------------------------------------------------------------------------
74// Data traits
75
76namespace traits {
77
79 template <class matrix_t>
80 struct layout_trait<
81 matrix_t,
82 typename std::enable_if<is_eigen_type<matrix_t> &&
83 (matrix_t::InnerStrideAtCompileTime == 1 ||
84 matrix_t::OuterStrideAtCompileTime == 1),
85 int>::type> {
86 static constexpr Layout value =
87 (matrix_t::IsVectorAtCompileTime)
88 ? Layout::Strided
89 : ((matrix_t::IsRowMajor) ? Layout::RowMajor
90 : Layout::ColMajor);
91 };
92
93 template <class matrix_t>
94 struct real_type_traits<
95 matrix_t,
96 typename std::enable_if<is_eigen_type<matrix_t>, int>::type> {
97 using type = Eigen::Matrix<real_type<typename matrix_t::Scalar>,
98 matrix_t::RowsAtCompileTime,
99 matrix_t::ColsAtCompileTime,
100 matrix_t::IsRowMajor ? Eigen::RowMajor
101 : Eigen::ColMajor,
102 matrix_t::MaxRowsAtCompileTime,
103 matrix_t::MaxColsAtCompileTime>;
104 };
105
106 template <class matrix_t>
107 struct complex_type_traits<
108 matrix_t,
109 typename std::enable_if<is_eigen_type<matrix_t>, int>::type> {
110 using type = Eigen::Matrix<complex_type<typename matrix_t::Scalar>,
111 matrix_t::RowsAtCompileTime,
112 matrix_t::ColsAtCompileTime,
113 matrix_t::IsRowMajor ? Eigen::RowMajor
114 : Eigen::ColMajor,
115 matrix_t::MaxRowsAtCompileTime,
116 matrix_t::MaxColsAtCompileTime>;
117 };
118
120 template <typename U>
121 struct CreateFunctor<U,
122 typename std::enable_if<is_eigen_type<U>, int>::type> {
123 static constexpr int Rows_ =
124 (U::RowsAtCompileTime == 1) ? 1 : Eigen::Dynamic;
125 static constexpr int Cols_ =
126 (U::ColsAtCompileTime == 1) ? 1 : Eigen::Dynamic;
127 static constexpr int Options_ =
128 (U::IsRowMajor) ? Eigen::RowMajor : Eigen::ColMajor;
129
130 template <typename T>
131 constexpr auto operator()(std::vector<T>& v,
132 Eigen::Index m,
133 Eigen::Index n = 1) const
134 {
135 assert(m >= 0 && n >= 0);
136 v.resize(0);
137 return Eigen::Matrix<T, Rows_, Cols_, Options_>(m, n);
138 }
139 };
140
142 template <typename U, int m, int n>
143 struct CreateStaticFunctor<
144 U,
145 m,
146 n,
147 typename std::enable_if<is_eigen_type<U>, int>::type> {
148 static constexpr int Options_ =
149 (U::IsRowMajor) ? Eigen::RowMajor : Eigen::ColMajor;
150 static_assert(m >= 0 && n >= -1);
151
152 template <typename T>
153 constexpr auto operator()(T* v) const
154 {
155 constexpr int Cols_ = (n == -1) ? 1 : n;
156 return Eigen::Matrix<T, m, Cols_, Options_, m, Cols_>();
157 }
158 };
159} // namespace traits
160
161// -----------------------------------------------------------------------------
162// Data descriptors for Eigen datatypes
163
164// Size
165template <
166 class T,
167 int Rows,
168 int Cols,
169 int Options,
170 int MaxRows,
171 int MaxCols
172#if __cplusplus >= 201703L
173 // Avoids conflict with std::size
174 ,
175 std::enable_if_t<!(traits::internal::isStdComplex<std::decay_t<T>>::value),
176 int> = 0
177#endif
178 >
179constexpr auto size(
180 const Eigen::Matrix<T, Rows, Cols, Options, MaxRows, MaxCols>& x) noexcept
181{
182 return x.size();
183}
184template <class Derived>
185constexpr auto size(const Eigen::EigenBase<Derived>& x) noexcept
186{
187 return x.size();
188}
189// Number of rows
190template <class T>
191constexpr auto nrows(const Eigen::EigenBase<T>& x) noexcept
192{
193 return x.rows();
194}
195// Number of columns
196template <class T>
197constexpr auto ncols(const Eigen::EigenBase<T>& x) noexcept
198{
199 return x.cols();
200}
201
202// -----------------------------------------------------------------------------
203/*
204 * It is important to separate block operations between:
205 * 1. Data types that do not derive from Eigen::Block
206 * 2. Data types that derive from Eigen::Block
207 *
208 * This is done to avoid the creation of classes like:
209 * Block< Block< Block< MaxtrixXd > > >
210 * which would increase the number of template specializations in <T>LAPACK.
211 *
212 * Moreover, recursive algorithms call themselves using blocks of thier
213 * matrices. This creates a dead lock at compile time if the sizes are
214 * given at runtime, and the compiler cannot deduce where to stop the
215 * recursion for the template specialization.
216 */
217
218#define isSlice(SliceSpec) !std::is_convertible<SliceSpec, Eigen::Index>::value
219
220// Block operations for Eigen::Dense that are not derived from Eigen::Block
221
222template <
223 class T,
224 class SliceSpecRow,
225 class SliceSpecCol,
226 typename std::enable_if<isSlice(SliceSpecRow) && isSlice(SliceSpecCol) &&
227 traits::is_eigen_type<T> &&
228 !traits::internal::is_eigen_block<T>,
229 int>::type = 0>
230constexpr auto slice(T& A, SliceSpecRow&& rows, SliceSpecCol&& cols) noexcept
231{
232 return A.block(rows.first, cols.first, rows.second - rows.first,
233 cols.second - cols.first);
234}
235
236template <class T,
237 typename SliceSpecCol,
238 typename std::enable_if<traits::is_eigen_type<T> &&
239 !traits::internal::is_eigen_block<T>,
240 int>::type = 0>
241constexpr auto slice(T& A, Eigen::Index rowIdx, SliceSpecCol&& cols) noexcept
242{
243 return A.row(rowIdx).segment(cols.first, cols.second - cols.first);
244}
245
246template <class T,
247 typename SliceSpecRow,
248 typename std::enable_if<traits::is_eigen_type<T> &&
249 !traits::internal::is_eigen_block<T>,
250 int>::type = 0>
251constexpr auto slice(T& A, SliceSpecRow&& rows, Eigen::Index colIdx) noexcept
252{
253 return A.col(colIdx).segment(rows.first, rows.second - rows.first);
254}
255
256template <class T,
257 typename SliceSpec,
258 typename std::enable_if<traits::is_eigen_type<T> &&
259 !traits::internal::is_eigen_block<T>,
260 int>::type = 0>
261constexpr auto slice(T& x, SliceSpec&& range) noexcept
262{
263 return x.segment(range.first, range.second - range.first);
264}
265
266template <class T,
267 typename SliceSpec,
268 typename std::enable_if<traits::is_eigen_type<T> &&
269 !traits::internal::is_eigen_block<T>,
270 int>::type = 0>
271constexpr auto rows(T& A, SliceSpec&& rows) noexcept
272{
273 return A.middleRows(rows.first, rows.second - rows.first);
274}
275
276template <class T,
277 typename std::enable_if<traits::is_eigen_type<T> &&
278 !traits::internal::is_eigen_block<T>,
279 int>::type = 0>
280constexpr auto row(T& A, Eigen::Index rowIdx) noexcept
281{
282 return A.row(rowIdx);
283}
284
285template <class T,
286 typename SliceSpec,
287 typename std::enable_if<traits::is_eigen_type<T> &&
288 !traits::internal::is_eigen_block<T>,
289 int>::type = 0>
290constexpr auto cols(T& A, SliceSpec&& cols) noexcept
291{
292 return A.middleCols(cols.first, cols.second - cols.first);
293}
294
295template <class T,
296 typename std::enable_if<traits::is_eigen_type<T> &&
297 !traits::internal::is_eigen_block<T>,
298 int>::type = 0>
299constexpr auto col(T& A, Eigen::Index colIdx) noexcept
300{
301 return A.col(colIdx);
302}
303
304// Block operations for Eigen::Dense that are derived from Eigen::Block
305
306template <
307 class XprType,
308 int BlockRows,
309 int BlockCols,
310 bool InnerPanel,
311 class SliceSpecRow,
312 class SliceSpecCol,
313 typename std::enable_if<isSlice(SliceSpecRow) && isSlice(SliceSpecCol),
314 int>::type = 0>
315constexpr auto slice(Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
316 SliceSpecRow&& rows,
317 SliceSpecCol&& cols) noexcept
318{
319 assert(rows.second <= A.rows());
320 assert(cols.second <= A.cols());
321
322 return Eigen::Block<XprType, Eigen::Dynamic, Eigen::Dynamic, InnerPanel>(
323 A.nestedExpression(), A.startRow() + rows.first,
324 A.startCol() + cols.first, rows.second - rows.first,
325 cols.second - cols.first);
326}
327
328template <class XprType,
329 int BlockRows,
330 int BlockCols,
331 bool InnerPanel,
332 typename SliceSpecCol>
333constexpr auto slice(Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
334 Eigen::Index rowIdx,
335 SliceSpecCol&& cols) noexcept
336{
337 assert(rowIdx < A.rows());
338 assert(cols.second <= A.cols());
339
340 return Eigen::Block<XprType, 1, Eigen::Dynamic, InnerPanel>(
341 A.nestedExpression(), A.startRow() + rowIdx, A.startCol() + cols.first,
342 1, cols.second - cols.first);
343}
344
345template <class XprType,
346 int BlockRows,
347 int BlockCols,
348 bool InnerPanel,
349 typename SliceSpecRow>
350constexpr auto slice(Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
351 SliceSpecRow&& rows,
352 Eigen::Index colIdx) noexcept
353{
354 assert(rows.second <= A.rows());
355 assert(colIdx < A.cols());
356
357 return Eigen::Block<XprType, Eigen::Dynamic, 1, InnerPanel>(
358 A.nestedExpression(), A.startRow() + rows.first, A.startCol() + colIdx,
359 rows.second - rows.first, 1);
360}
361
362template <class XprType, int BlockRows, bool InnerPanel, typename SliceSpec>
363constexpr auto slice(Eigen::Block<XprType, BlockRows, 1, InnerPanel>& x,
364 SliceSpec&& range) noexcept
365{
366 assert(range.second <= x.size());
367
368 return Eigen::Block<XprType, Eigen::Dynamic, 1, InnerPanel>(
369 x.nestedExpression(), x.startRow() + range.first, x.startCol(),
370 range.second - range.first, 1);
371}
372
373template <class XprType, int BlockCols, bool InnerPanel, typename SliceSpec>
374constexpr auto slice(Eigen::Block<XprType, 1, BlockCols, InnerPanel>& x,
375 SliceSpec&& range) noexcept
376{
377 assert(range.second <= x.size());
378
379 return Eigen::Block<XprType, 1, Eigen::Dynamic, InnerPanel>(
380 x.nestedExpression(), x.startRow(), x.startCol() + range.first, 1,
381 range.second - range.first);
382}
383
384template <class XprType,
385 int BlockRows,
386 int BlockCols,
387 bool InnerPanel,
388 typename SliceSpec>
389constexpr auto rows(Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
390 SliceSpec&& rows) noexcept
391{
392 assert(rows.second <= A.rows());
393
394 return Eigen::Block<XprType, Eigen::Dynamic, BlockCols, InnerPanel>(
395 A.nestedExpression(), A.startRow() + rows.first, A.startCol(),
396 rows.second - rows.first, A.cols());
397}
398
399template <class XprType, int BlockRows, int BlockCols, bool InnerPanel>
400constexpr auto row(Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
401 Eigen::Index rowIdx) noexcept
402{
403 assert(rowIdx < A.rows());
404
405 return Eigen::Block<XprType, 1, BlockCols, InnerPanel>(
406 A.nestedExpression(), A.startRow() + rowIdx, A.startCol(), 1, A.cols());
407}
408
409template <class XprType,
410 int BlockRows,
411 int BlockCols,
412 bool InnerPanel,
413 typename SliceSpec>
414constexpr auto cols(Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
415 SliceSpec&& cols) noexcept
416{
417 assert(cols.second <= A.cols());
418
419 return Eigen::Block<XprType, BlockRows, Eigen::Dynamic, InnerPanel>(
420 A.nestedExpression(), A.startRow(), A.startCol() + cols.first, A.rows(),
421 cols.second - cols.first);
422}
423
424template <class XprType, int BlockRows, int BlockCols, bool InnerPanel>
425constexpr auto col(Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
426 Eigen::Index colIdx) noexcept
427{
428 assert(colIdx < A.cols());
429
430 return Eigen::Block<XprType, BlockRows, 1, InnerPanel>(
431 A.nestedExpression(), A.startRow(), A.startCol() + colIdx, A.rows(), 1);
432}
433
434template <
435 class XprType,
436 int BlockRows,
437 int BlockCols,
438 bool InnerPanel,
439 class SliceSpecRow,
440 class SliceSpecCol,
441 typename std::enable_if<isSlice(SliceSpecRow) && isSlice(SliceSpecCol),
442 int>::type = 0>
443constexpr auto slice(
444 const Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
445 SliceSpecRow&& rows,
446 SliceSpecCol&& cols) noexcept
447{
448 assert(rows.second <= A.rows());
449 assert(cols.second <= A.cols());
450
451 return Eigen::Block<const XprType, Eigen::Dynamic, Eigen::Dynamic,
452 InnerPanel>(
453 A.nestedExpression(), A.startRow() + rows.first,
454 A.startCol() + cols.first, rows.second - rows.first,
455 cols.second - cols.first);
456}
457
458template <class XprType,
459 int BlockRows,
460 int BlockCols,
461 bool InnerPanel,
462 typename SliceSpecCol>
463constexpr auto slice(
464 const Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
465 Eigen::Index rowIdx,
466 SliceSpecCol&& cols) noexcept
467{
468 assert(rowIdx < A.rows());
469 assert(cols.second <= A.cols());
470
471 return Eigen::Block<const XprType, 1, Eigen::Dynamic, InnerPanel>(
472 A.nestedExpression(), A.startRow() + rowIdx, A.startCol() + cols.first,
473 1, cols.second - cols.first);
474}
475
476template <class XprType,
477 int BlockRows,
478 int BlockCols,
479 bool InnerPanel,
480 typename SliceSpecRow>
481constexpr auto slice(
482 const Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
483 SliceSpecRow&& rows,
484 Eigen::Index colIdx) noexcept
485{
486 assert(rows.second <= A.rows());
487 assert(colIdx < A.cols());
488
489 return Eigen::Block<const XprType, Eigen::Dynamic, 1, InnerPanel>(
490 A.nestedExpression(), A.startRow() + rows.first, A.startCol() + colIdx,
491 rows.second - rows.first, 1);
492}
493
494template <class XprType, int BlockRows, bool InnerPanel, typename SliceSpec>
495constexpr auto slice(const Eigen::Block<XprType, BlockRows, 1, InnerPanel>& x,
496 SliceSpec&& range) noexcept
497{
498 assert(range.second <= x.size());
499
500 return Eigen::Block<const XprType, Eigen::Dynamic, 1, InnerPanel>(
501 x.nestedExpression(), x.startRow() + range.first, x.startCol(),
502 range.second - range.first, 1);
503}
504
505template <class XprType, int BlockCols, bool InnerPanel, typename SliceSpec>
506constexpr auto slice(const Eigen::Block<XprType, 1, BlockCols, InnerPanel>& x,
507 SliceSpec&& range) noexcept
508{
509 assert(range.second <= x.size());
510
511 return Eigen::Block<const XprType, 1, Eigen::Dynamic, InnerPanel>(
512 x.nestedExpression(), x.startRow(), x.startCol() + range.first, 1,
513 range.second - range.first);
514}
515
516template <class XprType,
517 int BlockRows,
518 int BlockCols,
519 bool InnerPanel,
520 typename SliceSpec>
521constexpr auto rows(
522 const Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
523 SliceSpec&& rows) noexcept
524{
525 assert(rows.second <= A.rows());
526
527 return Eigen::Block<const XprType, Eigen::Dynamic, BlockCols, InnerPanel>(
528 A.nestedExpression(), A.startRow() + rows.first, A.startCol(),
529 rows.second - rows.first, A.cols());
530}
531
532template <class XprType, int BlockRows, int BlockCols, bool InnerPanel>
533constexpr auto row(
534 const Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
535 Eigen::Index rowIdx) noexcept
536{
537 assert(rowIdx < A.rows());
538
539 return Eigen::Block<const XprType, 1, BlockCols, InnerPanel>(
540 A.nestedExpression(), A.startRow() + rowIdx, A.startCol(), 1, A.cols());
541}
542
543template <class XprType,
544 int BlockRows,
545 int BlockCols,
546 bool InnerPanel,
547 typename SliceSpec>
548constexpr auto cols(
549 const Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
550 SliceSpec&& cols) noexcept
551{
552 assert(cols.second <= A.cols());
553
554 return Eigen::Block<const XprType, BlockRows, Eigen::Dynamic, InnerPanel>(
555 A.nestedExpression(), A.startRow(), A.startCol() + cols.first, A.rows(),
556 cols.second - cols.first);
557}
558
559template <class XprType, int BlockRows, int BlockCols, bool InnerPanel>
560constexpr auto col(
561 const Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel>& A,
562 Eigen::Index colIdx) noexcept
563{
564 assert(colIdx < A.cols());
565
566 return Eigen::Block<const XprType, BlockRows, 1, InnerPanel>(
567 A.nestedExpression(), A.startRow(), A.startCol() + colIdx, A.rows(), 1);
568}
569
570#undef isSlice
571
573template <class T,
574 typename std::enable_if<traits::internal::is_eigen_matrix<T>,
575 int>::type = 0>
576constexpr auto diag(T& A, int diagIdx = 0) noexcept
577{
578 return A.diagonal(diagIdx);
579}
580
581// Transpose view
582template <class matrix_t,
583 typename std::enable_if<(traits::is_eigen_type<matrix_t> &&
584 !matrix_t::IsVectorAtCompileTime),
585 int>::type = 0>
586constexpr auto transpose_view(matrix_t& A) noexcept
587{
588 using T = typename matrix_t::Scalar;
589 using Stride = Eigen::OuterStride<>;
590 assert(A.innerStride() == 1);
591
592 constexpr int Rows_ = matrix_t::ColsAtCompileTime;
593 constexpr int Cols_ = matrix_t::RowsAtCompileTime;
594
595 using transpose_t = Eigen::Matrix<
596 T, Rows_, Cols_,
597 (matrix_t::IsRowMajor) ? Eigen::ColMajor : Eigen::RowMajor,
598 matrix_t::MaxColsAtCompileTime, matrix_t::MaxRowsAtCompileTime>;
599
600 using map_t = Eigen::Map<transpose_t, Eigen::Unaligned, Stride>;
601
602 return map_t((T*)A.data(), A.cols(), A.rows(), A.outerStride());
603}
604
605// Reshape view into matrices
606template <class matrix_t,
607 typename std::enable_if<(traits::is_eigen_type<matrix_t> &&
608 !matrix_t::IsVectorAtCompileTime),
609 int>::type = 0>
610auto reshape(matrix_t& A, Eigen::Index m, Eigen::Index n)
611{
612 using idx_t = Eigen::Index;
613 using T = typename matrix_t::Scalar;
614 using newm_t =
615 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic,
616 matrix_t::IsRowMajor ? Eigen::RowMajor : Eigen::ColMajor>;
617 using map_t = Eigen::Map<newm_t, Eigen::Unaligned, Eigen::OuterStride<>>;
618
619 // constants
620 const idx_t size = A.size();
621 const idx_t new_size = m * n;
622 const idx_t stride = A.outerStride();
623 const bool is_contiguous =
624 (size <= 1) ||
625 (matrix_t::IsRowMajor ? (stride == A.cols() || A.rows() <= 1)
626 : (stride == A.rows() || A.cols() <= 1));
627
628 // Check arguments
629 if (new_size > size)
630 throw std::domain_error("New size is larger than current size");
631 if (A.innerStride() != 1)
632 throw std::domain_error(
633 "Reshaping is not available for matrices with inner stride "
634 "different than 1.");
635
636 if (is_contiguous) {
637 const idx_t s = size - new_size;
638 if constexpr (matrix_t::IsRowMajor)
639 return std::make_pair(map_t((T*)A.data(), m, n, n),
640 map_t((T*)A.data() + new_size, 1, s, s));
641 else
642 return std::make_pair(map_t((T*)A.data(), m, n, m),
643 map_t((T*)A.data() + new_size, s, 1, s));
644 }
645 else {
646 if (m == A.rows() || n == 0) {
647 return std::make_pair(
648 map_t((T*)A.data(), m, n, stride),
649 map_t((T*)A.data() + (matrix_t::IsRowMajor ? n : n * stride),
650 A.rows(), A.cols() - n, stride));
651 }
652 else if (n == A.cols() || m == 0) {
653 return std::make_pair(
654 map_t((T*)A.data(), m, n, stride),
655 map_t((T*)A.data() + (matrix_t::IsRowMajor ? m * stride : m),
656 A.rows() - m, A.cols(), stride));
657 }
658 else {
659 throw std::domain_error(
660 "Cannot reshape to non-contiguous matrix if both the number of "
661 "rows and columns are different.");
662 }
663 }
664}
665template <class vector_t,
666 typename std::enable_if<(traits::is_eigen_type<vector_t> &&
667 vector_t::IsVectorAtCompileTime),
668 int>::type = 0>
669auto reshape(vector_t& v, Eigen::Index m, Eigen::Index n)
670{
671 using idx_t = Eigen::Index;
672 constexpr idx_t Rows_ = vector_t::IsRowMajor ? 1 : Eigen::Dynamic;
673 constexpr idx_t Cols_ = vector_t::IsRowMajor ? Eigen::Dynamic : 1;
674
675 using T = typename vector_t::Scalar;
676 using newm_t = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic,
677 (vector_t::IsRowMajor) ? Eigen::RowMajor
678 : Eigen::ColMajor>;
679 using newv_t = Eigen::Matrix<T, Rows_, Cols_,
680 (vector_t::IsRowMajor) ? Eigen::RowMajor
681 : Eigen::ColMajor>;
682 using map0_t = Eigen::Map<newm_t, Eigen::Unaligned, Eigen::OuterStride<>>;
683 using map1_t = Eigen::Map<newv_t, Eigen::Unaligned, Eigen::InnerStride<>>;
684
685 // constants
686 const idx_t size = v.size();
687 const idx_t new_size = m * n;
688 const idx_t s = size - new_size;
689 const idx_t stride = v.innerStride();
690 const bool is_contiguous = (size <= 1 || stride == 1);
691
692 // Check arguments
693 if (new_size > size)
694 throw std::domain_error("New size is larger than current size");
695 if (!is_contiguous && (newm_t::IsRowMajor ? (n > 1) : (m > 1)))
696 throw std::domain_error(
697 "New sizes are not compatible with the current vector.");
698
699 if (is_contiguous) {
700 return std::make_pair(
701 map0_t((T*)v.data(), m, n, (newm_t::IsRowMajor) ? n : m),
702 map1_t((T*)v.data() + new_size, s, Eigen::InnerStride<>(1)));
703 }
704 else {
705 return std::make_pair(map0_t((T*)v.data(), m, n, stride),
706 map1_t((T*)v.data() + new_size * stride, s,
707 Eigen::InnerStride<>(stride)));
708 }
709}
710
711// Reshape view into vectors
712template <class matrix_t,
713 typename std::enable_if<(traits::is_eigen_type<matrix_t> &&
714 !matrix_t::IsVectorAtCompileTime),
715 int>::type = 0>
716auto reshape(matrix_t& A, Eigen::Index n)
717{
718 using idx_t = Eigen::Index;
719 constexpr idx_t Rows_ = matrix_t::IsRowMajor ? 1 : Eigen::Dynamic;
720 constexpr idx_t Cols_ = matrix_t::IsRowMajor ? Eigen::Dynamic : 1;
721
722 using T = typename matrix_t::Scalar;
723 using newm_t =
724 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic,
725 matrix_t::IsRowMajor ? Eigen::RowMajor : Eigen::ColMajor>;
726 using newv_t =
727 Eigen::Matrix<T, Rows_, Cols_,
728 matrix_t::IsRowMajor ? Eigen::RowMajor : Eigen::ColMajor>;
729 using map0_t = Eigen::Map<newv_t, Eigen::Unaligned, Eigen::InnerStride<>>;
730 using map1_t = Eigen::Map<newm_t, Eigen::Unaligned, Eigen::OuterStride<>>;
731
732 // constants
733 const idx_t size = A.size();
734 const idx_t stride = A.outerStride();
735 const bool is_contiguous =
736 (size <= 1) ||
737 (matrix_t::IsRowMajor ? (stride == A.cols() || A.rows() <= 1)
738 : (stride == A.rows() || A.cols() <= 1));
739
740 // Check arguments
741 if (n > size)
742 throw std::domain_error("New size is larger than current size");
743 if (A.innerStride() != 1)
744 throw std::domain_error(
745 "Reshaping is not available for matrices with inner stride "
746 "different than 1.");
747
748 if (is_contiguous) {
749 const idx_t s = size - n;
750 return std::make_pair(map0_t((T*)A.data(), n, Eigen::InnerStride<>(1)),
751 (matrix_t::IsRowMajor)
752 ? map1_t((T*)A.data() + n, 1, s, s)
753 : map1_t((T*)A.data() + n, s, 1, s));
754 }
755 else {
756 if (n == 0) {
757 return std::make_pair(
758 map0_t((T*)A.data(), 0, Eigen::InnerStride<>(1)),
759 map1_t((T*)A.data(), A.rows(), A.cols(), stride));
760 }
761 else if (n == A.rows()) {
762 if constexpr (matrix_t::IsRowMajor)
763 return std::make_pair(
764 map0_t((T*)A.data(), n, Eigen::InnerStride<>(stride)),
765 map1_t((T*)A.data() + 1, A.rows(), A.cols() - 1, stride));
766 else
767 return std::make_pair(
768 map0_t((T*)A.data(), n, Eigen::InnerStride<>(1)),
769 map1_t((T*)A.data() + stride, A.rows(), A.cols() - 1,
770 stride));
771 }
772 else if (n == A.cols()) {
773 if constexpr (matrix_t::IsRowMajor)
774 return std::make_pair(
775 map0_t((T*)A.data(), n, Eigen::InnerStride<>(1)),
776 map1_t((T*)A.data() + stride, A.rows() - 1, A.cols(),
777 stride));
778 else
779 return std::make_pair(
780 map0_t((T*)A.data(), n, Eigen::InnerStride<>(stride)),
781 map1_t((T*)A.data() + 1, A.rows() - 1, A.cols(), stride));
782 }
783 else {
784 throw std::domain_error(
785 "Cannot reshape to non-contiguous matrix into a vector if the "
786 "new size is different from the number of rows and columns.");
787 }
788 }
789}
790template <class vector_t,
791 typename std::enable_if<(traits::is_eigen_type<vector_t> &&
792 vector_t::IsVectorAtCompileTime),
793 int>::type = 0>
794auto reshape(vector_t& v, Eigen::Index n)
795{
796 using idx_t = Eigen::Index;
797 constexpr idx_t Rows_ = vector_t::IsRowMajor ? 1 : Eigen::Dynamic;
798 constexpr idx_t Cols_ = vector_t::IsRowMajor ? Eigen::Dynamic : 1;
799
800 using T = typename vector_t::Scalar;
801 using newv_t =
802 Eigen::Matrix<T, Rows_, Cols_,
803 vector_t::IsRowMajor ? Eigen::RowMajor : Eigen::ColMajor>;
804 using map_t = Eigen::Map<newv_t, Eigen::Unaligned, Eigen::InnerStride<>>;
805
806 // constants
807 const idx_t stride = v.innerStride();
808
809 // Check arguments
810 if (n > v.size())
811 throw std::domain_error("New size is larger than current size");
812
813 return std::make_pair(map_t((T*)v.data(), n, Eigen::InnerStride<>(stride)),
814 map_t((T*)v.data() + n * stride, v.size() - n,
815 Eigen::InnerStride<>(stride)));
816}
817
818// -----------------------------------------------------------------------------
819// Deduce matrix and vector type from two provided ones
820
821namespace traits {
822
823 template <typename T>
824 constexpr bool cast_to_eigen_type = is_eigen_type<T> || is_stdvector_type<T>
825#ifdef TLAPACK_LEGACYARRAY_HH
826 || is_legacy_type<T>
827#endif
828#ifdef TLAPACK_MDSPAN_HH
829 || is_mdspan_type<T>
830#endif
831 ;
832
833 // for two types
834 // should be especialized for every new matrix class
835 template <typename matrixA_t, typename matrixB_t>
836 struct matrix_type_traits<
837 matrixA_t,
838 matrixB_t,
839 typename std::enable_if<
840 ((is_eigen_type<matrixA_t> ||
841 is_eigen_type<matrixB_t>)&&cast_to_eigen_type<matrixA_t> &&
842 cast_to_eigen_type<matrixB_t>),
843 int>::type> {
844 using T = scalar_type<type_t<matrixA_t>, type_t<matrixB_t>>;
845
846 static constexpr Layout LA = layout<matrixA_t>;
847 static constexpr Layout LB = layout<matrixB_t>;
848 static constexpr int L =
849 ((LA == Layout::RowMajor) && (LB == Layout::RowMajor))
850 ? Eigen::RowMajor
851 : Eigen::ColMajor;
852
853 using type = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, L>;
854 };
855
856 // for two types
857 // should be especialized for every new vector class
858 template <typename vecA_t, typename vecB_t>
859 struct vector_type_traits<
860 vecA_t,
861 vecB_t,
862 typename std::enable_if<
863 ((is_eigen_type<vecA_t> ||
864 is_eigen_type<vecB_t>)&&cast_to_eigen_type<vecA_t> &&
865 cast_to_eigen_type<vecB_t>),
866 int>::type> {
867 using T = scalar_type<type_t<vecA_t>, type_t<vecB_t>>;
868 using type = Eigen::Matrix<T, Eigen::Dynamic, 1>;
869 };
870
871#if !defined(TLAPACK_MDSPAN_HH) && !defined(TLAPACK_LEGACYARRAY_HH)
872 template <class vecA_t, class vecB_t>
873 struct matrix_type_traits<
874 vecA_t,
875 vecB_t,
876 std::enable_if_t<traits::is_stdvector_type<vecA_t> &&
877 traits::is_stdvector_type<vecB_t>,
878 int>> {
879 using T = scalar_type<type_t<vecA_t>, type_t<vecB_t>>;
880 using type = Eigen::Matrix<T, Eigen::Dynamic, 1, Eigen::ColMajor>;
881 };
882
883 template <class vecA_t, class vecB_t>
884 struct vector_type_traits<
885 vecA_t,
886 vecB_t,
887 std::enable_if_t<traits::is_stdvector_type<vecA_t> &&
888 traits::is_stdvector_type<vecB_t>,
889 int>> {
890 using T = scalar_type<type_t<vecA_t>, type_t<vecB_t>>;
891 using type = Eigen::Matrix<T, Eigen::Dynamic, 1, Eigen::ColMajor>;
892 };
893#endif
894
895} // namespace traits
896
897// -----------------------------------------------------------------------------
898// Cast to Legacy arrays
899
900template <class T, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
901constexpr auto legacy_matrix(
902 const Eigen::Matrix<T, Rows, Cols, Options, MaxRows, MaxCols>& A) noexcept
903{
904 using matrix_t = Eigen::Matrix<T, Rows, Cols, Options, MaxRows, MaxCols>;
905 using idx_t = Eigen::Index;
906
907 if constexpr (matrix_t::IsVectorAtCompileTime)
908 return legacy::Matrix<T, idx_t>{Layout::ColMajor, 1, A.size(),
909 (T*)A.data(), A.innerStride()};
910 else {
911 constexpr Layout L = layout<matrix_t>;
912 return legacy::Matrix<T, idx_t>{L, A.rows(), A.cols(), (T*)A.data(),
913 A.outerStride()};
914 }
915}
916
917template <class Derived>
918constexpr auto legacy_matrix(
919 const Eigen::MapBase<Derived, Eigen::ReadOnlyAccessors>& A) noexcept
920{
921 using T = typename Derived::Scalar;
922 using idx_t = Eigen::Index;
923
924 if constexpr (Derived::IsVectorAtCompileTime) {
925 assert(A.outerStride() == 1 ||
926 (A.innerStride() == 1 &&
927 traits::internal::is_eigen_block<Derived>));
928 return legacy::Matrix<T, idx_t>{Layout::ColMajor, 1, A.size(),
929 (T*)A.data(), A.innerStride()};
930 }
931 else {
932 assert(A.innerStride() == 1);
933 constexpr Layout L = layout<Derived>;
934 return legacy::Matrix<T, idx_t>{L, A.rows(), A.cols(), (T*)A.data(),
935 A.outerStride()};
936 }
937}
938
939template <class T, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
940constexpr auto legacy_vector(
941 const Eigen::Matrix<T, Rows, Cols, Options, MaxRows, MaxCols>& A) noexcept
942{
943 using matrix_t = Eigen::Matrix<T, Rows, Cols, Options, MaxRows, MaxCols>;
944 using idx_t = Eigen::Index;
945
946 if constexpr (matrix_t::IsVectorAtCompileTime)
947 return legacy::Vector<T, idx_t>{A.size(), (T*)A.data(),
948 A.innerStride()};
949 else {
950 assert(A.rows() == 1 || A.cols() == 1);
951 return legacy::Vector<T, idx_t>{A.size(), (T*)A.data(),
952 A.outerStride()};
953 }
954}
955
956template <class Derived>
957constexpr auto legacy_vector(
958 const Eigen::MapBase<Derived, Eigen::ReadOnlyAccessors>& A) noexcept
959{
960 using T = typename Derived::Scalar;
961 using idx_t = Eigen::Index;
962
963 if constexpr (Derived::IsVectorAtCompileTime) {
964 assert(A.outerStride() == 1 ||
965 (A.innerStride() == 1 &&
966 traits::internal::is_eigen_block<Derived>));
967 return legacy::Vector<T, idx_t>{A.size(), (T*)A.data(),
968 A.innerStride()};
969 }
970 else {
971 assert(A.innerStride() == 1);
972 assert(A.rows() == 1 || A.cols() == 1);
973 return legacy::Vector<T, idx_t>{A.size(), (T*)A.data(),
974 A.outerStride()};
975 }
976}
977
978} // namespace tlapack
979
980#endif // TLAPACK_EIGEN_HH
tlapack::Layout
Layout
Definition types.hpp:24
tlapack::Layout::ColMajor
@ ColMajor
Column-major layout.
tlapack::traits::is_eigen_type
constexpr bool is_eigen_type
True if T is derived from Eigen::DenseBase.
Definition eigen.hpp:69
tlapack::traits::internal::is_eigen_matrix
constexpr bool is_eigen_matrix
True if T is derived from Eigen::EigenMatrix<T>
Definition eigen.hpp:50
tlapack::diag
constexpr auto diag(T &A, int diagIdx=0) noexcept
Get the Diagonal of an Eigen Matrix.
Definition eigen.hpp:576
tlapack::traits::internal::is_eigen_block
constexpr bool is_eigen_block
True if T is derived from Eigen::EigenBlock.
Definition eigen.hpp:56
tlapack::traits::internal::is_eigen_dense
constexpr bool is_eigen_dense
True if T is derived from Eigen::EigenDense<T>
Definition eigen.hpp:44
tlapack::real_type
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
tlapack::legacy::Matrix
Describes a row- or column-major matrix.
Definition types.hpp:443
tlapack::traits::CreateFunctor
Functor for data creation.
Definition arrayTraits.hpp:89
tlapack::traits::CreateFunctor::operator()
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
tlapack::traits::CreateStaticFunctor
Functor for data creation with static size.
Definition arrayTraits.hpp:141
tlapack::traits::CreateStaticFunctor::operator()
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
tlapack::traits::complex_type_traits
Complex type traits for the list of types Types.
Definition scalar_type_traits.hpp:145
tlapack::traits::layout_trait
Trait to determine the layout of a given data structure.
Definition arrayTraits.hpp:75
tlapack::traits::matrix_type_traits
Matrix type deduction.
Definition arrayTraits.hpp:176
tlapack::traits::real_type_traits
Real type traits for the list of types Types.
Definition scalar_type_traits.hpp:71
tlapack::traits::vector_type_traits
Vector type deduction.
Definition arrayTraits.hpp:203