1#ifndef DIAGONAL_MATRIX_H_
2#define DIAGONAL_MATRIX_H_
17template <
int n,
typename T>
22 "DiagonalMatrix requires Complex or arithmetic type");
25 using base_type = hila::arithmetic_type<T>;
26 using argument_type = T;
36 template <
typename S, std::enable_if_t<(hila::is_assignable<T &, S>::value),
int> = 0>
38 for (
int i = 0; i < n; i++)
45 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
47 assert(rhs.size() == n &&
"Matrix/Vector initializer list size must match variable size");
49 for (
auto it = rhs.begin(); it != rhs.end(); it++, i++) {
57 static constexpr int rows() {
60 static constexpr int columns() {
63 static constexpr int size() {
72 inline T
e(
const int i)
const {
77 inline T &
e(
const int i) const_function {
82 T
e(
const int i,
const int j)
const {
117 for (
int i = 0; i < n; i++) {
133 template <
typename S>
135 for (
int i = 0; i < n; i++) {
136 if (
e(i) != rhs.
e(i))
145 template <
typename S,
typename Mtype>
147 for (
int i = 0; i < n; i++) {
148 for (
int j = 0; j < n; j++) {
150 if (
e(i) != rhs.
e(i, i))
153 if (rhs.
e(i, j) != 0)
165 template <
typename S>
167 return !(*
this == rhs);
176#pragma hila loop_function
177 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
180 for (
int i = 0; i < n; i++) {
187#pragma hila loop_function
188 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
191 for (
int i = 0; i < n; i++) {
198#pragma hila loop_function
199 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
201 assert(rhs.size() == n &&
"Initializer list has a wrong size in assignment");
203 for (
auto it = rhs.begin(); it != rhs.end(); it++, i++) {
211#pragma hila loop_function
212 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
214 for (
int i = 0; i < n; i++) {
220#pragma hila loop_function
221 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
223 for (
int i = 0; i < n; i++) {
230#pragma hila loop_function
231 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
233 for (
int i = 0; i < n; i++) {
239#pragma hila loop_function
240 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
242 for (
int i = 0; i < n; i++) {
249#pragma hila loop_function
250 template <
typename S,
251 std::enable_if_t<hila::is_assignable<T &, hila::type_mul<T, S>>::value,
int> = 0>
253 for (
int i = 0; i < n; i++) {
259#pragma hila loop_function
260 template <
typename S,
261 std::enable_if_t<hila::is_assignable<T &, hila::type_mul<T, S>>::value,
int> = 0>
263 for (
int i = 0; i < n; i++) {
270#pragma hila loop_function
271 template <
typename S,
272 std::enable_if_t<hila::is_assignable<T &, hila::type_div<T, S>>::value,
int> = 0>
274 for (
int i = 0; i < n; i++) {
280#pragma hila loop_function
281 template <
typename S,
282 std::enable_if_t<hila::is_assignable<T &, hila::type_div<T, S>>::value,
int> = 0>
284 for (
int i = 0; i < n; i++) {
294 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value,
int> = 0>
296 for (
int i = 0; i < n; i++)
313 for (
int i = 0; i < n; i++)
337 for (
int i = 0; i < n; i++) {
338 res.c[i] =
::abs(c[i]);
350 for (
int i = 0; i < n; i++) {
358 for (
int i = 0; i < n; i++) {
359 res.c[i] = ::imag(c[i]);
368 template <typename S = T, std::enable_if_t<hila::is_arithmetic<S>::value,
int> = 0>
371 for (
int i = 1; i < n; i++) {
378 template <typename S = T, std::enable_if_t<hila::is_arithmetic<S>::value,
int> = 0>
381 for (
int i = 1; i < n; i++) {
390 for (
int i = 1; i < n; i++)
397 for (
int i = 1; i < n; i++)
402 auto squarenorm()
const {
403 hila::arithmetic_type<T> res(0);
404 for (
int i = 0; i < n; i++)
405 res += ::squarenorm(c[i]);
409 hila::arithmetic_type<T> norm()
const {
410 return sqrt(squarenorm());
419 static_assert(hila::is_floating_point<hila::arithmetic_type<T>>::value,
420 "DiagonalMatrix random() requires non-integral type elements");
421 for (
int i = 0; i < n; i++) {
429 static_assert(hila::is_floating_point<hila::arithmetic_type<T>>::value,
430 "DiagonaMatrix gaussian_random() requires non-integral type elements");
431 for (
int i = 0; i < n; i++) {
441 std::string
str(
int prec = 8,
char separator =
' ')
const {
442 return this->asArray().str(prec, separator);
450 for (
int i = 0; i < n; i++) {
451 for (
int j = 0; j < n; j++) {
477 return *(
reinterpret_cast<const Vector<n,T> *
>(
this));
484 hila::sort order = hila::sort::ascending)
const {
485 return this->asArray().sort(permutation, order).asDiagonalMatrix();
490 return this->asArray().sort(order).asDiagonalMatrix();
496template <
int n,
typename T>
501template <
int n,
typename T>
506template <
int n,
typename T>
511template <
int n,
typename T>
513 return arg.adjoint();
516template <
int n,
typename T>
521template <
int n,
typename T>
526template <
int n,
typename T>
531template <
int n,
typename T>
536template <
int n,
typename T>
538 return arg.squarenorm();
541template <
int n,
typename T>
546template <
int n,
typename T>
562template <
typename Mt,
typename S,
typename Enable =
void>
563struct diagonalmatrix_scalar_op_s {
569template <
typename Mt,
typename S>
570struct diagonalmatrix_scalar_op_s<
572 typename
std::enable_if_t<std::is_convertible<hila::type_plus<hila::number_type<Mt>, S>,
573 hila::number_type<Mt>>::value>> {
575 using type =
typename std::conditional<
576 hila::is_floating_point<hila::arithmetic_type<Mt>>::value, Mt,
581template <
typename Mt,
typename S>
582using diagonalmatrix_scalar_type =
typename diagonalmatrix_scalar_op_s<Mt, S>::type;
592template <
int n,
typename T,
typename S,
593 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value,
int> = 0,
594 typename Rtype = hila::diagonalmatrix_scalar_type<DiagonalMatrix<n, T>, S>>
597 for (
int i = 0; i < n; i++)
598 res.e(i) = a.
e(i) + b;
603template <
int n,
typename T,
typename S,
604 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value,
int> = 0,
605 typename Rtype = hila::diagonalmatrix_scalar_type<DiagonalMatrix<n, T>, S>>
611template <
int n,
typename T,
typename S,
612 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value,
int> = 0,
613 typename Rtype = hila::diagonalmatrix_scalar_type<DiagonalMatrix<n, T>, S>>
616 for (
int i = 0; i < n; i++)
617 res.e(i) = a.
e(i) - b;
622template <
int n,
typename T,
typename S,
623 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value,
int> = 0,
624 typename Rtype = hila::diagonalmatrix_scalar_type<DiagonalMatrix<n, T>, S>>
627 for (
int i = 0; i < n; i++)
628 res.e(i) = b - a.
e(i);
633template <
int n,
typename T,
typename S,
634 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value,
int> = 0,
635 typename Rtype = hila::diagonalmatrix_scalar_type<DiagonalMatrix<n, T>, S>>
638 for (
int i = 0; i < n; i++)
639 res.e(i) = a.
e(i) * b;
644template <
int n,
typename T,
typename S,
645 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value,
int> = 0,
646 typename Rtype = hila::diagonalmatrix_scalar_type<DiagonalMatrix<n, T>, S>>
652template <
int n,
typename T,
typename S,
653 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value,
int> = 0,
654 typename Rtype = hila::diagonalmatrix_scalar_type<DiagonalMatrix<n, T>, S>>
657 for (
int i = 0; i < n; i++)
658 res.e(i) = a.
e(i) / b;
663template <
int n,
typename T,
typename S,
664 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value,
int> = 0,
665 typename Rtype = hila::diagonalmatrix_scalar_type<DiagonalMatrix<n, T>, S>>
668 for (
int i = 0; i < n; i++)
669 res.e(i) = b / a.
e(i);
675template <
int n,
typename A,
typename B,
typename R = hila::type_plus<A, B>>
678 for (
int i = 0; i < n; i++)
679 res.
e(i) = a.
e(i) + b.
e(i);
683template <
int n,
typename A,
typename B,
typename R = hila::type_minus<A, B>>
686 for (
int i = 0; i < n; i++)
687 res.
e(i) = a.
e(i) - b.
e(i);
691template <
int n,
typename A,
typename B,
typename R = hila::type_mul<A, B>>
694 for (
int i = 0; i < n; i++)
695 res.
e(i) = a.
e(i) * b.
e(i);
699template <
int n,
typename A,
typename B,
typename R = hila::type_div<A, B>>
702 for (
int i = 0; i < n; i++)
703 res.
e(i) = a.
e(i) / b.
e(i);
709template <
int n,
typename T,
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0,
710 typename Rtype = hila::mat_x_mat_type<Matrix<n, n, T>, Mtype>>
716 static_assert(mc == n && mr == n,
"Matrix sizes do not match");
720 for (
int i = 0; i < n; i++)
725template <
int n,
typename T,
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0,
726 typename Rtype = hila::mat_x_mat_type<Matrix<n, n, T>, Mtype>>
731template <
int n,
typename T,
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0,
732 typename Rtype = hila::mat_x_mat_type<Matrix<n, n, T>, Mtype>>
738 static_assert(mc == n && mr == n,
"Matrix sizes do not match");
741 for (
int i = 0; i < n; i++)
746template <
int n,
typename T,
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0,
747 typename Rtype = hila::mat_x_mat_type<Matrix<n, n, T>, Mtype>>
752 static_assert(mc == n && mr == n,
"Matrix sizes do not match");
755 for (
int i = 0; i < n; i++)
763template <
int n,
typename T,
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0,
764 typename Rtype = hila::mat_x_mat_type<Matrix<n, n, T>, Mtype>>
770 static_assert(mr == n,
"Matrix sizes do not match");
773 for (
int i = 0; i < n; i++)
774 for (
int j = 0; j < mc; j++)
775 r.e(i, j) = a.
e(i) * b.
e(i, j);
779template <
int n,
typename T,
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0,
780 typename Rtype = hila::mat_x_mat_type<Matrix<n, n, T>, Mtype>>
786 static_assert(mc == n,
"Matrix sizes do not match");
789 for (
int i = 0; i < mr; i++)
790 for (
int j = 0; j < n; j++)
791 r.e(i, j) = b.
e(i, j) * a.
e(j);
796template <
int n,
typename T,
typename Mtype, std::enable_if_t<Mtype::is_matrix(),
int> = 0,
797 typename Rtype = hila::mat_x_mat_type<Matrix<n, n, T>, Mtype>>
803 static_assert(mc == n,
"Matrix sizes do not match");
806 for (
int i = 0; i < mr; i++)
807 for (
int j = 0; j < n; j++)
808 r.e(i, j) = b.
e(i, j) / a.
e(j);
816template <
int n,
typename T>
818 for (
int i = 0; i < n; i++)
819 a.c[i] = sqrt(a.c[i]);
823template <
int n,
typename T>
825 for (
int i = 0; i < n; i++)
826 a.c[i] = cbrt(a.c[i]);
830template <
int n,
typename T>
832 for (
int i = 0; i < n; i++)
833 a.c[i] = exp(a.c[i]);
837template <
int n,
typename T>
839 for (
int i = 0; i < n; i++)
840 a.c[i] =
log(a.c[i]);
844template <
int n,
typename T>
846 for (
int i = 0; i < n; i++)
847 a.c[i] = sin(a.c[i]);
851template <
int n,
typename T>
853 for (
int i = 0; i < n; i++)
854 a.c[i] = cos(a.c[i]);
858template <
int n,
typename T>
860 for (
int i = 0; i < n; i++)
861 a.c[i] = tan(a.c[i]);
865template <
int n,
typename T>
867 for (
int i = 0; i < n; i++)
868 a.c[i] = asin(a.c[i]);
872template <
int n,
typename T>
874 for (
int i = 0; i < n; i++)
875 a.c[i] = acos(a.c[i]);
879template <
int n,
typename T>
881 for (
int i = 0; i < n; i++)
882 a.c[i] = atan(a.c[i]);
886template <
int n,
typename T>
888 for (
int i = 0; i < n; i++)
889 a.c[i] = sinh(a.c[i]);
893template <
int n,
typename T>
895 for (
int i = 0; i < n; i++)
896 a.c[i] = cosh(a.c[i]);
900template <
int n,
typename T>
902 for (
int i = 0; i < n; i++)
903 a.c[i] = tanh(a.c[i]);
907template <
int n,
typename T>
909 for (
int i = 0; i < n; i++)
910 a.c[i] = asinh(a.c[i]);
914template <
int n,
typename T>
916 for (
int i = 0; i < n; i++)
917 a.c[i] = acosh(a.c[i]);
921template <
int n,
typename T>
923 for (
int i = 0; i < n; i++)
924 a.c[i] = atanh(a.c[i]);
931 int n,
typename T,
typename S,
932 std::enable_if_t<hila::is_arithmetic<S>::value || hila::contains_complex<T>::value,
int> = 0>
934 for (
int i = 0; i < n; i++)
935 a.c[i] = pow(a.c[i], p);
940template <
int n,
typename T,
typename S,
941 std::enable_if_t<!hila::contains_complex<T>::value,
int> = 0>
944 for (
int i = 0; i < n; i++)
945 res.
e(i) = pow(a.
e(i), p);
954template <
typename Ntype,
typename T,
int n,
955 std::enable_if_t<hila::is_arithmetic<T>::value,
int> = 0>
958 for (
int i = 0; i < n; i++)
963template <typename Ntype, typename T, int n, std::enable_if_t<hila::is_complex<T>::value,
int> = 0>
966 for (
int i = 0; i < n; i++)
967 res.c[i] = cast_to<Ntype>(mat.c[i]);
972template <
int n,
typename T>
979template <
int n,
typename T>
981 return to_string(A.asArray(), prec, separator);
984template <
int n,
typename T>
986 return prettyprint(A.
toMatrix(), prec);
Array< n, m, T > conj(const Array< n, m, T > &arg)
Return conjugate Array.
std::ostream & operator<<(std::ostream &strm, const Array< n, m, T > &A)
Stream operator.
Array< n, m, hila::arithmetic_type< T > > imag(const Array< n, m, T > &arg)
Return imaginary part of Array.
hila::arithmetic_type< T > squarenorm(const Array< n, m, T > &rhs)
Return square norm of Array.
Array< n, m, hila::arithmetic_type< T > > real(const Array< n, m, T > &arg)
Return real part of Array.
Define type DiagonalMatrix<n,T>
bool operator==(const Matrix_t< n, n, S, Mtype > &rhs) const
Boolean operator == to compare with square matrix.
Matrix< n, n, T > toMatrix() const
convert to generic matrix
const DiagonalMatrix & transpose() const
transpose - leaves diagonal matrix as is
DiagonalMatrix & operator=(const DiagonalMatrix< n, S > &rhs)
bool operator!=(const S &rhs) const
Boolean operator != to check if matrices are exactly different.
Vector< n, T > column(int i) const
Returns column vector i.
T e(const int i) const
Element access - e(i) gives diagonal element i.
DiagonalMatrix conj() const
conj - conjugate elements
DiagonalMatrix & random()
Fills Matrix with random elements.
std::string str(int prec=8, char separator=' ') const
convert to string for printing
RowVector< n, T > row(int i) const
Return row from Diagonal matrix.
DiagonalMatrix dagger() const
dagger - conjugate elements
DiagonalMatrix()=default
Define default constructors to ensure std::is_trivial.
DiagonalMatrix< n, T > operator-() const
Unary - operator.
DiagonalMatrix & fill(const S rhs)
fill with constant value - same as assignment
T max() const
Find max or min value - only for arithmetic types.
static constexpr int rows()
Returns matrix size - all give same result.
const auto & operator+() const
Unary + operator.
DiagonalMatrix< n, T > sort(Vector< N, int > &permutation, hila::sort order=hila::sort::ascending) const
implement sort as casting to matrix
bool operator==(const DiagonalMatrix< n, S > &rhs) const
Boolean operator == to determine if two matrices are exactly the same.
auto real() const
real and imaginary parts of diagonal matrix
auto abs() const
absolute value of all elements
The main matrix type template Matrix_t. This is a base class type for "useful" types which are deriv...
static constexpr int rows()
Define constant methods rows(), columns() - may be useful in template code.
static constexpr int columns()
Returns column length.
T e(const int i, const int j) const
Standard array indexing operation for matrices.
Matrix class which defines matrix operations.
Complex< T > dagger(const Complex< T > &val)
Return dagger of Complex number.
T abs(const Complex< T > &a)
Return absolute value of Complex number.
T arg(const Complex< T > &a)
Return argument of Complex number.
Definition of Matrix types.
Invert diagonal + const. matrix using Sherman-Morrison formula.
T gaussian_random()
Template function T hila::gaussian_random<T>(),generates gaussian random value of type T,...
logger_class log
Now declare the logger.
double random()
Real valued uniform random number generator.
decltype(std::declval< A >()+std::declval< B >()) type_plus
std::string to_string(const Array< n, m, T > &A, int prec=8, char separator=' ')
Converts Array object to string.
hila::is_complex_or_arithmetic<T>::value