HILA
Loading...
Searching...
No Matches
matrix.h
Go to the documentation of this file.
1#ifndef HILA_MATRIX_H_
2#define HILA_MATRIX_H_
3
4/**
5 * @file matrix.h
6 * @brief Definition of Matrix types
7 * @details This file contains base matrix type Matrix_t which defines all general matrix type
8 * operations Matrix types are Matrix, #Vector, #RowVector, #SquareMatrix of which Matrix is defined
9 * as a class and the rest are special cases of the Matrix class.
10 *
11 */
12
13#include <type_traits>
14#include <sstream>
15#include "plumbing/defs.h"
16#include "datatypes/cmplx.h"
17
18// forward definitions of needed classes
19template <const int n, const int m, typename T, typename Mtype>
20class Matrix_t;
21
22template <const int n, const int m, typename T = double>
23class Array;
24
25template <int n, int m, typename T>
26class Matrix;
27
28template <int n, typename T>
29class DiagonalMatrix;
30
31/**
32 * @brief Vector is defined as 1-column Matrix
33 */
34template <int n, typename T>
36
37/**
38 * @brief RowVector is a 1-row Matrix
39 */
40template <int n, typename T>
42
43/**
44 * @brief Square matrix is defined as alias with special case of Matrix<n,n,T>
45 */
46template <int n, typename T>
48
49// template <const int n, const int m, typename T>
50// class DaggerMatrix;
51
52// Special case - m.column(), column of a matrix (not used now)
53// #include "matrix_column.h"
54
55
56/// @brief type to store the return combo of svd:
57/// {U, D, V} where U and V are nxn unitary / orthogonal,
58/// and D is real diagonal singular value matrices.
59/// @tparam M - type of input matrix
60template <typename M>
61struct svd_result {
62 static_assert(M::is_matrix() && M::rows() == M::columns(), "SVD only for square matrix");
63 M U;
64 DiagonalMatrix<M::size(), hila::arithmetic_type<M>> singularvalues;
65 M V;
66};
67
68/// @brief type to store the return value of eigen_hermitean():
69/// {E, U} where E is nxn DiagonalMatrix containing eigenvalues and
70/// U nxn unitary matrix, with eigenvector columns
71/// @tparam M - type of input matrix
72template <typename M>
74 static_assert(M::is_matrix() && M::rows() == M::columns(),
75 "Eigenvalues only for square matrix");
76 DiagonalMatrix<M::size(), hila::arithmetic_type<M>> eigenvalues;
77 M eigenvectors;
78};
79
80/**
81 * @brief The main \f$ n \times m \f$ matrix type template Matrix_t. This is a base class type for
82 * "useful" types which are derived from this.
83 *
84 * @details Uses curiously recurring template pattern (CRTP), where the last template parameter is
85 * the template itself
86 *
87 * Example: the Matrix type below is defined as
88 * @code{.cpp}
89 * template <int n, int m, typename T>
90 * class Matrix : public Matrix_t<n, m, T, Matrix<n, m, T>> { .. }
91 * @endcode
92 *
93 * This pattern is used because stupid c++ makes it complicated to write generic code, in this case
94 * derived functions to return derived type
95 *
96 * @tparam n Number of rows
97 * @tparam m Number of columns
98 * @tparam T Matrix element type
99 * @tparam Mtype Specific "Matrix" type for CRTP
100 */
101
102// Helper struct for getting the floating point number epsilons without
103// having to use std::numeric_limits .
104
105template <const int n, const int m, typename T, typename Mtype>
106class Matrix_t {
107
108 public:
109 /// The data as a one dimensional array
110 T c[n * m];
111
112 public:
113 static_assert(hila::is_complex_or_arithmetic<T>::value || hila::is_extended<T>::value,
114 "Matrix requires Complex or arithmetic type");
115
116 // std incantation for field types
117 using base_type = hila::arithmetic_type<T>;
118 using argument_type = T;
119
120 // help for templates, can use T::is_matrix()
121 // Not very useful outside template parameters
122 static constexpr bool is_matrix() {
123 return true;
124 }
125
126 /**
127 * @brief Returns true if Matrix is a vector
128 *
129 * @return true
130 * @return false
131 */
132 static constexpr bool is_vector() {
133 return (n == 1 || m == 1);
134 }
135
136 /**
137 * @brief Returns true if matrix is a square matrix
138 *
139 * @return true
140 * @return false
141 */
142 static constexpr bool is_square() {
143 return (n == m);
144 }
145
146 /// Define default constructors to ensure std::is_trivial
147 // Move constructor not needed
148 Matrix_t() = default;
149 ~Matrix_t() = default;
150 Matrix_t(const Matrix_t &v) = default;
151
152 // constructor from scalar -- keep it explicit! Not good for auto use
153 // NOTE: I forgot why this should be kept explicit
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++) {
159 if (i == j)
160 e(i, j) = rhs;
161 else
162 e(i, j) = 0;
163 }
164 }
165
166 // Construct from a different type matrix
167 template <typename S, typename MT,
168 std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
169 Matrix_t(const Matrix_t<n, m, S, MT> &rhs) out_only {
170 for (int i = 0; i < n * m; i++) {
171 c[i] = rhs.c[i];
172 }
173 }
174
175 // construct from 0
176 inline Matrix_t(const std::nullptr_t &z) {
177 for (int i = 0; i < n * m; i++) {
178 c[i] = 0;
179 }
180 }
181
182 // Construct matrix automatically from right-size initializer list
183 // This does not seem to be dangerous, so keep non-explicit
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");
188 int i = 0;
189 for (auto it = rhs.begin(); it != rhs.end(); it++, i++) {
190 c[i] = *it;
191 }
192 }
193
194 // cast to curious type
195
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);
200 }
201
202
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);
207 }
208
209 // automatically cast to generic matrix
210
211 inline operator Matrix<n, m, T> &() {
212 return *reinterpret_cast<Matrix<n, m, T> *>(this);
213 }
214
215 inline operator const Matrix<n, m, T> &() const {
216 return *reinterpret_cast<const Matrix<n, m, T> *>(this);
217 }
218
219 /// Define constant methods rows(), columns() - may be useful in template code
220
221 /**
222 * @brief Returns row length
223 *
224 * @return constexpr int
225 */
226 static constexpr int rows() {
227 return n;
228 }
229 /**
230 * @brief Returns column length
231 *
232 * @return constexpr int
233 */
234 static constexpr int columns() {
235 return m;
236 }
237
238
239 /**
240 * @brief Returns size of #Vector or square Matrix
241 *
242 * @tparam q row size n
243 * @tparam p column size m
244 * @return constexpr int
245 */
246 // size for row vector
247 template <int q = n, int p = m, std::enable_if_t<q == 1, int> = 0>
248 static constexpr int size() {
249 return p;
250 }
251 // size for column vector
252 template <int q = n, int p = m, std::enable_if_t<p == 1, int> = 0>
253 static constexpr int size() {
254 return q;
255 }
256 // size for square matrix
257 template <int q = n, int p = m, std::enable_if_t<q == p, int> = 0>
258 static constexpr int size() {
259 return q;
260 }
261
262 /**
263 * @brief Standard array indexing operation for matrices
264 * @details Accessing singular elements is insufficient, but Matrix elements are often quite
265 * small.
266 *
267 * Exammple for matrix:
268 * \code
269 * Matrix<n,m,MyType> M;
270 * MyType a = M.e(i,j); \\ i <= n, j <= m
271 * \endcode
272 *
273 *
274 * @param i Row index
275 * @param j Column index
276 * @return T matrix element type
277 */
278
279 inline T e(const int i, const int j) const {
280 // return elem[i][j];
281 return c[i * m + j];
282 }
283
284 /// @internal const_function implementation. See const_function for details
285 inline T &e(const int i, const int j) const_function {
286 // return elem[i][j];
287 return c[i * m + j];
288 }
289 // declare single e here too in case we have a vector
290 // (n || m == 1)
291 /**
292 * @brief Standard array indexing operation for vectors
293 * @details Accessing singular elements is insufficient, but Vector elements are often quite
294 * small.
295 * \code {.cpp}
296 * Vector<n,MyType> V;
297 * MyType a = V.e(i) \\ i <= n
298 * \endcode
299 * @tparam q Number of rows
300 * @tparam p Number of columns
301 * @param i Index
302 * @return T
303 */
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 {
306 return c[i];
307 }
308
309 /// @internal const_function implementation. See const_function for details
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 {
312 return c[i];
313 }
314
315 /**
316 * @brief Indexing operation [] defined only for vectors.
317 *
318 * @details Example:
319 *
320 * \code {.cpp}
321 * Vector<n,MyType> V;
322 * MyType a = V[i] \\ i <= n
323 * \endcode
324 *
325 * @tparam q row size n
326 * @tparam p column size m
327 * @param i row or vector index depending on which is being indexed
328 * @return T
329 */
330
331 template <int q = n, int p = m, std::enable_if_t<(q == 1 || p == 1), int> = 0>
332 inline T operator[](const int i) const {
333 return c[i];
334 }
335 /// @internal const_function implementation. See const_function for details
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 {
338 return c[i];
339 }
340
341 /**
342 * @brief Return reference to row in a matrix
343 * @details Since the Matrix data is ordered in a row major order accessing a row returns a
344 * reference to the row.
345 *
346 * \code{.cpp}
347 * Matrix<n,n,T> M;
348 * M.random();
349 * RowVector<n,T> V = M.row(i);
350 * \endcode
351 *
352 * @param r index of row to be referenced
353 * @return const RowVector<m, T>&
354 */
355 RowVector<m, T> row(int r) const {
357 for (int i = 0; i < m; i++)
358 v[i] = e(r, i);
359 return v;
360 }
361
362 /**
363 * @brief Set row of Matrix with #RowVector if types are assignable
364 * @details
365 * \code{.cpp}
366 * RowVector<n,T> V;
367 * V.random();
368 * Matrix<n,n,T> M;
369 * M.set_row(i,V);
370 * \endcode
371 * @tparam S RowVector type
372 * @param r Index of row to be set
373 * @param v RowVector to be set
374 */
375 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
376 void set_row(int r, const RowVector<m, S> &v) {
377 for (int i = 0; i < m; i++)
378 e(r, i) = v[i];
379 }
380
381 /**
382 * @brief Returns column vector as value at index c
383 * \code{.cpp}
384 * Matrix<n,n,T> M;
385 * M.random();
386 * Vector<n,T> V = M.column(i);
387 * \endcode
388 * @param c index of column vector to be returned
389 * @return const Vector<n, T>
390 */
391 Vector<n, T> column(int c) const {
392 Vector<n, T> v;
393 for (int i = 0; i < n; i++)
394 v[i] = e(i, c);
395 return v;
396 }
397
398
399 // matrix_col_t<n,T,Mtype> column(int c) {
400
401 // return matrix_col_t<n,T,Mtype>(*this,c);
402 // }
403
404
405 /// get column of a matrix
406 // hila_matrix_column_t<n, T, Mtype> column(int c) {
407 // return hila_matrix_column_t<n, T, Mtype>(*this, c);
408 // }
409
410 /**
411 * @brief Set column of Matrix with #Vector if types are assignable
412 * @details
413 * \code{.cpp}
414 * Vector<n,T> V;
415 * V.random();
416 * Matrix<n,n,T> M;
417 * M.set_column(i,V);
418 * \endcode
419 * @tparam S Vector type
420 * @param c Index of column to be set
421 * @param v #Vector to be set
422 */
423 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
424 void set_column(int c, const Vector<n, S> &v) {
425 for (int i = 0; i < n; i++)
426 e(i, c) = v[i];
427 }
428
429 /**
430 * @brief Return diagonal of square matrix
431 * @details If called for non square matrix the program will throw an error.
432 *
433 * \code{.cpp}
434 * Matrix<n,n,T> M = 1;
435 * DiagonalMatrix<n,T> D = M.diagonal();
436 * \endcode
437 * @return Vector<n, T> returned vector.
438 */
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);
444 return res;
445 }
446
447 /**
448 * @brief Set diagonal of square matrix to #Vector which is passed to the method
449 * @details If called for non square matrix the program will throw an error.
450 *
451 * Example:
452 *
453 * \code {.cpp}
454 * SquareMatrix<n,MyType> S = 0; \\ Zero matrix
455 * Vector<n,MyType> V = 1; \\ Vector assigned to 1 at all elements
456 * S.set_diagonal(V); \\ Results in Identity matrix of size n
457 * \endcode
458 *
459 * @tparam S type vector to assign values to
460 * @param v Vector to assign to diagonal
461 */
462 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
463 void set_diagonal(const Vector<n, S> &v) {
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);
467 }
468
469
470 /**
471 * @brief Cast Matrix to Array
472 * @details used for array operations
473 *
474 * @return Array<n, m, T>&
475 */
476 const Array<n, m, T> &asArray() const {
477 return *reinterpret_cast<const Array<n, m, T> *>(this);
478 }
479 // Same as above but with const_function, see const_function for details
480 Array<n, m, T> &asArray() const_function {
481 return *reinterpret_cast<Array<n, m, T> *>(this);
482 }
483
484 /**
485 * @brief Cast Vector to DiagonalMatrix
486 *
487 * @return DiagonalMatrix<n,T>
488 */
489 #pragma hila loop_function
490 template <int mm = m, std::enable_if_t<mm == 1, int> = 0>
492 return *reinterpret_cast<const DiagonalMatrix<n, T> *>(this);
493 }
494 // Same as above but with const_function, see const_function for details
495 template <int mm = m, std::enable_if_t<mm == 1, int> = 0>
496 DiagonalMatrix<n, T> &asDiagonalMatrix() const_function {
497 return *reinterpret_cast<DiagonalMatrix<n, T> *>(this);
498 }
499
500
501 // casting from one Matrix (number) type to another: do not do this automatically.
502 // but require an explicit cast operator. This makes it easier to write code.
503 // or should it be automatic? keep/remove explicit?
504 // TODO: CHECK AVX CONVERSIONS
505
506 // template <typename S, typename Rtype,
507 // std::enable_if_t<
508 // Rtype::is_matrix() && Rtype::rows() == n && Rtype::columns() == m
509 // &&
510 // hila::is_assignable<typename (Rtype::argument_type) &,
511 // T>::value,
512 // int> = 0>
513 // explicit operator Rtype() const {
514 // Rtype res;
515 // for (int i = 0; i < n * m; i++)
516 // res.c[i] = c[i];
517 // return res;
518 // }
519
520 /**
521 * @brief Unary operator
522 * @details
523 *
524 * Negation operation
525 *
526 * \code {.cpp}
527 * M == -M;
528 * \endcode
529 *
530 */
531 inline Mtype operator-() const {
532 Mtype res;
533 for (int i = 0; i < n * m; i++) {
534 res.c[i] = -c[i];
535 }
536 return res;
537 }
538
539 /**
540 * @brief Addition operator
541 * @memberof Matrix_t
542 * @details
543 *
544 * Identity operation
545 *
546 * \code {.cpp}
547 * M == +M;
548 * \endcode
549 *
550 * @tparam Mtype1 Matrix type for a
551 * @tparam Mtype2 Matrix type for b
552 * @param a Left matrix
553 * @param b Right matrix
554 * @return Rtype
555 */
556 inline const auto &operator+() const {
557 return *this;
558 }
559
560 /**
561 * @brief Boolean operator == to determine if two matrices are exactly the same
562 * @details Tolerance for equivalence is zero, meaning that the matrices must be exactly the
563 * same.
564 *
565 * @tparam S Type for Matrix which is being compared to
566 * @param rhs right hand side Matrix which we are comparing
567 * @return true
568 * @return false
569 */
570 template <typename S, int n2, int m2>
571 bool operator==(const Matrix<n2, m2, S> &rhs) const {
572 if constexpr (n != n2 || m != m2)
573 return false;
574
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))
578 return false;
579 }
580 return true;
581 }
582
583 /**
584 * @brief Boolean operator != to check if matrices are exactly different
585 * @details if matrices are exactly the same then this will return false
586 *
587 * @tparam S Type for MAtrix which is being compared to
588 * @param rhs right hand side Matrix which we are comparing
589 * @return true
590 * @return false
591 */
592 template <typename S>
593 bool operator!=(const Matrix<n, m, S> &rhs) const {
594 return !(*this == rhs);
595 }
596
597 /**
598 * @brief Copy matrix assignment
599 * @details
600 * \code {.cpp}
601 * Matrix<n,m,MyType> M_0;
602 * .
603 * . M_0 has values assigned to it
604 * .
605 * Matrix<n,m,MyType> M; \\ undefined matrix
606 * M = M_0; \\ Assignment from M_0
607 * \endcode
608 * @param rhs Matrix to assign
609 * @return
610 */
611
612 template <typename S, typename MT,
613 std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
614 inline auto &operator=(const Matrix_t<n, m, S, MT> &rhs) out_only & {
615 for (int i = 0; i < n * m; i++) {
616 c[i] = rhs.c[i];
617 }
618 return *this;
619 }
620
621 // assign from 0
622 // #pragma hila loop_function
623 /**
624 * @brief Zero assignment
625 *
626 * \code {.cpp}
627 * Matrix<n,m,MyType> M;
628 * M = 0; Zero matrix;
629 * \endcode
630 *
631 * @param z 0
632 * @return Mtype&
633 */
634 inline Mtype &operator=(const std::nullptr_t &z) out_only {
635 for (int i = 0; i < n * m; i++) {
636 c[i] = 0;
637 }
638 return *this;
639 }
640
641 // Assign from "scalar" for square matrix
642 // #pragma hila loop_function
643 /**
644 * @brief Assignment from scalar
645 * @details Assigns the scalar to the diagonal elements as \f$ M = I\cdot a\f$
646 *
647 * \code {.cpp}
648 * MyType a = hila::random();
649 * Matrix<n,m,MyType> M;
650 * M = a; M = I*a
651 * \endcode
652
653 * @tparam S Scalar type to assign
654 * @param rhs Scalar to assign
655 * @return Mtype&
656 */
657 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value && n == m, int> = 0>
658 inline auto &operator=(const S rhs) out_only & {
659
660 for (int i = 0; i < n; i++)
661 for (int j = 0; j < n; j++) {
662 if (i == j)
663 e(i, j) = rhs;
664 else
665 e(i, j) = 0;
666 }
667 return *this;
668 }
669
670 // Assign from diagonal matrix
671 // #pragma hila loop_function
672
673 /**
674 * @brief Assignment from diagonal matrix
675 *
676 * @tparam S Element type of Diagonal matrix
677 * @param rhs Diagonal matrix to assign
678 * @return Mtype&
679 */
680 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
681 inline auto &operator=(const DiagonalMatrix<n, S> &rhs) out_only & {
682 static_assert(n == m,
683 "Assigning DiagonalMatrix to Matrix possible only for square matrices");
684
685 for (int i = 0; i < n; i++)
686 for (int j = 0; j < n; j++) {
687 if (i == j)
688 e(i, j) = rhs.e(i);
689 else
690 e(i, j) = 0;
691 }
692 return *this;
693 }
694
695 /**
696 * @brief Initializer list assignment
697 * @details
698 * \code{.cpp}
699 * Matrix<2,2,int> M ;
700 * M = {1, 0
701 * 0, 1};
702 * \endcode
703 * @tparam S Element type of initializer list
704 * @param rhs Initializer list to assign
705 * @return Mtype&
706 */
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");
710 int i = 0;
711 for (auto it = rhs.begin(); it != rhs.end(); it++, i++) {
712 c[i] = *it;
713 }
714 return *this;
715 }
716
717 // Delete the rvalue-assign op
718 template <typename S>
719 Matrix_t &operator=(const S &s) && = delete;
720
721
722 /**
723 * @brief Add assign operator Matrix to Matrix
724 * @details
725 *
726 * \code {.cpp}
727 * Matrix<n,m,MyType> M,N;
728 * M = 1;
729 * N = 1;
730 * M += N; \\M = 2*I
731 * \endcode
732 * @tparam S Element type of rhs
733 * @tparam MT Matrix type of rhs
734 * @param rhs Matrix to add
735 * @return Mtype&
736 */
737
738 template <typename S, typename MT,
739 std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
740 auto &operator+=(const Matrix_t<n, m, S, MT> &rhs) & {
741 for (int i = 0; i < n * m; i++) {
742 c[i] += rhs.c[i];
743 }
744 return *this;
745 }
746
747 /**
748 * @brief Subtract assign operator Matrix to MAtrix
749 * @details
750 *
751 * \code {.cpp}
752 * Matrix<n,m,MyType> M,N;
753 * M = 3;
754 * N = 1;
755 * M -= N; \\M = 2*I
756 * \endcode
757 *
758 * @param rhs Matrix to subtract with
759 * @return Mtype&
760 */
761
762 template <typename S, typename MT,
763 std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
764 auto &operator-=(const Matrix_t<n, m, S, MT> &rhs) & {
765 for (int i = 0; i < n * m; i++) {
766 c[i] -= rhs.c[i];
767 }
768 return *this;
769 }
770
771 // add assign a scalar to square matrix
772 // #pragma hila loop_function
773 /**
774 * @brief Add assign operator scalar to Matrix
775 * @details
776 * Adds scalar \f$ a \f$ to __square__ matrix as \f$ M + a\cdot\mathbb{1} \f$
777 *
778 * \code {.cpp}
779 * Matrix<n,m,MyType> M = 1;
780 * M += 1 ; \\ M = 2*I
781 * \endcode
782 *
783 * @tparam S Element type of scalar rhs
784 * @param rhs Sclar to add
785 * @return Mtype&
786 */
787 template <typename S,
788 std::enable_if_t<hila::is_assignable<T &, hila::type_plus<T, S>>::value, int> = 0>
789 auto &operator+=(const S &rhs) & {
790
791 static_assert(n == m, "rows != columns : scalar addition possible for square matrix only!");
792
793 for (int i = 0; i < n; i++) {
794 e(i, i) += rhs;
795 }
796 return *this;
797 }
798
799 // subtract assign type T and convertible
800 // #pragma hila loop_function
801 /**
802 * @brief Subtract assign operator scalar to Matrix
803 * @details Subtract scalar \f$ a \f$ to __square__ matrix as \f$ M - a\cdot\mathbb{1} \f$
804 *
805 * \code {.cpp}
806 * Matrix<n,m,MyType> M = 3;
807 * M -= 1 ; \\ M = 2*I
808 * \endcode
809 * @tparam S Element type of scalar rhs
810 * @param rhs scalar to add
811 * @return Mtype&
812 */
813 template <typename S,
814 std::enable_if_t<hila::is_assignable<T &, hila::type_minus<T, S>>::value, int> = 0>
815 auto &operator-=(const S rhs) & {
816 static_assert(n == m,
817 "rows != columns : scalar subtraction possible for square matrix only!");
818 for (int i = 0; i < n; i++) {
819 e(i, i) -= rhs;
820 }
821 return *this;
822 }
823
824 /**
825 * @brief Multiply assign scalar or matrix
826 * @details Multiplication works as defined for matrix multiplication
827 *
828 * Matrix multiply assign only defined for square matrices, since the matrix dimensions would
829 * change otherwise.
830 *
831 * \code {.cpp}
832 * Matrix<n,m,MyType> M,N;
833 * .
834 * . Fill matrices M and N
835 * .
836 * M *= N; \\ M = M*N
837 * \endcode
838 *
839 * @param rhs Matrix to multiply with
840 * @return template <int p, typename S, typename MT,
841 * std::enable_if_t<hila::is_assignable<T &, hila::type_mul<T, S>>::value, int> = 0>&
842 */
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>
845 auto &operator*=(const Matrix_t<m, p, S, MT> &rhs) & {
846 static_assert(m == p, "can't assign result of *= to lhs Matrix, because doing so "
847 "would change it's dimensions");
848 *this = *this * rhs;
849 return *this;
850 }
851
852
853 /*
854 // same type square matrices:
855 template <int p,typename S,typename MT,
856 std::enable_if_t<hila::is_assignable<T&,hila::type_mul<T,S>>::value, int> = 0>
857 Mtype& operator*=(const Matrix_t<m,p,S,MT>& rhs)& {
858 static_assert(m==p,"can't assign result of *= to lhs Matrix, because doing so "
859 "would change it's dimensions");
860
861 S tmp_row[m];
862 int i,j,k;
863 for(i=0; i<m; ++i) {
864 for(j=0; j<m; ++j) {
865 tmp_row[j]=e(i,j);
866 }
867 for(j=0; j<m; ++j) {
868 e(i,j)=tmp_row[0]*rhs.e(0,j);
869 for(k=1; k<m; ++k) {
870 e(i,j)+=tmp_row[k]*rhs.e(k,j);
871 }
872 }
873 }
874 return *this;
875 }
876 */
877
878 // multiply assign with scalar
879 // #pragma hila loop_function
880 /**
881 * @brief Multiply assign operator scalar to Matrix
882 * @details Multiply Matrix uniformly with scalar
883 *
884 * \code {.cpp}
885 * Matrix<n,m,MyType> M;
886 * .
887 * . Fill whole matrix with 1
888 * .
889 * M *= 2 ; \\ M is filled with 2
890 * \endcode
891 * @tparam S Scalar type of rhs
892 * @param rhs Scalar to multiply
893 * @return Mtype&
894 */
895 template <typename S,
896 std::enable_if_t<hila::is_assignable<T &, hila::type_mul<T, S>>::value, int> = 0>
897 auto &operator*=(const S rhs) & {
898 for (int i = 0; i < n * m; i++) {
899 c[i] *= rhs;
900 }
901 return *this;
902 }
903
904 /**
905 * @brief Divide assign oprator scalar with matrix
906 * @details Divide works as defined for scalar matrix division.
907 *
908 *
909 * \code {.cpp}
910 * Matrix<n,m,MyType> M;
911 * .
912 * . Fill whole matrix with 2
913 * .
914 * M /= 2 ; \\ M is filled with 1
915 * \endcode
916 * @tparam S Scalar type of rhs
917 * @param rhs Scalar to divide with
918 * @return
919 */
920
921 template <typename S,
922 std::enable_if_t<hila::is_assignable<T &, hila::type_div<T, S>>::value, int> = 0>
923 auto &operator/=(const S rhs) & {
924 for (int i = 0; i < n * m; i++) {
925 c[i] /= rhs;
926 }
927 return *this;
928 }
929
930 /**
931 * @brief Addition assign operator for DiagonalMatrix to Matrix
932 * @details Works as long as Matrix to assign to is square
933 * @tparam S Element type of rhs
934 * @param rhs DiagonalMatrix to add
935 * @return Mtype&
936 */
937 template <typename S,
938 std::enable_if_t<hila::is_assignable<T &, hila::type_plus<T, S>>::value, int> = 0>
939 auto &operator+=(const DiagonalMatrix<n, S> &rhs) & {
940 static_assert(n == m, "Assigning DiagonalMatrix possible only for square matrix");
941
942 for (int i = 0; i < n; i++)
943 e(i, i) += rhs.e(i);
944 return *this;
945 }
946
947 /**
948 * @brief Subtract assign operator for DiagonalMatrix to Matrix
949 * @details Works as long as Matrix to assign to is square
950 * @tparam S Element type of rhs
951 * @param rhs DiagonalMatrix to add
952 * @return Mtype&
953 */
954 template <typename S,
955 std::enable_if_t<hila::is_assignable<T &, hila::type_plus<T, S>>::value, int> = 0>
956 auto &operator-=(const DiagonalMatrix<n, S> &rhs) & {
957 static_assert(n == m, "Assigning DiagonalMatrix possible only for square matrix");
958
959 for (int i = 0; i < n; i++)
960 e(i, i) -= rhs.e(i);
961 return *this;
962 }
963
964 /**
965 * @brief Multiply assign operator for DiagonalMatrix to Matrix
966 * @details Simply defined as matrix multiplication, but since rhs is guaranteed to be diagonal
967 * the method is optimized to skip most of the elements.
968 *
969 * @tparam S Element type of rhs
970 * @param rhs DiagonalMatrix to multiply
971 * @return Mtype&
972 */
973 template <typename S,
974 std::enable_if_t<hila::is_assignable<T &, hila::type_mul<T, S>>::value, int> = 0>
975 auto &operator*=(const DiagonalMatrix<m, S> &rhs) & {
976
977 for (int i = 0; i < n; i++)
978 for (int j = 0; j < m; j++)
979 e(i, j) *= rhs.e(j);
980
981 return *this;
982 }
983
984 /**
985 * @brief Divide assign operator for DiagonalMatrix to Matrix
986 * @details Well defined since rhs is guaranteed to be Diagonal.
987 *
988 * Let M be the Matrix which we divide and D be the DiagonalMatrix which we divide with then the
989 * operation is defined as \f$ M \cdot D^{-1} \f$.
990 * @tparam S Element type of rhs
991 * @param rhs DiagonalMatrix to divide
992 * @return Mtype&
993 */
994 template <typename S,
995 std::enable_if_t<hila::is_assignable<T &, hila::type_div<T, S>>::value, int> = 0>
996 auto &operator/=(const DiagonalMatrix<m, S> &rhs) & {
997
998 for (int i = 0; i < n; i++)
999 for (int j = 0; j < m; j++)
1000 e(i, j) /= rhs.e(j);
1001
1002 return *this;
1003 }
1004
1005
1006 /**
1007 * @brief Matrix fill
1008 * @details Fills the matrix with element if it is assignable to matrix type T
1009 *
1010 * Works as follows:
1011 *
1012 * \code {.cpp}
1013 * Matrix<n,m,MyType> M;
1014 * M.fill(2) \\ Matrix is filled with 2
1015 * \endcode
1016 *
1017 * @tparam S Scalar type to of rhs
1018 * @param rhs Element to fill matrix with
1019 * @return const Mtype&
1020 */
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++)
1024 c[i] = rhs;
1025 return *this;
1026 }
1027
1028 /**
1029 * @brief Transpose of matrix
1030 * @details Return type for square matrix is same input, for non square return type is
1031 * Matrix<n,m,MyType>
1032 * @tparam mm
1033 * @tparam std::conditional<n, Mtype, Matrix<m, n, T>>::type
1034 * @return const Rtype
1035 */
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>
1039 inline Rtype transpose() const {
1040 Rtype res;
1041 for (int i = 0; i < n; i++)
1042 for (int j = 0; j < m; j++) {
1043 res.e(j, i) = e(i, j);
1044 }
1045 return res;
1046 }
1047
1048 /**
1049 * @brief Transpose of vector
1050 * @details Returns reference
1051 *
1052 * @tparam mm
1053 * @return const RowVector<n, T>&
1054 */
1055 template <int mm = m, std::enable_if_t<mm == 1, int> = 0>
1056 inline const RowVector<n, T> &transpose() const {
1057 return *reinterpret_cast<const RowVector<n, T> *>(this);
1058 }
1059
1060 /**
1061 * @brief Transpose of RowVector
1062 * @details Returns reference
1063 *
1064 * @tparam mm
1065 * @return const RowVector<n, T>&
1066 */
1067 template <int nn = n, std::enable_if_t<nn == 1 && m != 1, int> = 0>
1068 inline const Vector<m, T> &transpose() const {
1069 return *reinterpret_cast<const Vector<m, T> *>(this);
1070 }
1071
1072
1073 /**
1074 * @brief Returns complex conjugate of Matrix
1075 *
1076 * @return const Mtype
1077 */
1078 inline Mtype conj() const {
1079 Mtype res;
1080 for (int i = 0; i < n * m; i++) {
1081 res.c[i] = ::conj(c[i]);
1082 }
1083 return res;
1084 }
1085
1086 /**
1087 * @brief Hermitian conjugate of matrix
1088 * @details for square matrix return type is same, for non square it is Matrix<m,n,MyType>
1089 *
1090 * @tparam std::conditional<n, Mtype, Matrix<m, n, T>>::type
1091 * @return const Rtype
1092 */
1093 template <typename Rtype = typename std::conditional<n == m, Mtype, Matrix<m, n, T>>::type>
1094 inline Rtype dagger() const {
1095 Rtype res;
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));
1099 }
1100 return res;
1101 }
1102
1103 /**
1104 * @brief Adjoint of matrix
1105 * @details Alias to dagger
1106 *
1107 * @tparam std::conditional<n, Mtype, Matrix<m, n, T>>::type
1108 * @return Rtype
1109 */
1110 template <typename Rtype = typename std::conditional<n == m, Mtype, Matrix<m, n, T>>::type>
1111 inline Rtype adjoint() const {
1112 return dagger();
1113 }
1114
1115 /**
1116 * @brief Returns absolute value of Matrix
1117 * @details For Matrix<n,m,Complex<T>> case type is changed to Matrix<n,m,T> as expected since
1118 * absolute value of a complex number is real
1119 *
1120 * @tparam M
1121 * @return Mtype
1122 */
1123 auto abs() const {
1125 for (int i = 0; i < n * m; i++) {
1126 res.c[i] = ::abs(c[i]);
1127 }
1128 return res;
1129 }
1130
1131
1132 // It seems that using special "Dagger" type makes the code slower!
1133 // Disable it now
1134 // inline const DaggerMatrix<m,n,T> & dagger() const {
1135 // return *reinterpret_cast<const DaggerMatrix<m,n,T>*>(this);
1136 // }
1137
1138 /**
1139 * @brief Returns real part of Matrix or #Vector
1140 *
1141 * @return Matrix<n, m, hila::arithmetic_type<T>>
1142 */
1145 for (int i = 0; i < m * n; i++) {
1146 res.c[i] = ::real(c[i]);
1147 }
1148 return res;
1149 }
1150
1151 /**
1152 * @brief Returns imaginary part of Matrix or #Vector
1153 *
1154 * @return Matrix<n, m, hila::arithmetic_type<T>>
1155 */
1158 for (int i = 0; i < m * n; i++) {
1159 res.c[i] = ::imag(c[i]);
1160 }
1161 return res;
1162 }
1163
1164 /**
1165 * @brief Computes Trace for Matrix
1166 * @details Not define for non square matrices. Static assert will stop code from compiling if
1167 * executed for non square matrices.
1168 *
1169 * @return T
1170 */
1171 T trace() const {
1172 static_assert(n == m, "trace not defined for non square matrices");
1173 T result(0);
1174 for (int i = 0; i < n; i++) {
1175 result += e(i, i);
1176 }
1177 return result;
1178 }
1179
1180 /**
1181 * @brief Multiply with given matrix and compute trace of result
1182 * @details Slightly cheaper operation than \code (M*N).trace() \endcode
1183 *
1184 * @tparam p
1185 * @tparam q
1186 * @tparam S
1187 * @tparam MT
1188 * @param rm
1189 * @return hila::type_mul<T, S>
1190 */
1191 template <int p, int q, typename S, typename MT>
1192 hila::type_mul<T, S> mul_trace(const Matrix_t<p, q, S, MT> &rm) const {
1193
1194 static_assert(p == m && q == n, "mul_trace(): argument matrix size mismatch");
1195
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);
1200 }
1201 return res;
1202 }
1203
1204 /**
1205 * @brief Calculate square norm - sum of squared elements
1206 *
1207 * @return hila::arithmetic_type<T>
1208 */
1209 hila::arithmetic_type<T> squarenorm() const {
1210 hila::arithmetic_type<T> result(0);
1211 for (int i = 0; i < n * m; i++) {
1212 result += ::squarenorm(c[i]);
1213 }
1214 return result;
1215 }
1216
1217 /**
1218 * @brief Calculate vector norm - sqrt of squarenorm
1219 *
1220 * @tparam S
1221 * @return hila::arithmetic_type<T>
1222 */
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 {
1226 return sqrt(squarenorm());
1227 }
1228
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()));
1233 }
1234
1235 /**
1236 * @brief Find max of Matrix only for arithmetic types
1237 */
1238 template <typename S = T, std::enable_if_t<hila::is_arithmetic<S>::value, int> = 0>
1239 T max() const {
1240 T res = c[0];
1241 for (int i = 1; i < n * m; i++) {
1242 if (res < c[i])
1243 res = c[i];
1244 }
1245 return res;
1246 }
1247
1248 /**
1249 * @brief Find min of Matrix only for arithmetic types
1250 */
1251 template <typename S = T, std::enable_if_t<hila::is_arithmetic<S>::value, int> = 0>
1252 T min() const {
1253 T res = c[0];
1254 for (int i = 1; i < n * m; i++) {
1255 if (res > c[i])
1256 res = c[i];
1257 }
1258 return res;
1259 }
1260
1261 auto max_abs() const {
1262 hila::arithmetic_type<T> tres, res = 0;
1263 for (int i = 0; i < n * m; i++) {
1264 tres = ::abs(c[i]);
1265 if (tres > res) {
1266 res = tres;
1267 }
1268 }
1269 return res;
1270 }
1271
1272 auto min_abs() const {
1273 hila::arithmetic_type<T> tres, res = ::abs(c[0]);
1274 for (int i = 1; i < n * m; i++) {
1275 tres = ::abs(c[i]);
1276 if (tres < res) {
1277 res = tres;
1278 }
1279 }
1280 return res;
1281 }
1282
1283 // dot product - (*this).dagger() * rhs
1284 // could be done as well by writing the operation as above!
1285 /**
1286 * @brief Dot product
1287 * @details Only works between two #Vector objects
1288 *
1289 * \code {.cpp}
1290 * Vector<m,MyType> V,W;
1291 * .
1292 * .
1293 * .
1294 * V.dot(W); // valid operation
1295 * \endcode
1296 *
1297 * \code {.cpp}
1298 * RowVector<m,MyType> V,W;
1299 * .
1300 * .
1301 * .
1302 * V.dot(W); // not valid operation
1303 * \endcode
1304 *
1305 *
1306 * @tparam p Row length for rhs
1307 * @tparam q Column length for rhs
1308 * @tparam S Type for rhs
1309 * @tparam R Gives resulting type of lhs and rhs multiplication
1310 * @param rhs Vector to compute dot product with
1311 * @return R Value of dot product
1312 */
1313 template <int p, int q, typename S, typename R = hila::type_mul<T, S>>
1314 inline R dot(const Matrix<p, q, S> &rhs) const {
1315 static_assert(m == 1 && q == 1 && p == n,
1316 "dot() product only for vectors of the same length");
1317
1318 R r = 0;
1319 for (int i = 0; i < n; i++) {
1320 r += ::conj(c[i]) * rhs.e(i);
1321 }
1322 return r;
1323 }
1324
1325 // outer product - (*this) * rhs.dagger(), sizes (n,1) and (p,1)
1326 // gives n * p matrix
1327 // could be done as well by the above operation!
1328
1329 /**
1330 * @brief Outer product
1331 * @details Only works between two #Vector objects
1332 *
1333 * \code{.cpp}
1334 * Vector<n,MyType> V,W;
1335 * .
1336 * .
1337 * .
1338 * V.outer_product(W); \\ Valid operation
1339 * \endcode
1340 *
1341 * \code {.cpp}
1342 * RowVector<m,MyType> V,W;
1343 * .
1344 * .
1345 * .
1346 * V.outer_product(W); // not valid operation
1347 * \endcode
1348 *
1349 * @tparam p Row length for rhs
1350 * @tparam q Column length for rhs
1351 * @tparam S Element type for rhs
1352 * @tparam R Type between lhs and rhs multiplication
1353 * @param rhs Vector to compute outer product with
1354 * @return Matrix<n, p, R>
1355 */
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");
1359
1360 Matrix<n, p, R> res;
1361 for (int i = 0; i < n; i++) {
1362 for (int j = 0; j < p; j++) {
1363 res.e(i, j) = c[i] * ::conj(rhs.e(j));
1364 }
1365 }
1366 return res;
1367 }
1368
1369
1370 /// dot with matrix - matrix
1371 // template <int p, std::enable_if_t<(p > 1 || m > 1), int> = 0>
1372 // inline Matrix<m, p, T> dot(const Matrix<n, p, T> &rhs) const {
1373 // Matrix<m, p, T> res;
1374 // for (int i = 0; i < m; i++)
1375 // for (int j = 0; j < p; j++) {
1376 // res.e(i, j) = 0;
1377 // for (int k = 0; k < n; j++)
1378 // res.e(i, j) += ::conj(e(k, i)) * rhs.e(k, j);
1379 // }
1380 // }
1381
1382 /// Generate random elements
1383
1384 /**
1385 * @brief Fills Matrix with random elements from uniform distribution
1386 * @details
1387 * \code{.cpp}
1388 * Matrix<n,m,T> M;
1389 * M.random();
1390 * \endcode
1391 * @return Mtype&
1392 */
1393 Mtype &random() out_only {
1394
1395 static_assert(hila::is_floating_point<hila::arithmetic_type<T>>::value,
1396 "Matrix/Vector random() requires non-integral type elements");
1397
1398 for (int i = 0; i < n * m; i++) {
1399 hila::random(c[i]);
1400 }
1401 return *this;
1402 }
1403
1404 /**
1405 * @brief Fills Matrix with gaussian random elements from gaussian distribution
1406 * @details
1407 * \code {.cpp}
1408 * Matrix<n,m,T> M;
1409 * M.gaussian_random();
1410 * \endcode
1411 *
1412 * @param width
1413 * @return Mtype&
1414 */
1415 Mtype &gaussian_random(base_type width = 1.0) out_only {
1416
1417 static_assert(hila::is_floating_point<hila::arithmetic_type<T>>::value,
1418 "Matrix/Vector gaussian_random() requires non-integral type elements");
1419
1420 // for Complex numbers gaussian_random fills re and im efficiently
1421 if constexpr (hila::is_complex<T>::value) {
1422 for (int i = 0; i < n * m; i++) {
1423 c[i].gaussian_random(width);
1424 }
1425 } else {
1426 // now not complex matrix
1427 // if n*m even, max i in loop below is n*m-2.
1428 // if n*m odd, max i is n*m-3
1429 double gr;
1430 for (int i = 0; i < n * m - 1; i += 2) {
1431 c[i] = hila::gaussrand2(gr) * width;
1432 c[i + 1] = gr * width;
1433 }
1434 if constexpr ((n * m) % 2 > 0) {
1435 c[n * m - 1] = hila::gaussrand() * width;
1436 }
1437 }
1438 return *this;
1439 }
1440
1441 /**
1442 * @brief permute columns of Matrix
1443 * @details Reordering is done based off of permutation vector
1444 *
1445 * @param permutation Vector of integers to permute columns
1446 * @return Mtype
1447 */
1448 Mtype permute_columns(const Vector<m, int> &permutation) const {
1449 Mtype res;
1450 for (int i = 0; i < m; i++)
1451 res.set_column(i, this->column(permutation[i]));
1452 return res;
1453 }
1454
1455 /**
1456 * @brief permute rows of Matrix
1457 * @details Reordering is done based off of permutation vector
1458 *
1459 * @param permutation Vector of integers to permute rows
1460 * @return Mtype
1461 */
1462 Mtype permute_rows(const Vector<n, int> &permutation) const {
1463 Mtype res;
1464 for (int i = 0; i < n; i++)
1465 res.set_row(i, this->row(permutation[i]));
1466 return res;
1467 }
1468
1469
1470 /**
1471 * @brief Permute elements of vector
1472 * @details Reordering is done based off of permutation vector
1473 *
1474 * @param permutation Vector of integers to permute rows
1475 * @return Mtype
1476 */
1477 template <int N>
1478 Mtype permute(const Vector<N, int> &permutation) const {
1479 static_assert(
1480 n == 1 || m == 1,
1481 "permute() only for vectors, use permute_rows() or permute_columns() for matrices");
1482 static_assert(N == Mtype::size(), "Incorrect size of permutation vector");
1483
1484 Mtype res;
1485 for (int i = 0; i < N; i++) {
1486 res[i] = (*this)[permutation[i]];
1487 }
1488 return res;
1489 }
1490
1491 /**
1492 * @brief Sort method for #Vector which returns permutation order
1493 * @details Order can be past as argument as either ascending (default) or descending
1494 *
1495 * __Ascending__
1496 * @code {.cpp}
1497 * Vector<n,MyType> V;
1498 * Vector<n,int> perm;
1499 * V.random();
1500 * V.sort(perm);
1501 * V.permute(perm);
1502 * @endcode
1503 *
1504 * __Descending__
1505 * @code {.cpp}
1506 * Vector<n,MyType> V;
1507 * Vector<n,int> perm;
1508 * V.random();
1509 * V.sort(perm,hila::sort::descending);
1510 * V.permute(perm);
1511 * @endcode
1512 * @tparam N
1513 * @param permutation
1514 * @param order
1515 * @return Mtype
1516 */
1517#pragma hila novector
1518 template <int N>
1519 Mtype sort(Vector<N, int> &permutation, hila::sort order = hila::sort::ascending) const {
1520
1521 static_assert(n == 1 || m == 1, "Sorting possible only for vectors");
1522 static_assert(hila::is_arithmetic<T>::value,
1523 "Sorting possible only for arithmetic vector elements");
1524 static_assert(N == Mtype::size(), "Incorrect size of permutation vector");
1525
1526 for (int i = 0; i < N; i++)
1527 permutation[i] = i;
1528 if (hila::sort::unsorted == order) {
1529 return *this;
1530 }
1531
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]);
1536 }
1537 } else {
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]);
1541 }
1542 }
1543
1544 return this->permute(permutation);
1545 }
1546
1547 /**
1548 * @brief Sort method for #Vector
1549 * @details Order can be past as argument as either ascending (default) or descending
1550 *
1551 * __Ascending__
1552 * @code {.cpp}
1553 * Vector<n,MyType> V;
1554 * V.random();
1555 * V.sort(); // V is sorted in ascending order
1556 * @endcode
1557 *
1558 * __Descending__
1559 * @code {.cpp}
1560 * V.sort(hila::sort::descending);
1561 * @endcode
1562 *
1563 * @param order Order to sort in
1564 * @return Mtype
1565 */
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");
1569
1570 Vector<Mtype::size(), int> permutation;
1571 return sort(permutation, order);
1572 }
1573
1574
1575 /**
1576 * @brief Multiply \f$ n \times m \f$-matrix from the left by \f$ n \times m \f$ matrix defined
1577 * by \f$ 2 \times 2 \f$ sub matrix
1578 * @details The multiplication is defined as follows, let \f$M\f$ as the \f$ 2 \times 2 \f$
1579 * input matrix and \f$B\f$ be `(this)` matrix, being the matrix stored in the object this
1580 * method is called for. Let \f$A = I\f$ be a \f$ n \times m \f$ unit matrix. We then set the
1581 * values of A to be: \f{align}{ A_{p,p} = M_{0,0}, \hspace{5px} A_{p,q} = M_{0,1}, \hspace{5px}
1582 * A_{q,p} = M_{1,0}, \hspace{5px} A_{q,q} = M_{1,1}. \f}
1583 *
1584 * Then the resulting matrix will be:
1585 *
1586 * \f{align}{ B = A \cdot B \f}
1587 *
1588 * @tparam R Element type of M
1589 * @tparam Mt Matrix type of M
1590 * @param p First row and column
1591 * @param q Second row and column
1592 * @param M \f$ 2 \times 2\f$ Matrix to multiply with
1593 */
1594 template <typename R, typename Mt>
1595 void mult_by_2x2_left(int p, int q, const Matrix_t<2, 2, R, Mt> &M) {
1596 static_assert(hila::contains_complex<T>::value || !hila::contains_complex<R>::value,
1597 "Cannot multiply real matrix with complex 2x2 block matrix");
1598
1599 Vector<2, T> a;
1600 for (int i = 0; i < m; i++) {
1601 a.e(0) = this->e(p, i);
1602 a.e(1) = this->e(q, i);
1603
1604 a = M * a;
1605
1606 this->e(p, i) = a.e(0);
1607 this->e(q, i) = a.e(1);
1608 }
1609 }
1610
1611 /**
1612 * @brief Multiply \f$ n \times m \f$-matrix from the right by \f$ n \times m \f$ matrix
1613 * defined by \f$ 2 \times 2 \f$ sub matrix
1614 * @details See Matrix::mult_by_2x2_left, only difference being that the multiplication is from
1615 * the right.
1616 *
1617 * @tparam R Element type of M
1618 * @tparam Mt Matrix type of M
1619 * @param p First row and column
1620 * @param q Second row and column
1621 * @param M \f$ 2 \times 2\f$ Matrix to multiply with
1622 */
1623 template <typename R, typename Mt>
1624 void mult_by_2x2_right(int p, int q, const Matrix_t<2, 2, R, Mt> &M) {
1625 static_assert(hila::contains_complex<T>::value || !hila::contains_complex<R>::value,
1626 "Cannot multiply real matrix with complex 2x2 block matrix");
1627
1629 for (int i = 0; i < n; i++) {
1630 a.e(0) = this->e(i, p);
1631 a.e(1) = this->e(i, q);
1632
1633 a = a * M;
1634
1635 this->e(i, p) = a.e(0);
1636 this->e(i, q) = a.e(1);
1637 }
1638 }
1639
1640
1641 // Following methods are defined in matrix_linalg.h
1642 //
1643
1644 /**
1645 * @brief following calculate the determinant of a square matrix
1646 * det() is the generic interface, using laplace for small matrices and LU for large
1647 */
1648
1649#pragma hila novector
1650 T det_lu() const;
1651#pragma hila novector
1652 T det_laplace() const;
1653#pragma hila novector
1654 T det() const;
1655
1656 /**
1657 * @brief Calculate eigenvalues and -vectors of an hermitean (or symmetric) matrix.
1658 * @details Returns the number of Jacobi iterations, or -1 if did not converge.
1659 * Arguments will contain eigenvalues and eigenvectors in columns of the "eigenvectors" matrix.
1660 * Computation is done in double precision despite the input matrix types.
1661 * @param eigenvaluevec Vector of computed eigenvalue
1662 * @param eigenvectors Eigenvectors as columns of \f$ n \times n \f$ Matrix
1663 *
1664 */
1665
1666#pragma hila novector
1667 template <typename Et, typename Mt, typename MT>
1668 int eigen_hermitean(out_only DiagonalMatrix<n, Et> &eigenvalues,
1669 out_only Matrix_t<n, n, Mt, MT> &eigenvectors,
1670 enum hila::sort sorted = hila::sort::unsorted) const;
1671
1672 eigen_result<Mtype> eigen_hermitean(enum hila::sort sorted = hila::sort::unsorted) const;
1673
1674 // Temporary interface to eigenvalues, to be removed
1675 template <typename Et, typename Mt, typename MT>
1676 int eigen_jacobi(out_only Vector<n, Et> &eigenvaluevec,
1677 out_only Matrix_t<n, n, Mt, MT> &eigenvectors,
1678 enum hila::sort sorted = hila::sort::unsorted) const {
1679
1681 int i = eigen_hermitean(d, eigenvectors, sorted);
1682 eigenvaluevec = d.asArray().asVector();
1683 return i;
1684 }
1685
1686
1687#pragma hila novector
1688 template <typename Et, typename Mt, typename MT>
1689 int svd(out_only Matrix_t<n, n, Mt, MT> &_U, out_only DiagonalMatrix<n, Et> &_D,
1690 out_only Matrix_t<n, n, Mt, MT> &_V,
1691 enum hila::sort sorted = hila::sort::unsorted) const;
1692
1693 svd_result<Mtype> svd(enum hila::sort sorted = hila::sort::unsorted) const;
1694
1695
1696#pragma hila novector
1697 template <typename Et, typename Mt, typename MT>
1698 int svd_pivot(out_only Matrix_t<n, n, Mt, MT> &_U, out_only DiagonalMatrix<n, Et> &_D,
1699 out_only Matrix_t<n, n, Mt, MT> &_V,
1700 enum hila::sort sorted = hila::sort::unsorted) const;
1701
1702 svd_result<Mtype> svd_pivot(enum hila::sort sorted = hila::sort::unsorted) const;
1703
1704
1705 //////// matrix_linalg.h
1706
1707 /// Convert to string for printing
1708 ///
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++) {
1714 text << e(i, j);
1715 if (i < n - 1 || j < m - 1)
1716 text << separator;
1717 }
1718 }
1719 return text.str();
1720 }
1721};
1722
1723/**
1724 * @brief \f$ n \times m \f$ Matrix class which defines matrix operations.
1725 * @details The Matrix class is a derived class of Matrix_t, which is the general definition of a
1726 * matrix class. See Matrix_t details section for reasoning.
1727 *
1728 * All mathematical operations for Matrix are inherited from Matrix_t and are visible on this page.
1729 * __NOTE__: Some method documentation is not being displayed, for example the assignment operator
1730 * documentation is not being inherited. See Matrix_t::operator= for use.
1731 *
1732 * To see all methods of initializing a matrix see constructor method #Matrix::Matrix
1733 *
1734 * __NOTE__: In the documentation examples n,m are integers and MyType is a HILA [standard
1735 * type](@ref standard) or Complex.
1736 *
1737 * @tparam n Number of rows
1738 * @tparam m Number of columns
1739 * @tparam T Data type Matrix
1740 */
1741template <int n, int m, typename T>
1742class Matrix : public Matrix_t<n, m, T, Matrix<n, m, T>> {
1743
1744 public:
1745 /// std incantation for field types
1746 using base_type = hila::arithmetic_type<T>;
1747 using argument_type = T;
1748
1749 /**
1750 * @brief Default constructo
1751 * @details
1752 * __NOTE__: n,m are integers and MyType is a HILA [standard type](@ref standard) or Complex.
1753 *
1754 * __Default constructor__:
1755 *
1756 * Allocates undefined \f$ n\times m\f$ Array.
1757 *
1758 * \code{.cpp}
1759 * Matrix<n,m,MyType> M;
1760 * \endcode
1761 */
1762 Matrix() = default;
1763 ~Matrix() = default;
1764 /// Default copy constructor
1765 Matrix(const Matrix &v) = default;
1766
1767
1768 // constructor from scalar -- keep it explicit! Not good for auto use
1769 /**
1770 * @brief Scalar constructor
1771 * @details
1772 *
1773 * Construct with given scalar at diagonal elements \f$ M = \mathbf{I}\cdot x\f$. Matrix must be
1774 * square \f$n == m\f$
1775 *
1776 * \code{.cpp}
1777 * MyType x = hila::random();
1778 * Matrix<n,m,MyType> M = x;
1779 * \endcode
1780
1781 * @tparam S Type for scalar
1782 * @tparam nn Number of rows
1783 * @tparam mm Numebr of columns
1784 * @param rhs Scalar element to assign
1785 */
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>
1788 explicit Matrix(const S rhs) out_only {
1789 for (int i = 0; i < n; i++)
1790 for (int j = 0; j < n; j++) {
1791 if (i == j)
1792 this->e(i, j) = rhs;
1793 else
1794 this->e(i, j) = 0;
1795 }
1796 }
1797
1798 /**
1799 * @brief Copy constructor
1800 *
1801 * Construction from previously defined matrix if types are compatible. For example the code
1802 * below only works if assignment from MyOtherType to MyType is defined.
1803 *
1804 * \code{.cpp}
1805 * Matrix<n,m,MyOtherType> M_0;
1806 * //
1807 * // M_0 is filled with content
1808 * //
1809 * Matrix<n,m,MyType> M = M_0;
1810 * \endcode
1811
1812 * @tparam S Element type for copied matrix
1813 * @tparam MT Matrix type for copied matrix
1814 * @param rhs Matrix to copy
1815 */
1816 template <typename S, typename MT,
1817 std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
1818 Matrix(const Matrix_t<n, m, S, MT> &rhs) out_only {
1819 for (int i = 0; i < n * m; i++) {
1820 this->c[i] = rhs.c[i];
1821 }
1822 }
1823
1824 /**
1825 * @brief Zero constructor
1826 * @details Constructing from 0 sets the whole Matrix to zero
1827 * @internal
1828 * @param z 0
1829 */
1830 Matrix(const std::nullptr_t &z) out_only {
1831 for (int i = 0; i < n * m; i++) {
1832 this->c[i] = 0;
1833 }
1834 }
1835
1836 // Construct matrix automatically from right-size initializer list
1837 // This does not seem to be dangerous, so keep non-explicit
1838 /**
1839 * @brief Initializer list constructor
1840 *
1841 * Construction from c++ initializer list.
1842 *
1843 * \code{.cpp}
1844 * Matrix<2,2,int> M = {1, 0
1845 * 0, 1};
1846 * \endcode
1847 * @tparam S Element type for initializer list
1848 * @param rhs
1849 */
1850 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
1851 Matrix(std::initializer_list<S> rhs) {
1852 assert(rhs.size() == n * m &&
1853 "Matrix/Vector initializer list size must match variable size");
1854 int i = 0;
1855 for (auto it = rhs.begin(); it != rhs.end(); it++, i++) {
1856 this->c[i] = *it;
1857 }
1858 }
1859
1860
1861 // use the Base::Base -trick to inherit assignments
1862 using Matrix_t<n, m, T, Matrix<n, m, T>>::operator=;
1863 Matrix &operator=(const Matrix &v) out_only & = default;
1864};
1865
1866namespace hila {
1867
1868//////////////////////////////////////////////////////////////////////////
1869// Tool to get "right" result type for matrix (T1) + (T2) -op, where
1870// T1 and T2 are either complex or arithmetic matrices
1871// Use: hila::mat_x_mat_type<T1,T2>
1872// - T1 == T2 returns T1
1873// - If T1 or T2 contains complex, returns Matrix<rows,cols,Complex<typeof T1+T2>>
1874// - otherwise returns the "larger" type
1875
1876template <typename T1, typename T2>
1877struct matrix_op_type_s {
1878 using type = Matrix<T1::rows(), T1::columns(), hila::ntype_op<T1, T2>>;
1879};
1880
1881// if types are the same
1882template <typename T>
1883struct matrix_op_type_s<T, T> {
1884 using type = T;
1885};
1886
1887// useful format mat_x_mat_type<T1,T2>
1888template <typename T1, typename T2>
1889using mat_x_mat_type = typename matrix_op_type_s<T1, T2>::type;
1890
1891////////////////////////////////////////////////////////////////////////////////
1892// Matrix + scalar result type:
1893// hila::mat_scalar_type<Mt,S>
1894// - if result is convertible to Mt, return Mt
1895// - if Mt is not complex and S is, return
1896// Matrix<Complex<type_sum(scalar_type(Mt),scalar_type(S))>>
1897// - otherwise return Matrix<type_sum>
1898
1899template <typename Mt, typename S, typename Enable = void>
1900struct matrix_scalar_op_s {
1901 using type =
1902 Matrix<Mt::rows(), Mt::columns(),
1903 Complex<hila::type_plus<hila::arithmetic_type<Mt>, hila::arithmetic_type<S>>>>;
1904};
1905
1906template <typename Mt, typename S>
1907struct matrix_scalar_op_s<
1908 Mt, S,
1909 typename std::enable_if_t<std::is_convertible<hila::type_plus<hila::number_type<Mt>, S>,
1910 hila::number_type<Mt>>::value>> {
1911 // using type = Mt;
1912 using type = typename std::conditional<
1913 hila::is_floating_point<hila::arithmetic_type<Mt>>::value, Mt,
1914 Matrix<Mt::rows(), Mt::columns(),
1915 hila::type_plus<hila::arithmetic_type<Mt>, hila::arithmetic_type<S>>>>::type;
1916};
1917
1918template <typename Mt, typename S>
1919using mat_scalar_type = typename matrix_scalar_op_s<Mt, S>::type;
1920
1921
1922} // namespace hila
1923
1924//////////////////////////////////////////////////////////////////////////
1925
1926/**
1927 * @brief Return transposed Matrix or #Vector
1928 * @details Wrapper around Matrix::transpose
1929 *
1930 * Can be called as:
1931 * \code {.cpp}
1932 * Matrix<n,m,MyType> M,M_transposed;
1933 * M_transposed = transpose(M);
1934 * \endcode
1935 *
1936 *
1937 * @tparam Mtype Matrix type
1938 * @param arg Object to be transposed
1939 * @return auto Mtype transposed matrix
1940 */
1941template <typename Mtype, std::enable_if_t<Mtype::is_matrix(), int> = 0>
1942inline auto transpose(const Mtype &arg) {
1943 return arg.transpose();
1944}
1945
1946/**
1947 * @brief Return conjugate Matrix or #Vector
1948 * @details Wrapper around Matrix::conj
1949 *
1950 * Can be called as:
1951 * \code {.cpp}
1952 * Matrix<n,m,MyType> M,M_conjugate;
1953 * M_conjugate = conj(M);
1954 * \endcode
1955 *
1956 *
1957 * @tparam Mtype Matrix type
1958 * @param arg Object to be conjugated
1959 * @return auto Mtype conjugated matrix
1960 */
1961template <typename Mtype, std::enable_if_t<Mtype::is_matrix(), int> = 0>
1962inline auto conj(const Mtype &arg) {
1963 return arg.conj();
1964}
1965
1966/**
1967 * @brief Return adjoint Matrix
1968 * @details Wrapper around Matrix::adjoint
1969 *
1970 * Can be called as:
1971 * \code {.cpp}
1972 * Matrix<n,m,MyType> M,M_adjoint;
1973 * M_conjugate = adjoint(M);
1974 * \endcode
1975 *
1976 *
1977 * @tparam Mtype Matrix type
1978 * @param arg Object to compute adjoint of
1979 * @return auto Mtype adjoint matrix
1980 */
1981template <typename Mtype, std::enable_if_t<Mtype::is_matrix(), int> = 0>
1982inline auto adjoint(const Mtype &arg) {
1983 return arg.adjoint();
1984}
1985
1986/**
1987 * @brief Return dagger of Matrix
1988 * @details Wrapper around Matrix::adjoint
1989 *
1990 * Same as adjoint
1991 *
1992 *
1993 * @tparam Mtype Matrix type
1994 * @param arg Object to compute dagger of
1995 * @return auto Mtype dagger matrix
1996 */
1997template <typename Mtype, std::enable_if_t<Mtype::is_matrix(), int> = 0>
1998inline auto dagger(const Mtype &arg) {
1999 return arg.adjoint();
2000}
2001
2002/**
2003 * @brief Return absolute value Matrix or #Vector
2004 * @details Wrapper around Matrix::abs
2005 *
2006 * Can be called as:
2007 * \code {.cpp}
2008 * Matrix<n,m,MyType> M,M_abs;
2009 * M_abs = abs(M);
2010 * \endcode
2011 *
2012 *
2013 * @tparam Mtype Matrix type
2014 * @param arg Object to compute absolute value of
2015 * @return auto Mtype absolute value Matrix or Vector
2016 */
2017template <typename Mtype, std::enable_if_t<Mtype::is_matrix(), int> = 0>
2018inline auto abs(const Mtype &arg) {
2019 return arg.abs();
2020}
2021
2022/**
2023 * @brief Return trace of square Matrix
2024 * @details Wrapper around Matrix::trace
2025 *
2026 * Can be called as:
2027 * \code {.cpp}
2028 * Matrix<n,m,MyType> M;
2029 * MyType T;
2030 * T = trace(M);
2031 * \endcode
2032 *
2033 *
2034 * @tparam Mtype Matrix type
2035 * @param arg Object to compute trace of
2036 * @return auto Trace of Matrix
2037 */
2038template <typename Mtype, std::enable_if_t<Mtype::is_matrix(), int> = 0>
2039inline auto trace(const Mtype &arg) {
2040 return arg.trace();
2041}
2042
2043/**
2044 * @brief Return real of Matrix or #Vector
2045 * @details Wrapper around Matrix::real
2046 *
2047 * Can be called as:
2048 * \code {.cpp}
2049 * Matrix<n,m,MyType> M,M_real;
2050 * M_real = real(M);
2051 * \endcode
2052 *
2053 *
2054 * @tparam Mtype Matrix type
2055 * @param arg Object to compute real part of
2056 * @return auto Mtype real part of arg
2057 */
2058template <typename Mtype, std::enable_if_t<Mtype::is_matrix(), int> = 0>
2059inline auto real(const Mtype &arg) {
2060 return arg.real();
2061}
2062
2063/**
2064 * @brief Return imaginary of Matrix or #Vector
2065 * @details Wrapper around Matrix::imag
2066 *
2067 * Can be called as:
2068 * \code {.cpp}
2069 * Matrix<n,m,MyType> M,M_imag;
2070 * M_imag = imag(M);
2071 * \endcode
2072 *
2073 *
2074 * @tparam Mtype Matrix type
2075 * @param arg Object to compute imag part of
2076 * @return auto Mtype imag part of arg
2077 */
2078template <typename Mtype, std::enable_if_t<Mtype::is_matrix(), int> = 0>
2079inline auto imag(const Mtype &arg) {
2080 return arg.imag();
2081}
2082
2083
2084///////////////////////////////////////////////////////////////////////////
2085// Now matrix additions: matrix + matrix
2086
2087
2088/**
2089 * @brief Addition operator Matrix + Matrix
2090 * @memberof Matrix_t
2091 * @details Addition operator between matrices is defined in the usual way (element wise).
2092 *
2093 * __NOTE__: Matrices must share same dimensionality.
2094 *
2095 * \code {.cpp}
2096 * Matrix<n,m,MyType> M, N, S;
2097 * M.fill(1);
2098 * N.fill(1);
2099 * S = M + N; // Resulting matrix is uniformly 2
2100 * \endcode
2101 * @tparam Mtype1 Matrix type of a
2102 * @tparam Mtype2 Matrix type of b
2103 * @param a Matrix to add
2104 * @param b Matrix to add
2105 * @return Rtype Return matrix of compatible type between Mtype1 and Mtype2
2106 */
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>
2111inline Rtype operator+(const Mtype1 &a, const Mtype2 &b) {
2112
2113 constexpr int n = Mtype1::rows();
2114 constexpr int m = Mtype1::columns();
2115
2116 static_assert(n == Mtype2::rows() && m == Mtype2::columns(), "Matrix sizes do not match");
2117
2118 Rtype r;
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);
2122 return r;
2123}
2124
2125/**
2126 * @brief Real micro-optimization Matrix + Matrix - no extra creation of variable and copy.
2127 * @internal
2128 * @tparam Mtype1
2129 * @tparam Mtype2
2130 * @tparam Mtype2,
2131 * std::enable_if_t<Mtype1::is_matrix() && Mtype2::is_matrix(), int,
2132 * typename Rtype, Mtype2>
2133 * @param a
2134 * @param b
2135 * @return Rtype Return matrix of compatible type between Mtype1 and Mtype2
2136 */
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>
2141inline Rtype operator+(Mtype1 a, const Mtype2 &b) {
2142
2143 constexpr int n = Mtype1::rows();
2144 constexpr int m = Mtype1::columns();
2145
2146 static_assert(n == Mtype2::rows() && m == Mtype2::columns(), "Matrix sizes do not match");
2147
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);
2151 return a;
2152}
2153
2154/**
2155 * @brief Subtraction operator Matrix - Matrix
2156 * @memberof Matrix_t
2157 * @details
2158 * Subtraction operator between matrices is defined in the usual way (element wise).
2159 *
2160 * __NOTE__: Matrices must share same dimensionality.
2161 *
2162 * \code {.cpp}
2163 * Matrix<n,m,MyType> M, N, S;
2164 * M.fill(1);
2165 * N.fill(1);
2166 * S = M - N; // Resulting matrix is uniformly 0
2167 * \endcode
2168 *
2169 * @tparam Mtype1 Matrix type of a
2170 * @tparam Mtype2 Matrix type of b
2171 * @param a Matrix to subtract from
2172 * @param b Matrix to subtract
2173 * @return Rtype Return matrix of compatible type between Mtype1 and Mtype2
2174 */
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>>
2178inline Rtype operator-(const Mtype1 &a, const Mtype2 &b) {
2179
2180 constexpr int n = Mtype1::rows();
2181 constexpr int m = Mtype1::columns();
2182
2183 static_assert(n == Mtype2::rows() && m == Mtype2::columns(), "Matrix sizes do not match");
2184
2185 Rtype r;
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);
2189 return r;
2190}
2191
2192/**
2193 * @brief Addition operator Matrix + scalar
2194 * @memberof Matrix_t
2195 * @details
2196 * Addition operator between matrix and scalar is defined in the usual way, where the scalar is
2197 * added to the diagonal elements.
2198 * *
2199 * \f$ M + 2 = M + 2\cdot\mathbb{1} \f$
2200 *
2201 * \code {.cpp}
2202 * Matrix<n,m,MyType> M,S;
2203 * M.fill(0);
2204 * S = M + 1; // Resulting matrix is identity matrix
2205 * \endcode
2206 *
2207 * @tparam Mtype Matrix type of b
2208 * @tparam S Type of scalar
2209 * @param a Matrix to add to
2210 * @param b Scalar to add
2211 * @return Rtype
2212 */
2213template <typename Mtype, typename S,
2214 std::enable_if_t<Mtype::is_matrix() && hila::is_complex_or_arithmetic<S>::value, int> = 0,
2215 typename Rtype = hila::mat_scalar_type<Mtype, S>>
2216inline Rtype operator+(const Mtype &a, const S &b) {
2217
2218 static_assert(Mtype::rows() == Mtype::columns(),
2219 "Matrix + scalar possible only for square matrix");
2220
2221 Rtype r;
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);
2225 if (i == j)
2226 r.e(i, j) += b;
2227 }
2228 return r;
2229}
2230
2231/**
2232 * @brief Addition operator scalar + Matrix
2233 * @memberof Matrix_t
2234 * @details
2235 * Addition operator between Scalar and Matrix is defined in the usual way, where the scalar is
2236 * treated as diagonal matrix which is then added to.
2237 * \f$ 2 + M = 2\cdot\mathbb{1} + M \f$
2238 *
2239 * \code {.cpp}
2240 * Matrix<n,m,MyType> M,R;
2241 * M = 0; // M = 0*I
2242 * R = 1 + M; // Resulting matrix is identity matrix
2243 * \endcode
2244 * @tparam Mtype Matrix type of a
2245 * @tparam S scalar type of b
2246 * @param b Matrix to add
2247 * @param a Scalar to add to
2248 * @return Rtype
2249 */
2250template <typename Mtype, typename S,
2251 std::enable_if_t<Mtype::is_matrix() && hila::is_complex_or_arithmetic<S>::value, int> = 0,
2252 typename Rtype = hila::mat_scalar_type<Mtype, S>>
2253inline Rtype operator+(const S &b, const Mtype &a) {
2254
2255 static_assert(Mtype::rows() == Mtype::columns(),
2256 "Matrix + scalar possible only for square matrix");
2257 return a + b;
2258}
2259
2260// matrix - scalar
2261/**
2262 * @brief Subtraction operator Matrix - scalar
2263 * @memberof Matrix_t
2264 * @details
2265 *
2266 * Subtraction operator between matrix and scalar is defined in the usual way, where the scalar
2267 * is subtracted from the diagonal elements.
2268 *
2269 * \f$ M - 2 = M - 2\cdot\mathbb{1} \f$
2270 *
2271 * \code {.cpp}
2272 * Matrix<n,m,MyType> M,R;
2273 * M = 2; // M = 2*I
2274 * R = M - 1; // Resulting matrix is identity matrix
2275 * \endcode
2276 * @tparam Mtype Matrix type of a
2277 * @tparam S Scalar type of b
2278 * @param a Matrix to subtract from
2279 * @param b Scalar to subtract
2280 * @return Rtype
2281 */
2282template <typename Mtype, typename S,
2283 std::enable_if_t<Mtype::is_matrix() && hila::is_complex_or_arithmetic<S>::value, int> = 0,
2284 typename Rtype = hila::mat_scalar_type<Mtype, S>>
2285inline Rtype operator-(const Mtype &a, const S &b) {
2286
2287 static_assert(Mtype::rows() == Mtype::columns(),
2288 "Matrix - scalar possible only for square matrix");
2289
2290 Rtype r;
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);
2294 if (i == j)
2295 r.e(i, j) -= b;
2296 }
2297 return r;
2298}
2299
2300/**
2301 * @brief Subtraction operator Scalar - Matrix
2302 * @memberof Matrix_t
2303 * @details
2304 * Subtraction operator between Scalar and Matrix is defined in the usual way, where the scalar is
2305 * treated as diagonal matrix which is then subtracted from.
2306 *
2307 * \f$ 2 - M = 2\cdot\mathbb{1} - M \f$
2308 *
2309 * \code {.cpp}
2310 * Matrix<n,m,MyType> M,R;
2311 * M = 1; // M = 1*I
2312 * R = 2 - M; // Resulting matrix is identity matrix
2313 * \endcode
2314 * @tparam Mtype Matrix type a
2315 * @tparam S Scalar type of b
2316 * @param b Scalar to subtract from
2317 * @param a Matrix to subtract
2318 * @return Rtype
2319 */
2320template <typename Mtype, typename S,
2321 std::enable_if_t<Mtype::is_matrix() && hila::is_complex_or_arithmetic<S>::value, int> = 0,
2322 typename Rtype = hila::mat_scalar_type<Mtype, S>>
2323inline Rtype operator-(const S &b, const Mtype &a) {
2324
2325 static_assert(Mtype::rows() == Mtype::columns(),
2326 "Scalar - matrix possible only for square matrix");
2327
2328 Rtype r;
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);
2332 if (i == j)
2333 r.e(i, j) += b;
2334 }
2335 return r;
2336}
2337
2338////////////////////////////////////////
2339// matrix * matrix is the most important bit
2340
2341/**
2342 * @brief Multiplication operator Square Matrix * Square Matrix
2343 * @internal
2344 * @memberof Matrix_t
2345 * @tparam Mt Matrix type
2346 * @param a First matrix to multiply
2347 * @param b Second matrix to multiply
2348 * @return Mt
2349 */
2350template <typename Mt, std::enable_if_t<Mt::is_matrix() && Mt::rows() == Mt::columns(), int> = 0>
2351inline Mt operator*(const Mt &a, const Mt &b) {
2352
2353 constexpr int n = Mt::rows();
2354
2355 Mt res;
2356
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);
2362 }
2363 }
2364 return res;
2365}
2366
2367/**
2368 * @brief Multiplication operator
2369 * @memberof Matrix_t
2370 * @details Multiplication type depends on the original types of the multiplied matrices. Defined
2371 * for the following operations.
2372 *
2373 * Matrix * Matrix
2374 *
2375 * [Standard matrix multiplication](https://en.wikipedia.org/wiki/Matrix_multiplication) where LHS
2376 * columns must match RHS rows
2377 *
2378 * \code {.cpp}
2379 * Matrix<2, 3, double> M;
2380 * Matrix<3, 2, double> N;
2381 * M.random();
2382 * N.random();
2383 * auto S = N * M; // Results in a 3 x 3 Matrix since N.rows() = 3 and M.columns = 3
2384 * \endcode
2385 *
2386 * #Vector * RowVector is same as outer product which is equivalent to a matrix
2387 * multiplication
2388 *
2389 * \code {.cpp}
2390 * auto S = W*V // Result in n x n Matrix of type MyType
2391 * \endcode
2392 *
2393 * @tparam Mt1 Matrix type for a
2394 * @tparam Mt2 Matrix type for b
2395 * @tparam n Number of rows
2396 * @tparam m Number of columns
2397 * @param a Left Matrix or Vector
2398 * @param b Right Matrix or RowVector
2399 * @return Matrix<n, m, R>
2400 */
2401template <typename Mt1, typename Mt2,
2402 std::enable_if_t<Mt1::is_matrix() && Mt2::is_matrix() && !std::is_same<Mt1, Mt2>::value,
2403 int> = 0,
2404 typename R = hila::type_mul<hila::number_type<Mt1>, hila::number_type<Mt2>>,
2405 int n = Mt1::rows(), int m = Mt2::columns()>
2406inline Matrix<n, m, R> operator*(const Mt1 &a, const Mt2 &b) {
2407
2408 constexpr int p = Mt1::columns();
2409 static_assert(p == Mt2::rows(), "Matrix size: LHS columns != RHS rows");
2410
2411 Matrix<n, m, R> res;
2412
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);
2419 }
2420 }
2421 } else if constexpr (m == 1) {
2422 // matrix x vector
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);
2427 }
2428 }
2429 } else if constexpr (n == 1) {
2430 // vector x matrix
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);
2435 }
2436 }
2437 }
2438
2439 return res;
2440}
2441
2442/**
2443 * @brief Multiplication operator RowVector * Vector
2444 * @memberof Matrix_t
2445 * @details
2446 * Defined as standard [dot product](https://en.wikipedia.org/wiki/Dot_product) between
2447 * vectors as long as vectors are of same length
2448 *
2449 * \code {.cpp}
2450 * RowVector<n,MyType> V;
2451 * Vector<n,MyType> W;
2452 * //
2453 * // Fill V and W
2454 * //
2455 * auto S = V*W; // Results in MyType scalar
2456 * \endcode
2457 *
2458 * @tparam m Dimensions of RowVector
2459 * @tparam n Dimension of Vector
2460 * @tparam T1 Element type of RowVector
2461 * @tparam T2 Element type of Vector
2462 * @tparam Rt Return type of T1 \f$ \cdot \f$ T2 product
2463 * @param A RowVector to multiply
2464 * @param B Vector to multiply
2465 * @return Rt
2466 */
2467template <int m, int n, typename T1, typename T2, typename Rt = hila::ntype_op<T1, T2>>
2469
2470 static_assert(m == n, "Vector lengths do not match");
2471
2472 Rt res;
2473 res = A.e(0) * B.e(0);
2474 for (int i = 1; i < m; i++) {
2475 res += A.e(i) * B.e(i);
2476 }
2477 return res;
2478}
2479
2480///////////////////////////////////////////////////////////////////////////
2481
2482/**
2483 * @brief Multiplication operator Matrix * scalar
2484 * @memberof Matrix_t
2485 * @details Multiplication operator between Matrix and Scalar are defined in the usual way (element
2486 * wise)
2487 *
2488 * \code {.cpp}
2489 * Matrix<n,n,MyType> M;
2490 * M.fill(1);
2491 * auto S = M*2; // Resulting Matrix is uniformly 2
2492 * \endcode
2493 *
2494 * @tparam Mt Matrix type of mat
2495 * @tparam s Scalar type rhs
2496 * @param mat Matrix to multiply
2497 * @param rhs Scalar to multiply
2498 * @return Rt
2499 */
2500template <typename Mt, typename S,
2501 std::enable_if_t<(Mt::is_matrix() && hila::is_complex_or_arithmetic<S>::value), int> = 0,
2502 typename Rt = hila::mat_scalar_type<Mt, S>>
2503inline Rt operator*(const Mt &mat, const S &rhs) {
2504 Rt res;
2505 for (int i = 0; i < Rt::rows() * Rt::columns(); i++) {
2506 res.c[i] = mat.c[i] * rhs;
2507 }
2508 return res;
2509}
2510
2511/**
2512 * @brief Multiplication operator Scalar * Matrix
2513 * @memberof Matrix_t
2514 * @details Multiplication operator between Matrix and Scalar are defined in the usual way (element
2515 * wise)
2516 *
2517 * \code {.cpp}
2518 * Matrix<n,n,MyType> M;
2519 * M.fill(1);
2520 * auto S = 2*M; // Resulting Matrix is uniformly 2
2521 * \endcode
2522 *
2523 * @tparam Mt Matrix type of mat
2524 * @tparam s Scalar type rhs
2525 * @param mat Matrix to multiply
2526 * @param rhs Scalar to multiply
2527 * @return Rt
2528 */
2529template <typename Mt, typename S,
2530 std::enable_if_t<(Mt::is_matrix() && hila::is_complex_or_arithmetic<S>::value), int> = 0,
2531 typename Rt = hila::mat_scalar_type<Mt, S>>
2532inline Rt operator*(const S &rhs, const Mt &mat) {
2533 return mat * rhs; // assume commutes
2534}
2535
2536// matrix / scalar
2537
2538/**
2539 * @brief Division operator
2540 * @memberof Matrix_t
2541 * @details Defined for following operations
2542 *
2543 * __Matrix / Scalar:__
2544 *
2545 * Division operator between Matrix and Scalar are defined in the usual way (element wise)
2546 *
2547 * \code {.cpp}
2548 * Matrix<n,m,MyType> M;
2549 * M.fill(2);
2550 * auto S = M/2; // Resulting matrix is uniformly 1
2551 * \endcode
2552 *
2553 * @tparam Mt Matrix type
2554 * @tparam S Scalar type
2555 * @param mat Matrix to divide scalar with
2556 * @param rhs Scalar to divide with
2557 * @return Rt Resulting Matrix
2558 */
2559template <typename Mt, typename S,
2560 std::enable_if_t<(Mt::is_matrix() && hila::is_complex_or_arithmetic<S>::value), int> = 0,
2561 typename Rt = hila::mat_scalar_type<Mt, S>>
2562inline Rt operator/(const Mt &mat, const S &rhs) {
2563 Rt res;
2564 for (int i = 0; i < Rt::rows() * Rt::columns(); i++) {
2565 res.c[i] = mat.c[i] / rhs;
2566 }
2567 return res;
2568}
2569
2570
2571/**
2572 * @brief Returns trace of product of two matrices
2573 * @details Wrapper around Matrix::mul_trace in the form
2574 * \code {.cpp}
2575 * mul_trace(a,b) // -> a.mul_trace(b)
2576 * \endcode
2577 *
2578 * @tparam Mtype1 Matrix type for a
2579 * @tparam Mtype2 Matrix type for b
2580 * @param a Left Matrix
2581 * @param b Right Matrix
2582 * @return auto Resulting trace
2583 */
2584template <typename Mtype1, typename Mtype2,
2585 std::enable_if_t<Mtype1::is_matrix() && Mtype2::is_matrix(), int> = 0>
2586inline auto mul_trace(const Mtype1 &a, const Mtype2 &b) {
2587
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);
2591}
2592
2593/**
2594 * @brief compute product of two matrices and write result to existing matrix
2595 *
2596 * @tparam Mt1, Mt2, Mt3 Matrix types
2597 * @param a Left Matrix of type Mt1
2598 * @param b Right Matrix of type Mt2
2599 * @param c Matrix of type Mt3 to which result gets written
2600 * @return void
2601 */
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();
2611 int i, j, k;
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);
2617 }
2618 }
2619 }
2620}
2621
2622/**
2623 * @brief compute hermitian conjugate of product of two matrices and write result to existing
2624 * matrix
2625 * @tparam Mt1, Mt2, Mt3 Matrix types
2626 * @param a Left Matrix of type Mt1
2627 * @param b Right Matrix of type Mt2
2628 * @param c Matrix of type Mt3 to which result gets written
2629 * @return void
2630 */
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();
2640 int i, j, k;
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));
2646 }
2647 }
2648 }
2649}
2650
2651/**
2652 * @brief compute product of a matrix and a scalar and write result to existing matrix
2653 *
2654 * @tparam Mt1 Matrix type1, S Scalar type, Mt2 Matrix type2
2655 * @param a Left Matrix of type Mt1
2656 * @param b Right Scalar of type S
2657 * @param c Matrix of type Mt2 to which result gets written
2658 * @return void
2659 */
2660template <typename Mt1, typename S, typename Mt2,
2661 std::enable_if_t<(Mt1::is_matrix() && Mt2::is_matrix() &&
2663 int> = 0>
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();
2669 int i, j, k;
2670 for (i = 0; i < n; ++i) {
2671 for (j = 0; j < m; ++j) {
2672 res.e(i, j) = a.e(i, j) * b;
2673 }
2674 }
2675}
2676
2677/**
2678 * @brief compute product of a scalar and a matrix and write result to existing matrix
2679 *
2680 * @tparam S Scalar type, Mt1 Matrix type1, Mt2 Matrix type2
2681 * @param a Left Scalar of type S
2682 * @param b Right Matrix of type Mt1
2683 * @param c Matrix of type Mt2 to which result gets written
2684 * @return void
2685 */
2686template <typename S, typename Mt1, typename Mt2,
2687 std::enable_if_t<(Mt1::is_matrix() && Mt2::is_matrix() &&
2689 int> = 0>
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();
2695 int i, j, k;
2696 for (i = 0; i < n; ++i) {
2697 for (j = 0; j < m; ++j) {
2698 res.e(i, j) = b.e(i, j) * a;
2699 }
2700 }
2701}
2702
2703/**
2704 * @brief compute product of two matrices and add result to existing matrix
2705 *
2706 * @tparam Mt1, Mt2, Mt3 Matrix types
2707 * @param a Left Matrix of type Mt1
2708 * @param b Right Matrix of type Mt2
2709 * @param c Matrix of type Mt3 to which result gets written
2710 * @return void
2711 */
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();
2721 int i, j, k;
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);
2726 }
2727 }
2728 }
2729}
2730
2731/**
2732 * @brief compute product of a matrix and a scalar and add result to existing matrix
2733 *
2734 * @tparam Mt1 Matrix type1, S Scalar type, Mt2 Matrix type2
2735 * @param a Left Matrix of type Mt1
2736 * @param b Right Scalar of type S
2737 * @param c Matrix of type Mt2 to which result gets written
2738 * @return void
2739 */
2740template <typename Mt1, typename S, typename Mt2,
2741 std::enable_if_t<(Mt1::is_matrix() && Mt2::is_matrix() &&
2743 int> = 0>
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();
2749 int i, j, k;
2750 for (i = 0; i < n; ++i) {
2751 for (j = 0; j < m; ++j) {
2752 res.e(i, j) += a.e(i, j) * b;
2753 }
2754 }
2755}
2756
2757/**
2758 * @brief compute product of a scalar and a matrix and add result to existing matrix
2759 *
2760 * @tparam S Scalar type, Mt1 Matrix type1, Mt2 Matrix type2
2761 * @param a Left Scalar of type S
2762 * @param b Right Matrix of type Mt1
2763 * @param c Matrix of type Mt2 to which result gets written
2764 * @return void
2765 */
2766template <typename S, typename Mt1, typename Mt2,
2767 std::enable_if_t<(Mt1::is_matrix() && Mt2::is_matrix() &&
2769 int> = 0>
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();
2775 int i, j, k;
2776 for (i = 0; i < n; ++i) {
2777 for (j = 0; j < m; ++j) {
2778 res.e(i, j) += b.e(i, j) * a;
2779 }
2780 }
2781}
2782
2783
2784/**
2785 * @brief compute product of two matrices and subtract result from existing matrix
2786 *
2787 * @tparam Mt1, Mt2, Mt3 Matrix types
2788 * @param a Left Matrix of type Mt1
2789 * @param b Right Matrix of type Mt2
2790 * @param c Matrix of type Mt3 to which result gets written
2791 * @return void
2792 */
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();
2802 int i, j, k;
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);
2807 }
2808 }
2809 }
2810}
2811
2812/**
2813 * @brief compute product of a matrix and a scalar and subtract result from existing matrix
2814 *
2815 * @tparam Mt1 Matrix type1, S Scalar type, Mt2 Matrix type2
2816 * @param a Left Matrix of type Mt1
2817 * @param b Right Scalar of type S
2818 * @param c Matrix of type Mt2 to which result gets written
2819 * @return void
2820 */
2821template <typename Mt1, typename S, typename Mt2,
2822 std::enable_if_t<(Mt1::is_matrix() && Mt2::is_matrix() &&
2824 int> = 0>
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();
2830 int i, j, k;
2831 for (i = 0; i < n; ++i) {
2832 for (j = 0; j < m; ++j) {
2833 res.e(i, j) -= a.e(i, j) * b;
2834 }
2835 }
2836}
2837
2838/**
2839 * @brief compute product of a scalar and a matrix and subtract result from existing matrix
2840 *
2841 * @tparam S Scalar type, Mt1 Matrix type1, Mt2 Matrix type2
2842 * @param a Left Scalar of type S
2843 * @param b Right Matrix of type Mt1
2844 * @param c Matrix of type Mt2 to which result gets written
2845 * @return void
2846 */
2847template <typename S, typename Mt1, typename Mt2,
2848 std::enable_if_t<(Mt1::is_matrix() && Mt2::is_matrix() &&
2850 int> = 0>
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();
2856 int i, j, k;
2857 for (i = 0; i < n; ++i) {
2858 for (j = 0; j < m; ++j) {
2859 res.e(i, j) -= b.e(i, j) * a;
2860 }
2861 }
2862}
2863
2864//////////////////////////////////////////////////////////////////////////////////
2865
2866
2867// Stream operator
2868/**
2869 * @brief Stream operator
2870 * @memberof Matrix
2871 * @details Naive Stream operator for formatting Matrix for printing. Simply puts elements one after
2872 * the other in row major order
2873 *
2874 * @tparam n
2875 * @tparam m
2876 * @tparam T
2877 * @tparam MT
2878 * @param strm
2879 * @param A
2880 * @return std::ostream&
2881 */
2882template <int n, int m, typename T, typename MT>
2883std::ostream &operator<<(std::ostream &strm, const Matrix_t<n, m, T, MT> &A) {
2884 for (int i = 0; i < n; i++)
2885 for (int j = 0; j < m; j++) {
2886 strm << A.e(i, j);
2887 if (i < n - 1 || j < m - 1)
2888 strm << ' ';
2889 }
2890 return strm;
2891}
2892
2893// Convert to string for "pretty" printing
2894//
2895
2896namespace hila {
2897
2898/**
2899 * @brief Change basic number type of Matrix/Vector
2900 *
2901 * hila::cast_to<double>(a);
2902 */
2903
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]);
2910 return res;
2911}
2912
2913template <typename Ntype, typename T, int n, int m,
2914 std::enable_if_t<hila::is_complex<T>::value, int> = 0>
2915Matrix<n, m, Complex<Ntype>> cast_to(const Matrix<n, m, T> &mat) {
2917 for (int i = 0; i < n * m; i++)
2918 res.c[i] = cast_to<Ntype>(mat.c[i]);
2919 return res;
2920}
2921
2922
2923
2924/**
2925 * @brief Converts Matrix_t object to string
2926 *
2927 * @tparam n Number of rows
2928 * @tparam m Number of columns
2929 * @tparam T Matrix element type
2930 * @tparam MT Matrix type
2931 * @param A Matrix to convert to string
2932 * @param prec Precision of T
2933 * @param separator Separator between elements
2934 * @return std::string
2935 */
2936template <int n, int m, typename T, typename MT>
2937std::string to_string(const Matrix_t<n, m, T, MT> &A, int prec = 8, char separator = ' ') {
2938 std::stringstream strm;
2939 strm.precision(prec);
2940
2941 for (int i = 0; i < n; i++)
2942 for (int j = 0; j < m; j++) {
2943 strm << hila::to_string(A.e(i, j), prec, separator);
2944 if (i < n - 1 || j < m - 1)
2945 strm << separator;
2946 }
2947 return strm.str();
2948}
2949
2950/**
2951 * @brief Formats Matrix_t object in a human readable way
2952 * @details Example 2 x 3 matrix is the following
2953 * \code {.cpp}
2954 * Matrix<2, 3, double> W;
2955 * W.random();
2956 * hila::out0 << hila::prettyprint(W ,4) << '\n';
2957 * // Result:
2958 * // [ 0.8555 0.006359 0.3193 ]
2959 * // [ 0.237 0.8871 0.8545 ]
2960 * \endcode
2961 *
2962 *
2963 * @tparam n Number of rows
2964 * @tparam m Number of columns
2965 * @tparam T Matrix element type
2966 * @tparam MT Matrix type
2967 * @param A Matrix to be formatted
2968 * @param prec Precision of Matrix element
2969 * @return std::string
2970 */
2971template <int n, int m, typename T, typename MT>
2972std::string prettyprint(const Matrix_t<n, m, T, MT> &A, int prec = 8) {
2973 std::stringstream strm;
2974 strm.precision(prec);
2975
2976 if constexpr (n == 1) {
2977 // print a vector, horizontally
2978 strm << '[';
2979 for (int i = 0; i < n * m; i++)
2980 strm << ' ' << prettyprint(A.e(i), prec);
2981 strm << " ]";
2982 // if (n > 1)
2983 // strm << "^T";
2984 } else {
2985 // now a matrix - split the matrix on lines.
2986 // do it so that columns are equal width...
2987 std::vector<std::string> lines, columns;
2988 lines.resize(n);
2989 columns.resize(n);
2990
2991 for (int i = 0; i < n; i++)
2992 lines[i] = "[ ";
2993 for (int j = 0; j < m; j++) {
2994 int size = 0;
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();
3001 }
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, ' ');
3006 }
3007 }
3008 for (int i = 0; i < n - 1; i++) {
3009 strm << lines[i] << "]\n";
3010 }
3011 strm << lines[n - 1] << "]";
3012 }
3013 return strm.str();
3014}
3015
3016} // namespace hila
3017
3018/**
3019 * @brief Returns square norm of Matrix
3020 * @details Wrapper around Matrix::squarenorm - sum of squared elements
3021 *
3022 * Can be called as:
3023 * \code {.cpp}
3024 * Matrix<n,m,MyType> M;
3025 * auto a = squarenorm(M);
3026 * \endcode
3027 *
3028 * @tparam Mt Matrix type
3029 * @param rhs Matrix to compute square norm of
3030 * @return auto
3031 */
3032template <typename Mt, std::enable_if_t<Mt::is_matrix(), int> = 0>
3033inline auto squarenorm(const Mt &rhs) {
3034 return rhs.squarenorm();
3035}
3036
3037// Vector norm - sqrt of squarenorm()
3038
3039/**
3040 * @brief Returns vector norm of Matrix
3041 * @details Wrapper around Matrix::norm - sqrt of sum of squared elements
3042 * @tparam Mt Matrix type
3043 * @param rhs Matrix to compute norm of
3044 * @return auto
3045 */
3046template <typename Mt, std::enable_if_t<Mt::is_matrix(), int> = 0>
3047inline auto norm(const Mt &rhs) {
3048 return rhs.norm();
3049}
3050
3051/**
3052 * @brief Calculate exp of a square matrix
3053 * @memberof Matrix
3054 * @details Computes up to certain order given as argument
3055 *
3056 * __Evaluation is done as__:
3057 * \f{align}{ \exp(H) &= 1 + H + \frac{H^2}{2!} + \frac{H^2}{2!} + \frac{H^3}{3!} \\
3058 * &= 1 + H\cdot(1 + (\frac{H}{2})\cdot (1 + (\frac{H}{3})\cdot(...))) \f}
3059 * Done backwards in order to reduce accumulation of errors
3060 * @tparam n Number of rowsMa
3061 * @tparam T Matrix element type
3062 * @tparam MT Matrix type
3063 * @param mat Matrix to compute exponential for
3064 * @param order order to compute exponential to
3065 * @return Matrix_t<n, m, T, MT>
3066 */
3067template <int n, int m, typename T, typename MT>
3068inline Matrix_t<n, m, T, MT> exp(const Matrix_t<n, m, T, MT> &mat, const int order = 20) {
3069 static_assert(n == m, "exp() only for square matrices");
3070 if(order>0) {
3071 hila::arithmetic_type<T> one = 1.0;
3072 Matrix_t<n, m, T, MT> r = mat;
3073
3074 // doing the loop step-by-step should reduce temporaries
3075 for (int k = order; k > 1; k--) {
3076 r *= (one / k);
3077 r += one;
3078 r *= mat;
3079 }
3080 r += one;
3081
3082 return r;
3083 } else {
3084 Matrix_t<n, m, T, MT> r(1.0);
3085 return r;
3086 }
3087}
3088
3089/**
3090 * @brief Calculate exp of a square matrix
3091 * @details Adds higher order terms till matrix norm converges within floating point accuracy and
3092 * returns required number of iterations
3093
3094 * @tparam n Number of rowsMa
3095 * @tparam T Matrix element type
3096 * @tparam MT Matrix type
3097 * @param mat Matrix to compute exponential for
3098 * @param niter (output) number of iteration performed till converges was reached
3099 * @return Matrix_t<n, m, T, MT>
3100 */
3101template <int n, int m, typename T, typename MT, typename atype = hila::arithmetic_type<T>>
3102inline Matrix_t<n, m, T, MT> altexp(const Matrix_t<n, m, T, MT> &mat, out_only int &niter) {
3103 static_assert(n == m, "exp() only for square matrices");
3104 atype one = 1.0;
3105 Matrix_t<n, m, T, MT> r = mat, mpow = mat;
3106 r += one;
3107 atype orsqnorm, rsqnorm = r.squarenorm();
3108 niter = 1;
3109 for (int k = 2; k < n * 15; ++k) {
3110 mpow *= mat;
3111 mpow *= (one / k);
3112 r += mpow;
3113 orsqnorm = rsqnorm;
3114 rsqnorm = r.squarenorm();
3115 if(rsqnorm == orsqnorm) {
3116 niter = k;
3117 break;
3118 }
3119 }
3120
3121 return r;
3122}
3123
3124/**
3125 * @brief Calculate exp of a square matrix
3126 * @details Adds hihger order terms till matrix norm converges within floating point accuracy
3127
3128 * @tparam n Number of rowsMa
3129 * @tparam T Matrix element type
3130 * @tparam MT Matrix type
3131 * @param mat Matrix to compute exponential for
3132 * @return Matrix_t<n, m, T, MT>
3133 */
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");
3137 atype one = 1.0;
3138 Matrix_t<n, m, T, MT> r = mat, mpow = mat;
3139 r += one;
3140 atype orsqnorm, rsqnorm = r.squarenorm();
3141 for (int k = 2; k < n * 15; ++k) {
3142 mpow *= mat;
3143 mpow *= (one / k);
3144 r += mpow;
3145 orsqnorm = rsqnorm;
3146 rsqnorm = r.squarenorm();
3147 if (rsqnorm == orsqnorm) {
3148 break;
3149 }
3150 }
3151 return r;
3152}
3153
3154/**
3155 * @brief Calculate mmat*exp(mat) and trace(mmat*dexp(mat))
3156 * @details exp and dexp computation is done using factorized, truncated
3157 * Taylor series method
3158
3159 * @tparam n Number of rowsMa
3160 * @tparam T Matrix element type
3161 * @tparam MT Matrix type
3162 * @param mat Matrix to compute exponential for
3163 * @param mmat Matrix to multiply with
3164 * @param r Matrix to which mmat*exp(mat) gets stored
3165 * @param dr matrix to which trace(mmat*dexp(mat)) gets stored
3166 * @param order = 20 integer input parameter giving the Taylor series truncation order
3167 * @return void
3168 */
3169template <int n, int m, typename T, typename MT>
3170inline void mult_exp(const Matrix_t<n, m, T, MT> &mat, const Matrix_t<n, m, T, MT> &mmat,
3171 out_only Matrix_t<n, m, T, MT> &r, out_only Matrix_t<n, m, T, MT> &dr,
3172 const int order = 20) {
3173 static_assert(n == m, "mult_exp() only for square matrices");
3174
3175 if (order > 0) {
3176 hila::arithmetic_type<T> one = 1.0;
3177 r = mat;
3178 dr = mmat;
3179 // doing the loop step-by-step should reduce temporaries
3180 for (int k = order; k > 1; k--) {
3181 r *= (one / k);
3182 r += one;
3183
3184 dr *= mat;
3185 dr *= (one / k);
3186 dr += r * mmat;
3187
3188 r *= mat;
3189 }
3190 r += one;
3191 r = mmat * r;
3192 } else {
3193 r = mmat;
3194 dr = 0;
3195 }
3196}
3197
3198
3199// Calculate exp of a square matrix
3200// using iterative Cayley-Hamilton described in arXiv:2404.07704
3201/**
3202 * @brief Calculate exp of a square matrix
3203 * @details Computation is done using iterative Cayley-Hamilton (cf. from arXiv:2404.07704)
3204
3205 * @tparam n Number of rowsMa
3206 * @tparam T Matrix element type
3207 * @tparam MT Matrix type
3208 * @param mat Matrix to compute exponential for
3209 * @param omat Matrix to which exponential of mat gets stored (optional)
3210 * @param pl array of n+1 temporary nxn Matrices (optional)
3211 * @return void (if omat is provided) or Matrix_t<n,m,T,MT>
3212 */
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>
3215inline int chexp(const Matrix_t<n, m, T, MT> &mat, out_only Matrix_t<n, m, T, MT> &omat,
3216 Mt(out_only &pl)[n]) {
3217 static_assert(n == m, "chexp() only for square matrices");
3218 // determine scaling factor:
3219 hila::arithmetic_type<T> sclim = sqrt(2.0), matnorm=norm(mat), sfac=1.0;
3220 int nb = 0;
3221 while (matnorm * sfac >= sclim) {
3222 sfac *= 0.5;
3223 ++nb;
3224 }
3225 // compute the first n matrix powers of mat and the corresponding traces :
3226 // the i-th matrix power of mat[][] is stored in pl[i][][]
3227 T trpl[n + 1]; // the trace of pl[i][][] is stored in trpl[i]
3228 trpl[0] = n;
3229 pl[1] = mat;
3230 pl[1] *= sfac;
3231 trpl[1] = trace(pl[1]);
3232 int i, j, k;
3233 for (i = 2; i < n; ++i) {
3234 j = i / 2;
3235 k = i % 2;
3236 mult(pl[j], pl[j + k], pl[i]);
3237 trpl[i] = trace(pl[i]);
3238 }
3239 j = n / 2;
3240 k = n % 2;
3241 trpl[n] = mul_trace(pl[j], pl[j + k]);
3242
3243 // compute the characteristic polynomial coefficients crpl[] from the traced powers trpl[] :
3244 T crpl[n + 1];
3245 crpl[n] = 1;
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];
3250 }
3251 crpl[n - j] /= j;
3252 }
3253
3254
3255 int mmax = 25 * n; // maximum number of Cayley-Hamilton iterations if no convergence is reached
3256 T al[n], pal[n]; // temp. Cayley-Hamilton coefficents
3257 hila::arithmetic_type<T>
3258 wpf = 1.0,
3259 twpf = 1.0; // initial values for power series coefficnet and its running sum
3260
3261 // set initial values for the n entries in al[] and pal[] :
3262 for (i = 0; i < n; ++i) {
3263 pal[i] = 0;
3264 al[i] = wpf;
3265 wpf /= (i + 1); // compute (i+1)-th power series coefficent from the i-th coefficient
3266 twpf += wpf;
3267 }
3268 pal[n - 1] = 1.0;
3269
3270 // next we iteratively add higher order power series terms to al[] till al[] stops changing
3271 // more precisely: the iteration will terminate as soon as twpf stops changing. Here twpf
3272 // is the sum \sum_{i=0}^{j} s_i/i!, with s_i referring to the magnitude the vector pal[]
3273 // would have at iteration i, if no renormalization were used.
3274 T cho; // temporary variables for iteration
3275 hila::arithmetic_type<T> ttwpf; // temporary variable for convergence check
3276 hila::arithmetic_type<T> s, rs = 1.0, rss; // temp variables used for renormalization of pal[]
3277 for (j = n; j < mmax; ++j) {
3278 s = 0;
3279 cho = pal[n - 1] * rs;
3280 for (i = n - 1; i > 0; --i) {
3281 pal[i] = pal[i - 1] * rs - cho * crpl[i];
3282 s += ::squarenorm(pal[i]);
3283 al[i] += wpf * pal[i];
3284 }
3285 pal[0] = -cho * crpl[0];
3286 s += ::squarenorm(pal[0]);
3287 al[0] += wpf * pal[0];
3288
3289 s = sqrt(s);
3290 if (s > 1.0) {
3291 // if s is bigger than 1, normalize pal[] by a factor rs=1.0/s in next itaration,
3292 // and multiply wpf by s to compensate
3293 wpf *= s / (j + 1);
3294 rs = 1.0 / s;
3295 rss = 1.0;
3296 } else {
3297 wpf /= (j + 1);
3298 rs = 1.0;
3299 rss = s;
3300 }
3301 ttwpf = twpf;
3302 twpf += wpf * rss;
3303 if (ttwpf == twpf) {
3304 // terminate iteration when numeric value of twpf stops changing
3305 break;
3306 }
3307 }
3308 // if(hila::myrank()==0) {
3309 // std::cout<<"chexp niter: "<<j<<" ("<<j-n<<")"<<std::endl;
3310 // }
3311
3312 // compute coefficients for the nb-times squared exponential of the scaled matrix:
3313 for (k = 0; k < nb; ++k) {
3314 cho = al[0];
3315 for(i=0; i<n; ++i) {
3316 trpl[i] = pal[i] = al[i];
3317 al[i] *= cho;
3318 }
3319 for (j = 1; j < n; ++j) {
3320 cho = pal[n - 1];
3321 for (i = n - 1; i > 0; --i) {
3322 pal[i] = pal[i - 1] - cho * crpl[i];
3323 al[i] += trpl[j] * pal[i];
3324 }
3325 pal[0] = -cho * crpl[0];
3326 al[0] += trpl[j] * pal[0];
3327 }
3328 }
3329
3330 // form output matrix:
3331 omat = al[0];
3332 for (k = 1; k < n;++k) {
3333 mult_add(al[k], pl[k], omat);
3334 }
3335
3336 return j;
3337}
3338
3339// overload wrapper for chexp where omat is not provided
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);
3346 return pl[0];
3347}
3348
3349// overload wrapper for chexp which creates the temporary matrix array pl[n+1] internally
3350template <int n, int m, typename T, typename MT>
3351inline int chexp(const Matrix_t<n, m, T, MT> &mat, out_only Matrix_t<n, m, T, MT> &omat) {
3352 static_assert(n == m, "chexp() only for square matrices");
3354 return chexp(mat, omat, pl);
3355}
3356
3357// overload wrapper for chexp where omat is not provided
3358// and which creates the temporary matrix array pl[n+1] internally
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);
3364 return pl[0];
3365}
3366
3367// overload wrapper for chexp where omat is not provided but niter
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>
3370inline Matrix_t<n, m, T, MT> chexp(const Matrix_t<n, m, T, MT> &mat, out_only int &niter, Mt(out_only &pl)[n]) {
3371 static_assert(n == m, "chexp() only for square matrices");
3372 niter = chexp(mat, pl[0], pl);
3373 return pl[0];
3374}
3375
3376// overload wrapper for chexp where omat is not provided but niter,
3377// and which creates the temporary matrix array pl[n+1] internally
3378template <int n, int m, typename T, typename MT>
3379inline Matrix_t<n, m, T, MT> chexp(const Matrix_t<n, m, T, MT> &mat, out_only int &niter) {
3380 static_assert(n == m, "chexp() only for square matrices");
3382 niter = chexp(mat, pl[0], pl);
3383 return pl[0];
3384}
3385
3386
3387// Calculate exp and dexp of a square matrix
3388// using iterative Cayley-Hamilton described in arXiv:2404.07704
3389/**
3390 * @brief Calculate exp and dexp of a square matrix
3391 * @details Computation is done using iterative Cayley-Hamilton (cf. from arXiv:2404.07704)
3392 * @tparam n Number of rowsMa
3393 * @tparam T Matrix element type
3394 * @tparam MT Matrix type
3395 * @param mat Matrix to compute exponential for
3396 * @param omat Matrix to which exponential of mat gets stored
3397 * @param domat matrix of Matrices to which the derivatives of the exponential
3398 * with respect to the components of mat gets stored
3399 * @return void
3400 */
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>
3403inline void chexp(const Matrix_t<n, m, T, MT> &mat, out_only Matrix_t<n, m, T, MT> &omat,
3404 Mt(out_only &domat)[n][m]) {
3405 static_assert(n == m, "chexp() only for square matrices");
3406 // determine scaling factor:
3407 hila::arithmetic_type<T> sclim = sqrt(2.0), matnorm = norm(mat), sfac = 1.0;
3408 int nb = 0;
3409 while (matnorm * sfac >= sclim) {
3410 sfac *= 0.5;
3411 ++nb;
3412 }
3413
3414 // compute the first n matrix powers of mat and the corresponding traces :
3415 Matrix_t<n, m, T, MT> pl[n]; // the i-th matrix power of mat[][] is stored in pl[i][][]
3416 T trpl[n + 1]; // the trace of pl[i][][] is stored in trpl[i]
3417 trpl[0] = (T)n;
3418 pl[1] = mat;
3419 pl[1] *= sfac;
3420 trpl[1] = trace(pl[1]);
3421 int i, j, k, l;
3422 for (i = 2; i < n; ++i) {
3423 j = i / 2;
3424 k = i % 2;
3425 mult(pl[j], pl[j + k], pl[i]);
3426 trpl[i] = trace(pl[i]);
3427 }
3428 j = n / 2;
3429 k = n % 2;
3430 trpl[n] = mul_trace(pl[j], pl[j + k]);
3431
3432 // compute the characteristic polynomial coefficients crpl[] from the traced powers trpl[] :
3433 T crpl[n + 1];
3434 crpl[n] = 1;
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];
3439 }
3440 crpl[n - j] /= j;
3441 }
3442
3443 const int mmax = 25 * n; // maximum number of iterations if no convergence is reached
3444 T al[n], pal[n]; // temp. Cayley-Hamilton coefficents
3445 hila::arithmetic_type<T>
3446 wpf = 1.0,
3447 twpf = 1.0; // initial values for power series coefficnet and its running sum
3448
3449 Matrix_t<n, m, T, MT> kmats = 0, kh = 0; // matrices required for derivative terms
3450
3451 // set initial values for the n entries in al[] and pal[], as well as for kmats and kh:
3452 // (note: since kmats and kh are symmetric, we operate only on their upper triangles)
3453 for(i = 0; i < n; ++i) {
3454 pal[i] = 0;
3455 al[i] = wpf;
3456 wpf /= (i + 1); // compute (i+1)-th power series coefficent from the i-th coefficient
3457 k = i / 2;
3458 for (j = i - k; j <= i; ++j) {
3459 kmats.e(i - j, j) = wpf;
3460 }
3461 twpf += wpf;
3462 }
3463 pal[n - 1] = 1.0;
3464 k = (n - 1) / 2;
3465 for (i = n - 1 - k; i < n; ++i) {
3466 kh.e(n - 1 - i, i) = 1;
3467 }
3468
3469 // next we iteratively add higher order power series terms to al[] till al[] stops changing
3470 // more precisely: the iteration will terminate as soon as twpf stops changing. Here twpf
3471 // is the sum \sum_{i=0}^{j} s_i/i!, with s_i referring to the magnitude the vector pal[]
3472 // would have at iteration i, if no renormalization were used.
3473 T cho; // temporary variables for iteration
3474 hila::arithmetic_type<T> ttwpf; // temporary variable for convergence check
3475 hila::arithmetic_type<T> s, rs = 1.0, rss; // temp variables used for renormalization of pal[]
3476 for (j = n; j < mmax; ++j) {
3477 s = 0;
3478 cho = pal[n - 1] * rs;
3479 for (i = n - 1; i > 0; --i) {
3480 pal[i] = pal[i - 1] * rs - cho * crpl[i];
3481 s += ::squarenorm(pal[i]);
3482 al[i] += wpf * pal[i];
3483 }
3484 pal[0] = -cho * crpl[0];
3485 s += ::squarenorm(pal[0]);
3486 al[0] += wpf * pal[0];
3487 s = sqrt(s);
3488 if (s > 1.0) {
3489 // if s is bigger than 1, normalize pal[] by a factor rs=1.0/s in next itaration,
3490 // and multiply wpf by s to compensate
3491 wpf *= s / (j + 1);
3492 rs = 1.0 / s;
3493 rss = 1.0;
3494 } else {
3495 wpf /= (j + 1);
3496 rs = 1.0;
3497 rss = s;
3498 }
3499 ttwpf = twpf;
3500 twpf += wpf * rss;
3501 if (ttwpf == twpf) {
3502 // terminate iteration when numeric value of twpf stops changing
3503 break;
3504 }
3505
3506 // add new terms to kmats and update kh :
3507 // (note: since kmats and kh are symmetric, we operate only on their upper triangles)
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);
3513 }
3514 if (i > 0) {
3515 kh.e(i, i) = kh.e(i - 1, i) * rs - cho * crpl[i];
3516 } else {
3517 kh.e(i, i) = pal[i] * rs - cho * crpl[i];
3518 }
3519 kmats.e(i, i) += wpf * kh.e(i, i);
3520 }
3521 }
3522
3523
3524 // compute coefficients for the nb-times squared exponential of the scaled matrix:
3525 for (k = 0; k < nb; ++k) {
3526 cho = al[0];
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); // define symmetric kh[][] to avoid case-distinctions
3534 // in computations below
3535 kmats.e(i, j) = kh.e(0, i) * al[j] + kh.e(0, j) * al[i];
3536 }
3537 al[i] *= cho;
3538 }
3539 for (l = 1; l < n; ++l) {
3540 cho = pal[n - 1];
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];
3546 }
3547 }
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];
3552 }
3553 }
3554 }
3555
3556
3557 // form output matrix omat = exp(mat) :
3558 omat = al[0];
3559 for (k = 1; k < n; ++k) {
3560 mult_add(al[k], pl[k], omat);
3561 }
3562
3563 // form output matrix of derivative matrices domat[ic1][ic2] = dexp(mat)/dmat^{ic2}_{ic1} :
3564 // ( i.e. domat[ic1][ic2].e(k, l) = dexp(mat)^k_l/dmat^{ic2}_{ic1}
3565 // = \sum_{i,j} kmats_{i,j} (mat^i)^{ic1}_l (mat^j)^k_{ic2} )
3566 // (note: in principle one could use symmetry domat[ic1][ic2].e(k, l) = domat[l][k].e(ic2, ic1),
3567 // which would amount to (i, j) <--> (j, i) with i = ic1 + n * ic2; j = l + n * k; )
3568
3569 // i=0: (treat i=0 case separately, since pl[0]=id is not used to avoid matrix-mult. by id)
3570 // j=0: (treat j=0 case separately, since pl[0]=id is not used)
3571 kh = kmats.e(0, 0);
3572 // j>0:
3573 for (j = 1; j < n; ++j) {
3574 mult_add(kmats.e(0, j), pl[j], kh);
3575 }
3576 int ic1, ic2;
3577 for (ic1 = 0; ic1 < n; ++ic1) {
3578 for (ic2 = 0; ic2 < n; ++ic2) {
3579 Mt &tdomat = domat[ic1][ic2];
3580 tdomat = 0;
3581 for (k = 0; k < n; ++k) {
3582 tdomat.e(k, ic1) += kh.e(k, ic2);
3583 }
3584 }
3585 }
3586 // i>0:
3587 for (i = 1; i < n; ++i) {
3588 // j=0: (treat j=0 case separately, since pl[0]=id is not used)
3589 kh = kmats.e(0, i);
3590 // j>0: (note: kmats is symmetric; only have upper triangle set)
3591 for (j = 1; j < i; ++j) {
3592 mult_add(kmats.e(j, i), pl[j], kh);
3593 }
3594 for (j = i; j < n; ++j) {
3595 mult_add(kmats.e(i, j), pl[j], kh);
3596 }
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);
3603 }
3604 }
3605 }
3606 }
3607 }
3608
3609}
3610
3611
3612/**
3613 * @brief Calculate exp(mat).dagger()*mmat*exp(mat) and trace(exp(mat).dagger*mmat*dexp(mat))
3614 * @details exp and dexp computation is done using iterative Cayley-Hamilton (arXiv:2404.07704)
3615 * Calculate omat[i][j] = (exp(mat).dagger() * mmat * exp(mat))[i][j]
3616 * and domat[i][j] = trace(exp(mat).dagger() * mmat * dexp(mat)/dmat[j][i]) for
3617 * given matrices mat and mmat, using iterative Cayley-Hamilton method (arXiv:2404.07704)
3618 * for computing exp(mat) and dexp(mat)/dmat[][]
3619 * @tparam n Number of rowsMa
3620 * @tparam T Matrix element type
3621 * @tparam MT Matrix type
3622 * @param mat Matrix to compute exponential for
3623 * @param mmat Matrix to multiply with
3624 * @param omat Matrix to which exp(mat).dagger()*mmat*exp(mat) gets stored
3625 * @param domat matrix to which trace(exp(mat).dagger()*mmat*dexp(mat)) gets stored
3626 * @return void
3627 */
3628template <int n, int m, typename T, typename MT>
3629inline void mult_chexp(const Matrix_t<n, m, T, MT> &mat, const Matrix_t<n, m, T, MT> &mmat,
3630 out_only Matrix_t<n, m, T, MT> &omat,
3631 out_only Matrix_t<n, m, T, MT> &domat) {
3632 static_assert(n == m, "mult_chexp() only for square matrices");
3633 // determine scaling factor:
3634 hila::arithmetic_type<T> sclim = sqrt(2.0), matnorm = norm(mat), sfac = 1.0;
3635 int nb = 0;
3636 while (matnorm * sfac >= sclim) {
3637 sfac *= 0.5;
3638 ++nb;
3639 }
3640
3641 // compute the first n matrix powers of mat and the corresponding traces :
3642 Matrix_t<n, m, T, MT> pl[n]; // the i-th matrix power of mat[][] is stored in pl[i][][]
3643 T trpl[n + 1]; // the trace of pl[i][][] is stored in trpl[i]
3644 Matrix_t<n, m, T, MT> &texp = pl[0]; // temp. storage for result of exponentiation
3645
3646 pl[1] = mat;
3647 pl[1] *= sfac;
3648 trpl[1] = trace(pl[1]);
3649 int i, j, k, l;
3650 for (i = 2; i < n; ++i) {
3651 j = i / 2;
3652 k = i % 2;
3653 mult(pl[j], pl[j + k], pl[i]);
3654 trpl[i] = trace(pl[i]);
3655 }
3656 // n-th power of mat
3657 j = n / 2;
3658 k = n % 2;
3659 trpl[n] = mul_trace(pl[j], pl[j + k]);
3660
3661 // compute the characteristic polynomial coefficients crpl[] from the traced powers trpl[] :
3662 T crpl[n + 1];
3663 crpl[n] = 1;
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];
3668 }
3669 crpl[n - j] /= j;
3670 }
3671
3672 const int mmax = 25 * n; // maximum number of iterations if no convergence is reached
3673 T al[n], pal[n]; // temp. Cayley-Hamilton coefficents
3674 hila::arithmetic_type<T>
3675 wpf = 1.0,
3676 twpf = 1.0; // initial values for power series coefficnet and its running sum
3677
3678 Matrix_t<n, m, T, MT> kmats = 0, kh = 0; // matrices required for derivative terms
3679
3680 // set initial values for the n entries in al[] and pal[], as well as for kmats and kh:
3681 // (note: since kmats and kh are symmetric, we operate only on their upper triangles)
3682 for (i = 0; i < n; ++i) {
3683 pal[i] = 0;
3684 al[i] = wpf;
3685 wpf /= (i + 1); // compute (i+1)-th power series coefficent from the i-th coefficient
3686 k = i / 2;
3687 for (j = i - k; j <= i; ++j) {
3688 kmats.e(i - j, j) = wpf;
3689 }
3690 twpf += wpf;
3691 }
3692 pal[n - 1] = 1.0;
3693 k = (n - 1) / 2;
3694 for (i = n - 1 - k; i < n; ++i) {
3695 kh.e(n - 1 - i, i) = 1.0;
3696 }
3697
3698 // next we iteratively add higher order power series terms to al[] till al[] stops changing
3699 // more precisely: the iteration will terminate as soon as twpf stops changing. Here twpf
3700 // is the sum \sum_{i=0}^{j} s_i/i!, with s_i referring to the magnitude the vector pal[]
3701 // would have at iteration i, if no renormalization were used.
3702 T cho; // temporary variable for iteration
3703 hila::arithmetic_type<T> ttwpf; // temporary variable for convergence check
3704 hila::arithmetic_type<T> s, rs = 1.0, rss; // temp variables used for renormalization of pal[]
3705 for (j = n; j < mmax; ++j) {
3706 s = 0;
3707 cho = pal[n - 1] * rs;
3708 for (i = n - 1; i > 0; --i) {
3709 pal[i] = pal[i - 1] * rs - cho * crpl[i];
3710 s += ::squarenorm(pal[i]);
3711 al[i] += wpf * pal[i];
3712 }
3713 pal[0] = -cho * crpl[0];
3714 s += ::squarenorm(pal[0]);
3715 al[0] += wpf * pal[0];
3716
3717 s = sqrt(s);
3718 if (s > 1.0) {
3719 // if s is bigger than 1, normalize pal[] by a factor rs=1.0/s in next itaration,
3720 // and multiply wpf by s to compensate
3721 wpf *= s / (j + 1);
3722 rs = 1.0 / s;
3723 rss = 1.0;
3724 } else {
3725 wpf /= (j + 1);
3726 rs = 1.0;
3727 rss = s;
3728 }
3729 ttwpf = twpf;
3730 twpf += wpf * rss;
3731 if (ttwpf == twpf) {
3732 // terminate iteration when numeric value of twpf stops changing
3733 break;
3734 }
3735
3736 // add new terms to kmats and update kh :
3737 // (note: since kmats and kh are symmetric, we operate only on their upper triangles)
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);
3743 }
3744 if(i>0) {
3745 kh.e(i, i) = kh.e(i - 1, i) * rs - cho * crpl[i];
3746 } else {
3747 kh.e(i, i) = pal[i] * rs - cho * crpl[i];
3748 }
3749 kmats.e(i, i) += wpf * kh.e(i, i);
3750 }
3751 }
3752
3753
3754 // compute coefficients for the nb-times squared exponential of the scaled matrix:
3755 for (k = 0; k < nb; ++k) {
3756 cho = al[0];
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); // define symmetric kh[][] to avoid case-distinctions
3764 // in computations below
3765 kmats.e(i, j) = kh.e(0, i) * al[j] + kh.e(0, j) * al[i];
3766 }
3767 al[i] *= cho;
3768 }
3769 for (l = 1; l < n; ++l) {
3770 cho = pal[n - 1];
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];
3776 }
3777 }
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];
3782 }
3783 }
3784 }
3785
3786
3787 // from matrix texp = exp(mat):
3788 texp = al[0];
3789 for (k = 1; k < n; ++k) {
3790 mult_add(al[k], pl[k], texp);
3791 }
3792
3793 Matrix_t<n, m, T, MT> tomat; // temp. storage for compuatation of derivative term
3794 // set tomat = exp(mat) * mmat
3795 // tomat = texp.dagger() * mmat;
3796 mult(texp.dagger(), mmat, tomat);
3797
3798 // computing domat[ic1][ic2] = tr(exp(mat).dagger() * mmat * dexp(mat)/dU^{ic2}_{ic1})
3799 // from kmats[][] and the matrix powers of mat[][]:
3800 // domat = \sum_{i=0}^{n-1} pl[i] * tomat * \sum_{j=0}^{n-1} kmats[i][j] * pl[j]
3801
3802 // i=0: (treat i=0 case separately, since pl[0]=id is not used to avoid matrix-mult. by id)
3803 // j=0: (treat j=0 case separately, since pl[0]=id is not used)
3804 kh = kmats.e(0, 0);
3805 // j>0:
3806 for (j = 1; j < n; ++j) {
3807 mult_add(kmats.e(0, j), pl[j], kh);
3808 }
3809
3810 mult(tomat, kh, domat);
3811 // i>0:
3812 for (i = 1; i < n; ++i) {
3813 // j=0: (treat j=0 case separately, since pl[0]=id is not used)
3814 kh = kmats.e(0, i);
3815 // j>0: (note: kmats is symmetric; only have upper triangle set)
3816 for (j = 1; j < i; ++j) {
3817 mult_add(kmats.e(j, i), pl[j], kh);
3818 }
3819 for (j = i; j < n; ++j) {
3820 mult_add(kmats.e(i, j), pl[j], kh);
3821 }
3822 mult(pl[i], tomat, omat);
3823 mult_add(omat, kh, domat);
3824 }
3825
3826 // setting omat = exp(mat).dagger() * mmat * exp(mat) = tomat * exp(mat) :
3827 // omat = tomat * texp;
3828 mult(tomat, texp, omat);
3829}
3830
3831
3832// Calculate exp and dexp of a square matrix
3833// using iterative Cayley-Hamilton described in arXiv:2404.07704
3834/**
3835 * @brief Calculate exp(mat) and the decomposition k_{i,j} of dexp in terms bilinears of powers
3836 * of mat
3837 * @details Computation is done using iterative Cayley-Hamilton (cf. from arXiv:2404.07704)
3838 * exp(mat)^a_b = \sum_{i=0}^{n-1} r_{i} (mat^i)^a_b
3839 * dexp(X)^a_b/dX^c_d|_X=mat = \sum_{i,j=0}^{n-1} k_{i,j} (mat^i)^d_b (mat^j)^a_c
3840 * @tparam n Number of rows of Matrix
3841 * @tparam m Number of columns of Matrix
3842 * @tparam T Matrix element type
3843 * @tparam MT Matrix type
3844 * @param mat Matrix of which to compute exponential
3845 * @param omat Matrix to which exponential of mat gets stored
3846 * @param kmats Matrix to which decomposition coefficients of dexp/dX|X=mat in terms
3847 * of bilinears of powers of mat get stroed
3848 * @return void
3849 */
3850template <int n, int m, typename T, typename MT>
3851inline void chexpk(const Matrix_t<n, m, T, MT> &mat, out_only Matrix_t<n, m, T, MT> &omat,
3852 out_only Matrix_t<n, m, T, MT> &kmats) {
3853 static_assert(n == m, "chexpk() only for square matrices");
3854 // determine scaling factor:
3855 hila::arithmetic_type<T> sclim = sqrt(2.0), matnorm = norm(mat), sfac = 1.0;
3856 int nb = 0;
3857 while (matnorm * sfac >= sclim) {
3858 sfac *= 0.5;
3859 ++nb;
3860 }
3861
3862 // compute the first n matrix powers of mat and the corresponding traces :
3863 Matrix_t<n, m, T, MT> pl[n]; // the i-th matrix power of mat[][] is stored in pl[i][][]
3864 T trpl[n + 1]; // the trace of pl[i][][] is stored in trpl[i]
3865 trpl[0] = (T)n;
3866 pl[1] = mat;
3867 pl[1] *= sfac;
3868 trpl[1] = trace(pl[1]);
3869 int i, j, k, l;
3870 for (i = 2; i < n; ++i) {
3871 j = i / 2;
3872 k = i % 2;
3873 mult(pl[j], pl[j + k], pl[i]);
3874 trpl[i] = trace(pl[i]);
3875 }
3876 j = n / 2;
3877 k = n % 2;
3878 trpl[n] = mul_trace(pl[j], pl[j + k]);
3879
3880 // compute the characteristic polynomial coefficients crpl[] from the traced powers trpl[] :
3881 T crpl[n + 1];
3882 crpl[n] = 1;
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];
3887 }
3888 crpl[n - j] /= j;
3889 }
3890
3891 const int mmax = 25 * n; // maximum number of iterations if no convergence is reached
3892 T al[n], pal[n]; // temp. Cayley-Hamilton coefficents
3893 hila::arithmetic_type<T>
3894 wpf = 1.0,
3895 twpf = 1.0; // initial values for power series coefficnet and its running sum
3896
3897 Matrix_t<n, m, T, MT> kh = 0; // temp. matrix required for derivative terms
3898 kmats = 0;
3899
3900 // set initial values for the n entries in al[] and pal[], as well as for kmats and kh:
3901 // (note: since kmats and kh are symmetric, we operate only on their upper triangles)
3902 for (i = 0; i < n; ++i) {
3903 pal[i] = 0;
3904 al[i] = wpf;
3905 wpf /= (i + 1); // compute (i+1)-th power series coefficent from the i-th coefficient
3906 k = i / 2;
3907 for (j = i - k; j <= i; ++j) {
3908 kmats.e(i - j, j) = wpf;
3909 }
3910 twpf += wpf;
3911 }
3912 pal[n - 1] = 1.0;
3913 k = (n - 1) / 2;
3914 for (i = n - 1 - k; i < n; ++i) {
3915 kh.e(n - 1 - i, i) = 1;
3916 }
3917
3918 // next we iteratively add higher order power series terms to al[] till al[] stops changing
3919 // more precisely: the iteration will terminate as soon as twpf stops changing. Here twpf
3920 // is the sum \sum_{i=0}^{j} s_i/i!, with s_i referring to the magnitude the vector pal[]
3921 // would have at iteration i, if no renormalization were used.
3922 T cho, tcrp; // temporary variables for iteration
3923 hila::arithmetic_type<T> ttwpf; // temporary variable for convergence check
3924 hila::arithmetic_type<T> s, rs = 1.0, rss; // temp variables used for renormalization of pal[]
3925 int jmax = mmax - 1;
3926 for (j = n; j < mmax; ++j) {
3927 s = 0;
3928 cho = pal[n - 1] * rs;
3929 for (i = n - 1; i > 0; --i) {
3930 pal[i] = pal[i - 1] * rs - cho * crpl[i];
3931 s += ::squarenorm(pal[i]);
3932 al[i] += wpf * pal[i];
3933 }
3934 pal[0] = -cho * crpl[0];
3935 s += ::squarenorm(pal[0]);
3936 al[0] += wpf * pal[0];
3937
3938 s = sqrt(s);
3939 if (s > 1.0) {
3940 // if s is bigger than 1, normalize pal[] by a factor rs=1.0/s in next itaration,
3941 // and multiply wpf by s to compensate
3942 wpf *= s / (j + 1);
3943 rs = 1.0 / s;
3944 rss = 1.0;
3945 } else {
3946 wpf /= (j + 1);
3947 rs = 1.0;
3948 rss = s;
3949 }
3950 ttwpf = twpf;
3951 twpf += wpf * rss;
3952 if (ttwpf == twpf) {
3953 // terminate iteration when numeric value of twpf stops changing
3954 jmax = j;
3955 break;
3956 }
3957
3958 // add new terms to kmats and update kh :
3959 // (note: since kmats and kh are symmetric, we operate only on their upper triangles)
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);
3965 }
3966 if (i > 0) {
3967 kh.e(i, i) = kh.e(i - 1, i) * rs - cho * crpl[i];
3968 } else {
3969 kh.e(i, i) = pal[i] * rs - cho * crpl[i];
3970 }
3971 kmats.e(i, i) += wpf * kh.e(i, i);
3972 }
3973 }
3974
3975 // compute coefficients for the nb-times squared exponential of the scaled matrix:
3976 for (k = 0; k < nb; ++k) {
3977 cho = al[0];
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); // define symmetric kh[][] to avoid case-distinctions
3985 // in computations below
3986 kmats.e(i, j) = kh.e(0, i) * al[j] + kh.e(0, j) * al[i];
3987 }
3988 al[i] *= cho;
3989 }
3990 for (l = 1; l < n; ++l) {
3991 cho = pal[n - 1];
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];
3997 }
3998 }
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];
4003 }
4004 }
4005 }
4006
4007
4008 // form output matrix omat :
4009 omat = al[0];
4010 for (k = 1; k < n; ++k) {
4011 mult_add(al[k], pl[k], omat);
4012 }
4013}
4014
4015// overload wrapper for chexpk where omat is not provided
4016template <int n, int m, typename T, typename MT>
4018 out_only Matrix_t<n, m, T, MT> &kmat) {
4019 static_assert(n == m, "chexpk() only for square matrices");
4021 chexpk(mat, omat, kmat);
4022 return omat;
4023}
4024
4025/**
4026 * @brief Calculate exp(mat).dagger()*mmat*exp(mat) and trace(exp(mat).dagger*mmat*dexp(mat))
4027 * from output of chexpk
4028 * @details exp is given and dexp computation is done using decompisition matrix kmats_{i,j},
4029 * as computed by chexpk using iterative Cayley-Hamilton (arXiv:2404.07704)
4030 * Calculate omat[i][j] = (exp(mat).dagger() * mmat * exp(mat))[i][j]
4031 * and domat[i][j] = trace(exp(mat).dagger() * mmat * dexp(mat)/dmat[j][i]) for
4032 * given matrices mat, exp(mat), kmats, and mmat
4033 * @tparam n Number of rowsMa
4034 * @tparam T Matrix element type
4035 * @tparam MT Matrix type
4036 * @param mat Matrix to compute exponential for
4037 * @param texp Matrix containing exp(mat)
4038 * @param kmats Matrix containing decomposition coefficients k_{i,j} of
4039 * dexp(X)^a_b/dX^c_d = k_{i,j} (X^i)^d_b (X^j)^a_c
4040 * @param mmat Matrix to multiply with
4041 * @param omat Matrix to which exp(mat).dagger()*mmat*exp(mat) gets stored
4042 * @param domat matrix to which trace(exp(mat).dagger()*mmat*dexp(mat)) gets stored
4043 * @return void
4044 */
4045template <int n, int m, typename T, typename MT>
4047 const Matrix_t<n, m, T, MT> &kmats, const Matrix_t<n, m, T, MT> &mmat,
4048 out_only Matrix_t<n, m, T, MT> &omat,
4049 out_only Matrix_t<n, m, T, MT> &domat) {
4050 static_assert(n == m, "mult_chexp() only for square matrices");
4051 // determine scaling factor:
4052 hila::arithmetic_type<T> sclim = sqrt(2.0), matnorm = norm(mat), sfac = 1.0;
4053 int nb = 0;
4054 while (matnorm * sfac >= sclim) {
4055 sfac *= 0.5;
4056 ++nb;
4057 }
4058
4059 // compute the first n matrix powers of mat and the corresponding traces :
4060 Matrix_t<n, m, T, MT> pl[n]; // the i-th matrix power of tU[][] is stored in pl[i][][]
4061 Matrix_t<n, m, T, MT> tomat, kh; // temp. storage for compuatation of derivative term
4062
4063 pl[1] = mat;
4064 pl[1] *= sfac;
4065 int i, j, k;
4066 for (i = 2; i < n; ++i) {
4067 j = i / 2;
4068 k = i % 2;
4069 mult(pl[j], pl[j + k], pl[i]);
4070 }
4071
4072 // tomat = texp.dagger() * mmat;
4073 mult(texp.dagger(), mmat, tomat);
4074
4075 // computing domat[ic1][ic2] = tr(exp(mat).dagger() * mmat * dexp(mat)/dU^{ic2}_{ic1})
4076 // from kmats[][] and the matrix powers of mat[][]:
4077 // domat = \sum_{i=0}^{n-1} pl[i] * tomat * \sum_{j=0}^{n-1} kmats[i][j] * pl[j]
4078
4079 // i=0: (treat i=0 case separately, since pl[0]=id is not used to avoid matrix-mult. by id)
4080 // j=0: (treat j=0 case separately, since pl[0]=id is not used)
4081 kh = kmats.e(0, 0);
4082 // j>0:
4083 for (j = 1; j < n; ++j) {
4084 mult_add(kmats.e(0, j), pl[j], kh);
4085 }
4086
4087 mult(tomat, kh, domat);
4088 // i>0:
4089 for (i = 1; i < n; ++i) {
4090 // j=0: (treat j=0 case separately, since pl[0]=id is not used)
4091 kh = kmats.e(0, i);
4092 // j>0: (note: kmats is symmetric; only have upper triangle set)
4093 for (j = 1; j < i; ++j) {
4094 mult_add(kmats.e(j, i), pl[j], kh);
4095 }
4096 for (j = i; j < n; ++j) {
4097 mult_add(kmats.e(i, j), pl[j], kh);
4098 }
4099 mult(pl[i], tomat, omat);
4100 mult_add(omat, kh, domat);
4101 }
4102
4103 // setting omat = exp(mat).dagger() * mmat * exp(mat) = tomat * exp(mat) :
4104 // omat = tomat * texp;
4105 mult(tomat, texp, omat);
4106}
4107
4108
4109// Calculate exp of a square matrix
4110// using iterative Cayley-Hamilton described in arXiv:2404.07704
4111/**
4112 * @brief Calculate exp of a square matrix
4113 * @details Computation is done using iterative Cayley-Hamilton (cf. from arXiv:2404.07704) with
4114 minimal temporary storage
4115
4116 * @tparam n Number of rowsMa
4117 * @tparam T Matrix element type
4118 * @tparam MT Matrix type
4119 * @param mat Matrix to compute exponential for
4120 * @return Matrix_t<n, m, T, MT>
4121 */
4122template <int n, int m, typename T, typename MT>
4124 static_assert(n == m, "chsexp() only for square matrices");
4125
4126 // compute the characteristic polynomial coefficients crpl[] with the Faddeev-LeVerrier
4127 // algorithm :
4128 int i, j, k;
4130 T crpl[n + 1];
4131 crpl[n] = 1.0;
4132 int ip = 0;
4133 tB[ip] = 1.;
4134 T tc = trace(mat);
4135 crpl[n - 1] = tc;
4136 tB[1 - ip] = mat;
4137 for (k = 2; k <= n; ++k) {
4138 tB[1 - ip] -= tc;
4139 mult(mat, tB[1 - ip], tB[ip]);
4140 tc = trace(tB[ip]) / k;
4141 crpl[n - k] = tc;
4142 ip = 1 - ip;
4143 }
4144
4145
4146 int mmax = 25 * n; // maximum number of Cayley-Hamilton iterations if no convergence is reached
4147 T al[n], pal[n]; // temp. Cayley-Hamilton coefficents
4148 hila::arithmetic_type<T> wpf = 1.0, twpf = 1.0, ttwpf; // leading coefficient of power series
4149 // set initial values for the n entries in al[] and pal[] :
4150 for (i = 0; i < n; ++i) {
4151 pal[i] = 0;
4152 al[i] = wpf;
4153 wpf /= (i + 1); // compute (i+1)-th power series coefficent from the i-th coefficient
4154 twpf += wpf;
4155 }
4156 pal[n - 1] = 1.0;
4157
4158 // next we iteratively add higher order power series terms to al[] till al[] stops changing
4159 // more precisely: the iteration will terminate as soon as twpf stops changing. Here twpf
4160 // is the sum \sum_{i=0}^{j} s_i/i!, with s_i referring to the magnitude the vector pal[]
4161 // would have at iteration i, if no renormalization were used.
4162 T ch, cho; // temporary variables for iteration
4163 hila::arithmetic_type<T> s, rs = 1.0, rss; // temp variables used for renormalization of pal[]
4164 for (j = n; j < mmax; ++j) {
4165 pal[n - 1] *= rs;
4166 ch = -pal[n - 1] * crpl[0];
4167 cho = pal[0] * rs;
4168 pal[0] = ch;
4169 s = ::squarenorm(ch);
4170 al[0] += wpf * ch;
4171 for (i = 1; i < n; ++i) {
4172 ch = cho - pal[n - 1] * crpl[i];
4173 cho = pal[i] * rs;
4174 pal[i] = ch;
4175 s += ::squarenorm(ch);
4176 al[i] += wpf * ch;
4177 }
4178
4179 s = sqrt(s);
4180 if (s > 1.0) {
4181 // if s is bigger than 1, normalize pal[] by a factor rs=1.0/s in next itaration,
4182 // and multiply wpf by s to compensate
4183 wpf *= s / (j + 1);
4184 rs = 1.0 / s;
4185 rss = 1.0;
4186 } else {
4187 wpf /= (j + 1);
4188 rs = 1.0;
4189 rss = s;
4190 }
4191 ttwpf = twpf;
4192 twpf += wpf * rss;
4193 if (ttwpf == twpf) {
4194 // terminate iteration
4195 break;
4196 }
4197 }
4198 // if(hila::myrank()==0) {
4199 // std::cout<<"chsexp niter: "<<j<<" ("<<j-n<<")"<<std::endl;
4200 // }
4201
4202 // form output matrix:
4203 ip = 0;
4204 tB[ip] = mat;
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];
4210 ip = 1 - ip;
4211 }
4212
4213 return tB[ip];
4214}
4215
4216
4217#include "datatypes/array.h"
4218
4219#include "datatypes/diagonal_matrix.h"
4220
4221#include "datatypes/matrix_linalg.h"
4222
4223// #include "datatypes/dagger.h"
4224
4225#endif
Definition of Array class.
Array type
Definition array.h:43
Complex definition.
Definition cmplx.h:58
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...
Definition matrix.h:106
void set_column(int c, const Vector< n, S > &v)
get column of a matrix
Definition matrix.h:424
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.
Definition matrix.h:1595
T trace() const
Computes Trace for Matrix.
Definition matrix.h:1171
Vector< n, T > column(int c) const
Returns column vector as value at index c.
Definition matrix.h:391
hila::arithmetic_type< T > norm() const
Calculate vector norm - sqrt of squarenorm.
Definition matrix.h:1225
void set_row(int r, const RowVector< m, S > &v)
Set row of Matrix with RowVector if types are assignable.
Definition matrix.h:376
auto & operator=(const Matrix_t< n, m, S, MT > &rhs) &
Copy matrix assignment.
Definition matrix.h:614
static constexpr int rows()
Define constant methods rows(), columns() - may be useful in template code.
Definition matrix.h:226
const auto & fill(const S rhs)
Matrix fill.
Definition matrix.h:1022
static constexpr int columns()
Returns column length.
Definition matrix.h:234
Mtype & random()
dot with matrix - matrix
Definition matrix.h:1393
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.
Definition matrix.h:1624
Rtype operator-(const S &b, const Mtype &a)
Subtraction operator Scalar - Matrix.
Definition matrix.h:2323
Mtype sort(Vector< N, int > &permutation, hila::sort order=hila::sort::ascending) const
Sort method for Vector which returns permutation order.
Definition matrix.h:1519
Matrix< n, m, hila::arithmetic_type< T > > imag() const
Returns imaginary part of Matrix or Vector.
Definition matrix.h:1156
void set_diagonal(const Vector< n, S > &v)
Set diagonal of square matrix to Vector which is passed to the method.
Definition matrix.h:463
auto & operator-=(const Matrix_t< n, m, S, MT > &rhs) &
Subtract assign operator Matrix to MAtrix.
Definition matrix.h:764
Mtype sort(hila::sort order=hila::sort::ascending) const
Sort method for Vector.
Definition matrix.h:1567
Rtype operator+(const Mtype &a, const S &b)
Addition operator Matrix + scalar.
Definition matrix.h:2216
auto & operator=(const S rhs) &
Assignment from scalar.
Definition matrix.h:658
const RowVector< n, T > & transpose() const
Transpose of vector.
Definition matrix.h:1056
auto & operator/=(const DiagonalMatrix< m, S > &rhs) &
Divide assign operator for DiagonalMatrix to Matrix.
Definition matrix.h:996
Matrix< n, m, R > operator*(const Mt1 &a, const Mt2 &b)
Multiplication operator.
Definition matrix.h:2406
static constexpr int size()
Returns size of Vector or square Matrix.
Definition matrix.h:248
Mtype permute(const Vector< N, int > &permutation) const
Permute elements of vector.
Definition matrix.h:1478
Rtype operator+(const Mtype1 &a, const Mtype2 &b)
Addition operator Matrix + Matrix.
Definition matrix.h:2111
Mtype & operator=(const std::nullptr_t &z)
Zero assignment.
Definition matrix.h:634
Matrix_t()=default
Define default constructors to ensure std::is_trivial.
Mtype conj() const
Returns complex conjugate of Matrix.
Definition matrix.h:1078
Mtype permute_rows(const Vector< n, int > &permutation) const
permute rows of Matrix
Definition matrix.h:1462
static constexpr bool is_vector()
Returns true if Matrix is a vector.
Definition matrix.h:132
T det() const
determinant function - if matrix size is < 5, use Laplace, otherwise LU
auto abs() const
Returns absolute value of Matrix.
Definition matrix.h:1123
Matrix< n, p, R > outer_product(const Matrix< p, q, S > &rhs) const
Outer product.
Definition matrix.h:1357
T max() const
Find max of Matrix only for arithmetic types.
Definition matrix.h:1239
Rtype dagger() const
Hermitian conjugate of matrix.
Definition matrix.h:1094
const auto & operator+() const
Addition operator.
Definition matrix.h:556
T e(const int i) const
Standard array indexing operation for vectors.
Definition matrix.h:305
Matrix< n, m, hila::arithmetic_type< T > > real() const
Returns real part of Matrix or Vector.
Definition matrix.h:1143
bool operator!=(const Matrix< n, m, S > &rhs) const
Boolean operator != to check if matrices are exactly different.
Definition matrix.h:593
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.
Definition matrix.h:2468
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.
Definition matrix.h:2253
const Vector< m, T > & transpose() const
Transpose of RowVector.
Definition matrix.h:1068
auto & operator+=(const S &rhs) &
Add assign operator scalar to Matrix.
Definition matrix.h:789
bool operator==(const Matrix< n2, m2, S > &rhs) const
Boolean operator == to determine if two matrices are exactly the same.
Definition matrix.h:571
T c[n *m]
The data as a one dimensional array.
Definition matrix.h:110
auto & operator+=(const Matrix_t< n, m, S, MT > &rhs) &
Add assign operator Matrix to Matrix.
Definition matrix.h:740
auto & operator-=(const DiagonalMatrix< n, S > &rhs) &
Subtract assign operator for DiagonalMatrix to Matrix.
Definition matrix.h:956
Mtype & gaussian_random(base_type width=1.0)
Fills Matrix with gaussian random elements from gaussian distribution.
Definition matrix.h:1415
T min() const
Find min of Matrix only for arithmetic types.
Definition matrix.h:1252
auto & operator-=(const S rhs) &
Subtract assign operator scalar to Matrix.
Definition matrix.h:815
T e(const int i, const int j) const
Standard array indexing operation for matrices.
Definition matrix.h:279
T det_laplace() const
Determinant of matrix, using Laplace method.
RowVector< m, T > row(int r) const
Return reference to row in a matrix.
Definition matrix.h:355
const DiagonalMatrix< n, T > & asDiagonalMatrix() const
Cast Vector to DiagonalMatrix.
Definition matrix.h:491
Rtype operator-(const Mtype &a, const S &b)
Subtraction operator Matrix - scalar.
Definition matrix.h:2285
Rtype transpose() const
Transpose of matrix.
Definition matrix.h:1039
auto & operator*=(const DiagonalMatrix< m, S > &rhs) &
Multiply assign operator for DiagonalMatrix to Matrix.
Definition matrix.h:975
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.
Definition matrix.h:1192
T operator[](const int i) const
Indexing operation [] defined only for vectors.
Definition matrix.h:332
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.
Definition matrix.h:2178
auto & operator+=(const DiagonalMatrix< n, S > &rhs) &
Addition assign operator for DiagonalMatrix to Matrix.
Definition matrix.h:939
auto & operator=(const DiagonalMatrix< n, S > &rhs) &
Assignment from diagonal matrix.
Definition matrix.h:681
Mtype permute_columns(const Vector< m, int > &permutation) const
permute columns of Matrix
Definition matrix.h:1448
R dot(const Matrix< p, q, S > &rhs) const
Dot product.
Definition matrix.h:1314
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.
Definition matrix.h:708
auto & operator*=(const S rhs) &
Multiply assign operator scalar to Matrix.
Definition matrix.h:897
const Array< n, m, T > & asArray() const
Cast Matrix to Array.
Definition matrix.h:476
Mtype operator-() const
Unary operator.
Definition matrix.h:531
auto & operator/=(const S rhs) &
Divide assign oprator scalar with matrix.
Definition matrix.h:923
DiagonalMatrix< n, T > diagonal()
Return diagonal of square matrix.
Definition matrix.h:439
std::string str(int prec=8, char separator=' ') const
Definition matrix.h:1709
Rtype adjoint() const
Adjoint of matrix.
Definition matrix.h:1111
auto & operator*=(const Matrix_t< m, p, S, MT > &rhs) &
Multiply assign scalar or matrix.
Definition matrix.h:845
hila::arithmetic_type< T > squarenorm() const
Calculate square norm - sum of squared elements.
Definition matrix.h:1209
static constexpr bool is_square()
Returns true if matrix is a square matrix.
Definition matrix.h:142
Matrix class which defines matrix operations.
Definition matrix.h:1742
Matrix(const std::nullptr_t &z)
Zero constructor.
Definition matrix.h:1830
Matrix(const Matrix &v)=default
Default copy constructor.
Matrix(const Matrix_t< n, m, S, MT > &rhs)
Copy constructor.
Definition matrix.h:1818
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.
Definition matrix.h:3068
Matrix(std::initializer_list< S > rhs)
Initializer list constructor.
Definition matrix.h:1851
hila::arithmetic_type< T > base_type
std incantation for field types
Definition matrix.h:1746
Matrix(const S rhs)
Scalar constructor.
Definition matrix.h:1788
Matrix()=default
Default constructo.
std::ostream & operator<<(std::ostream &strm, const Matrix_t< n, m, T, MT > &A)
Stream operator.
Definition matrix.h:2883
Definition of Complex types.
T arg(const Complex< T > &a)
Return argument of Complex number.
Definition cmplx.h:1278
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.
Definition matrix.h:4123
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))
Definition matrix.h:3629
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.
Definition matrix.h:3215
auto dagger(const Mtype &arg)
Return dagger of Matrix.
Definition matrix.h:1998
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))
Definition matrix.h:3170
Mt operator*(const Mt &a, const Mt &b)
Multiplication operator Square Matrix * Square Matrix.
Definition matrix.h:2351
Rtype operator+(Mtype1 a, const Mtype2 &b)
Real micro-optimization Matrix + Matrix - no extra creation of variable and copy.
Definition matrix.h:2141
auto abs(const Mtype &arg)
Return absolute value Matrix or Vector.
Definition matrix.h:2018
auto norm(const Mt &rhs)
Returns vector norm of Matrix.
Definition matrix.h:3047
void mult_add(const Mt1 &a, const Mt2 &b, Mt3 &res)
compute product of two matrices and add result to existing matrix
Definition matrix.h:2714
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...
Definition matrix.h:4046
auto conj(const Mtype &arg)
Return conjugate Matrix or Vector.
Definition matrix.h:1962
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.
Definition matrix.h:3851
auto squarenorm(const Mt &rhs)
Returns square norm of Matrix.
Definition matrix.h:3033
void mult_sub(const Mt1 &a, const Mt2 &b, Mt3 &res)
compute product of two matrices and subtract result from existing matrix
Definition matrix.h:2795
auto transpose(const Mtype &arg)
Return transposed Matrix or Vector.
Definition matrix.h:1942
auto real(const Mtype &arg)
Return real of Matrix or Vector.
Definition matrix.h:2059
auto imag(const Mtype &arg)
Return imaginary of Matrix or Vector.
Definition matrix.h:2079
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
Definition matrix.h:2633
void mult(const Mt1 &a, const Mt2 &b, Mt3 &res)
compute product of two matrices and write result to existing matrix
Definition matrix.h:2604
auto trace(const Mtype &arg)
Return trace of square Matrix.
Definition matrix.h:2039
auto mul_trace(const Mtype1 &a, const Mtype2 &b)
Returns trace of product of two matrices.
Definition matrix.h:2586
auto adjoint(const Mtype &arg)
Return adjoint Matrix.
Definition matrix.h:1982
Matrix_t< n, m, T, MT > altexp(const Matrix_t< n, m, T, MT > &mat, int &niter)
Calculate exp of a square matrix.
Definition matrix.h:3102
Implement hila::swap for gauge fields.
Definition array.h:982
double random()
Real valued uniform random number generator.
Definition hila_gpu.cpp:121
Array< n, m, Ntype > cast_to(const Array< n, m, T > &mat)
Array casting operation.
Definition array.h:1292
double gaussrand()
Gaussian random generation routine.
Definition random.cpp:240
double gaussrand2(double &out2)
hila::gaussrand2 returns 2 Gaussian distributed random numbers with variance .
Definition random.cpp:188
decltype(std::declval< A >()+std::declval< B >()) type_plus
Definition type_tools.h:103
std::string to_string(const Array< n, m, T > &A, int prec=8, char separator=' ')
Converts Array object to string.
Definition array.h:996
type to store the return value of eigen_hermitean(): {E, U} where E is nxn DiagonalMatrix containing ...
Definition matrix.h:73
hila::is_complex_or_arithmetic<T>::value
Definition cmplx.h:708
type to store the return combo of svd: {U, D, V} where U and V are nxn unitary / orthogonal,...
Definition matrix.h:61