19template <const
int n, const
int m,
typename T,
typename Mtype>
22template <const
int n, const
int m,
typename T =
double>
25template <
int n,
int m,
typename T>
28template <
int n,
typename T>
34template <
int n,
typename T>
40template <
int n,
typename T>
46template <
int n,
typename T>
62 static_assert(M::is_matrix() && M::rows() == M::columns(),
"SVD only for square matrix");
64 DiagonalMatrix<M::size(), hila::arithmetic_type<M>> singularvalues;
74 static_assert(M::is_matrix() && M::rows() == M::columns(),
75 "Eigenvalues only for square matrix");
101template <const
int n, const
int m,
typename T,
typename Mtype>
110 "Matrix requires Complex or arithmetic type");
113 using base_type = hila::arithmetic_type<T>;
114 using argument_type = T;
118 static constexpr bool is_matrix() {
129 return (n == 1 || m == 1);
148 template <
typename S,
int nn = n,
int mm = m,
149 std::enable_if_t<(hila::is_assignable<T &, S>::value && nn == mm),
int> = 0>
150 explicit inline Matrix_t(
const S rhs) {
151 for (
int i = 0; i < n; i++)
152 for (
int j = 0; j < n; j++) {
170 inline Matrix_t(
const std::nullptr_t &z) {
171 for (
int i = 0; i < n * m; i++) {
178 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
179 inline Matrix_t(std::initializer_list<S> rhs) {
180 assert(rhs.size() == n * m &&
181 "Matrix/Vector initializer list size must match variable size");
183 for (
auto it = rhs.begin(); it != rhs.end(); it++, i++) {
190 template <
typename Tm =
Mtype,
191 std::enable_if_t<!std::is_same<Tm, Matrix<n, m, T>>::value,
int> = 0>
192 inline operator Mtype &() {
193 return *
reinterpret_cast<Mtype *
>(
this);
197 template <
typename Tm =
Mtype,
198 std::enable_if_t<!std::is_same<Tm, Matrix<n, m, T>>::value,
int> = 0>
199 inline operator const Mtype &()
const {
200 return *
reinterpret_cast<const Mtype *
>(
this);
241 template <
int q = n,
int p = m, std::enable_if_t<q == 1,
int> = 0>
246 template <
int q = n,
int p = m, std::enable_if_t<p == 1,
int> = 0>
247 static constexpr int size() {
251 template <
int q = n,
int p = m, std::enable_if_t<q == p,
int> = 0>
252 static constexpr int size() {
272 inline T
e(
const int i,
const int j)
const {
278 inline T &
e(
const int i,
const int j) const_function {
297 template <
int q = n,
int p = m, std::enable_if_t<(q == 1 || p == 1),
int> = 0>
298 inline T
e(
const int i)
const {
303 template <
int q = n,
int p = m, std::enable_if_t<(q == 1 || p == 1),
int> = 0>
304 inline T &
e(
const int i) const_function {
323 template <
int q = n,
int p = m, std::enable_if_t<(q == 1 || p == 1),
int> = 0>
328 template <
int q = n,
int p = m, std::enable_if_t<(q == 1 || p == 1),
int> = 0>
329 inline T &
operator[](
const int i) const_function {
373 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
375 for (
int i = 0; i < m; i++)
391 for (
int i = 0; i < n; i++)
414 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
416 for (
int i = 0; i < n; i++)
431 static_assert(n == m,
"diagonal() method defined only for square matrices");
433 for (
int i = 0; i < n; i++)
434 res.
e(i) = (*this).e(i, i);
453 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
455 static_assert(n == m,
"set_diagonal() method defined only for square matrices");
456 for (
int i = 0; i < n; i++)
457 (*this).e(i, i) = v.
e(i);
480 template <
int mm = m, std::enable_if_t<mm == 1,
int> = 0>
485 template <
int mm = m, std::enable_if_t<mm == 1,
int> = 0>
523 for (
int i = 0; i < n * m; i++) {
560 template <
typename S>
562 for (
int i = 0; i < n; i++)
563 for (
int j = 0; j < m; j++) {
564 if (
e(i, j) != rhs.
e(i, j))
579 template <
typename S>
581 return !(*
this == rhs);
599 template <
typename S,
typename MT,
600 std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
602 for (
int i = 0; i < n * m; i++) {
622 for (
int i = 0; i < n * m; i++) {
644 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value && n == m,
int> = 0>
647 for (
int i = 0; i < n; i++)
648 for (
int j = 0; j < n; j++) {
667 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
669 static_assert(n == m,
670 "Assigning DiagonalMatrix to Matrix possible only for square matrices");
672 for (
int i = 0; i < n; i++)
673 for (
int j = 0; j < n; j++) {
694 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
696 assert(rhs.size() == n * m &&
"Initializer list has a wrong size in assignment");
698 for (
auto it = rhs.begin(); it != rhs.end(); it++, i++) {
721 template <
typename S,
typename MT,
722 std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
724 for (
int i = 0; i < n * m; i++) {
745 template <
typename S,
typename MT,
746 std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
748 for (
int i = 0; i < n * m; i++) {
770 template <
typename S,
771 std::enable_if_t<hila::is_assignable<T &, hila::type_plus<T, S>>::value,
int> = 0>
774 static_assert(n == m,
"rows != columns : scalar addition possible for square matrix only!");
776 for (
int i = 0; i < n; i++) {
796 template <
typename S,
797 std::enable_if_t<hila::is_assignable<T &, hila::type_minus<T, S>>::value,
int> = 0>
799 static_assert(n == m,
800 "rows != columns : scalar subtraction possible for square matrix only!");
801 for (
int i = 0; i < n; i++) {
827 template <
int p,
typename S,
typename MT,
828 std::enable_if_t<hila::is_assignable<T &, hila::type_mul<T, S>>::value,
int> = 0>
830 static_assert(m == p,
"can't assign result of *= to lhs Matrix, because doing so "
831 "would change it's dimensions");
853 template <
typename S,
854 std::enable_if_t<hila::is_assignable<T &, hila::type_mul<T, S>>::value,
int> = 0>
856 for (
int i = 0; i < n * m; i++) {
879 template <
typename S,
880 std::enable_if_t<hila::is_assignable<T &, hila::type_div<T, S>>::value,
int> = 0>
882 for (
int i = 0; i < n * m; i++) {
895 template <
typename S,
896 std::enable_if_t<hila::is_assignable<T &, hila::type_plus<T, S>>::value,
int> = 0>
898 static_assert(n == m,
"Assigning DiagonalMatrix possible only for square matrix");
900 for (
int i = 0; i < n; i++)
912 template <
typename S,
913 std::enable_if_t<hila::is_assignable<T &, hila::type_plus<T, S>>::value,
int> = 0>
915 static_assert(n == m,
"Assigning DiagonalMatrix possible only for square matrix");
917 for (
int i = 0; i < n; i++)
931 template <
typename S,
932 std::enable_if_t<hila::is_assignable<T &, hila::type_mul<T, S>>::value,
int> = 0>
935 for (
int i = 0; i < n; i++)
936 for (
int j = 0; j < m; j++)
952 template <
typename S,
953 std::enable_if_t<hila::is_assignable<T &, hila::type_div<T, S>>::value,
int> = 0>
956 for (
int i = 0; i < n; i++)
957 for (
int j = 0; j < m; j++)
979 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
981 for (
int i = 0; i < n * m; i++)
994 template <
int mm = m,
995 typename Rtype =
typename std::conditional<n == m, Mtype, Matrix<m, n, T>>::type,
996 std::enable_if_t<(mm != 1),
int> = 0>
999 for (
int i = 0; i < n; i++)
1000 for (
int j = 0; j < m; j++) {
1001 res.e(j, i) =
e(i, j);
1013 template <
int mm = m, std::enable_if_t<mm == 1,
int> = 0>
1025 template <
int nn = n, std::enable_if_t<nn == 1 && m != 1,
int> = 0>
1038 for (
int i = 0; i < n * m; i++) {
1051 template <
typename Rtype =
typename std::conditional<n == m, Mtype, Matrix<m, n, T>>::type>
1054 for (
int i = 0; i < n; i++)
1055 for (
int j = 0; j < m; j++) {
1056 res.e(j, i) =
::conj(
e(i, j));
1068 template <
typename Rtype =
typename std::conditional<n == m, Mtype, Matrix<m, n, T>>::type>
1083 for (
int i = 0; i < n * m; i++) {
1103 for (
int i = 0; i < m * n; i++) {
1116 for (
int i = 0; i < m * n; i++) {
1130 static_assert(n == m,
"trace not defined for non square matrices");
1132 for (
int i = 0; i < n; i++) {
1149 template <
int p,
int q,
typename S,
typename MT>
1152 static_assert(p == m && q == n,
"mul_trace(): argument matrix size mismatch");
1154 hila::type_mul<T, S> res = 0;
1155 for (
int i = 0; i < n; i++)
1156 for (
int j = 0; j < m; j++) {
1157 res +=
e(i, j) * rm.
e(j, i);
1168 hila::arithmetic_type<T> result(0);
1169 for (
int i = 0; i < n * m; i++) {
1181 template <
typename S = T,
1182 std::enable_if_t<hila::is_floating_point<hila::arithmetic_type<S>>::value,
int> = 0>
1183 hila::arithmetic_type<T>
norm()
const {
1187 template <
typename S = T,
1188 std::enable_if_t<!hila::is_floating_point<hila::arithmetic_type<S>>::value,
int> = 0>
1189 double norm()
const {
1190 return sqrt(
static_cast<double>(
squarenorm()));
1196 template <typename S = T, std::enable_if_t<hila::is_arithmetic<S>::value,
int> = 0>
1199 for (
int i = 1; i < n * m; i++) {
1209 template <typename S = T, std::enable_if_t<hila::is_arithmetic<S>::value,
int> = 0>
1212 for (
int i = 1; i < n * m; i++) {
1250 template <
int p,
int q,
typename S,
typename R = hila::type_mul<T, S>>
1252 static_assert(m == 1 && q == 1 && p == n,
1253 "dot() product only for vectors of the same length");
1256 for (
int i = 0; i < n; i++) {
1293 template <
int p,
int q,
typename S,
typename R = hila::type_mul<T, S>>
1295 static_assert(m == 1 && q == 1,
"outer_product() only for vectors");
1298 for (
int i = 0; i < n; i++) {
1299 for (
int j = 0; j < p; j++) {
1332 static_assert(hila::is_floating_point<hila::arithmetic_type<T>>::value,
1333 "Matrix/Vector random() requires non-integral type elements");
1335 for (
int i = 0; i < n * m; i++) {
1354 static_assert(hila::is_floating_point<hila::arithmetic_type<T>>::value,
1355 "Matrix/Vector gaussian_random() requires non-integral type elements");
1358 if constexpr (hila::is_complex<T>::value) {
1359 for (
int i = 0; i < n * m; i++) {
1360 c[i].gaussian_random(width);
1367 for (
int i = 0; i < n * m - 1; i += 2) {
1369 c[i + 1] = gr * width;
1371 if constexpr ((n * m) % 2 > 0) {
1387 for (
int i = 0; i < m; i++)
1388 res.
set_column(i, this->column(permutation[i]));
1401 for (
int i = 0; i < n; i++)
1402 res.
set_row(i, this->row(permutation[i]));
1418 "permute() only for vectors, use permute_rows() or permute_columns() for matrices");
1419 static_assert(N ==
Mtype::size(),
"Incorrect size of permutation vector");
1422 for (
int i = 0; i < N; i++) {
1423 res[i] = (*this)[permutation[i]];
1454#pragma hila novector
1458 static_assert(n == 1 || m == 1,
"Sorting possible only for vectors");
1460 "Sorting possible only for arithmetic vector elements");
1461 static_assert(N ==
Mtype::size(),
"Incorrect size of permutation vector");
1463 for (
int i = 0; i < N; i++)
1465 if (hila::sort::unsorted == order) {
1469 if (hila::sort::ascending == order) {
1470 for (
int i = 1; i < N; i++) {
1471 for (
int j = i; j > 0 &&
c[permutation[j - 1]] >
c[permutation[j]]; j--)
1472 hila::swap(permutation[j], permutation[j - 1]);
1475 for (
int i = 1; i < N; i++) {
1476 for (
int j = i; j > 0 &&
c[permutation[j - 1]] <
c[permutation[j]]; j--)
1477 hila::swap(permutation[j], permutation[j - 1]);
1481 return this->
permute(permutation);
1503#pragma hila novector
1504 Mtype sort(hila::sort order = hila::sort::ascending)
const {
1505 static_assert(n == 1 || m == 1,
"Sorting possible only for vectors");
1508 return sort(permutation, order);
1531 template <
typename R,
typename Mt>
1533 static_assert(hila::contains_complex<T>::value || !hila::contains_complex<R>::value,
1534 "Cannot multiply real matrix with complex 2x2 block matrix");
1537 for (
int i = 0; i < m; i++) {
1538 a.
e(0) = this->
e(p, i);
1539 a.
e(1) = this->
e(q, i);
1543 this->
e(p, i) = a.e(0);
1544 this->
e(q, i) = a.e(1);
1560 template <
typename R,
typename Mt>
1562 static_assert(hila::contains_complex<T>::value || !hila::contains_complex<R>::value,
1563 "Cannot multiply real matrix with complex 2x2 block matrix");
1566 for (
int i = 0; i < n; i++) {
1567 a.
e(0) = this->
e(i, p);
1568 a.
e(1) = this->
e(i, q);
1572 this->
e(i, p) = a.
e(0);
1573 this->
e(i, q) = a.
e(1);
1586#pragma hila novector
1588#pragma hila novector
1590#pragma hila novector
1603#pragma hila novector
1604 template <
typename Et,
typename Mt,
typename MT>
1607 enum hila::sort sorted = hila::sort::unsorted)
const;
1612 template <
typename Et,
typename Mt,
typename MT>
1615 enum hila::sort sorted = hila::sort::unsorted)
const {
1619 eigenvaluevec = d.asArray().asVector();
1624#pragma hila novector
1625 template <
typename Et,
typename Mt,
typename MT>
1628 enum hila::sort sorted = hila::sort::unsorted)
const;
1633#pragma hila novector
1634 template <
typename Et,
typename Mt,
typename MT>
1637 enum hila::sort sorted = hila::sort::unsorted)
const;
1646 std::string
str(
int prec = 8,
char separator =
' ')
const {
1647 std::stringstream text;
1648 text.precision(prec);
1649 for (
int i = 0; i < n; i++) {
1650 for (
int j = 0; j < m; j++) {
1652 if (i < n - 1 || j < m - 1)
1678template <
int n,
int m,
typename T>
1684 using argument_type = T;
1723 template <
typename S,
int nn = n,
int mm = m,
1724 std::enable_if_t<(hila::is_assignable<T &, S>::value && nn == mm),
int> = 0>
1726 for (
int i = 0; i < n; i++)
1727 for (
int j = 0; j < n; j++) {
1729 this->
e(i, j) = rhs;
1753 template <
typename S,
typename MT,
1754 std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
1756 for (
int i = 0; i < n * m; i++) {
1757 this->c[i] = rhs.c[i];
1768 for (
int i = 0; i < n * m; i++) {
1787 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
1789 assert(rhs.size() == n * m &&
1790 "Matrix/Vector initializer list size must match variable size");
1792 for (
auto it = rhs.begin(); it != rhs.end(); it++, i++) {
1818template <
typename T1,
typename T2>
1819struct matrix_op_type_s {
1820 using type =
Matrix<T1::rows(), T1::columns(), hila::ntype_op<T1, T2>>;
1824template <
typename T>
1825struct matrix_op_type_s<T, T> {
1830template <
typename T1,
typename T2>
1831using mat_x_mat_type =
typename matrix_op_type_s<T1, T2>::type;
1841template <
typename Mt,
typename S,
typename Enable =
void>
1842struct matrix_scalar_op_s {
1844 Matrix<Mt::rows(), Mt::columns(),
1848template <
typename Mt,
typename S>
1849struct matrix_scalar_op_s<
1851 typename
std::enable_if_t<std::is_convertible<hila::type_plus<hila::number_type<Mt>, S>,
1852 hila::number_type<Mt>>::value>> {
1854 using type =
typename std::conditional<
1855 hila::is_floating_point<hila::arithmetic_type<Mt>>::value, Mt,
1856 Matrix<Mt::rows(), Mt::columns(),
1860template <
typename Mt,
typename S>
1861using mat_scalar_type =
typename matrix_scalar_op_s<Mt, S>::type;
1883template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
1885 return arg.transpose();
1903template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
1923template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
1925 return arg.adjoint();
1939template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
1941 return arg.adjoint();
1959template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
1980template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
2000template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
2020template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
2049template <
typename Mtype1,
typename Mtype2,
2050 std::enable_if_t<Mtype1::is_matrix() && Mtype2::is_matrix(),
int> = 0,
2051 typename Rtype = hila::mat_x_mat_type<Mtype1, Mtype2>,
2052 std::enable_if_t<!std::is_same<Mtype1, Rtype>::value,
int> = 0>
2055 constexpr int n = Mtype1::rows();
2056 constexpr int m = Mtype1::columns();
2058 static_assert(n == Mtype2::rows() && m == Mtype2::columns(),
"Matrix sizes do not match");
2061 for (
int i = 0; i < n; i++)
2062 for (
int j = 0; j < m; j++)
2063 r.e(i, j) = a.e(i, j) + b.e(i, j);
2079template <
typename Mtype1,
typename Mtype2,
2080 std::enable_if_t<Mtype1::is_matrix() && Mtype2::is_matrix(),
int> = 0,
2081 typename Rtype = hila::mat_x_mat_type<Mtype1, Mtype2>,
2082 std::enable_if_t<std::is_same<Mtype1, Rtype>::value,
int> = 0>
2085 constexpr int n = Mtype1::rows();
2086 constexpr int m = Mtype1::columns();
2088 static_assert(n == Mtype2::rows() && m == Mtype2::columns(),
"Matrix sizes do not match");
2090 for (
int i = 0; i < n; i++)
2091 for (
int j = 0; j < m; j++)
2092 a.e(i, j) += b.e(i, j);
2117template <
typename Mtype1,
typename Mtype2,
2118 std::enable_if_t<Mtype1::is_matrix() && Mtype2::is_matrix(),
int> = 0,
2119 typename Rtype = hila::mat_x_mat_type<Mtype1, Mtype2>>
2122 constexpr int n = Mtype1::rows();
2123 constexpr int m = Mtype1::columns();
2125 static_assert(n == Mtype2::rows() && m == Mtype2::columns(),
"Matrix sizes do not match");
2128 for (
int i = 0; i < n; i++)
2129 for (
int j = 0; j < m; j++)
2130 r.e(i, j) = a.e(i, j) - b.e(i, j);
2155template <
typename Mtype,
typename S,
2157 typename Rtype = hila::mat_scalar_type<Mtype, S>>
2161 "Matrix + scalar possible only for square matrix");
2164 for (
int i = 0; i < Rtype::rows(); i++)
2165 for (
int j = 0; j < Rtype::columns(); j++) {
2166 r.e(i, j) = a.
e(i, j);
2192template <
typename Mtype,
typename S,
2194 typename Rtype = hila::mat_scalar_type<Mtype, S>>
2198 "Matrix + scalar possible only for square matrix");
2224template <
typename Mtype,
typename S,
2226 typename Rtype = hila::mat_scalar_type<Mtype, S>>
2230 "Matrix - scalar possible only for square matrix");
2233 for (
int i = 0; i < Rtype::rows(); i++)
2234 for (
int j = 0; j < Rtype::columns(); j++) {
2235 r.e(i, j) = a.
e(i, j);
2262template <
typename Mtype,
typename S,
2264 typename Rtype = hila::mat_scalar_type<Mtype, S>>
2268 "Scalar - matrix possible only for square matrix");
2271 for (
int i = 0; i < Rtype::rows(); i++)
2272 for (
int j = 0; j < Rtype::columns(); j++) {
2273 r.e(i, j) = -a.
e(i, j);
2292template <
typename Mt, std::enable_if_t<Mt::is_matrix() && Mt::rows() == Mt::columns(),
int> = 0>
2295 constexpr int n = Mt::rows();
2299 for (
int i = 0; i < n; i++)
2300 for (
int j = 0; j < n; j++) {
2301 res.e(i, j) = a.e(i, 0) * b.e(0, j);
2302 for (
int k = 1; k < n; k++) {
2303 res.e(i, j) += a.e(i, k) * b.e(k, j);
2343template <
typename Mt1,
typename Mt2,
2344 std::enable_if_t<Mt1::is_matrix() && Mt2::is_matrix() && !std::is_same<Mt1, Mt2>::value,
2346 typename R = hila::ntype_op<hila::number_type<Mt1>, hila::number_type<Mt2>>,
2347 int n = Mt1::rows(),
int m = Mt2::columns()>
2350 constexpr int p = Mt1::columns();
2351 static_assert(p == Mt2::rows(),
"Matrix size: LHS columns != RHS rows");
2355 if constexpr (n > 1 && m > 1) {
2356 for (
int i = 0; i < n; i++)
2357 for (
int j = 0; j < m; j++) {
2358 res.
e(i, j) = a.e(i, 0) * b.e(0, j);
2359 for (
int k = 1; k < p; k++) {
2360 res.
e(i, j) += a.e(i, k) * b.e(k, j);
2363 }
else if constexpr (m == 1) {
2365 for (
int i = 0; i < n; i++) {
2366 res.
e(i) = a.e(i, 0) * b.e(0);
2367 for (
int k = 1; k < p; k++) {
2368 res.
e(i) += a.e(i, k) * b.e(k);
2371 }
else if constexpr (n == 1) {
2373 for (
int j = 0; j < m; j++) {
2374 res.
e(j) = a.e(0) * b.e(0, j);
2375 for (
int k = 1; k < p; k++) {
2376 res.
e(j) += a.e(k) * b.e(k, j);
2409template <
int m,
int n,
typename T1,
typename T2,
typename Rt = hila::ntype_op<T1, T2>>
2412 static_assert(m == n,
"Vector lengths do not match");
2415 res = A.
e(0) * B.
e(0);
2416 for (
int i = 1; i < m; i++) {
2417 res += A.
e(i) * B.
e(i);
2442template <
typename Mt,
typename S,
2444 typename Rt = hila::mat_scalar_type<Mt, S>>
2445inline Rt
operator*(
const Mt &mat,
const S &rhs) {
2447 for (
int i = 0; i < Rt::rows() * Rt::columns(); i++) {
2448 res.
c[i] = mat.c[i] * rhs;
2471template <
typename Mt,
typename S,
2473 typename Rt = hila::mat_scalar_type<Mt, S>>
2474inline Rt
operator*(
const S &rhs,
const Mt &mat) {
2501template <
typename Mt,
typename S,
2503 typename Rt = hila::mat_scalar_type<Mt, S>>
2504inline Rt
operator/(
const Mt &mat,
const S &rhs) {
2506 for (
int i = 0; i < Rt::rows() * Rt::columns(); i++) {
2507 res.
c[i] = mat.c[i] / rhs;
2526template <
typename Mtype1,
typename Mtype2,
2527 std::enable_if_t<Mtype1::is_matrix() && Mtype2::is_matrix(),
int> = 0>
2530 static_assert(Mtype1::columns() == Mtype2::rows() && Mtype1::rows() == Mtype2::columns(),
2531 "mul_trace(a,b): matrix sizes are not compatible");
2532 return a.mul_trace(b);
2553template <
int n,
int m,
typename T,
typename MT>
2555 for (
int i = 0; i < n; i++)
2556 for (
int j = 0; j < m; j++) {
2558 if (i < n - 1 || j < m - 1)
2581template <
int n,
int m,
typename T,
typename MT>
2583 std::stringstream strm;
2584 strm.precision(prec);
2586 for (
int i = 0; i < n; i++)
2587 for (
int j = 0; j < m; j++) {
2589 if (i < n - 1 || j < m - 1)
2616template <
int n,
int m,
typename T,
typename MT>
2618 std::stringstream strm;
2619 strm.precision(prec);
2621 if constexpr (n == 1) {
2624 for (
int i = 0; i < n * m; i++)
2625 strm <<
' ' << prettyprint(A.
e(i), prec);
2632 std::vector<std::string> lines, columns;
2636 for (
int i = 0; i < n; i++)
2638 for (
int j = 0; j < m; j++) {
2640 for (
int i = 0; i < n; i++) {
2641 std::stringstream item;
2642 item << prettyprint(A.
e(i, j), prec);
2643 columns[i] = item.str();
2644 if (columns[i].size() > size)
2645 size = columns[i].size();
2647 for (
int i = 0; i < n; i++) {
2648 lines[i].append(size - columns[i].size(),
' ');
2649 lines[i].append(columns[i]);
2650 lines[i].append(1,
' ');
2653 for (
int i = 0; i < n - 1; i++) {
2654 strm << lines[i] <<
"]\n";
2656 strm << lines[n - 1] <<
"]";
2677template <
typename Mt, std::enable_if_t<Mt::is_matrix(),
int> = 0>
2679 return rhs.squarenorm();
2691template <
typename Mt, std::enable_if_t<Mt::is_matrix(),
int> = 0>
2700template <
typename Ntype,
typename T,
int n,
int m,
2701 std::enable_if_t<hila::is_arithmetic<T>::value,
int> = 0>
2704 for (
int i = 0; i < n * m; i++)
2705 res.
c[i] = mat.
c[i];
2709template <
typename Ntype,
typename T,
int n,
int m,
2710 std::enable_if_t<hila::is_complex<T>::value,
int> = 0>
2713 for (
int i = 0; i < n * m; i++)
2714 res.
c[i] = cast_to<Ntype>(mat.
c[i]);
2742template <
int n,
int m,
typename T,
typename MT>
2744 static_assert(n == m,
"exp() only for square matrices");
2747 hila::arithmetic_type<T> one = 1.0;
2749 r = mat * (one / order) + one;
2752 for (
int k = order - 1; k > 1; k--) {
2765#include "datatypes/diagonal_matrix.h"
2767#include "datatypes/matrix_linalg.h"
Definition of Array class.
Define type DiagonalMatrix<n,T>
T e(const int i) const
Element access - e(i) gives diagonal element i.
The main matrix type template Matrix_t. This is a base class type for "useful" types which are deriv...
Mtype & operator+=(const DiagonalMatrix< n, S > &rhs)
Addition assign operator for DiagonalMatrix to Matrix.
void set_column(int c, const Vector< n, S > &v)
get column of a matrix
void mult_by_2x2_left(int p, int q, const Matrix_t< 2, 2, R, Mt > &M)
Multiply -matrix from the left by matrix defined by sub matrix.
T trace() const
Computes Trace for Matrix.
Mtype & operator=(const Matrix_t< n, m, S, MT > &rhs)
Copy matrix assignment.
Vector< n, T > column(int c) const
Returns column vector as value at index c.
hila::arithmetic_type< T > norm() const
Calculate vector norm - sqrt of squarenorm.
void set_row(int r, const RowVector< m, S > &v)
Set row of Matrix with RowVector if types are assignable.
static constexpr int rows()
Define constant methods rows(), columns() - may be useful in template code.
static constexpr int columns()
Returns column length.
Mtype & random()
dot with matrix - matrix
const Mtype & fill(const S rhs)
Matrix fill.
void mult_by_2x2_right(int p, int q, const Matrix_t< 2, 2, R, Mt > &M)
Multiply -matrix from the right by matrix defined by sub matrix.
Mtype & operator-=(const Matrix_t< n, m, S, MT > &rhs)
Subtract assign operator Matrix to MAtrix.
Rtype operator-(const S &b, const Mtype &a)
Subtraction operator Scalar - Matrix.
Mtype sort(Vector< N, int > &permutation, hila::sort order=hila::sort::ascending) const
Sort method for Vector which returns permutation order.
Matrix< n, m, hila::arithmetic_type< T > > imag() const
Returns imaginary part of Matrix or Vector.
void set_diagonal(const Vector< n, S > &v)
Set diagonal of square matrix to Vector which is passed to the method.
Mtype & operator-=(const DiagonalMatrix< n, S > &rhs)
Subtract assign operator for DiagonalMatrix to Matrix.
Mtype sort(hila::sort order=hila::sort::ascending) const
Sort method for Vector.
Rtype operator+(const Mtype &a, const S &b)
Addition operator Matrix + scalar.
const RowVector< n, T > & transpose() const
Transpose of vector.
Mtype & gaussian_random(double width=1.0)
Fills Matrix with gaussian random elements from gaussian distribution.
const Vector< n, T > & transpose() const
Transpose of RowVector.
Matrix< n, m, R > operator*(const Mt1 &a, const Mt2 &b)
Multiplication operator.
static constexpr int size()
Returns size of Vector or square Matrix.
Mtype permute(const Vector< N, int > &permutation) const
Permute elements of vector.
Rtype operator+(const Mtype1 &a, const Mtype2 &b)
Addition operator Matrix + Matrix.
Mtype & operator=(const std::nullptr_t &z)
Zero assignment.
bool operator==(const Matrix< n, m, S > &rhs) const
Boolean operator == to determine if two matrices are exactly the same.
Matrix_t()=default
Define default constructors to ensure std::is_trivial.
Mtype & operator/=(const DiagonalMatrix< m, S > &rhs)
Divide assign operator for DiagonalMatrix to Matrix.
Mtype conj() const
Returns complex conjugate of Matrix.
Mtype permute_rows(const Vector< n, int > &permutation) const
permute rows of Matrix
static constexpr bool is_vector()
Returns true if Matrix is a vector.
T det() const
determinant function - if matrix size is < 5, use Laplace, otherwise LU
auto abs() const
Returns absolute value of Matrix.
const RowVector< m, T > & row(int r) const
Return reference to row in a matrix.
Mtype & operator=(const DiagonalMatrix< n, S > &rhs)
Assignment from diagonal matrix.
Matrix< n, p, R > outer_product(const Matrix< p, q, S > &rhs) const
Outer product.
T max() const
Find max of Matrix only for arithmetic types.
Mtype & operator*=(const S rhs)
Multiply assign operator scalar to Matrix.
Rtype dagger() const
Hermitian conjugate of matrix.
T e(const int i) const
Standard array indexing operation for vectors.
Matrix< n, m, hila::arithmetic_type< T > > real() const
Returns real part of Matrix or Vector.
bool operator!=(const Matrix< n, m, S > &rhs) const
Boolean operator != to check if matrices are exactly different.
int eigen_hermitean(DiagonalMatrix< n, Et > &eigenvalues, Matrix_t< n, n, Mt, MT > &eigenvectors, enum hila::sort sorted=hila::sort::unsorted) const
Calculate eigenvalues and -vectors of an hermitean (or symmetric) matrix.
Mtype & operator=(std::initializer_list< S > rhs)
Initializer list assignment.
const Mtype & operator+() const
Addition operator.
Mtype & operator+=(const Matrix_t< n, m, S, MT > &rhs)
Add assign operator Matrix to Matrix.
Rt operator*(const Matrix< 1, m, T1 > &A, const Matrix< n, 1, T2 > &B)
Multiplication operator RowVector * Vector.
int svd_pivot(Matrix_t< n, n, Mt, MT > &_U, DiagonalMatrix< n, Et > &_D, Matrix_t< n, n, Mt, MT > &_V, enum hila::sort sorted=hila::sort::unsorted) const
Singular value decomposition: divide matrix A = U S V*, where U,V unitary and S diagonal matrix of re...
Rtype operator+(const S &b, const Mtype &a)
Addition operator scalar + Matrix.
Mtype & operator+=(const S &rhs)
Add assign operator scalar to Matrix.
T c[n *m]
The data as a one dimensional array.
T min() const
Find min of Matrix only for arithmetic types.
T e(const int i, const int j) const
Standard array indexing operation for matrices.
T det_laplace() const
Determinant of matrix, using Laplace method.
Mtype & operator=(const S rhs)
Assignment from scalar.
Mtype & operator*=(const Matrix_t< m, p, S, MT > &rhs)
Multiply assign scalar or matrix.
Mtype & operator-=(const S rhs)
Subtract assign operator scalar to Matrix.
Mtype & operator/=(const S rhs)
Divide assign oprator scalar with matrix.
const DiagonalMatrix< n, T > & asDiagonalMatrix() const
Cast Vector to DiagonalMatrix.
Rtype operator-(const Mtype &a, const S &b)
Subtraction operator Matrix - scalar.
Rtype transpose() const
Transpose of matrix.
hila::type_mul< T, S > mul_trace(const Matrix_t< p, q, S, MT > &rm) const
Multiply with given matrix and compute trace of result.
T operator[](const int i) const
Indexing operation [] defined only for vectors.
int svd(Matrix_t< n, n, Mt, MT > &_U, DiagonalMatrix< n, Et > &_D, Matrix_t< n, n, Mt, MT > &_V, enum hila::sort sorted=hila::sort::unsorted) const
Singular value decomposition: divide matrix A = U S V*, where U,V unitary and S diagonal matrix of re...
Rtype operator-(const Mtype1 &a, const Mtype2 &b)
Subtraction operator Matrix - Matrix.
Mtype & operator*=(const DiagonalMatrix< m, S > &rhs)
Multiply assign operator for DiagonalMatrix to Matrix.
Mtype permute_columns(const Vector< m, int > &permutation) const
permute columns of Matrix
R dot(const Matrix< p, q, S > &rhs) const
Dot product.
T det_lu() const
following calculate the determinant of a square matrix det() is the generic interface,...
const Array< n, m, T > & asArray() const
Cast Matrix to Array.
Mtype operator-() const
Unary operator.
DiagonalMatrix< n, T > diagonal()
Return diagonal of square matrix.
std::string str(int prec=8, char separator=' ') const
Rtype adjoint() const
Adjoint of matrix.
hila::arithmetic_type< T > squarenorm() const
Calculate square norm - sum of squared elements.
static constexpr bool is_square()
Returns true if matrix is a square matrix.
Matrix class which defines matrix operations.
Matrix(const std::nullptr_t &z)
Zero constructor.
Matrix(const Matrix &v)=default
Default copy constructor.
Matrix(const Matrix_t< n, m, S, MT > &rhs)
Copy constructor.
Matrix_t< n, m, T, MT > exp(const Matrix_t< n, m, T, MT > &mat, const int order=20)
Calculate exp of a square matrix.
Matrix(std::initializer_list< S > rhs)
Initializer list constructor.
hila::arithmetic_type< T > base_type
std incantation for field types
Matrix(const S rhs)
Scalar constructor.
Matrix()=default
Default constructo.
std::ostream & operator<<(std::ostream &strm, const Matrix_t< n, m, T, MT > &A)
Stream operator.
Definition of Complex types.
T arg(const Complex< T > &a)
Return argument of Complex number.
This file defines all includes for HILA.
auto dagger(const Mtype &arg)
Return dagger of Matrix.
Mt operator*(const Mt &a, const Mt &b)
Multiplication operator Square Matrix * Square Matrix.
Rtype operator+(Mtype1 a, const Mtype2 &b)
Real micro-optimization Matrix + Matrix - no extra creation of variable and copy.
auto abs(const Mtype &arg)
Return absolute value Matrix or Vector.
auto norm(const Mt &rhs)
Returns vector norm of Matrix.
auto conj(const Mtype &arg)
Return conjugate Matrix or Vector.
auto squarenorm(const Mt &rhs)
Returns square norm of Matrix.
auto transpose(const Mtype &arg)
Return transposed Matrix or Vector.
auto real(const Mtype &arg)
Return real of Matrix or Vector.
auto imag(const Mtype &arg)
Return imaginary of Matrix or Vector.
auto trace(const Mtype &arg)
Return trace of square Matrix.
auto mul_trace(const Mtype1 &a, const Mtype2 &b)
Returns trace of product of two matrices.
auto adjoint(const Mtype &arg)
Return adjoint Matrix.
Invert diagonal + const. matrix using Sherman-Morrison formula.
double random()
Real valued uniform random number generator.
double gaussrand()
Gaussian random generation routine.
double gaussrand2(double &out2)
hila::gaussrand2 returns 2 Gaussian distributed random numbers with variance .
decltype(std::declval< A >()+std::declval< B >()) type_plus
std::string to_string(const Array< n, m, T > &A, int prec=8, char separator=' ')
Converts Array object to string.
type to store the return value of eigen_hermitean(): {E, U} where E is nxn DiagonalMatrix containing ...
hila::is_complex_or_arithmetic<T>::value
type to store the return combo of svd: {U, D, V} where U and V are nxn unitary / orthogonal,...