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");
105template <const
int n, const
int m,
typename T,
typename Mtype>
114 "Matrix requires Complex or arithmetic type");
117 using base_type = hila::arithmetic_type<T>;
118 using argument_type = T;
122 static constexpr bool is_matrix() {
133 return (n == 1 || m == 1);
154 template <
typename S,
int nn = n,
int mm = m,
155 std::enable_if_t<(hila::is_assignable<T &, S>::value && nn == mm),
int> = 0>
156 explicit inline Matrix_t(
const S rhs) {
157 for (
int i = 0; i < n; i++)
158 for (
int j = 0; j < n; j++) {
167 template <
typename S,
typename MT,
168 std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
170 for (
int i = 0; i < n * m; i++) {
176 inline Matrix_t(
const std::nullptr_t &z) {
177 for (
int i = 0; i < n * m; i++) {
184 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
185 inline Matrix_t(std::initializer_list<S> rhs) {
186 assert(rhs.size() == n * m &&
187 "Matrix/Vector initializer list size must match variable size");
189 for (
auto it = rhs.begin(); it != rhs.end(); it++, i++) {
196 template <
typename Tm =
Mtype,
197 std::enable_if_t<!std::is_same<Tm, Matrix<n, m, T>>::value,
int> = 0>
198 inline operator Mtype &() {
199 return *
reinterpret_cast<Mtype *
>(
this);
203 template <
typename Tm =
Mtype,
204 std::enable_if_t<!std::is_same<Tm, Matrix<n, m, T>>::value,
int> = 0>
205 inline operator const Mtype &()
const {
206 return *
reinterpret_cast<const Mtype *
>(
this);
247 template <
int q = n,
int p = m, std::enable_if_t<q == 1,
int> = 0>
252 template <
int q = n,
int p = m, std::enable_if_t<p == 1,
int> = 0>
253 static constexpr int size() {
257 template <
int q = n,
int p = m, std::enable_if_t<q == p,
int> = 0>
258 static constexpr int size() {
279 inline T
e(
const int i,
const int j)
const {
285 inline T &
e(
const int i,
const int j) const_function {
304 template <
int q = n,
int p = m, std::enable_if_t<(q == 1 || p == 1),
int> = 0>
305 inline T
e(
const int i)
const {
310 template <
int q = n,
int p = m, std::enable_if_t<(q == 1 || p == 1),
int> = 0>
311 inline T &
e(
const int i) const_function {
331 template <
int q = n,
int p = m, std::enable_if_t<(q == 1 || p == 1),
int> = 0>
336 template <
int q = n,
int p = m, std::enable_if_t<(q == 1 || p == 1),
int> = 0>
337 inline T &
operator[](
const int i) const_function {
357 for (
int i = 0; i < m; i++)
375 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
377 for (
int i = 0; i < m; i++)
393 for (
int i = 0; i < n; i++)
423 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
425 for (
int i = 0; i < n; i++)
440 static_assert(n == m,
"diagonal() method defined only for square matrices");
442 for (
int i = 0; i < n; i++)
443 res.
e(i) = (*this).e(i, i);
462 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
464 static_assert(n == m,
"set_diagonal() method defined only for square matrices");
465 for (
int i = 0; i < n; i++)
466 (*this).e(i, i) = v.
e(i);
489 #pragma hila loop_function
490 template <
int mm = m, std::enable_if_t<mm == 1,
int> = 0>
495 template <
int mm = m, std::enable_if_t<mm == 1,
int> = 0>
533 for (
int i = 0; i < n * m; i++) {
570 template <
typename S,
int n2,
int m2>
572 if constexpr (n != n2 || m != m2)
575 for (
int i = 0; i < n; i++)
576 for (
int j = 0; j < m; j++) {
577 if (
e(i, j) != rhs.
e(i, j))
592 template <
typename S>
594 return !(*
this == rhs);
612 template <
typename S,
typename MT,
613 std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
615 for (
int i = 0; i < n * m; i++) {
635 for (
int i = 0; i < n * m; i++) {
657 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value && n == m,
int> = 0>
660 for (
int i = 0; i < n; i++)
661 for (
int j = 0; j < n; j++) {
680 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
682 static_assert(n == m,
683 "Assigning DiagonalMatrix to Matrix possible only for square matrices");
685 for (
int i = 0; i < n; i++)
686 for (
int j = 0; j < n; j++) {
707 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
708 auto &
operator=(std::initializer_list<S> rhs) out_only & {
709 assert(rhs.size() == n * m &&
"Initializer list has a wrong size in assignment");
711 for (
auto it = rhs.begin(); it != rhs.end(); it++, i++) {
718 template <
typename S>
738 template <
typename S,
typename MT,
739 std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
741 for (
int i = 0; i < n * m; i++) {
762 template <
typename S,
typename MT,
763 std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
765 for (
int i = 0; i < n * m; i++) {
787 template <
typename S,
788 std::enable_if_t<hila::is_assignable<T &, hila::type_plus<T, S>>::value,
int> = 0>
791 static_assert(n == m,
"rows != columns : scalar addition possible for square matrix only!");
793 for (
int i = 0; i < n; i++) {
813 template <
typename S,
814 std::enable_if_t<hila::is_assignable<T &, hila::type_minus<T, S>>::value,
int> = 0>
816 static_assert(n == m,
817 "rows != columns : scalar subtraction possible for square matrix only!");
818 for (
int i = 0; i < n; i++) {
843 template <
int p,
typename S,
typename MT,
844 std::enable_if_t<hila::is_assignable<T &, hila::type_mul<T, S>>::value,
int> = 0>
846 static_assert(m == p,
"can't assign result of *= to lhs Matrix, because doing so "
847 "would change it's dimensions");
895 template <
typename S,
896 std::enable_if_t<hila::is_assignable<T &, hila::type_mul<T, S>>::value,
int> = 0>
898 for (
int i = 0; i < n * m; i++) {
921 template <
typename S,
922 std::enable_if_t<hila::is_assignable<T &, hila::type_div<T, S>>::value,
int> = 0>
924 for (
int i = 0; i < n * m; i++) {
937 template <
typename S,
938 std::enable_if_t<hila::is_assignable<T &, hila::type_plus<T, S>>::value,
int> = 0>
940 static_assert(n == m,
"Assigning DiagonalMatrix possible only for square matrix");
942 for (
int i = 0; i < n; i++)
954 template <
typename S,
955 std::enable_if_t<hila::is_assignable<T &, hila::type_plus<T, S>>::value,
int> = 0>
957 static_assert(n == m,
"Assigning DiagonalMatrix possible only for square matrix");
959 for (
int i = 0; i < n; i++)
973 template <
typename S,
974 std::enable_if_t<hila::is_assignable<T &, hila::type_mul<T, S>>::value,
int> = 0>
977 for (
int i = 0; i < n; i++)
978 for (
int j = 0; j < m; j++)
994 template <
typename S,
995 std::enable_if_t<hila::is_assignable<T &, hila::type_div<T, S>>::value,
int> = 0>
998 for (
int i = 0; i < n; i++)
999 for (
int j = 0; j < m; j++)
1000 e(i, j) /= rhs.e(j);
1021 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
1022 const auto &
fill(
const S rhs) out_only {
1023 for (
int i = 0; i < n * m; i++)
1036 template <
int mm = m,
1037 typename Rtype =
typename std::conditional<n == m, Matrix_t, Matrix<m, n, T>>::type,
1038 std::enable_if_t<(mm != 1 && n != 1),
int> = 0>
1041 for (
int i = 0; i < n; i++)
1042 for (
int j = 0; j < m; j++) {
1043 res.e(j, i) =
e(i, j);
1055 template <
int mm = m, std::enable_if_t<mm == 1,
int> = 0>
1067 template <
int nn = n, std::enable_if_t<nn == 1 && m != 1,
int> = 0>
1080 for (
int i = 0; i < n * m; i++) {
1093 template <
typename Rtype =
typename std::conditional<n == m, Mtype, Matrix<m, n, T>>::type>
1096 for (
int i = 0; i < n; i++)
1097 for (
int j = 0; j < m; j++) {
1098 res.e(j, i) =
::conj(
e(i, j));
1110 template <
typename Rtype =
typename std::conditional<n == m, Mtype, Matrix<m, n, T>>::type>
1125 for (
int i = 0; i < n * m; i++) {
1145 for (
int i = 0; i < m * n; i++) {
1158 for (
int i = 0; i < m * n; i++) {
1172 static_assert(n == m,
"trace not defined for non square matrices");
1174 for (
int i = 0; i < n; i++) {
1191 template <
int p,
int q,
typename S,
typename MT>
1194 static_assert(p == m && q == n,
"mul_trace(): argument matrix size mismatch");
1196 hila::type_mul<T, S> res = 0;
1197 for (
int i = 0; i < n; i++)
1198 for (
int j = 0; j < m; j++) {
1199 res +=
e(i, j) * rm.
e(j, i);
1210 hila::arithmetic_type<T> result(0);
1211 for (
int i = 0; i < n * m; i++) {
1223 template <
typename S = T,
1224 std::enable_if_t<hila::is_floating_point<hila::arithmetic_type<S>>::value,
int> = 0>
1225 hila::arithmetic_type<T>
norm()
const {
1229 template <
typename S = T,
1230 std::enable_if_t<!hila::is_floating_point<hila::arithmetic_type<S>>::value,
int> = 0>
1231 double norm()
const {
1232 return sqrt(
static_cast<double>(
squarenorm()));
1238 template <typename S = T, std::enable_if_t<hila::is_arithmetic<S>::value,
int> = 0>
1241 for (
int i = 1; i < n * m; i++) {
1251 template <typename S = T, std::enable_if_t<hila::is_arithmetic<S>::value,
int> = 0>
1254 for (
int i = 1; i < n * m; i++) {
1261 auto max_abs()
const {
1262 hila::arithmetic_type<T> tres, res = 0;
1263 for (
int i = 0; i < n * m; i++) {
1272 auto min_abs()
const {
1273 hila::arithmetic_type<T> tres, res =
::abs(
c[0]);
1274 for (
int i = 1; i < n * m; i++) {
1313 template <
int p,
int q,
typename S,
typename R = hila::type_mul<T, S>>
1315 static_assert(m == 1 && q == 1 && p == n,
1316 "dot() product only for vectors of the same length");
1319 for (
int i = 0; i < n; i++) {
1356 template <
int p,
int q,
typename S,
typename R = hila::type_mul<T, S>>
1358 static_assert(m == 1 && q == 1,
"outer_product() only for vectors");
1361 for (
int i = 0; i < n; i++) {
1362 for (
int j = 0; j < p; j++) {
1395 static_assert(hila::is_floating_point<hila::arithmetic_type<T>>::value,
1396 "Matrix/Vector random() requires non-integral type elements");
1398 for (
int i = 0; i < n * m; i++) {
1417 static_assert(hila::is_floating_point<hila::arithmetic_type<T>>::value,
1418 "Matrix/Vector gaussian_random() requires non-integral type elements");
1421 if constexpr (hila::is_complex<T>::value) {
1422 for (
int i = 0; i < n * m; i++) {
1423 c[i].gaussian_random(width);
1430 for (
int i = 0; i < n * m - 1; i += 2) {
1432 c[i + 1] = gr * width;
1434 if constexpr ((n * m) % 2 > 0) {
1450 for (
int i = 0; i < m; i++)
1451 res.
set_column(i, this->column(permutation[i]));
1464 for (
int i = 0; i < n; i++)
1465 res.
set_row(i, this->row(permutation[i]));
1481 "permute() only for vectors, use permute_rows() or permute_columns() for matrices");
1482 static_assert(N ==
Mtype::size(),
"Incorrect size of permutation vector");
1485 for (
int i = 0; i < N; i++) {
1486 res[i] = (*this)[permutation[i]];
1517#pragma hila novector
1521 static_assert(n == 1 || m == 1,
"Sorting possible only for vectors");
1523 "Sorting possible only for arithmetic vector elements");
1524 static_assert(N ==
Mtype::size(),
"Incorrect size of permutation vector");
1526 for (
int i = 0; i < N; i++)
1528 if (hila::sort::unsorted == order) {
1532 if (hila::sort::ascending == order) {
1533 for (
int i = 1; i < N; i++) {
1534 for (
int j = i; j > 0 &&
c[permutation[j - 1]] >
c[permutation[j]]; j--)
1535 hila::swap(permutation[j], permutation[j - 1]);
1538 for (
int i = 1; i < N; i++) {
1539 for (
int j = i; j > 0 &&
c[permutation[j - 1]] <
c[permutation[j]]; j--)
1540 hila::swap(permutation[j], permutation[j - 1]);
1544 return this->
permute(permutation);
1566#pragma hila novector
1567 Mtype sort(hila::sort order = hila::sort::ascending)
const {
1568 static_assert(n == 1 || m == 1,
"Sorting possible only for vectors");
1571 return sort(permutation, order);
1594 template <
typename R,
typename Mt>
1596 static_assert(hila::contains_complex<T>::value || !hila::contains_complex<R>::value,
1597 "Cannot multiply real matrix with complex 2x2 block matrix");
1600 for (
int i = 0; i < m; i++) {
1601 a.
e(0) = this->
e(p, i);
1602 a.
e(1) = this->
e(q, i);
1606 this->
e(p, i) = a.e(0);
1607 this->
e(q, i) = a.e(1);
1623 template <
typename R,
typename Mt>
1625 static_assert(hila::contains_complex<T>::value || !hila::contains_complex<R>::value,
1626 "Cannot multiply real matrix with complex 2x2 block matrix");
1629 for (
int i = 0; i < n; i++) {
1630 a.
e(0) = this->
e(i, p);
1631 a.
e(1) = this->
e(i, q);
1635 this->
e(i, p) = a.
e(0);
1636 this->
e(i, q) = a.
e(1);
1649#pragma hila novector
1651#pragma hila novector
1653#pragma hila novector
1666#pragma hila novector
1667 template <
typename Et,
typename Mt,
typename MT>
1670 enum hila::sort sorted = hila::sort::unsorted)
const;
1675 template <
typename Et,
typename Mt,
typename MT>
1678 enum hila::sort sorted = hila::sort::unsorted)
const {
1682 eigenvaluevec = d.asArray().asVector();
1687#pragma hila novector
1688 template <
typename Et,
typename Mt,
typename MT>
1691 enum hila::sort sorted = hila::sort::unsorted)
const;
1696#pragma hila novector
1697 template <
typename Et,
typename Mt,
typename MT>
1700 enum hila::sort sorted = hila::sort::unsorted)
const;
1709 std::string
str(
int prec = 8,
char separator =
' ')
const {
1710 std::stringstream text;
1711 text.precision(prec);
1712 for (
int i = 0; i < n; i++) {
1713 for (
int j = 0; j < m; j++) {
1715 if (i < n - 1 || j < m - 1)
1741template <
int n,
int m,
typename T>
1747 using argument_type = T;
1786 template <
typename S,
int nn = n,
int mm = m,
1787 std::enable_if_t<(hila::is_assignable<T &, S>::value && nn == mm),
int> = 0>
1789 for (
int i = 0; i < n; i++)
1790 for (
int j = 0; j < n; j++) {
1792 this->
e(i, j) = rhs;
1816 template <
typename S,
typename MT,
1817 std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
1819 for (
int i = 0; i < n * m; i++) {
1820 this->c[i] = rhs.c[i];
1831 for (
int i = 0; i < n * m; i++) {
1850 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
1852 assert(rhs.size() == n * m &&
1853 "Matrix/Vector initializer list size must match variable size");
1855 for (
auto it = rhs.begin(); it != rhs.end(); it++, i++) {
1863 Matrix &operator=(
const Matrix &v) out_only & =
default;
1876template <
typename T1,
typename T2>
1877struct matrix_op_type_s {
1878 using type =
Matrix<T1::rows(), T1::columns(), hila::ntype_op<T1, T2>>;
1882template <
typename T>
1883struct matrix_op_type_s<T, T> {
1888template <
typename T1,
typename T2>
1889using mat_x_mat_type =
typename matrix_op_type_s<T1, T2>::type;
1899template <
typename Mt,
typename S,
typename Enable =
void>
1900struct matrix_scalar_op_s {
1902 Matrix<Mt::rows(), Mt::columns(),
1906template <
typename Mt,
typename S>
1907struct matrix_scalar_op_s<
1909 typename std::enable_if_t<std::is_convertible<hila::type_plus<hila::number_type<Mt>, S>,
1910 hila::number_type<Mt>>::value>> {
1912 using type =
typename std::conditional<
1913 hila::is_floating_point<hila::arithmetic_type<Mt>>::value, Mt,
1914 Matrix<Mt::rows(), Mt::columns(),
1918template <
typename Mt,
typename S>
1919using mat_scalar_type =
typename matrix_scalar_op_s<Mt, S>::type;
1941template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
1943 return arg.transpose();
1961template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
1981template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
1983 return arg.adjoint();
1997template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
1999 return arg.adjoint();
2017template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
2038template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
2058template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
2078template <
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0>
2107template <
typename Mtype1,
typename Mtype2,
2108 std::enable_if_t<Mtype1::is_matrix() && Mtype2::is_matrix(),
int> = 0,
2109 typename Rtype = hila::mat_x_mat_type<Mtype1, Mtype2>,
2110 std::enable_if_t<!std::is_same<Mtype1, Rtype>::value,
int> = 0>
2113 constexpr int n = Mtype1::rows();
2114 constexpr int m = Mtype1::columns();
2116 static_assert(n == Mtype2::rows() && m == Mtype2::columns(),
"Matrix sizes do not match");
2119 for (
int i = 0; i < n; i++)
2120 for (
int j = 0; j < m; j++)
2121 r.e(i, j) = a.e(i, j) + b.e(i, j);
2137template <
typename Mtype1,
typename Mtype2,
2138 std::enable_if_t<Mtype1::is_matrix() && Mtype2::is_matrix(),
int> = 0,
2139 typename Rtype = hila::mat_x_mat_type<Mtype1, Mtype2>,
2140 std::enable_if_t<std::is_same<Mtype1, Rtype>::value,
int> = 0>
2143 constexpr int n = Mtype1::rows();
2144 constexpr int m = Mtype1::columns();
2146 static_assert(n == Mtype2::rows() && m == Mtype2::columns(),
"Matrix sizes do not match");
2148 for (
int i = 0; i < n; i++)
2149 for (
int j = 0; j < m; j++)
2150 a.e(i, j) += b.e(i, j);
2175template <
typename Mtype1,
typename Mtype2,
2176 std::enable_if_t<Mtype1::is_matrix() && Mtype2::is_matrix(),
int> = 0,
2177 typename Rtype = hila::mat_x_mat_type<Mtype1, Mtype2>>
2180 constexpr int n = Mtype1::rows();
2181 constexpr int m = Mtype1::columns();
2183 static_assert(n == Mtype2::rows() && m == Mtype2::columns(),
"Matrix sizes do not match");
2186 for (
int i = 0; i < n; i++)
2187 for (
int j = 0; j < m; j++)
2188 r.e(i, j) = a.e(i, j) - b.e(i, j);
2213template <
typename Mtype,
typename S,
2215 typename Rtype = hila::mat_scalar_type<Mtype, S>>
2219 "Matrix + scalar possible only for square matrix");
2222 for (
int i = 0; i < Rtype::rows(); i++)
2223 for (
int j = 0; j < Rtype::columns(); j++) {
2224 r.e(i, j) = a.
e(i, j);
2250template <
typename Mtype,
typename S,
2252 typename Rtype = hila::mat_scalar_type<Mtype, S>>
2256 "Matrix + scalar possible only for square matrix");
2282template <
typename Mtype,
typename S,
2284 typename Rtype = hila::mat_scalar_type<Mtype, S>>
2288 "Matrix - scalar possible only for square matrix");
2291 for (
int i = 0; i < Rtype::rows(); i++)
2292 for (
int j = 0; j < Rtype::columns(); j++) {
2293 r.e(i, j) = a.
e(i, j);
2320template <
typename Mtype,
typename S,
2322 typename Rtype = hila::mat_scalar_type<Mtype, S>>
2326 "Scalar - matrix possible only for square matrix");
2329 for (
int i = 0; i < Rtype::rows(); i++)
2330 for (
int j = 0; j < Rtype::columns(); j++) {
2331 r.e(i, j) = -a.
e(i, j);
2350template <
typename Mt, std::enable_if_t<Mt::is_matrix() && Mt::rows() == Mt::columns(),
int> = 0>
2353 constexpr int n = Mt::rows();
2357 for (
int i = 0; i < n; i++)
2358 for (
int j = 0; j < n; j++) {
2359 res.e(i, j) = a.e(i, 0) * b.e(0, j);
2360 for (
int k = 1; k < n; k++) {
2361 res.e(i, j) += a.e(i, k) * b.e(k, j);
2401template <
typename Mt1,
typename Mt2,
2402 std::enable_if_t<Mt1::is_matrix() && Mt2::is_matrix() && !std::is_same<Mt1, Mt2>::value,
2404 typename R = hila::type_mul<hila::number_type<Mt1>, hila::number_type<Mt2>>,
2405 int n = Mt1::rows(),
int m = Mt2::columns()>
2408 constexpr int p = Mt1::columns();
2409 static_assert(p == Mt2::rows(),
"Matrix size: LHS columns != RHS rows");
2413 if constexpr (n > 1 && m > 1) {
2414 for (
int i = 0; i < n; i++)
2415 for (
int j = 0; j < m; j++) {
2416 res.
e(i, j) = a.e(i, 0) * b.e(0, j);
2417 for (
int k = 1; k < p; k++) {
2418 res.
e(i, j) += a.e(i, k) * b.e(k, j);
2421 }
else if constexpr (m == 1) {
2423 for (
int i = 0; i < n; i++) {
2424 res.
e(i) = a.e(i, 0) * b.e(0);
2425 for (
int k = 1; k < p; k++) {
2426 res.
e(i) += a.e(i, k) * b.e(k);
2429 }
else if constexpr (n == 1) {
2431 for (
int j = 0; j < m; j++) {
2432 res.
e(j) = a.e(0) * b.e(0, j);
2433 for (
int k = 1; k < p; k++) {
2434 res.
e(j) += a.e(k) * b.e(k, j);
2467template <
int m,
int n,
typename T1,
typename T2,
typename Rt = hila::ntype_op<T1, T2>>
2470 static_assert(m == n,
"Vector lengths do not match");
2473 res = A.
e(0) * B.
e(0);
2474 for (
int i = 1; i < m; i++) {
2475 res += A.
e(i) * B.
e(i);
2500template <
typename Mt,
typename S,
2502 typename Rt = hila::mat_scalar_type<Mt, S>>
2503inline Rt
operator*(
const Mt &mat,
const S &rhs) {
2505 for (
int i = 0; i < Rt::rows() * Rt::columns(); i++) {
2506 res.
c[i] = mat.c[i] * rhs;
2529template <
typename Mt,
typename S,
2531 typename Rt = hila::mat_scalar_type<Mt, S>>
2532inline Rt
operator*(
const S &rhs,
const Mt &mat) {
2559template <
typename Mt,
typename S,
2561 typename Rt = hila::mat_scalar_type<Mt, S>>
2562inline Rt
operator/(
const Mt &mat,
const S &rhs) {
2564 for (
int i = 0; i < Rt::rows() * Rt::columns(); i++) {
2565 res.
c[i] = mat.c[i] / rhs;
2584template <
typename Mtype1,
typename Mtype2,
2585 std::enable_if_t<Mtype1::is_matrix() && Mtype2::is_matrix(),
int> = 0>
2588 static_assert(Mtype1::columns() == Mtype2::rows() && Mtype1::rows() == Mtype2::columns(),
2589 "mul_trace(a,b): matrix sizes are not compatible");
2590 return a.mul_trace(b);
2602template <
typename Mt1,
typename Mt2,
typename Mt3,
2603 std::enable_if_t<Mt1::is_matrix() && Mt2::is_matrix() && Mt3::is_matrix(),
int> = 0>
2604inline void mult(
const Mt1 &a,
const Mt2 &b, out_only Mt3 &res) {
2605 static_assert(Mt1::columns() == Mt2::rows() && Mt1::rows() == Mt3::rows() &&
2606 Mt2::columns() == Mt3::columns(),
2607 "mult(a,b,c): matrix sizes are not compatible");
2608 constexpr int n = Mt1::rows();
2609 constexpr int m = Mt2::columns();
2610 constexpr int l = Mt2::rows();
2612 for (i = 0; i < n; ++i) {
2613 for (j = 0; j < m; ++j) {
2614 res.e(i, j) = a.e(i, 0) * b.e(0, j);
2615 for (k = 1; k < l; ++k) {
2616 res.e(i, j) += a.e(i, k) * b.e(k, j);
2631template <
typename Mt1,
typename Mt2,
typename Mt3,
2632 std::enable_if_t<Mt1::is_matrix() && Mt2::is_matrix() && Mt3::is_matrix(),
int> = 0>
2633inline void mult_aa(
const Mt1 &a,
const Mt2 &b, out_only Mt3 &res) {
2634 static_assert(Mt1::columns() == Mt2::rows() && Mt1::rows() == Mt3::columns() &&
2635 Mt2::columns() == Mt3::rows(),
2636 "mult_aa(a,b,c): matrix sizes are not compatible");
2637 constexpr int n = Mt1::rows();
2638 constexpr int m = Mt2::columns();
2639 constexpr int l = Mt2::rows();
2641 for (i = 0; i < n; ++i) {
2642 for (j = 0; j < m; ++j) {
2643 res.e(j, i) =
conj(a.e(i, 0) * b.e(0, j));
2644 for (k = 1; k < l; ++k) {
2645 res.e(j, i) +=
conj(a.e(i, k) * b.e(k, j));
2660template <
typename Mt1,
typename S,
typename Mt2,
2661 std::enable_if_t<(Mt1::is_matrix() && Mt2::is_matrix() &&
2664inline void mult(
const Mt1 &a,
const S &b, out_only Mt2 &res) {
2665 static_assert(Mt1::columns() == Mt2::columns() && Mt1::rows() == Mt2::rows(),
2666 "mult(a,b,c): matrix sizes are not compatible");
2667 constexpr int n = Mt1::rows();
2668 constexpr int m = Mt1::columns();
2670 for (i = 0; i < n; ++i) {
2671 for (j = 0; j < m; ++j) {
2672 res.e(i, j) = a.e(i, j) * b;
2686template <
typename S,
typename Mt1,
typename Mt2,
2687 std::enable_if_t<(Mt1::is_matrix() && Mt2::is_matrix() &&
2690inline void mult(
const S &a,
const Mt1 &b, out_only Mt2 &res) {
2691 static_assert(Mt1::columns() == Mt2::columns() && Mt1::rows() == Mt2::rows(),
2692 "mult(a,b,c): matrix sizes are not compatible");
2693 constexpr int n = Mt1::rows();
2694 constexpr int m = Mt1::columns();
2696 for (i = 0; i < n; ++i) {
2697 for (j = 0; j < m; ++j) {
2698 res.e(i, j) = b.e(i, j) * a;
2712template <
typename Mt1,
typename Mt2,
typename Mt3,
2713 std::enable_if_t<Mt1::is_matrix() && Mt2::is_matrix() && Mt3::is_matrix(),
int> = 0>
2714inline void mult_add(
const Mt1 &a,
const Mt2 &b, Mt3 &res) {
2715 static_assert(Mt1::columns() == Mt2::rows() && Mt1::rows() == Mt3::rows() &&
2716 Mt2::columns() == Mt3::columns(),
2717 "mult_add(a,b,c): matrix sizes are not compatible");
2718 constexpr int n = Mt1::rows();
2719 constexpr int m = Mt2::columns();
2720 constexpr int l = Mt2::rows();
2722 for (i = 0; i < n; ++i) {
2723 for (j = 0; j < m; ++j) {
2724 for (k = 0; k < l; ++k) {
2725 res.e(i, j) += a.e(i, k) * b.e(k, j);
2740template <
typename Mt1,
typename S,
typename Mt2,
2741 std::enable_if_t<(Mt1::is_matrix() && Mt2::is_matrix() &&
2744inline void mult_add(
const Mt1 &a,
const S &b, Mt2 &res) {
2745 static_assert(Mt1::columns() == Mt2::columns() && Mt1::rows() == Mt2::rows(),
2746 "mult_add(a,b,c): matrix sizes are not compatible");
2747 constexpr int n = Mt1::rows();
2748 constexpr int m = Mt1::columns();
2750 for (i = 0; i < n; ++i) {
2751 for (j = 0; j < m; ++j) {
2752 res.e(i, j) += a.e(i, j) * b;
2766template <
typename S,
typename Mt1,
typename Mt2,
2767 std::enable_if_t<(Mt1::is_matrix() && Mt2::is_matrix() &&
2770inline void mult_add(
const S &a,
const Mt1 &b, Mt2 &res) {
2771 static_assert(Mt1::columns() == Mt2::columns() && Mt1::rows() == Mt2::rows(),
2772 "mult_add(a,b,c): matrix sizes are not compatible");
2773 constexpr int n = Mt1::rows();
2774 constexpr int m = Mt1::columns();
2776 for (i = 0; i < n; ++i) {
2777 for (j = 0; j < m; ++j) {
2778 res.e(i, j) += b.e(i, j) * a;
2793template <
typename Mt1,
typename Mt2,
typename Mt3,
2794 std::enable_if_t<Mt1::is_matrix() && Mt2::is_matrix() && Mt3::is_matrix(),
int> = 0>
2795inline void mult_sub(
const Mt1 &a,
const Mt2 &b, Mt3 &res) {
2796 static_assert(Mt1::columns() == Mt2::rows() && Mt1::rows() == Mt3::rows() &&
2797 Mt2::columns() == Mt3::columns(),
2798 "mult_sub(a,b,c): matrix sizes are not compatible");
2799 constexpr int n = Mt1::rows();
2800 constexpr int m = Mt2::columns();
2801 constexpr int l = Mt2::rows();
2803 for (i = 0; i < n; ++i) {
2804 for (j = 0; j < m; ++j) {
2805 for (k = 0; k < l; ++k) {
2806 res.e(i, j) -= a.e(i, k) * b.e(k, j);
2821template <
typename Mt1,
typename S,
typename Mt2,
2822 std::enable_if_t<(Mt1::is_matrix() && Mt2::is_matrix() &&
2825inline void mult_sub(
const Mt1 &a,
const S &b, Mt2 &res) {
2826 static_assert(Mt1::columns() == Mt2::columns() && Mt1::rows() == Mt2::rows(),
2827 "mult_sub(a,b,c): matrix sizes are not compatible");
2828 constexpr int n = Mt1::rows();
2829 constexpr int m = Mt1::columns();
2831 for (i = 0; i < n; ++i) {
2832 for (j = 0; j < m; ++j) {
2833 res.e(i, j) -= a.e(i, j) * b;
2847template <
typename S,
typename Mt1,
typename Mt2,
2848 std::enable_if_t<(Mt1::is_matrix() && Mt2::is_matrix() &&
2851inline void mult_sub(
const S &a,
const Mt1 &b, Mt2 &res) {
2852 static_assert(Mt1::columns() == Mt2::columns() && Mt1::rows() == Mt2::rows(),
2853 "mult_sub(a,b,c): matrix sizes are not compatible");
2854 constexpr int n = Mt1::rows();
2855 constexpr int m = Mt1::columns();
2857 for (i = 0; i < n; ++i) {
2858 for (j = 0; j < m; ++j) {
2859 res.e(i, j) -= b.e(i, j) * a;
2882template <
int n,
int m,
typename T,
typename MT>
2884 for (
int i = 0; i < n; i++)
2885 for (
int j = 0; j < m; j++) {
2887 if (i < n - 1 || j < m - 1)
2904template <
typename Ntype,
typename T,
int n,
int m,
2905 std::enable_if_t<hila::is_arithmetic_or_extended<T>::value,
int> = 0>
2908 for (
int i = 0; i < n * m; i++)
2909 res.
c[i] = cast_to<Ntype>(mat.
c[i]);
2913template <
typename Ntype,
typename T,
int n,
int m,
2914 std::enable_if_t<hila::is_complex<T>::value,
int> = 0>
2917 for (
int i = 0; i < n * m; i++)
2918 res.
c[i] = cast_to<Ntype>(mat.
c[i]);
2936template <
int n,
int m,
typename T,
typename MT>
2938 std::stringstream strm;
2939 strm.precision(prec);
2941 for (
int i = 0; i < n; i++)
2942 for (
int j = 0; j < m; j++) {
2944 if (i < n - 1 || j < m - 1)
2971template <
int n,
int m,
typename T,
typename MT>
2973 std::stringstream strm;
2974 strm.precision(prec);
2976 if constexpr (n == 1) {
2979 for (
int i = 0; i < n * m; i++)
2980 strm <<
' ' << prettyprint(A.
e(i), prec);
2987 std::vector<std::string> lines, columns;
2991 for (
int i = 0; i < n; i++)
2993 for (
int j = 0; j < m; j++) {
2995 for (
int i = 0; i < n; i++) {
2996 std::stringstream item;
2997 item << prettyprint(A.
e(i, j), prec);
2998 columns[i] = item.str();
2999 if (columns[i].size() > size)
3000 size = columns[i].size();
3002 for (
int i = 0; i < n; i++) {
3003 lines[i].append(size - columns[i].size(),
' ');
3004 lines[i].append(columns[i]);
3005 lines[i].append(1,
' ');
3008 for (
int i = 0; i < n - 1; i++) {
3009 strm << lines[i] <<
"]\n";
3011 strm << lines[n - 1] <<
"]";
3032template <
typename Mt, std::enable_if_t<Mt::is_matrix(),
int> = 0>
3034 return rhs.squarenorm();
3046template <
typename Mt, std::enable_if_t<Mt::is_matrix(),
int> = 0>
3067template <
int n,
int m,
typename T,
typename MT>
3069 static_assert(n == m,
"exp() only for square matrices");
3071 hila::arithmetic_type<T> one = 1.0;
3075 for (
int k = order; k > 1; k--) {
3101template <
int n,
int m,
typename T,
typename MT,
typename atype = hila::arithmetic_type<T>>
3103 static_assert(n == m,
"exp() only for square matrices");
3109 for (
int k = 2; k < n * 15; ++k) {
3115 if(rsqnorm == orsqnorm) {
3134template <
int n,
int m,
typename T,
typename MT,
typename atype = hila::arithmetic_type<T>>
3136 static_assert(n == m,
"exp() only for square matrices");
3141 for (
int k = 2; k < n * 15; ++k) {
3147 if (rsqnorm == orsqnorm) {
3169template <
int n,
int m,
typename T,
typename MT>
3172 const int order = 20) {
3173 static_assert(n == m,
"mult_exp() only for square matrices");
3176 hila::arithmetic_type<T> one = 1.0;
3180 for (
int k = order; k > 1; k--) {
3213template <
int n,
int m,
typename T,
typename MT,
typename Mt,
3214 std::enable_if_t<Mt::is_matrix() && Mt::is_square() && Mt::rows() == n,
int> = 0>
3216 Mt(out_only &pl)[n]) {
3217 static_assert(n == m,
"chexp() only for square matrices");
3219 hila::arithmetic_type<T> sclim = sqrt(2.0), matnorm=
norm(mat), sfac=1.0;
3221 while (matnorm * sfac >= sclim) {
3231 trpl[1] =
trace(pl[1]);
3233 for (i = 2; i < n; ++i) {
3236 mult(pl[j], pl[j + k], pl[i]);
3237 trpl[i] =
trace(pl[i]);
3246 for (j = 1; j <= n; ++j) {
3247 crpl[n - j] = -trpl[j];
3248 for (i = 1; i < j; ++i) {
3249 crpl[n - j] -= crpl[n - (j - i)] * trpl[i];
3257 hila::arithmetic_type<T>
3262 for (i = 0; i < n; ++i) {
3275 hila::arithmetic_type<T> ttwpf;
3276 hila::arithmetic_type<T> s, rs = 1.0, rss;
3277 for (j = n; j < mmax; ++j) {
3279 cho = pal[n - 1] * rs;
3280 for (i = n - 1; i > 0; --i) {
3281 pal[i] = pal[i - 1] * rs - cho * crpl[i];
3283 al[i] += wpf * pal[i];
3285 pal[0] = -cho * crpl[0];
3287 al[0] += wpf * pal[0];
3303 if (ttwpf == twpf) {
3313 for (k = 0; k < nb; ++k) {
3315 for(i=0; i<n; ++i) {
3316 trpl[i] = pal[i] = al[i];
3319 for (j = 1; j < n; ++j) {
3321 for (i = n - 1; i > 0; --i) {
3322 pal[i] = pal[i - 1] - cho * crpl[i];
3323 al[i] += trpl[j] * pal[i];
3325 pal[0] = -cho * crpl[0];
3326 al[0] += trpl[j] * pal[0];
3332 for (k = 1; k < n;++k) {
3340template <
int n,
int m,
typename T,
typename MT,
typename Mt,
3341 std::enable_if_t<Mt::is_matrix() && Mt::is_square() && Mt::rows()==n,
int> = 0>
3343 Mt(out_only &pl)[n]) {
3344 static_assert(n == m,
"chexp() only for square matrices");
3345 chexp(mat, pl[0], pl);
3350template <
int n,
int m,
typename T,
typename MT>
3352 static_assert(n == m,
"chexp() only for square matrices");
3354 return chexp(mat, omat, pl);
3359template <
int n,
int m,
typename T,
typename MT>
3361 static_assert(n == m,
"chexp() only for square matrices");
3363 chexp(mat, pl[0], pl);
3368template <
int n,
int m,
typename T,
typename MT,
typename Mt,
3369 std::enable_if_t<Mt::is_matrix() && Mt::is_square() && Mt::rows() == n,
int> = 0>
3371 static_assert(n == m,
"chexp() only for square matrices");
3372 niter =
chexp(mat, pl[0], pl);
3378template <
int n,
int m,
typename T,
typename MT>
3380 static_assert(n == m,
"chexp() only for square matrices");
3382 niter =
chexp(mat, pl[0], pl);
3401template <
int n,
int m,
typename T,
typename MT,
typename Mt,
3402 std::enable_if_t<Mt::is_matrix() && Mt::is_square() && Mt::rows() == n,
int> = 0>
3404 Mt(out_only &domat)[n][m]) {
3405 static_assert(n == m,
"chexp() only for square matrices");
3407 hila::arithmetic_type<T> sclim = sqrt(2.0), matnorm =
norm(mat), sfac = 1.0;
3409 while (matnorm * sfac >= sclim) {
3420 trpl[1] =
trace(pl[1]);
3422 for (i = 2; i < n; ++i) {
3425 mult(pl[j], pl[j + k], pl[i]);
3426 trpl[i] =
trace(pl[i]);
3435 for (j = 1; j <= n; ++j) {
3436 crpl[n - j] = -trpl[j];
3437 for (i = 1; i < j; ++i) {
3438 crpl[n - j] -= crpl[n - (j - i)] * trpl[i];
3443 const int mmax = 25 * n;
3445 hila::arithmetic_type<T>
3453 for(i = 0; i < n; ++i) {
3458 for (j = i - k; j <= i; ++j) {
3459 kmats.
e(i - j, j) = wpf;
3465 for (i = n - 1 - k; i < n; ++i) {
3466 kh.e(n - 1 - i, i) = 1;
3474 hila::arithmetic_type<T> ttwpf;
3475 hila::arithmetic_type<T> s, rs = 1.0, rss;
3476 for (j = n; j < mmax; ++j) {
3478 cho = pal[n - 1] * rs;
3479 for (i = n - 1; i > 0; --i) {
3480 pal[i] = pal[i - 1] * rs - cho * crpl[i];
3482 al[i] += wpf * pal[i];
3484 pal[0] = -cho * crpl[0];
3486 al[0] += wpf * pal[0];
3501 if (ttwpf == twpf) {
3508 for (i = n - 1; i >= 0; --i) {
3509 cho = kh.e(i, n - 1) * rs;
3510 for (k = n - 1; k > i; --k) {
3511 kh.e(i, k) = kh.e(i, k - 1) * rs - cho * crpl[k];
3512 kmats.
e(i, k) += wpf * kh.e(i, k);
3515 kh.e(i, i) = kh.e(i - 1, i) * rs - cho * crpl[i];
3517 kh.e(i, i) = pal[i] * rs - cho * crpl[i];
3519 kmats.
e(i, i) += wpf * kh.e(i, i);
3525 for (k = 0; k < nb; ++k) {
3527 for (i = 0; i < n; ++i) {
3528 trpl[i] = pal[i] = al[i];
3529 kh.e(i, i) = 0.5 * kmats.
e(i, i);
3530 kmats.
e(i, i) = 2.0 * kh.e(0, i) * al[i];
3531 for (j = i + 1; j < n; ++j) {
3532 kh.e(j, i) = kh.e(i, j) =
3533 0.5 * kmats.
e(i, j);
3535 kmats.
e(i, j) = kh.e(0, i) * al[j] + kh.e(0, j) * al[i];
3539 for (l = 1; l < n; ++l) {
3541 for (i = n - 1; i > 0; --i) {
3542 pal[i] = pal[i - 1] - cho * crpl[i];
3543 al[i] += trpl[l] * pal[i];
3544 for (j = i; j < n; ++j) {
3545 kmats.
e(i, j) += kh.e(i, l) * pal[j] + kh.e(l, j) * pal[i];
3548 pal[0] = -cho * crpl[0];
3549 al[0] += trpl[l] * pal[0];
3550 for (j = 0; j < n; ++j) {
3551 kmats.
e(0, j) += kh.e(0, l) * pal[j] + kh.e(l, j) * pal[0];
3559 for (k = 1; k < n; ++k) {
3573 for (j = 1; j < n; ++j) {
3577 for (ic1 = 0; ic1 < n; ++ic1) {
3578 for (ic2 = 0; ic2 < n; ++ic2) {
3579 Mt &tdomat = domat[ic1][ic2];
3581 for (k = 0; k < n; ++k) {
3582 tdomat.e(k, ic1) += kh.e(k, ic2);
3587 for (i = 1; i < n; ++i) {
3591 for (j = 1; j < i; ++j) {
3594 for (j = i; j < n; ++j) {
3597 for (ic1 = 0; ic1 < n; ++ic1) {
3598 for (ic2 = 0; ic2 < n; ++ic2) {
3599 Mt &tdomat = domat[ic1][ic2];
3600 for (k = 0; k < n; ++k) {
3601 for (l = 0; l < n; ++l) {
3602 tdomat.e(k, l) += pl[i].
e(ic1, l) * kh.e(k, ic2);
3628template <
int n,
int m,
typename T,
typename MT>
3632 static_assert(n == m,
"mult_chexp() only for square matrices");
3634 hila::arithmetic_type<T> sclim = sqrt(2.0), matnorm =
norm(mat), sfac = 1.0;
3636 while (matnorm * sfac >= sclim) {
3648 trpl[1] =
trace(pl[1]);
3650 for (i = 2; i < n; ++i) {
3653 mult(pl[j], pl[j + k], pl[i]);
3654 trpl[i] =
trace(pl[i]);
3664 for (j = 1; j <= n; ++j) {
3665 crpl[n - j] = -trpl[j];
3666 for (i = 1; i < j; ++i) {
3667 crpl[n - j] -= crpl[n - (j - i)] * trpl[i];
3672 const int mmax = 25 * n;
3674 hila::arithmetic_type<T>
3682 for (i = 0; i < n; ++i) {
3687 for (j = i - k; j <= i; ++j) {
3688 kmats.
e(i - j, j) = wpf;
3694 for (i = n - 1 - k; i < n; ++i) {
3695 kh.e(n - 1 - i, i) = 1.0;
3703 hila::arithmetic_type<T> ttwpf;
3704 hila::arithmetic_type<T> s, rs = 1.0, rss;
3705 for (j = n; j < mmax; ++j) {
3707 cho = pal[n - 1] * rs;
3708 for (i = n - 1; i > 0; --i) {
3709 pal[i] = pal[i - 1] * rs - cho * crpl[i];
3711 al[i] += wpf * pal[i];
3713 pal[0] = -cho * crpl[0];
3715 al[0] += wpf * pal[0];
3731 if (ttwpf == twpf) {
3738 for (i = n-1; i >= 0; --i) {
3739 cho = kh.e(i, n - 1) * rs;
3740 for (k = n - 1; k > i; --k) {
3741 kh.e(i, k) = kh.e(i,k-1) * rs - cho * crpl[k];
3742 kmats.
e(i, k) += wpf * kh.e(i, k);
3745 kh.e(i, i) = kh.e(i - 1, i) * rs - cho * crpl[i];
3747 kh.e(i, i) = pal[i] * rs - cho * crpl[i];
3749 kmats.
e(i, i) += wpf * kh.e(i, i);
3755 for (k = 0; k < nb; ++k) {
3757 for (i = 0; i < n; ++i) {
3758 trpl[i] = pal[i] = al[i];
3759 kh.e(i, i) = 0.5 * kmats.
e(i, i);
3760 kmats.
e(i, i) = 2.0 * kh.e(0, i) * al[i];
3761 for (j = i + 1; j < n; ++j) {
3762 kh.e(j, i) = kh.e(i, j) =
3763 0.5 * kmats.
e(i, j);
3765 kmats.
e(i, j) = kh.e(0, i) * al[j] + kh.e(0, j) * al[i];
3769 for (l = 1; l < n; ++l) {
3771 for (i = n - 1; i > 0; --i) {
3772 pal[i] = pal[i - 1] - cho * crpl[i];
3773 al[i] += trpl[l] * pal[i];
3774 for (j = i; j < n; ++j) {
3775 kmats.
e(i, j) += kh.e(i, l) * pal[j] + kh.e(l, j) * pal[i];
3778 pal[0] = -cho * crpl[0];
3779 al[0] += trpl[l] * pal[0];
3780 for (j = 0; j < n; ++j) {
3781 kmats.
e(0, j) += kh.e(0, l) * pal[j] + kh.e(l, j) * pal[0];
3789 for (k = 1; k < n; ++k) {
3806 for (j = 1; j < n; ++j) {
3810 mult(tomat, kh, domat);
3812 for (i = 1; i < n; ++i) {
3816 for (j = 1; j < i; ++j) {
3819 for (j = i; j < n; ++j) {
3822 mult(pl[i], tomat, omat);
3828 mult(tomat, texp, omat);
3850template <
int n,
int m,
typename T,
typename MT>
3853 static_assert(n == m,
"chexpk() only for square matrices");
3855 hila::arithmetic_type<T> sclim = sqrt(2.0), matnorm =
norm(mat), sfac = 1.0;
3857 while (matnorm * sfac >= sclim) {
3868 trpl[1] =
trace(pl[1]);
3870 for (i = 2; i < n; ++i) {
3873 mult(pl[j], pl[j + k], pl[i]);
3874 trpl[i] =
trace(pl[i]);
3883 for (j = 1; j <= n; ++j) {
3884 crpl[n - j] = -trpl[j];
3885 for (i = 1; i < j; ++i) {
3886 crpl[n - j] -= crpl[n - (j - i)] * trpl[i];
3891 const int mmax = 25 * n;
3893 hila::arithmetic_type<T>
3902 for (i = 0; i < n; ++i) {
3907 for (j = i - k; j <= i; ++j) {
3908 kmats.
e(i - j, j) = wpf;
3914 for (i = n - 1 - k; i < n; ++i) {
3915 kh.
e(n - 1 - i, i) = 1;
3923 hila::arithmetic_type<T> ttwpf;
3924 hila::arithmetic_type<T> s, rs = 1.0, rss;
3925 int jmax = mmax - 1;
3926 for (j = n; j < mmax; ++j) {
3928 cho = pal[n - 1] * rs;
3929 for (i = n - 1; i > 0; --i) {
3930 pal[i] = pal[i - 1] * rs - cho * crpl[i];
3932 al[i] += wpf * pal[i];
3934 pal[0] = -cho * crpl[0];
3936 al[0] += wpf * pal[0];
3952 if (ttwpf == twpf) {
3960 for (i = n - 1; i >= 0; --i) {
3961 cho = kh.
e(i, n - 1) * rs;
3962 for (k = n - 1; k > i; --k) {
3963 kh.
e(i, k) = kh.
e(i, k - 1) * rs - cho * crpl[k];
3964 kmats.e(i, k) += wpf * kh.
e(i, k);
3967 kh.
e(i, i) = kh.
e(i - 1, i) * rs - cho * crpl[i];
3969 kh.
e(i, i) = pal[i] * rs - cho * crpl[i];
3971 kmats.e(i, i) += wpf * kh.
e(i, i);
3976 for (k = 0; k < nb; ++k) {
3978 for (i = 0; i < n; ++i) {
3979 trpl[i] = pal[i] = al[i];
3980 kh.
e(i, i) = 0.5 * kmats.e(i, i);
3981 kmats.e(i, i) = 2.0 * kh.
e(0, i) * al[i];
3982 for (j = i + 1; j < n; ++j) {
3983 kh.
e(j, i) = kh.
e(i, j) =
3984 0.5 * kmats.e(i, j);
3986 kmats.e(i, j) = kh.
e(0, i) * al[j] + kh.
e(0, j) * al[i];
3990 for (l = 1; l < n; ++l) {
3992 for (i = n - 1; i > 0; --i) {
3993 pal[i] = pal[i - 1] - cho * crpl[i];
3994 al[i] += trpl[l] * pal[i];
3995 for (j = i; j < n; ++j) {
3996 kmats.e(i, j) += kh.
e(i, l) * pal[j] + kh.
e(l, j) * pal[i];
3999 pal[0] = -cho * crpl[0];
4000 al[0] += trpl[l] * pal[0];
4001 for (j = 0; j < n; ++j) {
4002 kmats.e(0, j) += kh.
e(0, l) * pal[j] + kh.
e(l, j) * pal[0];
4010 for (k = 1; k < n; ++k) {
4016template <
int n,
int m,
typename T,
typename MT>
4019 static_assert(n == m,
"chexpk() only for square matrices");
4045template <
int n,
int m,
typename T,
typename MT>
4050 static_assert(n == m,
"mult_chexp() only for square matrices");
4052 hila::arithmetic_type<T> sclim = sqrt(2.0), matnorm =
norm(mat), sfac = 1.0;
4054 while (matnorm * sfac >= sclim) {
4066 for (i = 2; i < n; ++i) {
4069 mult(pl[j], pl[j + k], pl[i]);
4083 for (j = 1; j < n; ++j) {
4087 mult(tomat, kh, domat);
4089 for (i = 1; i < n; ++i) {
4093 for (j = 1; j < i; ++j) {
4096 for (j = i; j < n; ++j) {
4099 mult(pl[i], tomat, omat);
4105 mult(tomat, texp, omat);
4122template <
int n,
int m,
typename T,
typename MT>
4124 static_assert(n == m,
"chsexp() only for square matrices");
4137 for (k = 2; k <= n; ++k) {
4139 mult(mat, tB[1 - ip], tB[ip]);
4140 tc =
trace(tB[ip]) / k;
4148 hila::arithmetic_type<T> wpf = 1.0, twpf = 1.0, ttwpf;
4150 for (i = 0; i < n; ++i) {
4163 hila::arithmetic_type<T> s, rs = 1.0, rss;
4164 for (j = n; j < mmax; ++j) {
4166 ch = -pal[n - 1] * crpl[0];
4171 for (i = 1; i < n; ++i) {
4172 ch = cho - pal[n - 1] * crpl[i];
4193 if (ttwpf == twpf) {
4205 tB[ip] *= al[n - 1];
4206 tB[ip] += al[n - 2];
4207 for (i = 2; i < n; ++i) {
4208 mult(tB[ip], mat, tB[1 - ip]);
4209 tB[1 - ip] += al[n - i - 1];
4219#include "datatypes/diagonal_matrix.h"
4221#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...
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.
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.
auto & operator=(const Matrix_t< n, m, S, MT > &rhs) &
Copy matrix assignment.
static constexpr int rows()
Define constant methods rows(), columns() - may be useful in template code.
const auto & fill(const S rhs)
Matrix fill.
static constexpr int columns()
Returns column length.
Mtype & random()
dot with matrix - matrix
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.
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.
auto & operator-=(const Matrix_t< n, m, S, MT > &rhs) &
Subtract assign operator Matrix 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.
auto & operator=(const S rhs) &
Assignment from scalar.
const RowVector< n, T > & transpose() const
Transpose of vector.
auto & operator/=(const DiagonalMatrix< m, S > &rhs) &
Divide assign operator for DiagonalMatrix to Matrix.
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.
Matrix_t()=default
Define default constructors to ensure std::is_trivial.
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.
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.
Rtype dagger() const
Hermitian conjugate of matrix.
const auto & operator+() const
Addition operator.
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.
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.
const Vector< m, T > & transpose() const
Transpose of RowVector.
auto & operator+=(const S &rhs) &
Add assign operator scalar to Matrix.
bool operator==(const Matrix< n2, m2, S > &rhs) const
Boolean operator == to determine if two matrices are exactly the same.
T c[n *m]
The data as a one dimensional array.
auto & operator+=(const Matrix_t< n, m, S, MT > &rhs) &
Add assign operator Matrix to Matrix.
auto & operator-=(const DiagonalMatrix< n, S > &rhs) &
Subtract assign operator for DiagonalMatrix to Matrix.
Mtype & gaussian_random(base_type width=1.0)
Fills Matrix with gaussian random elements from gaussian distribution.
T min() const
Find min of Matrix only for arithmetic types.
auto & operator-=(const S rhs) &
Subtract assign operator scalar to Matrix.
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.
RowVector< m, T > row(int r) const
Return reference to row in a 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.
auto & operator*=(const DiagonalMatrix< m, S > &rhs) &
Multiply assign operator for DiagonalMatrix to 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.
auto & operator+=(const DiagonalMatrix< n, S > &rhs) &
Addition assign operator for DiagonalMatrix to Matrix.
auto & operator=(const DiagonalMatrix< n, S > &rhs) &
Assignment from diagonal 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,...
auto & operator=(std::initializer_list< S > rhs) &
Initializer list assignment.
auto & operator*=(const S rhs) &
Multiply assign operator scalar to Matrix.
const Array< n, m, T > & asArray() const
Cast Matrix to Array.
Mtype operator-() const
Unary operator.
auto & operator/=(const S rhs) &
Divide assign oprator scalar with matrix.
DiagonalMatrix< n, T > diagonal()
Return diagonal of square matrix.
std::string str(int prec=8, char separator=' ') const
Rtype adjoint() const
Adjoint of matrix.
auto & operator*=(const Matrix_t< m, p, S, MT > &rhs) &
Multiply assign scalar or 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.
Matrix_t< n, m, T, MT > chsexp(const Matrix_t< n, m, T, MT > &mat)
Calculate exp of a square matrix.
void mult_chexp(const Matrix_t< n, m, T, MT > &mat, const Matrix_t< n, m, T, MT > &mmat, Matrix_t< n, m, T, MT > &omat, Matrix_t< n, m, T, MT > &domat)
Calculate exp(mat).dagger()*mmat*exp(mat) and trace(exp(mat).dagger*mmat*dexp(mat))
int chexp(const Matrix_t< n, m, T, MT > &mat, Matrix_t< n, m, T, MT > &omat, Mt(&pl)[n])
Calculate exp of a square matrix.
auto dagger(const Mtype &arg)
Return dagger of Matrix.
void mult_exp(const Matrix_t< n, m, T, MT > &mat, const Matrix_t< n, m, T, MT > &mmat, Matrix_t< n, m, T, MT > &r, Matrix_t< n, m, T, MT > &dr, const int order=20)
Calculate mmat*exp(mat) and trace(mmat*dexp(mat))
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.
void mult_add(const Mt1 &a, const Mt2 &b, Mt3 &res)
compute product of two matrices and add result to existing matrix
void mult_chexpk_fast(const Matrix_t< n, m, T, MT > &mat, const Matrix_t< n, m, T, MT > &texp, const Matrix_t< n, m, T, MT > &kmats, const Matrix_t< n, m, T, MT > &mmat, Matrix_t< n, m, T, MT > &omat, Matrix_t< n, m, T, MT > &domat)
Calculate exp(mat).dagger()*mmat*exp(mat) and trace(exp(mat).dagger*mmat*dexp(mat)) from output of ch...
auto conj(const Mtype &arg)
Return conjugate Matrix or Vector.
void chexpk(const Matrix_t< n, m, T, MT > &mat, Matrix_t< n, m, T, MT > &omat, Matrix_t< n, m, T, MT > &kmats)
Calculate exp(mat) and the decomposition k_{i,j} of dexp in terms bilinears of powers of mat.
auto squarenorm(const Mt &rhs)
Returns square norm of Matrix.
void mult_sub(const Mt1 &a, const Mt2 &b, Mt3 &res)
compute product of two matrices and subtract result from existing 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.
void mult_aa(const Mt1 &a, const Mt2 &b, Mt3 &res)
compute hermitian conjugate of product of two matrices and write result to existing matrix
void mult(const Mt1 &a, const Mt2 &b, Mt3 &res)
compute product of two matrices and write result to existing matrix
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.
Matrix_t< n, m, T, MT > altexp(const Matrix_t< n, m, T, MT > &mat, int &niter)
Calculate exp of a square matrix.
Implement hila::swap for gauge fields.
double random()
Real valued uniform random number generator.
Array< n, m, Ntype > cast_to(const Array< n, m, T > &mat)
Array casting operation.
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,...