HILA
Loading...
Searching...
No Matches
array.h
Go to the documentation of this file.
1#ifndef HILA_ARRAY_H_
2#define HILA_ARRAY_H_
3
4/**
5 * @file array.h
6 * @brief Definition of Array class
7 * @details This file contains the definitions of Array class and utility functions related to it.
8 */
9
10#include "datatypes/matrix.h"
11
12/**
13 * @brief \f$ n\times m \f$ Array type
14 * @details Acts as array class which stores data in a simple C style array.
15 *
16 * Main functionality which the Array class offers is to supplement the fall backs of storing
17 * information in a Matrix data structure.
18 *
19 * For example assigning a value to each element with the Matrix class is not directly possible
20 * using the assignment operator=. This is because assignment with matrices is defined as
21 * \f$ M = a = a*I \f$ where M is a general matrix, a is a scalar and I an identity matrix. The
22 * result would only assign a to the diagonal elements.
23 *
24 * Unlike the Matrix object, the Array object assigns element wise. This allows filling the Array
25 * with the assignment operator. Additionally element wise operations are useable as long as they
26 * are defined for the Array type. For example the operation:
27 *
28 * \code
29 * Array<n,m,double> m.random();
30 * sin(m);
31 * \endcode
32 *
33 * is defined, since the \f$sin\f$ function is defined for doubles.
34 *
35 * The above operation would not work for matrices, but with casting operations. Matrix::asArray and
36 * Array::asMatrix one can take full advantage of element wise operations.
37 *
38 * @tparam n Number of rows
39 * @tparam m Number of columns
40 * @tparam T Array element type
41 */
42template <const int n, const int m, typename T>
43class Array {
44 public:
45 static_assert(hila::is_complex_or_arithmetic<T>::value || hila::is_extended<T>::value,
46 "Array requires Complex or arithmetic type");
47
48 // Store Array contents in one dimensional C style array
49 T c[n * m];
50
51 // std incantation for field types
52 using base_type = hila::arithmetic_type<T>;
53 using argument_type = T;
54
55 static constexpr bool is_array() {
56 return true;
57 }
58
59 /**
60 * @brief Default constructor
61 * @details The following ways of constructing a matrix are:
62 *
63 * Allocates undefined \f$ n\times m\f$ Array.
64 *
65 * \code{.cpp}
66 * Array<n,m,MyType> A;
67 * \endcode
68 *
69 */
70 Array() = default;
71 ~Array() = default;
72 /**
73 * @brief Copy constructor
74 * @details
75 * Construction from previously defined Array.
76 *
77 * \code{.cpp}
78 * Array<n,m,MyOtherType> A_0;
79 * //
80 * // A_0 is filled with content
81 * //
82 * Array<n,m,MyType> A = A_0;
83 * \endcode
84 *
85 * @param v Array to assign
86 */
87 Array(const Array<n, m, T> &v) = default;
88
89 /**
90 * @brief Scalar constructor
91 * @details
92 *
93 * Construct with given scalar which is assigned to all elements in array
94 *
95 * \code{.cpp}
96 * MyType x = hila::random();
97 * Array<n,m,MyType> A = x;
98 * \endcode
99 *
100 * @tparam S Type for scalar
101 * @param rhs Scalar to assign
102 */
103 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
104 explicit inline Array(const S rhs) {
105 for (int i = 0; i < n * m; i++) {
106 this->c[i] = rhs;
107 }
108 }
109
110 // and make non-explicit constructor from 0
111 inline Array(const std::nullptr_t &z) {
112 for (int i = 0; i < n * m; i++)
113 c[i] = 0;
114 }
115
116 // Construct array automatically from right-size initializer list
117 // This does not seem to be dangerous, so keep non-explicit
118
119 /**
120 * @brief Initializer list constructor
121 *
122 * Construction from c++ initializer list.
123 *
124 * \code{.cpp}
125 * Array<2,2,int> A = {1, 0
126 * 0, 1};
127 * \endcode
128 * @tparam S Element type of initializer list
129 * @param rhs Initializer list to assign
130 */
131 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
132 inline Array(std::initializer_list<S> rhs) {
133 assert(rhs.size() == n * m && "Array initializer list size must match variable size");
134 int i = 0;
135 for (auto it = rhs.begin(); it != rhs.end(); it++, i++) {
136 c[i] = *it;
137 }
138 }
139
140 /**
141 * @brief Returns number of rows
142 *
143 * @return constexpr int
144 */
145 constexpr int rows() const {
146 return n;
147 }
148 /**
149 * @brief Returns number of columns
150 *
151 * @return constexpr int
152 */
153 constexpr int columns() const {
154 return m;
155 }
156
157 /**
158 * @brief Returns size of #Vector Array or square Array
159 *
160 * @tparam q Number of rows n
161 * @tparam p Number of columns m
162 * @return constexpr int
163 */
164 template <int q = n, int p = m, std::enable_if_t<q == 1, int> = 0>
165 constexpr int size() const {
166 return p;
167 }
168
169 template <int q = n, int p = m, std::enable_if_t<p == 1, int> = 0>
170 constexpr int size() const {
171 return q;
172 }
173
174 template <int q = n, int p = m, std::enable_if_t<q == p, int> = 0>
175 constexpr int size() const {
176 return q;
177 }
178
179 /**
180 * @brief Standard array indexing operation for 2D Array
181 *
182 * @details Accessing singular elements is insufficient, but matrix elements are often quite
183 * small.
184 *
185 * Example for 2D Array:
186 * \code
187 * Array<n,m,MyType> A;
188 * MyType a = A.e(i,j); // i <= n, j <= m
189 * \endcode
190 *
191 * @param i row index
192 * @param j column index
193 * @return T matrix element type
194 */
195 inline T e(const int i, const int j) const {
196 return c[i * m + j];
197 }
198
199 /// @brief Const overload
200 /// @internal
201 inline T &e(const int i, const int j) const_function {
202 return c[i * m + j];
203 }
204
205 /// @brief Const overload
206 /// @internal
207 template <int q = n, int p = m, std::enable_if_t<(q == 1 || p == 1), int> = 0>
208 inline T e(const int i) const {
209 return c[i];
210 }
211
212 /**
213 * @brief Standard array indexing operation for 1D Array
214 * @details Example for Vector Array:
215 * \code {.cpp}
216 * Array1d<n,MyType> A;
217 * MyType a = A.e(i); // i <= n
218 * \endcode
219 * @param i Vector index type
220 * @return T matrix element type
221 */
222 template <int q = n, int p = m, std::enable_if_t<(q == 1 || p == 1), int> = 0>
223 inline T &e(const int i) const_function {
224 return c[i];
225 }
226
227 /**
228 * @brief Cast Array to Matrix
229 *
230 * @return Matrix<n, m, T>&
231 */
232 Matrix<n, m, T> &asMatrix() const_function {
233 return *reinterpret_cast<Matrix<n, m, T> *>(this);
234 }
235
236//#pragma hila loop_function
237 const Matrix<n, m, T> &asMatrix() const {
238 return *reinterpret_cast<const Matrix<n, m, T> *>(this);
239 }
240
241 /**
242 * @brief Cast Array1D to Vector
243 * @details asMatrix will perform the same operation.
244 *
245 * @return Vector<n, T>&
246 */
247 Vector<n, T> &asVector() const_function {
248 static_assert(1 == m, "asVector() only for column arrays");
249 return *reinterpret_cast<Vector<n, T> *>(this);
250 }
251
252//#pragma hila loop_function
253 const Vector<n, T> &asVector() const {
254 static_assert(1 == m, "asVector() only for column arrays");
255 return *reinterpret_cast<const Vector<n, T> *>(this);
256 }
257
258 DiagonalMatrix<n, T> &asDiagonalMatrix() const_function {
259 static_assert(1 == m, "asDiagonalMatrix() only for column arrays");
260 return *reinterpret_cast<DiagonalMatrix<n, T> *>(this);
261 }
262
263//#pragma hila loop_function
264 const DiagonalMatrix<n, T> &asDiagonalMatrix() const {
265 static_assert(1 == m, "asDiagonalMatrix() only for column arrays");
266 return *reinterpret_cast<const DiagonalMatrix<n, T> *>(this);
267 }
268
269
270 // casting from one Array (number) type to another
271 // TODO: CHECK AVX CONVERSIONS
272 template <typename S, std::enable_if_t<hila::is_assignable<S &, T>::value, int> = 0>
273 operator Array<n, m, S>() {
274 Array<n, m, S> res;
275 for (int i = 0; i < n * m; i++) {
276 res.c[i] = c[i];
277 }
278 return res;
279 }
280
281 /**
282 * @brief Unary - operator
283 * @details Returns Array with the signs of all the elements in the Arrat flipped.
284 *
285 * @return Array<n,m,T>
286 */
287 inline Array<n, m, T> operator-() const {
288 Array<n, m, T> res;
289 for (int i = 0; i < n * m; i++) {
290 res.c[i] = -c[i];
291 }
292 return res;
293 }
294
295 /**
296 * @brief Unary + operator
297 * @details Equivalent to identity operator meaning that Arrat stays as is.
298 *
299 * @return const Array<n, m, T>
300 */
301 inline Array<n, m, T> operator+() const {
302 return *this;
303 }
304
305 /**
306 * @brief Scalar assignment operator
307 * @details The following ways to assign an Array are:
308 *
309 *
310 * __Assignment from Array__:
311 *
312 * \code {.cpp}
313 * Array<n,m,MyType> A_0;
314 * .
315 * . A_0 has values assigned to it
316 * .
317 * Array<n,m,MyType> A; \\ undefined matrix
318 * A = A_0; \\ Assignment from A_0
319 * \endcode
320 *
321 * __Assignment from scalar__:
322 *
323 * Assignment from scalar assigns the scalar uniformly to the array
324 *
325 * \code {.cpp}
326 * MyType a = hila::random;
327 * Array<n,m,MyType> A;
328 * A = a;
329 * \endcode
330 *
331 */
332
333 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
334 inline Array<n, m, T> &operator=(const S rhs) out_only & {
335 for (int i = 0; i < n * m; i++) {
336 c[i] = rhs;
337 }
338 return *this;
339 }
340
341 template <typename S, std::enable_if_t<hila::is_assignable<T &, S>::value, int> = 0>
342 inline Array<n, m, T> &operator=(const Array<n, m, S> &rhs) out_only & {
343 for (int i = 0; i < n; i++)
344 for (int j = 0; j < m; j++) {
345 e(i, j) = rhs.e(i, j);
346 }
347 return *this;
348 }
349
350
351 /**
352 * @brief Compare equality of Arrays
353 *
354 * Two Arrays are equal iff Arrays are of same dimension and all elements compare to equal
355 * Note: Complex == scalar if arithmetic value is equal
356 *
357 * @tparam S
358 * @param rhs
359 * @return true if equal
360 */
361
362 template <typename S, int n1, int m1>
363 bool operator==(const Array<n1, m1, S> &rhs) const {
364 if constexpr (n != n1 || m != m1)
365 return false;
366
367 for (int i = 0; i < n; i++)
368 for (int j = 0; j < m; j++) {
369 if (e(i, j) != rhs.e(i, j))
370 return false;
371 }
372 return true;
373 }
374
375 /**
376 * @brief Compare non-equality of two Arrays
377 *
378 * Negation of operator==()
379 *
380 * @tparam S
381 * @tparam n1
382 * @tparam m1
383 * @param rhs
384 * @return true
385 * @return false
386 */
387
388 template <typename S, int n1, int m1>
389 bool operator!=(const Array<n1, m1, S> &rhs) const {
390 return !(*this == rhs);
391 }
392
393 /**
394 * @brief Add assign operator with Array or scalar
395 * @details Add assign operator can be used in the following ways
396 *
397 * __Add assign Array__:
398 *
399 * Requires that Arrays have same dimensions
400 *
401 * \code {.cpp}
402 * Array<n,m,MyType> A,B;
403 * A = 1;
404 * B = 1;
405 * A += B; \\ A is uniformly 2
406 * \endcode
407 *
408 * __Add assign scalar__:
409 *
410 * Adds scalar \f$ a \f$ to Array uniformly
411 *
412 * \code {.cpp}
413 * Array<n,m,MyType> A = 1;
414 * A += 1 ; \\ A is uniformly 2
415 * \endcode
416 *
417 * @tparam S Element type of rhs
418 * @param rhs Array to add
419 * @return Array<n, m, T>
420 */
421
422 template <typename S, std::enable_if_t<std::is_convertible<S, T>::value, int> = 0>
424 for (int i = 0; i < n; i++)
425 for (int j = 0; j < m; j++) {
426 e(i, j) += rhs.e(i, j);
427 }
428 return *this;
429 }
430
431 /**
432 * @brief Subtract assign operator with Array or scalar
433 * @details Subtract assign operator can be used in the following ways
434 *
435 * __Subtract assign Array__:
436 *
437 * \code {.cpp}
438 * Array<n,m,MyType> A,B;
439 * A = 3;
440 * B = 1;
441 * A -= B; \\ A is uniformly 2
442 * \endcode
443 *
444 * __Subtract assign scalar__:
445 *
446 * Subtract scalar uniformly from square matrix
447 *
448 * \code {.cpp}
449 * A<n,m,MyType> A = 3;
450 * A -= 1 ; \\ A is uniformly 2
451 * \endcode
452 * @tparam S
453 * @param rhs
454 * @return Array<n, m, T>&
455 */
456 template <typename S, std::enable_if_t<std::is_convertible<S, T>::value, int> = 0>
458 for (int i = 0; i < n; i++)
459 for (int j = 0; j < m; j++) {
460 e(i, j) -= rhs.e(i, j);
461 }
462 return *this;
463 }
464
465 // add assign type T and convertible
466 template <typename S,
467 std::enable_if_t<std::is_convertible<hila::type_plus<T, S>, T>::value, int> = 0>
468 Array<n, m, T> &operator+=(const S rhs) & {
469 for (int i = 0; i < n * m; i++) {
470 c[i] += rhs;
471 }
472 return *this;
473 }
474
475 // subtract assign type T and convertible
476 template <typename S,
477 std::enable_if_t<std::is_convertible<hila::type_minus<T, S>, T>::value, int> = 0>
478 Array<n, m, T> &operator-=(const S rhs) & {
479 for (int i = 0; i < n * m; i++) {
480 c[i] -= rhs;
481 }
482 return *this;
483 }
484
485 /**
486 * @brief Multiply assign scalar or array
487 * @details Multiplication works element wise
488 *
489 * Multiply assign operator can be used in the following ways
490 *
491 * __Multiply assign Array__:
492 *
493 * \code {.cpp}
494 * Array<n,m,MyType> A,B;
495 * A = 2;
496 * B = 2;
497 * A *= B; \\ A is uniformly 4
498 * \endcode
499 *
500 * __Multiply assign scalar__:
501 *
502 * \code {.cpp}
503 * Array<n,m,MyType> A = 1;
504
505 * A *= 2 ; \\ A is uniformly 2
506 * \endcode
507 * @tparam S
508 * @param rhs
509 * @return Array<n, m, T>&
510 */
511 template <typename S,
512 std::enable_if_t<std::is_convertible<hila::type_mul<T, S>, T>::value, int> = 0>
514 for (int i = 0; i < n; i++)
515 for (int j = 0; j < m; j++) {
516 e(i, j) *= rhs.e(i, j);
517 }
518 return *this;
519 }
520
521 /// multiply assign with scalar
522 template <typename S,
523 std::enable_if_t<std::is_convertible<hila::type_mul<T, S>, T>::value, int> = 0>
524 Array<n, m, T> &operator*=(const S rhs) & {
525 for (int i = 0; i < n * m; i++) {
526 c[i] *= rhs;
527 }
528 return *this;
529 }
530
531 /**
532 * @brief Division assign with array or scalar
533 * @details Division works element wise
534 *
535 * Division assign operator can be used in the following ways
536 *
537 * __Division assign Array__:
538 *
539 * \code {.cpp}
540 * Array<n,m,MyType> A,B;
541 * A = 2;
542 * B = 2;
543 * A /= B; \\ A is uniformly 1
544 * \endcode
545 *
546 * __Division assign scalar__:
547 *
548 * \code {.cpp}
549 * Array<n,m,MyType> A = 2;
550
551 * A /= 2 ; \\ A is uniformly 1
552 * \endcode
553 * @tparam S
554 * @param rhs
555 * @return Array<n, m, T>&
556 */
557 template <typename S,
558 std::enable_if_t<std::is_convertible<hila::type_div<T, S>, T>::value, int> = 0>
560 for (int i = 0; i < n; i++)
561 for (int j = 0; j < m; j++) {
562 e(i, j) /= rhs.e(i, j);
563 }
564 return *this;
565 }
566
567 // divide assign with scalar
568 template <typename S,
569 std::enable_if_t<std::is_convertible<hila::type_div<T, S>, T>::value, int> = 0>
570 Array<n, m, T> &operator/=(const S rhs) & {
571 for (int i = 0; i < n * m; i++) {
572 c[i] /= rhs;
573 }
574 return *this;
575 }
576
577 /**
578 * @brief Returns element wise Complex conjugate of Array
579 *
580 * @return Array<n, m, T>
581 */
582 inline Array<n, m, T> conj() const {
583 Array<n, m, T> res;
584 for (int i = 0; i < n * m; i++) {
585 res.c[i] = ::conj(c[i]);
586 }
587 return res;
588 }
589 /**
590 * @brief Returns real part of Array
591 *
592 * @return Array<n, m, T>
593 */
596 for (int i = 0; i < m * n; i++) {
597 res.c[i] = ::real(c[i]);
598 }
599 return res;
600 }
601
602 /// return imaginary part
605 for (int i = 0; i < m * n; i++) {
606 res.c[i] = ::imag(c[i]);
607 }
608 return res;
609 }
610
611 /// calculate square norm - sum of squared elements
612 hila::arithmetic_type<T> squarenorm() const {
613 hila::arithmetic_type<T> result = 0;
614 for (int i = 0; i < n * m; i++) {
615 result += ::squarenorm(c[i]);
616 }
617 return result;
618 }
619
620 /**
621 * @brief Fill Array with random elements
622 *
623 * @return Array<n, m, T>&
624 */
625 Array<n, m, T> &random() out_only {
626 for (int i = 0; i < n * m; i++) {
627 hila::random(c[i]);
628 }
629 return *this;
630 }
631
632 /**
633 * @brief Fill Array with Gaussian random elements
634 *
635 * @param width
636 * @return Array<n, m, T>&
637 */
638 Array<n, m, T> &gaussian_random(double width = 1.0) out_only {
639 for (int i = 0; i < n * m; i++) {
640 hila::gaussian_random(c[i], width);
641 }
642 return *this;
643 }
644
645 /// Convert to string for printing
646 std::string str(int prec = 8, char separator = ' ') const {
647 return this->asMatrix().str(prec, separator);
648 }
649
650 /// implement sort as casting to array
651#pragma hila novector
652 template <int N>
654 hila::sort order = hila::sort::ascending) const {
655 return this->asMatrix().sort(permutation, order).asArray();
656 }
657
658#pragma hila novector
659 Array<n, m, T> sort(hila::sort order = hila::sort::ascending) const {
660 return this->asMatrix().sort(order).asArray();
661 }
662};
663
664/**
665 * @brief Return conjugate Array
666 * @details Wrapper around Array::conj
667 *
668 * @tparam n Number of rows
669 * @tparam m Number of columns
670 * @tparam T Array element type
671 * @param arg Array to be conjugated
672 * @return Array<n, m, T>
673 */
674template <const int n, const int m, typename T>
676 return arg.conj();
677}
678/**
679 * @brief Return real part of Array
680 * @details Wrapper around Array::real
681 *
682 * @tparam n Number of rows
683 * @tparam m Number of columns
684 * @tparam T Array element type
685 * @param arg Array to return real part of
686 * @return Array<n, m, T>
687 */
688template <const int n, const int m, typename T>
690 return arg.real();
691}
692/**
693 * @brief Return imaginary part of Array
694 * @details Wrapper around Array::imag
695 *
696 * @tparam n Number of rows
697 * @tparam m Number of columns
698 * @tparam T Array element type
699 * @param arg Array to return real part of
700 * @return Array<n, m, T>
701 */
702template <const int n, const int m, typename T>
704 return arg.imag();
705}
706
707/**
708 * @brief Addition operator
709 * @details Defined for the following operations
710 *
711 * __Array + Array:__
712 *
713 * __NOTE__: Arrays must share same dimensionality
714 *
715 * \code {.cpp}
716 * Array<n,m,MyType> A, B, C;
717 * A = 1;
718 * B = 1;
719 * C = A + B; // C is uniformly 2
720 * \endcode
721 *
722 *
723 * __Scalar + Array / Array + Scalar:__
724 *
725 * __NOTE__: Exact definition exist in overloaded functions that can be viewed in source code.
726 * *
727 * \code {.cpp}
728 * Array<n,m,MyType> A,B;
729 * A = 1;
730 * B = A + 1; // B is uniformly 2
731 * \endcode
732 * @memberof Array
733 * @tparam n Number of rows
734 * @tparam m Number of columns
735 * @tparam T Array element type
736 * @param a
737 * @param b
738 * @return Array<n, m, T>
739 */
740template <int n, int m, typename A, typename B>
741inline auto operator+(const Array<n, m, A> &a, const Array<n, m, B> &b) {
743 for (int i = 0; i < n * m; i++)
744 res.c[i] = a.c[i] * b.c[i];
745 return res;
746}
747
748/**
749 * @brief Subtraction operator
750 * @details Defined for the following operations
751 *
752 * __Array - Array:__
753 *
754 * __NOTE__: Arrays must share same dimensionality
755 *
756 * \code {.cpp}
757 * Array<n,m,MyType> A, B, C;
758 * A = 2;
759 * B = 1;
760 * C = A - B; // C is uniformly 2
761 * \endcode
762 *
763 *
764 * __Scalar - Array / Array - Scalar:__
765 *
766 * __NOTE__: Exact definition exist in overloaded functions that can be viewed in source code.
767 * *
768 * \code {.cpp}
769 * Array<n,m,MyType> A,B;
770 * A = 2;
771 * B = A - 1; // B is uniformly 2
772 * \endcode
773 * @memberof Array
774 * @tparam n Number of rows
775 * @tparam m Number of columns
776 * @tparam T Array element type
777 * @param a
778 * @param b
779 * @return Array<n, m, T>
780 */
781template <int n, int m, typename A, typename B>
782inline auto operator-(const Array<n, m, A> &a, const Array<n, m, B> &b) {
784 for (int i = 0; i < n * m; i++)
785 res.c[i] = a.c[i] - b.c[i];
786 return res;
787}
788
789// Array + scalar
790template <int n, int m, typename T, typename S,
791 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value, int> = 0>
792inline auto operator+(const Array<n, m, T> &a, const S b) {
794 for (int i = 0; i < n * m; i++)
795 res.c[i] = a.c[i] + b;
796 return res;
797}
798
799// scalar + Array
800template <int n, int m, typename T, typename S,
801 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value, int> = 0>
802inline auto operator+(const S b, const Array<n, m, T> &a) {
803 return a + b;
804}
805
806// Array - scalar
807template <int n, int m, typename T, typename S,
808 std::enable_if_t<std::is_convertible<hila::type_minus<T, S>, T>::value, int> = 0>
809inline auto operator-(const Array<n, m, T> &a, const S b) {
811 for (int i = 0; i < n * m; i++)
812 res.c[i] = a.c[i] - b;
813 return res;
814}
815
816// scalar - Array
817template <int n, int m, typename T, typename S,
818 std::enable_if_t<std::is_convertible<hila::type_minus<T, S>, T>::value, int> = 0>
819inline auto operator-(const S b, const Array<n, m, T> &a) {
821 for (int i = 0; i < n * m; i++)
822 res.c[i] = b - a.c[i];
823 return res;
824}
825
826
827/**
828 * @brief Multiplication operator
829 * @details Defined for the following operations
830 *
831 * __Array*Array:__
832 *
833 * __NOTE__: Arrays must share same dimensionality
834 *
835 * \code {.cpp}
836 * Array<n,m,MyType> A, B, C;
837 * A = 2;
838 * B = 3;
839 * C = A*B; // C is uniformly 6
840 * \endcode
841 *
842 *
843 * __Scalar * Array / Array * Scalar:__
844 *
845 * __NOTE__: Exact definition exist in overloaded functions that can be viewed in source code.
846 * *
847 * \code {.cpp}
848 * Array<n,m,MyType> A,B;
849 * A = 2;
850 * B = A*3; // B is uniformly 6
851 * \endcode
852 * @memberof Array
853 * @tparam n Number of rows
854 * @tparam m Number of columns
855 * @tparam T Array element type
856 * @param a
857 * @param b
858 * @return Array<n, m, T>
859 */
860
861template <int n, int m, typename T>
863 a *= b;
864 return a;
865}
866
867// mult with different type arrays
868
869template <int n, int m, typename A, typename B,
870 std::enable_if_t<!std::is_same<A, B>::value, int> = 0>
871inline auto operator*(const Array<n, m, A> &a, const Array<n, m, B> &b) {
872
873 using R = hila::type_mul<A, B>;
874
875 Array<n, m, R> res;
876 for (int i = 0; i < n; i++)
877 for (int j = 0; j < m; j++)
878 res.e(i, j) = a.e(i, j) * b.e(i, j);
879
880 return res;
881}
882
883/**
884 * @brief Division operator
885 * @details Defined for the following operations
886 *
887 * __Array/Array:__
888 *
889 * __NOTE__: Arrays must share same dimensionality
890 *
891 * \code {.cpp}
892 * Array<n,m,MyType> A, B, C;
893 * A = 4;
894 * B = 2;
895 * C = A/B; // C is uniformly 2
896 * \endcode
897 *
898 *
899 * __Scalar / Array / Array / Scalar:__
900 *
901 * __NOTE__: Exact definition exist in overloaded functions that can be viewed in source code.
902 * *
903 * \code {.cpp}
904 * Array<n,m,MyType> A,B;
905 * A = 4;
906 * B = A/2; // B is uniformly 2
907 * \endcode
908 * @memberof Array
909 * @tparam n Number of rows
910 * @tparam m Number of columns
911 * @tparam T Array element type
912 * @param a
913 * @param b
914 * @return Array<n, m, T>
915 */
916template <int n, int m, typename A, typename B>
917inline auto operator/(const Array<n, m, A> &a, const Array<n, m, B> &b) {
919 for (int i = 0; i < n; i++)
920 for (int j = 0; j < m; j++)
921 res.e(i, j) = a.e(i, j) / b.e(i, j);
922 return res;
923}
924
925
926// Array * scalar
927template <int n, int m, typename T, typename S,
928 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value, int> = 0>
929inline auto operator*(const Array<n, m, T> &a, const S b) {
931 for (int i = 0; i < n * m; i++)
932 res.c[i] = a.c[i] * b;
933 return res;
934}
935
936// scalar * Array
937template <int n, int m, typename T, typename S,
938 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value, int> = 0>
939inline auto operator*(const S b, const Array<n, m, T> &a) {
940 return a * b;
941}
942
943
944// Array / scalar
945template <int n, int m, typename T, typename S,
946 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value, int> = 0>
947inline auto operator/(const Array<n, m, T> &a, const S b) {
949 for (int i = 0; i < n * m; i++)
950 res.c[i] = a.c[i] / b;
951 return res;
952}
953
954// scalar / Array
955template <int n, int m, typename T, typename S,
956 std::enable_if_t<hila::is_complex_or_arithmetic<S>::value, int> = 0>
957inline auto operator/(const S b, const Array<n, m, T> &a) {
959 for (int i = 0; i < n * m; i++)
960 res.c[i] = b / a.c[i];
961 return res;
962}
963
964
965/**
966 * @brief Stream operator
967 * @details Naive Stream operator for formatting Matrix for printing. Simply puts elements one after
968 * the other in row major order
969 *
970 * @tparam n Number of rows
971 * @tparam m Number of columns
972 * @tparam T Array element type
973 * @param strm OS stream
974 * @param A Array to write
975 * @return std::ostream&
976 */
977template <int n, int m, typename T>
978std::ostream &operator<<(std::ostream &strm, const Array<n, m, T> &A) {
979 return operator<<(strm, A.asMatrix());
980}
981
982namespace hila {
983
984/**
985 * @brief Converts Array object to string
986 *
987 * @tparam n Number of rows
988 * @tparam m Number of columns
989 * @tparam T Array element type
990 * @param A Array to convert to string
991 * @param prec Precision of T
992 * @param separator Separator between elements
993 * @return std::string
994 */
995template <int n, int m, typename T>
996std::string to_string(const Array<n, m, T> &A, int prec = 8, char separator = ' ') {
997 return to_string(A.asMatrix(), prec, separator);
998}
999
1000template <int n, int m, typename T>
1001std::string prettyprint(const Array<n, m, T> &A, int prec = 8) {
1002 return prettyprint(A.asMatrix(), prec);
1003}
1004
1005} // namespace hila
1006
1007
1008/**
1009 * @brief Return square norm of Array
1010 * @details Wrapper around Array::squarenrom
1011 *
1012 * @tparam n Number of rows
1013 * @tparam m Number of columns
1014 * @tparam T Array element type
1015 * @param rhs Array to compute squarenorm of
1016 * @return hila::arithmetic_type<T>
1017 */
1018template <int n, int m, typename T>
1019inline hila::arithmetic_type<T> squarenorm(const Array<n, m, T> &rhs) {
1020 return rhs.squarenorm();
1021}
1022
1023/** @name Arithmetic operations
1024 * @details all operations are applied linearly to the Array A
1025 *
1026 * Most functions are self explanitory, but if necessary function will have a detailed section with
1027 * additional information.
1028 * @memberof Array
1029 * @tparam n Number of rows
1030 * @tparam m Number of columns
1031 * @tparam T Array element type
1032 * @param a Input array
1033 * @return Array<n, m, T>
1034 * @{
1035 */
1036
1037/**
1038 * @brief Square root
1039 */
1040template <int n, int m, typename T>
1042 for (int i = 0; i < n * m; i++)
1043 a.c[i] = sqrt(a.c[i]);
1044 return a;
1045}
1046
1047/**
1048 * @brief Cuberoot
1049 */
1050template <int n, int m, typename T>
1052 for (int i = 0; i < n * m; i++)
1053 a.c[i] = cbrt(a.c[i]);
1054 return a;
1055}
1056
1057/**
1058 * @brief Exponential
1059 */
1060template <int n, int m, typename T>
1062 for (int i = 0; i < n * m; i++)
1063 a.c[i] = exp(a.c[i]);
1064 return a;
1065}
1066
1067/**
1068 * @brief Logarithm
1069 */
1070template <int n, int m, typename T>
1072 for (int i = 0; i < n * m; i++)
1073 a.c[i] = log(a.c[i]);
1074 return a;
1075}
1076
1077/**
1078 * @brief Sine
1079 */
1080template <int n, int m, typename T>
1082 for (int i = 0; i < n * m; i++)
1083 a.c[i] = sin(a.c[i]);
1084 return a;
1085}
1086
1087/**
1088 * @brief Cosine
1089 */
1090template <int n, int m, typename T>
1092 for (int i = 0; i < n * m; i++)
1093 a.c[i] = cos(a.c[i]);
1094 return a;
1095}
1096
1097/**
1098 * @brief Tangent
1099 */
1100template <int n, int m, typename T>
1102 for (int i = 0; i < n * m; i++)
1103 a.c[i] = tan(a.c[i]);
1104 return a;
1105}
1106
1107/**
1108 * @brief Inverse Sine
1109 */
1110template <int n, int m, typename T>
1112 for (int i = 0; i < n * m; i++)
1113 a.c[i] = asin(a.c[i]);
1114 return a;
1115}
1116
1117/**
1118 * @brief Inverse Cosine
1119 */
1120template <int n, int m, typename T>
1122 for (int i = 0; i < n * m; i++)
1123 a.c[i] = acos(a.c[i]);
1124 return a;
1125}
1126
1127/**
1128 * @brief Inverse Tangent
1129 */
1130template <int n, int m, typename T>
1132 for (int i = 0; i < n * m; i++)
1133 a.c[i] = atan(a.c[i]);
1134 return a;
1135}
1136
1137/**
1138 * @brief Hyperbolic Sine
1139 */
1140template <int n, int m, typename T>
1142 for (int i = 0; i < n * m; i++)
1143 a.c[i] = sinh(a.c[i]);
1144 return a;
1145}
1146
1147/**
1148 * @brief Hyperbolic Cosine
1149 */
1150template <int n, int m, typename T>
1152 for (int i = 0; i < n * m; i++)
1153 a.c[i] = cosh(a.c[i]);
1154 return a;
1155}
1156
1157/**
1158 * @brief Hyperbolic tangent
1159 */
1160template <int n, int m, typename T>
1162 for (int i = 0; i < n * m; i++)
1163 a.c[i] = tanh(a.c[i]);
1164 return a;
1165}
1166
1167/**
1168 * @brief Inverse Hyperbolic Sine
1169 */
1170template <int n, int m, typename T>
1172 for (int i = 0; i < n * m; i++)
1173 a.c[i] = asinh(a.c[i]);
1174 return a;
1175}
1176
1177/**
1178 * @brief Inverse Hyperbolic Cosine
1179 */
1180template <int n, int m, typename T>
1182 for (int i = 0; i < n * m; i++)
1183 a.c[i] = acosh(a.c[i]);
1184 return a;
1185}
1186
1187/**
1188 * @brief Inverse Hyperbolic Tangent
1189 */
1190template <int n, int m, typename T>
1192 for (int i = 0; i < n * m; i++)
1193 a.c[i] = atanh(a.c[i]);
1194 return a;
1195}
1196
1197/**
1198 * @brief Power
1199 *
1200 * @details Array can be raised homogeneously to the power of a integer or scalar b or Array can be
1201 * raised to the power of another Array b element wise.
1202 *
1203 * @param b Array, Integer or Real scalar to raise to the power of
1204 */
1205template <int n, int m, typename T>
1207 for (int i = 0; i < n * m; i++)
1208 a.c[i] = pow(a.c[i], b);
1209 return a;
1210}
1211
1212template <int n, int m, typename T>
1213inline Array<n, m, T> pow(Array<n, m, T> a, T b) {
1214 for (int i = 0; i < n * m; i++)
1215 a.c[i] = pow(a.c[i], b);
1216 return a;
1217}
1218
1219template <int n, int m, typename T>
1220inline Array<n, m, T> pow(Array<n, m, T> a, const Array<n, m, T> &b) {
1221 for (int i = 0; i < n * m; i++)
1222 a.c[i] = pow(a.c[i], b.c[i]);
1223 return a;
1224}
1225
1226/**
1227 * @brief Rounding
1228 * @details Rounding function if T is arithmetic type
1229 * @tparam T Array element type, must be arithmetic type
1230 */
1231template <int n, int m, typename T, std::enable_if_t<hila::is_arithmetic<T>::value, int> = 0>
1233 for (int i = 0; i < n * m; i++)
1234 a.c[i] = round(a.c[i]);
1235 return a;
1236}
1237
1238/**
1239 * @brief Floor
1240 * @details Flooring function if T is arithmetic type
1241 * @tparam T Array element type, must be arithmetic type
1242 */
1243template <int n, int m, typename T, std::enable_if_t<hila::is_arithmetic<T>::value, int> = 0>
1245 for (int i = 0; i < n * m; i++)
1246 a.c[i] = floor(a.c[i]);
1247 return a;
1248}
1249
1250/**
1251 * @brief Ceiling
1252 * @details Ceiling function if T is arithmetic type
1253 * @tparam T Array element type, must be arithmetic type
1254 */
1255template <int n, int m, typename T, std::enable_if_t<hila::is_arithmetic<T>::value, int> = 0>
1257 for (int i = 0; i < n * m; i++)
1258 a.c[i] = ceil(a.c[i]);
1259 return a;
1260}
1261
1262/**
1263 * @brief Truncation
1264 * @details Truncation function if T is arithmetic type
1265 * @tparam T Array element type, must be arithmetic type
1266 */
1267template <int n, int m, typename T, std::enable_if_t<hila::is_arithmetic<T>::value, int> = 0>
1269 for (int i = 0; i < n * m; i++)
1270 a.c[i] = trunc(a.c[i]);
1271 return a;
1272}
1273/** @} */
1274
1275namespace hila {
1276/**
1277 * @brief Array casting operation
1278 * @details Cast array to different number type or Complex type
1279 *
1280 * Allowed casting: number->number, number->Complex, Complex->Complex
1281 *
1282 * Not allowed casting: Complex->number
1283 * @tparam Ntype Type to cast to
1284 * @tparam T Input Array type
1285 * @tparam n Number of rows
1286 * @tparam m Number of columns
1287 * @param mat Array to cast
1288 * @return Array<n, m, Ntype>
1289 */
1290template <typename Ntype, typename T, int n, int m,
1291 std::enable_if_t<hila::is_arithmetic_or_extended<T>::value, int> = 0>
1294 for (int i = 0; i < n * m; i++)
1295 res.c[i] = mat.c[i];
1296 return res;
1297}
1298
1299template <typename Ntype, typename T, int n, int m,
1300 std::enable_if_t<hila::is_complex<T>::value, int> = 0>
1301auto cast_to(const Array<n, m, T> &mat) {
1303 for (int i = 0; i < n * m; i++)
1304 res.c[i] = cast_to<Ntype>(mat.c[i]);
1305 return res;
1306}
1307
1308} // namespace hila
1309
1310/// @brief Array1d is alias to Array where m = 1
1311template <int n, typename T = double>
1313
1314/// @brief Array2d is simply alias of Array
1315template <int n, int m, typename T = double>
1317
1318#endif
Array< n, m, T > conj(const Array< n, m, T > &arg)
Return conjugate Array.
Definition array.h:675
std::ostream & operator<<(std::ostream &strm, const Array< n, m, T > &A)
Stream operator.
Definition array.h:978
Array< n, m, hila::arithmetic_type< T > > imag(const Array< n, m, T > &arg)
Return imaginary part of Array.
Definition array.h:703
hila::arithmetic_type< T > squarenorm(const Array< n, m, T > &rhs)
Return square norm of Array.
Definition array.h:1019
Array< n, m, hila::arithmetic_type< T > > real(const Array< n, m, T > &arg)
Return real part of Array.
Definition array.h:689
Array type
Definition array.h:43
Array(const S rhs)
Scalar constructor.
Definition array.h:104
Array< n, m, T > & operator+=(const Array< n, m, S > &rhs) &
Add assign operator with Array or scalar.
Definition array.h:423
Array(std::initializer_list< S > rhs)
Initializer list constructor.
Definition array.h:132
auto operator/(const Array< n, m, A > &a, const Array< n, m, B > &b)
Division operator.
Definition array.h:917
Array< n, m, T > floor(Array< n, m, T > a)
Floor.
Definition array.h:1244
constexpr int rows() const
Returns number of rows.
Definition array.h:145
constexpr int size() const
Returns size of Vector Array or square Array.
Definition array.h:165
T & e(const int i, const int j)
Const overload.
Definition array.h:201
Array< n, m, T > & operator=(const S rhs) &
Scalar assignment operator.
Definition array.h:334
std::string str(int prec=8, char separator=' ') const
Convert to string for printing.
Definition array.h:646
Array< n, m, T > round(Array< n, m, T > a)
Rounding.
Definition array.h:1232
Array< n, m, T > operator+() const
Unary + operator.
Definition array.h:301
Array< n, m, T > asin(Array< n, m, T > a)
Inverse Sine.
Definition array.h:1111
auto operator-(const Array< n, m, A > &a, const Array< n, m, B > &b)
Subtraction operator.
Definition array.h:782
bool operator==(const Array< n1, m1, S > &rhs) const
Compare equality of Arrays.
Definition array.h:363
Array< n, m, T > & operator/=(const Array< n, m, S > &rhs) &
Division assign with array or scalar.
Definition array.h:559
Array< n, m, T > exp(Array< n, m, T > a)
Exponential.
Definition array.h:1061
T e(const int i) const
Const overload.
Definition array.h:208
Array< n, m, T > tanh(Array< n, m, T > a)
Hyperbolic tangent.
Definition array.h:1161
Array< n, m, T > trunc(Array< n, m, T > a)
Truncation.
Definition array.h:1268
T & e(const int i)
Standard array indexing operation for 1D Array.
Definition array.h:223
Array< n, m, T > sort(Vector< N, int > &permutation, hila::sort order=hila::sort::ascending) const
implement sort as casting to array
Definition array.h:653
Array< n, m, hila::arithmetic_type< T > > real() const
Returns real part of Array.
Definition array.h:594
Array< n, m, T > cosh(Array< n, m, T > a)
Hyperbolic Cosine.
Definition array.h:1151
Array< n, m, T > tan(Array< n, m, T > a)
Tangent.
Definition array.h:1101
auto operator+(const Array< n, m, A > &a, const Array< n, m, B > &b)
Addition operator.
Definition array.h:741
Array()=default
Default constructor.
Array< n, m, T > & operator*=(const S rhs) &
multiply assign with scalar
Definition array.h:524
Array< n, m, hila::arithmetic_type< T > > imag() const
return imaginary part
Definition array.h:603
Array< n, m, T > operator-() const
Unary - operator.
Definition array.h:287
Array< n, m, T > & operator-=(const Array< n, m, S > &rhs) &
Subtract assign operator with Array or scalar.
Definition array.h:457
constexpr int columns() const
Returns number of columns.
Definition array.h:153
T e(const int i, const int j) const
Standard array indexing operation for 2D Array.
Definition array.h:195
Array< n, m, T > cos(Array< n, m, T > a)
Cosine.
Definition array.h:1091
bool operator!=(const Array< n1, m1, S > &rhs) const
Compare non-equality of two Arrays.
Definition array.h:389
Array< n, m, T > cbrt(Array< n, m, T > a)
Cuberoot.
Definition array.h:1051
hila::arithmetic_type< T > squarenorm() const
calculate square norm - sum of squared elements
Definition array.h:612
Vector< n, T > & asVector()
Cast Array1D to Vector.
Definition array.h:247
Array< n, m, T > acos(Array< n, m, T > a)
Inverse Cosine.
Definition array.h:1121
Array< n, m, T > acosh(Array< n, m, T > a)
Inverse Hyperbolic Cosine.
Definition array.h:1181
Array< n, m, T > sqrt(Array< n, m, T > a)
Square root.
Definition array.h:1041
Array< n, m, T > conj() const
Returns element wise Complex conjugate of Array.
Definition array.h:582
Array< n, m, T > asinh(Array< n, m, T > a)
Inverse Hyperbolic Sine.
Definition array.h:1171
Matrix< n, m, T > & asMatrix()
Cast Array to Matrix.
Definition array.h:232
Array< n, m, T > & random()
Fill Array with random elements.
Definition array.h:625
Array< n, m, T > ceil(Array< n, m, T > a)
Ceiling.
Definition array.h:1256
Array< n, m, T > atan(Array< n, m, T > a)
Inverse Tangent.
Definition array.h:1131
Array< n, m, T > pow(Array< n, m, T > a, int b)
Power.
Definition array.h:1206
Array< n, m, T > & operator*=(const Array< n, m, S > &rhs) &
Multiply assign scalar or array.
Definition array.h:513
Array(const Array< n, m, T > &v)=default
Copy constructor.
Array< n, m, T > operator*(Array< n, m, T > a, const Array< n, m, T > &b)
Multiplication operator.
Definition array.h:862
Array< n, m, T > sin(Array< n, m, T > a)
Sine.
Definition array.h:1081
Array< n, m, T > sinh(Array< n, m, T > a)
Hyperbolic Sine.
Definition array.h:1141
Array< n, m, T > & gaussian_random(double width=1.0)
Fill Array with Gaussian random elements.
Definition array.h:638
Array< n, m, T > atanh(Array< n, m, T > a)
Inverse Hyperbolic Tangent.
Definition array.h:1191
Array< n, m, T > log(Array< n, m, T > a)
Logarithm.
Definition array.h:1071
Define type DiagonalMatrix<n,T>
Matrix class which defines matrix operations.
Definition matrix.h:1742
T arg(const Complex< T > &a)
Return argument of Complex number.
Definition cmplx.h:1278
Definition of Matrix types.
Implement hila::swap for gauge fields.
Definition array.h:982
T gaussian_random()
Template function T hila::gaussian_random<T>(),generates gaussian random value of type T,...
Definition random.h:179
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
std::string to_string(const Array< n, m, T > &A, int prec=8, char separator=' ')
Converts Array object to string.
Definition array.h:996
hila::is_complex_or_arithmetic<T>::value
Definition cmplx.h:708