14template <
int n,
typename Mtype>
18 for (
int i = 0; i < n - 1; i++) {
19 for (
int j = i + 1; j < n; j++) {
28 return ::sqrt(abs_mpq);
42template <
typename Dtype>
55 auto t0 = c * v.
e(0) + s * v.
e(1);
62 auto t0 = v.
e(0) * c -
::conj(s) * v.
e(1);
63 v.
e(1) = v.
e(0) * s + v.
e(1) * c;
68 template <
typename Mtype>
69 void mult_by_Givens_left(
Mtype &M,
int p,
int q)
const {
71 for (
int i = 0; i < M.columns(); i++) {
75 a = this->mult_vector_left(a);
82 template <
typename Mtype>
83 void mult_by_Givens_right(
Mtype &M,
int p,
int q)
const {
85 for (
int i = 0; i < M.rows(); i++) {
89 a = this->mult_row_right(a);
137template <
typename Dtype>
138GivensMatrix<Dtype> diagonalize_2x2(
const double mpp,
const double mqq,
const Dtype mpq) {
140 GivensMatrix<Dtype> res;
144 double a = (mqq - mpp);
147 double t = 2.0 / (
abs(a) + sqrt(a * a + 4 * mpq2));
150 res.c = 1.0 / sqrt(mpq2 * t * t + 1.0);
151 res.s = mpq * (t * res.c);
198template <
int n,
int m,
typename T,
typename Mtype>
199template <
typename Et,
typename Mt,
typename MT>
202 enum hila::sort sorted)
const {
204 static_assert(!hila::contains_complex<T>::value || hila::contains_complex<Mt>::value,
205 "Eigenvector matrix must be complex with complex original matrix");
207 static_assert(n == m,
"Eigensystem can be solved only for square matrices");
210 typename std::conditional<hila::contains_complex<T>::value,
Complex<double>,
double>::type;
214 constexpr double eps = 2.22e-16;
226 eigenvalues = M.diagonal().real();
228 for (rot = 0; rot < 100 * n * n; rot++) {
232 double abs_mpq = hila::linalg::find_largest_offdiag(M, p, q);
236 if (abs_mpq <= eps * sqrt(
::abs(eigenvalues.
e(p)) *
::abs(eigenvalues.
e(q)))) {
241 auto P = hila::linalg::diagonalize_2x2(eigenvalues.
e(p), eigenvalues.
e(q), M.e(p, q));
243 P.dagger().mult_by_Givens_left(M, p, q);
244 P.mult_by_Givens_right(M, p, q);
246 eigenvalues.
e(p) =
::real(M.e(p, p));
247 eigenvalues.
e(q) =
::real(M.e(q, q));
262 P.mult_by_Givens_right(V, p, q);
266 if (sorted == hila::sort::unsorted) {
275 E = eigenvalues.
sort(perm, sorted);
305template <
int n,
int m,
typename T,
typename Mtype>
309 this->eigen_hermitean(res.eigenvalues, res.eigenvectors, sorted);
342template <
int n,
int m,
typename T,
typename Mtype>
343template <
typename Et,
typename Mt,
typename MT>
347 enum hila::sort sorted)
const {
349 static_assert(!hila::contains_complex<T>::value || hila::contains_complex<Mt>::value,
350 "SVD: diagonalizing matrix must be complex with complex original matrix");
352 static_assert(n == m,
"SVD can be solved only for square matrices");
355 typename std::conditional<hila::contains_complex<T>::value,
Complex<double>,
double>::type;
359 constexpr double eps = 2.22e-16;
371 B = (M.dagger() * M);
373 for (rot = 0; rot < 100 *
sqr(n); rot++) {
377 double abs_pq = hila::linalg::find_largest_offdiag(B, p, q);
383 auto P = hila::linalg::diagonalize_2x2(
::real(B.
e(p, p)),
::real(B.
e(q, q)), B.
e(p, q));
386 P.mult_by_Givens_right(M, p, q);
387 P.mult_by_Givens_right(V, p, q);
391 auto col = M.column(p);
392 for (
int i = 0; i < n; i++) {
394 B.
e(p, i) = col.dot(M.column(i));
399 for (
int i = 0; i < n; i++) {
401 B.
e(q, i) = col.dot(M.column(i));
411 for (
int i = 0; i < n; i++) {
412 auto col = M.column(i);
414 M.set_column(i, col / S.
e(i));
417 if (sorted == hila::sort::unsorted) {
427 _S = S.
sort(perm, sorted);
429 _U = M.permute_columns(perm);
467template <
int n,
int m,
typename T,
typename Mtype>
468template <
typename Et,
typename Mt,
typename MT>
472 enum hila::sort sorted)
const {
474 static_assert(!hila::contains_complex<T>::value || hila::contains_complex<Mt>::value,
475 "SVD: diagonalizing matrix must be complex with complex original matrix");
477 static_assert(n == m,
"SVD can be solved only for square matrices");
480 typename std::conditional<hila::contains_complex<T>::value,
Complex<double>,
double>::type;
484 constexpr double eps = 2.22e-16;
497 for (rot = 0; cont && rot < 100 *
sqr(n);) {
500 for (
int p = 0; p < n - 1; p++) {
503 for (
int q = p + 1; q < n; q++) {
512 auto colq = M.column(q);
513 double Bqq = colq.squarenorm();
515 Dtype Bpq = colp.
dot(colq);
520 if (
::abs(Bpq) > eps * sqrt(Bpp * Bqq)) {
524 auto P = hila::linalg::diagonalize_2x2(Bpp, Bqq, Bpq);
527 P.mult_by_Givens_right(M, p, q);
528 P.mult_by_Givens_right(V, p, q);
536 for (
int i = 0; i < n; i++) {
537 auto col = M.column(i);
539 M.set_column(i, col / S.
e(i));
542 if (sorted == hila::sort::unsorted) {
552 _S = S.
sort(perm, sorted);
554 _U = M.permute_columns(perm);
576template <
int n,
int m,
typename T,
typename Mtype>
580 this->svd(res.U, res.singularvalues, res.V, sorted);
584template <
int n,
int m,
typename T,
typename Mtype>
588 this->svd_pivot(res.U, res.singularvalues, res.V, sorted);
598 typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0,
600Rtype Minor(
const Mtype &bigger,
int row,
int col) {
606 for (
int i = 0; i < n; i++)
609 for (
int j = 0; j < m; j++) {
611 result.e(ri, ci) = bigger.
e(i, j);
636template <
int n,
int m,
typename T,
typename Mtype>
639 static_assert(n == m,
"determinants defined only for square matrices");
641 if constexpr (n == 1) {
644 if constexpr (n == 2) {
645 return this->e(0, 0) * this->e(1, 1) - this->e(1, 0) * this->e(0, 1);
647 if constexpr (n == 3) {
648 return this->e(0, 0) * (this->e(1, 1) * this->e(2, 2) - this->e(2, 1) * this->e(1, 2)) -
649 this->e(0, 1) * (this->e(1, 0) * this->e(2, 2) - this->e(1, 2) * this->e(2, 0)) +
650 this->e(0, 2) * (this->e(1, 0) * this->e(2, 1) - this->e(1, 1) * this->e(2, 0));
653 if constexpr (n > 3) {
655 hila::arithmetic_type<T> parity = 1, opposite = -1;
656 for (
int i = 0; i < n; i++) {
657 Matrix<n - 1, m - 1, T> minor = hila::linalg::Minor(*
this, 0, i);
658 result += parity * minor.det_laplace() * (*this).e(0, i);
676template <
int n,
int m,
typename T,
typename Mtype>
679 static_assert(n == m,
"Determinant is defined only for square matrix");
684 hila::arithmetic_type<T> d = 1, big, tmp;
688 for (
int i = 0; i < n; i++) {
690 for (
int j = 0; j < n; j++) {
700 vv[i] = 1.0 / sqrt(big);
704 for (
int j = 0; j < n; j++) {
707 for (
int i = 0; i < j; i++) {
708 for (
int k = 0; k < i; k++) {
709 a.
e(i, j) -= a.
e(i, k) * a.
e(k, j);
715 for (
int i = j; i < n; i++) {
716 auto csum = a.
e(i, j);
717 for (
int k = 0; k < j; k++) {
718 csum -= a.
e(i, k) * a.
e(k, j);
721 auto dum = vv[i] *
::abs(csum);
729 for (
int k = 0; k < n; k++) {
730 hila::swap(a.
e(imax, k), a.
e(j, k));
737 if (
::abs(a.
e(j, j)) == 0.0)
742 auto dum = 1 / a.
e(j, j);
743 for (
int i = j + 1; i < n; i++) {
751 for (
int j = 0; j < n; j++) {
764template <
int n,
int m,
typename T,
typename Mtype>
766 static_assert(n == m,
"Determinant only for square matrix");
768 return this->det_laplace();
770 return this->det_lu();
778template <
int n,
int m,
typename T,
typename Mtype>
783template <
int n,
int m,
typename T,
typename Mtype>
788template <
int n,
int m,
typename T,
typename Mtype>
817template <
int N,
typename T,
typename C,
818 std::enable_if_t<hila::is_complex_or_arithmetic<C>::value,
int> = 0>
822 auto tmul = c / (1 + c * trace(B));
824 return B - tmul * B.asVector() * B.asVector().transpose();
Array< n, m, T > conj(const Array< n, m, T > &arg)
Return conjugate Array.
hila::arithmetic_type< T > squarenorm(const Array< n, m, T > &rhs)
Return square norm of Array.
Array< n, m, hila::arithmetic_type< T > > real(const Array< n, m, T > &arg)
Return real part of Array.
Define type DiagonalMatrix<n,T>
T e(const int i) const
Element access - e(i) gives diagonal element i.
DiagonalMatrix< n, T > sort(Vector< N, int > &permutation, hila::sort order=hila::sort::ascending) const
implement sort as casting to matrix
The main matrix type template Matrix_t. This is a base class type for "useful" types which are deriv...
static constexpr int rows()
Define constant methods rows(), columns() - may be useful in template code.
static constexpr int columns()
Returns column length.
T det() const
determinant function - if matrix size is < 5, use Laplace, otherwise LU
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.
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...
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.
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...
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,...
hila::arithmetic_type< T > squarenorm() const
Calculate square norm - sum of squared elements.
Matrix class which defines matrix operations.
Complex< T > dagger(const Complex< T > &val)
Return dagger of Complex number.
T abs(const Complex< T > &a)
Return absolute value of Complex number.
constexpr T sqr(const T &arg)
Define convenience function sqr(), returning square of argument.
Definition of Matrix types.
Invert diagonal + const. matrix using Sherman-Morrison formula.
type to store the return value of eigen_hermitean(): {E, U} where E is nxn DiagonalMatrix containing ...
type to store the return combo of svd: {U, D, V} where U and V are nxn unitary / orthogonal,...