18
19
20
24#include "./det_and_inverse.hpp"
25#include "../basic_array.hpp"
26#include "../declarations.hpp"
27#include "../exceptions.hpp"
28#include "../lapack/interface/cxx_interface.hpp"
29#include "../layout/policies.hpp"
30#include "../macros.hpp"
31#include "../traits.hpp"
36namespace nda::linalg {
42 auto _eigen_element_impl(M &&m,
char compz) {
43 using value_type =
typename std::decay_t<M>::value_type;
47 EXPECTS(is_matrix_square(m,
true));
48 EXPECTS(m.indexmap().is_contiguous());
51 int dim = m.extent(0);
53 array<
double, 1> ev(dim);
54 array<value_type, 1> work(lwork);
55 array<
double, 1> work2(is_complex_v<value_type> ? lwork : 0);
57#if defined(__has_feature
)
58#if __has_feature
(memory_sanitizer)
67 if constexpr (
not is_complex_v<value_type>) {
68 lapack::f77::syev(compz,
'U', dim, m.data(), dim, ev
.data(), work.data(), lwork, info);
70 lapack::f77::heev(compz,
'U', dim, m.data(), dim, ev
.data(), work.data(), lwork, work2
.data(), info);
72 if (info)
NDA_RUNTIME_ERROR <<
"Error in nda::linalg::detail::_eigen_element_impl: Diagonalization error";
79
80
81
84
85
86
87
88
89
90
91
92
93
94
95
98 auto m_copy = matrix<
typename M::value_type,
F_layout>(m);
99 auto ev = detail::_eigen_element_impl(m_copy,
'V');
100 return std::pair<array<
double, 1>,
typename M::regular_type>{ev, m_copy};
104
105
106
107
108
109
110
111
112
113
114
115 template <
typename M>
117 auto m_copy = matrix<
typename M::value_type,
F_layout>(m);
118 return detail::_eigen_element_impl(m_copy,
'N');
122
123
124
125
126
127
128
129
130
131
132
133 template <
typename M>
135 return detail::_eigen_element_impl(m,
'N');
ValueType * data() noexcept
Get a pointer to the actual data (in general this is not the beginning of thr memory block for a view...
#define NDA_RUNTIME_ERROR
Contiguous layout policy with Fortran-order (column-major order).