Block-Structured AMR Software Framework
AMReX_SmallMatrix.H
Go to the documentation of this file.
1 #ifndef AMREX_SMALL_MATRIX_H_
2 #define AMREX_SMALL_MATRIX_H_
3 #include <AMReX_Config.H>
4 
5 #include <AMReX_Algorithm.H>
6 #include <AMReX_BLassert.H>
7 #include <AMReX_Extension.H>
8 #include <AMReX_GpuQualifiers.H>
9 #include <AMReX_ConstexprFor.H>
10 
11 #include <algorithm>
12 #include <initializer_list>
13 #include <iostream>
14 #include <tuple>
15 #include <type_traits>
16 
17 namespace amrex {
18 
19  enum struct Order { C, F, RowMajor=C, ColumnMajor=F };
20 
33  template <class T, int NRows, int NCols, Order ORDER = Order::F, int StartIndex = 0>
34  struct SmallMatrix
35  {
36  using value_type = T;
37  using reference_type = T&;
38  static constexpr int row_size = NRows;
39  static constexpr int column_size = NCols;
40  static constexpr Order ordering = ORDER;
41  static constexpr int starting_index = StartIndex;
42 
50  constexpr SmallMatrix () = default;
51 
59  template <typename... Ts, int MM=NRows, int NN=NCols,
60  std::enable_if_t<MM==1 || NN==1, int> = 0>
62  constexpr explicit SmallMatrix (Ts... vs)
63  : m_mat{vs...}
64  {
65  static_assert(sizeof...(vs) <= std::max(NRows,NCols));
66  }
67 
80  explicit SmallMatrix (std::initializer_list<std::initializer_list<T>> const& init)
81  {
82  AMREX_ASSERT(NRows == init.size());
83  int i = StartIndex;
84  for (auto const& row : init) {
85  AMREX_ASSERT(NCols == row.size());
86  int j = StartIndex;
87  for (auto const& x : row) {
88  (*this)(i,j) = x;
89  ++j;
90  }
91  ++i;
92  }
93  }
94 
97  const T& operator() (int i, int j) const noexcept {
98  static_assert(StartIndex == 0 || StartIndex == 1);
99  if constexpr (StartIndex == 1) {
100  --i;
101  --j;
102  }
103  AMREX_ASSERT(i < NRows && j < NCols);
104  if constexpr (ORDER == Order::F) {
105  return m_mat[i+j*NRows];
106  } else {
107  return m_mat[j+i*NCols];
108  }
109  }
110 
113  T& operator() (int i, int j) noexcept {
114  static_assert(StartIndex == 0 || StartIndex == 1);
115  if constexpr (StartIndex == 1) {
116  --i;
117  --j;
118  }
119  AMREX_ASSERT(i < NRows && j < NCols);
120  if constexpr (ORDER == Order::F) {
121  return m_mat[i+j*NRows];
122  } else {
123  return m_mat[j+i*NCols];
124  }
125  }
126 
128  template <int MM=NRows, int NN=NCols, std::enable_if_t<(MM==1 || NN==1), int> = 0>
130  const T& operator() (int i) const noexcept {
131  static_assert(StartIndex == 0 || StartIndex == 1);
132  if constexpr (StartIndex == 1) {
133  --i;
134  }
135  AMREX_ASSERT(i < NRows*NCols);
136  return m_mat[i];
137  }
138 
140  template <int MM=NRows, int NN=NCols, std::enable_if_t<(MM==1 || NN==1), int> = 0>
142  T& operator() (int i) noexcept {
143  static_assert(StartIndex == 0 || StartIndex == 1);
144  if constexpr (StartIndex == 1) {
145  --i;
146  }
147  AMREX_ASSERT(i < NRows*NCols);
148  return m_mat[i];
149  }
150 
152  template <int MM=NRows, int NN=NCols, std::enable_if_t<(MM==1 || NN==1), int> = 0>
154  const T& operator[] (int i) const noexcept {
155  static_assert(StartIndex == 0 || StartIndex == 1);
156  if constexpr (StartIndex == 1) {
157  --i;
158  }
159  AMREX_ASSERT(i < NRows*NCols);
160  return m_mat[i];
161  }
162 
164  template <int MM=NRows, int NN=NCols, std::enable_if_t<(MM==1 || NN==1), int> = 0>
166  T& operator[] (int i) noexcept {
167  static_assert(StartIndex == 0 || StartIndex == 1);
168  if constexpr (StartIndex == 1) {
169  --i;
170  }
171  AMREX_ASSERT(i < NRows*NCols);
172  return m_mat[i];
173  }
174 
180  const T* begin () const noexcept { return m_mat; }
181 
187  const T* end () const noexcept { return m_mat + NRows*NCols; }
188 
194  T* begin () noexcept { return m_mat; }
195 
201  T* end () noexcept { return m_mat + NRows*NCols; }
202 
206  setVal (T val)
207  {
208  for (auto& x : m_mat) { x = val; }
209  return *this;
210  }
211 
213  template <int MM=NRows, int NN=NCols, std::enable_if_t<MM==NN, int> = 0>
214  static constexpr
217  Identity () noexcept {
218  static_assert(StartIndex == 0 || StartIndex == 1);
220  constexpr_for<StartIndex,NRows+StartIndex>(
221  [&] (int i) { I(i,i) = T(1); });
222  return I;
223  }
224 
226  static constexpr
229  Zero () noexcept {
231  return Z;
232  }
233 
237  transpose () const
238  {
240  for (int j = StartIndex; j < NRows+StartIndex; ++j) {
241  for (int i = StartIndex; i < NCols+StartIndex; ++i) {
242  r(i,j) = (*this)(j,i);
243  }
244  }
245  return r;
246  }
247 
249  template <int MM=NRows, int NN=NCols, std::enable_if_t<MM==NN,int> = 0>
253  {
254  static_assert(StartIndex == 0 || StartIndex == 1);
255  for (int j = 1+StartIndex; j < NCols+StartIndex; ++j) {
256  for (int i = StartIndex; i < j; ++i) {
257  amrex::Swap((*this)(i,j), (*this)(j,i));
258  }
259  }
260  return *this;
261  }
262 
265  T product () const
266  {
267  T p = 1;
268  for (auto const& x : m_mat) {
269  p *= x;
270  }
271  return p;
272  }
273 
276  T sum () const
277  {
278  T s = 0;
279  for (auto const& x : m_mat) {
280  s += x;
281  }
282  return s;
283  }
284 
286  template <int MM=NRows, int NN=NCols, std::enable_if_t<MM==NN,int> = 0>
288  T trace () const
289  {
290  T t = 0;
291  constexpr_for<StartIndex,MM+StartIndex>([&] (int i) { t += (*this)(i,i); });
292  return t;
293  }
294 
299  {
300  for (int n = 0; n < NRows*NCols; ++n) {
301  m_mat[n] += rhs.m_mat[n];
302  }
303  return *this;
304  }
305 
311  {
312  lhs += rhs;
313  return lhs;
314  }
315 
320  {
321  for (int n = 0; n < NRows*NCols; ++n) {
322  m_mat[n] -= rhs.m_mat[n];
323  }
324  return *this;
325  }
326 
332  {
333  lhs -= rhs;
334  return lhs;
335  }
336 
340  operator- () const
341  {
342  return (*this) * T(-1);
343  }
344 
349  {
350  for (auto& x : m_mat) {
351  x *= a;
352  }
353  return *this;
354  }
355 
360  {
361  m *= a;
362  return m;
363  }
364 
369  {
370  m *= a;
371  return m;
372  }
373 
375  template <class U, int N1, int N2, int N3, Order Ord, int SI>
379  SmallMatrix<U,N2,N3,Ord,SI> const& rhs);
380 
384  {
385  T r = 0;
386  for (int n = 0; n < NRows*NCols; ++n) {
387  r += m_mat[n] * rhs.m_mat[n];
388  }
389  return r;
390  }
391 
392  template <int N, std::enable_if_t<(N<NRows*NCols),int> = 0>
393  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
394  constexpr T const& get () const { return m_mat[N]; }
395 
396  template <int N, std::enable_if_t<(N<NRows*NCols),int> = 0>
397  [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
398  constexpr T& get () { return m_mat[N]; }
399 
400  private:
401  T m_mat[NRows*NCols];
402  };
403 
404  template <class U, int N1, int N2, int N3, Order Ord, int SI>
408  SmallMatrix<U,N2,N3,Ord,SI> const& rhs)
409  {
410  static_assert(SI == 0 || SI == 1);
412  if constexpr (Ord == Order::F) {
413  for (int j = SI; j < N3+SI; ++j) {
414  constexpr_for<SI,N1+SI>([&] (int i) { r(i,j) = U(0); });
415  for (int k = SI; k < N2+SI; ++k) {
416  auto b = rhs(k,j);
417  constexpr_for<SI,N1+SI>([&] (int i)
418  {
419  r(i,j) += lhs(i,k) * b;
420  });
421  }
422  }
423  } else {
424  for (int i = SI; i < N1+SI; ++i) {
425  constexpr_for<SI,N3+SI>([&] (int j) { r(i,j) = U(0); });
426  for (int k = SI; k < N2+SI; ++k) {
427  auto a = lhs(i,k);
428  constexpr_for<SI,N3+SI>([&] (int j)
429  {
430  r(i,j) += a * rhs(k,j);
431  });
432  }
433  }
434  }
435  return r;
436  }
437 
438  template <class T, int NRows, int NCols, Order ORDER, int SI>
439  std::ostream& operator<< (std::ostream& os,
441  {
442  for (int i = SI; i < NRows+SI; ++i) {
443  os << mat(i,SI);
444  for (int j = 1+SI; j < NCols+SI; ++j) {
445  os << " " << mat(i,j);
446  }
447  os << "\n";
448  }
449  return os;
450  }
451 
452  template <class T, int N, int StartIndex = 0>
454 
455  template <class T, int N, int StartIndex = 0>
457 }
458 
459 template <class T, int NRows, int NCols, amrex::Order ORDER, int StartIndex>
460 struct std::tuple_size<amrex::SmallMatrix<T,NRows,NCols,ORDER,StartIndex> >
461  : std::integral_constant<std::size_t,NRows*NCols> {};
462 
463 template <std::size_t N, class T, int NRows, int NCols, amrex::Order ORDER, int StartIndex>
464 struct std::tuple_element<N, amrex::SmallMatrix<T,NRows,NCols,ORDER,StartIndex> >
465 {
466  using type = T;
467 };
468 
469 #endif
470 
471 /*
472  * Notes on why SmallMatrix matrix{} is zero initialized.
473  *
474  * SmallMatrix is not an aggregate, because it has a user declared default
475  * constructor. The rule is that, for `SmallMatrix matrix{}` with an empty
476  * brace-enclosed initializer list, value-initialization is performed. The
477  * effects of value-initialization of SmallMatrix (which has a user-declared
478  * but not user-provided default constructor) are that the matrix object is
479  * first zero-initialized and then the object's default constructor is
480  * applied. Since the default constructor does nothing, the final result is
481  * the object is zero-initialized.
482  *
483  * Why is SmallMatrix's default constructor user-declared not user-provided?
484  * It's because we first declare it with `SmallMatrix () = default`.
485  *
486  * Reference:
487  * https://en.cppreference.com/w/cpp/language/list_initialization
488  * https://en.cppreference.com/w/cpp/language/value_initialization
489  * https://en.cppreference.com/w/cpp/language/zero_initialization
490  */
#define AMREX_ASSERT(EX)
Definition: AMReX_BLassert.H:38
#define AMREX_FORCE_INLINE
Definition: AMReX_Extension.H:119
#define AMREX_GPU_HOST_DEVICE
Definition: AMReX_GpuQualifiers.H:20
@ max
Definition: AMReX_ParallelReduce.H:17
Definition: AMReX_Amr.cpp:49
Order
Definition: AMReX_SmallMatrix.H:19
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Swap(T &t1, T &t2) noexcept
Definition: AMReX_Algorithm.H:75
std::ostream & operator<<(std::ostream &os, AmrMesh const &amr_mesh)
Definition: AMReX_AmrMesh.cpp:1236
Matrix class with compile-time size.
Definition: AMReX_SmallMatrix.H:35
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SmallMatrix< T, NCols, NRows, ORDER, StartIndex > transpose() const
Returns transposed matrix.
Definition: AMReX_SmallMatrix.H:237
static constexpr int row_size
Definition: AMReX_SmallMatrix.H:38
static constexpr AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SmallMatrix< T, NRows, NCols, ORDER, StartIndex > Identity() noexcept
Returns an identity matrix.
Definition: AMReX_SmallMatrix.H:217
T m_mat[NRows *NCols]
Definition: AMReX_SmallMatrix.H:401
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SmallMatrix< T, NRows, NCols, ORDER, StartIndex > & operator+=(SmallMatrix< T, NRows, NCols, ORDER, StartIndex > const &rhs)
Operator += performing matrix addition as in (*this) += rhs.
Definition: AMReX_SmallMatrix.H:298
static constexpr Order ordering
Definition: AMReX_SmallMatrix.H:40
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE friend SmallMatrix< T, NRows, NCols, ORDER, StartIndex > operator*(SmallMatrix< T, NRows, NCols, ORDER, StartIndex > m, T a)
Returns the product of a matrix and a scalar.
Definition: AMReX_SmallMatrix.H:359
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SmallMatrix< T, NRows, NCols, ORDER, StartIndex > operator-() const
Unary minus operator.
Definition: AMReX_SmallMatrix.H:340
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SmallMatrix< T, NRows, NCols, ORDER, StartIndex > & operator*=(T a)
Operator *= that scales this matrix in place by a scalar.
Definition: AMReX_SmallMatrix.H:348
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE SmallMatrix()=default
Default constructor.
static constexpr int starting_index
Definition: AMReX_SmallMatrix.H:41
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const T & operator()(int i, int j) const noexcept
Returns a const reference to the element at row i and column j.
Definition: AMReX_SmallMatrix.H:97
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SmallMatrix< T, NRows, NCols, ORDER, StartIndex > & setVal(T val)
Set all elements in the matrix to the given value.
Definition: AMReX_SmallMatrix.H:206
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T sum() const
Returns the sum of all elements in the matrix.
Definition: AMReX_SmallMatrix.H:276
static constexpr AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SmallMatrix< T, NRows, NCols, ORDER, StartIndex > Zero() noexcept
Returns a matrix initialized with zeros.
Definition: AMReX_SmallMatrix.H:229
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE SmallMatrix(Ts... vs)
Constructs column- or row-vector.
Definition: AMReX_SmallMatrix.H:62
T value_type
Definition: AMReX_SmallMatrix.H:36
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T dot(SmallMatrix< T, NRows, NCols, ORDER, StartIndex > const &rhs) const
Returns the dot product of two vectors.
Definition: AMReX_SmallMatrix.H:383
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE friend SmallMatrix< T, NRows, NCols, ORDER, StartIndex > operator+(SmallMatrix< T, NRows, NCols, ORDER, StartIndex > lhs, SmallMatrix< T, NRows, NCols, ORDER, StartIndex > const &rhs)
Binary operator + returning the result of maxtrix addition, lhs+rhs.
Definition: AMReX_SmallMatrix.H:309
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const T & operator[](int i) const noexcept
Returns a const reference to element i of a vector.
Definition: AMReX_SmallMatrix.H:154
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SmallMatrix(std::initializer_list< std::initializer_list< T >> const &init)
Constructs SmallMatrix with nested std::initializer_list.
Definition: AMReX_SmallMatrix.H:80
static constexpr int column_size
Definition: AMReX_SmallMatrix.H:39
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const T * begin() const noexcept
Definition: AMReX_SmallMatrix.H:180
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T * begin() noexcept
Definition: AMReX_SmallMatrix.H:194
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T product() const
Returns the product of all elements in the matrix.
Definition: AMReX_SmallMatrix.H:265
T & reference_type
Definition: AMReX_SmallMatrix.H:37
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T * end() noexcept
Definition: AMReX_SmallMatrix.H:201
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T trace() const
Returns the trace of a square matrix.
Definition: AMReX_SmallMatrix.H:288
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SmallMatrix< T, NRows, NCols, ORDER, StartIndex > & operator-=(SmallMatrix< T, NRows, NCols, ORDER, StartIndex > const &rhs)
Operator -= performing matrix subtraction as in (*this) -= rhs.
Definition: AMReX_SmallMatrix.H:319
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const T * end() const noexcept
Definition: AMReX_SmallMatrix.H:187
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SmallMatrix< T, NRows, NCols, ORDER, StartIndex > & transposeInPlace()
Transposes a square matrix in-place.
Definition: AMReX_SmallMatrix.H:252