Block-Structured AMR Software Framework
 
Loading...
Searching...
No Matches
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>
10
11#include <algorithm>
12#include <initializer_list>
13#include <iostream>
14#include <tuple>
15#include <type_traits>
16
17namespace 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>
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 >= 0 && j >= 0);
104 AMREX_ASSERT(i < NRows && j < NCols);
105 if constexpr (ORDER == Order::F) {
106 return m_mat[i+j*NRows];
107 } else {
108 return m_mat[j+i*NCols];
109 }
110 }
111
114 T& operator() (int i, int j) noexcept {
115 static_assert(StartIndex == 0 || StartIndex == 1);
116 if constexpr (StartIndex == 1) {
117 --i;
118 --j;
119 }
120 AMREX_ASSERT(i >= 0 && j >= 0);
121 AMREX_ASSERT(i < NRows && j < NCols);
122 if constexpr (ORDER == Order::F) {
123 return m_mat[i+j*NRows];
124 } else {
125 return m_mat[j+i*NCols];
126 }
127 }
128
130 template <int MM=NRows, int NN=NCols, std::enable_if_t<(MM==1 || NN==1), int> = 0>
132 const T& operator() (int i) const noexcept {
133 static_assert(StartIndex == 0 || StartIndex == 1);
134 if constexpr (StartIndex == 1) {
135 --i;
136 }
137 AMREX_ASSERT(i >= 0);
138 AMREX_ASSERT(i < NRows*NCols);
139 return m_mat[i];
140 }
141
143 template <int MM=NRows, int NN=NCols, std::enable_if_t<(MM==1 || NN==1), int> = 0>
145 T& operator() (int i) noexcept {
146 static_assert(StartIndex == 0 || StartIndex == 1);
147 if constexpr (StartIndex == 1) {
148 --i;
149 }
150 AMREX_ASSERT(i >= 0);
151 AMREX_ASSERT(i < NRows*NCols);
152 return m_mat[i];
153 }
154
156 template <int MM=NRows, int NN=NCols, std::enable_if_t<(MM==1 || NN==1), int> = 0>
158 const T& operator[] (int i) const noexcept {
159 static_assert(StartIndex == 0 || StartIndex == 1);
160 if constexpr (StartIndex == 1) {
161 --i;
162 }
163 AMREX_ASSERT(i >= 0);
164 AMREX_ASSERT(i < NRows*NCols);
165 return m_mat[i];
166 }
167
169 template <int MM=NRows, int NN=NCols, std::enable_if_t<(MM==1 || NN==1), int> = 0>
171 T& operator[] (int i) noexcept {
172 static_assert(StartIndex == 0 || StartIndex == 1);
173 if constexpr (StartIndex == 1) {
174 --i;
175 }
176 AMREX_ASSERT(i >= 0);
177 AMREX_ASSERT(i < NRows*NCols);
178 return m_mat[i];
179 }
180
186 const T* begin () const noexcept { return m_mat; }
187
193 const T* end () const noexcept { return m_mat + NRows*NCols; }
194
200 T* begin () noexcept { return m_mat; }
201
207 T* end () noexcept { return m_mat + NRows*NCols; }
208
212 setVal (T val)
213 {
214 for (auto& x : m_mat) { x = val; }
215 return *this;
216 }
217
219 template <int MM=NRows, int NN=NCols, std::enable_if_t<MM==NN, int> = 0>
220 static constexpr
223 Identity () noexcept {
224 static_assert(StartIndex == 0 || StartIndex == 1);
226 constexpr_for<StartIndex,NRows+StartIndex>(
227 [&] (int i) { I(i,i) = T(1); });
228 return I;
229 }
230
232 static constexpr
235 Zero () noexcept {
237 return Z;
238 }
239
243 transpose () const
244 {
246 for (int j = StartIndex; j < NRows+StartIndex; ++j) {
247 for (int i = StartIndex; i < NCols+StartIndex; ++i) {
248 r(i,j) = (*this)(j,i);
249 }
250 }
251 return r;
252 }
253
255 template <int MM=NRows, int NN=NCols, std::enable_if_t<MM==NN,int> = 0>
259 {
260 static_assert(StartIndex == 0 || StartIndex == 1);
261 for (int j = 1+StartIndex; j < NCols+StartIndex; ++j) {
262 for (int i = StartIndex; i < j; ++i) {
263 amrex::Swap((*this)(i,j), (*this)(j,i));
264 }
265 }
266 return *this;
267 }
268
271 T product () const
272 {
273 T p = 1;
274 for (auto const& x : m_mat) {
275 p *= x;
276 }
277 return p;
278 }
279
282 T sum () const
283 {
284 T s = 0;
285 for (auto const& x : m_mat) {
286 s += x;
287 }
288 return s;
289 }
290
292 template <int MM=NRows, int NN=NCols, std::enable_if_t<MM==NN,int> = 0>
294 T trace () const
295 {
296 T t = 0;
297 constexpr_for<StartIndex,MM+StartIndex>([&] (int i) { t += (*this)(i,i); });
298 return t;
299 }
300
305 {
306 for (int n = 0; n < NRows*NCols; ++n) {
307 m_mat[n] += rhs.m_mat[n];
308 }
309 return *this;
310 }
311
317 {
318 lhs += rhs;
319 return lhs;
320 }
321
326 {
327 for (int n = 0; n < NRows*NCols; ++n) {
328 m_mat[n] -= rhs.m_mat[n];
329 }
330 return *this;
331 }
332
338 {
339 lhs -= rhs;
340 return lhs;
341 }
342
347 {
348 return (*this) * T(-1);
349 }
350
355 {
356 for (auto& x : m_mat) {
357 x *= a;
358 }
359 return *this;
360 }
361
366 {
367 m *= a;
368 return m;
369 }
370
375 {
376 m *= a;
377 return m;
378 }
379
381 template <class U, int N1, int N2, int N3, Order Ord, int SI>
385 SmallMatrix<U,N2,N3,Ord,SI> const& rhs);
386
390 {
391 T r = 0;
392 for (int n = 0; n < NRows*NCols; ++n) {
393 r += m_mat[n] * rhs.m_mat[n];
394 }
395 return r;
396 }
397
398 template <int N, std::enable_if_t<(N<NRows*NCols),int> = 0>
399 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
400 constexpr T const& get () const { return m_mat[N]; }
401
402 template <int N, std::enable_if_t<(N<NRows*NCols),int> = 0>
403 [[nodiscard]] AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
404 constexpr T& get () { return m_mat[N]; }
405
406 private:
407 T m_mat[NRows*NCols];
408 };
409
410 template <class U, int N1, int N2, int N3, Order Ord, int SI>
415 {
416 static_assert(SI == 0 || SI == 1);
418 if constexpr (Ord == Order::F) {
419 for (int j = SI; j < N3+SI; ++j) {
420 constexpr_for<SI,N1+SI>([&] (int i) { r(i,j) = U(0); });
421 for (int k = SI; k < N2+SI; ++k) {
422 auto b = rhs(k,j);
423 constexpr_for<SI,N1+SI>([&] (int i)
424 {
425 r(i,j) += lhs(i,k) * b;
426 });
427 }
428 }
429 } else {
430 for (int i = SI; i < N1+SI; ++i) {
431 constexpr_for<SI,N3+SI>([&] (int j) { r(i,j) = U(0); });
432 for (int k = SI; k < N2+SI; ++k) {
433 auto a = lhs(i,k);
434 constexpr_for<SI,N3+SI>([&] (int j)
435 {
436 r(i,j) += a * rhs(k,j);
437 });
438 }
439 }
440 }
441 return r;
442 }
443
444 template <class T, int NRows, int NCols, Order ORDER, int SI>
445 std::ostream& operator<< (std::ostream& os,
447 {
448 for (int i = SI; i < NRows+SI; ++i) {
449 os << mat(i,SI);
450 for (int j = 1+SI; j < NCols+SI; ++j) {
451 os << " " << mat(i,j);
452 }
453 os << "\n";
454 }
455 return os;
456 }
457
458 template <class T, int N, int StartIndex = 0>
460
461 template <class T, int N, int StartIndex = 0>
463}
464
465template <class T, int NRows, int NCols, amrex::Order ORDER, int StartIndex>
466struct std::tuple_size<amrex::SmallMatrix<T,NRows,NCols,ORDER,StartIndex> >
467 : std::integral_constant<std::size_t,NRows*NCols> {};
468
469template <std::size_t N, class T, int NRows, int NCols, amrex::Order ORDER, int StartIndex>
470struct std::tuple_element<N, amrex::SmallMatrix<T,NRows,NCols,ORDER,StartIndex> >
471{
472 using type = T;
473};
474
475#endif
476
477/*
478 * Notes on why SmallMatrix matrix{} is zero initialized.
479 *
480 * SmallMatrix is not an aggregate, because it has a user declared default
481 * constructor. The rule is that, for `SmallMatrix matrix{}` with an empty
482 * brace-enclosed initializer list, value-initialization is performed. The
483 * effects of value-initialization of SmallMatrix (which has a user-declared
484 * but not user-provided default constructor) are that the matrix object is
485 * first zero-initialized and then the object's default constructor is
486 * applied. Since the default constructor does nothing, the final result is
487 * the object is zero-initialized.
488 *
489 * Why is SmallMatrix's default constructor user-declared not user-provided?
490 * It's because we first declare it with `SmallMatrix () = default`.
491 *
492 * Reference:
493 * https://en.cppreference.com/w/cpp/language/list_initialization
494 * https://en.cppreference.com/w/cpp/language/value_initialization
495 * https://en.cppreference.com/w/cpp/language/zero_initialization
496 */
#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
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:243
static constexpr int row_size
Definition AMReX_SmallMatrix.H:38
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:212
T m_mat[NRows *NCols]
Definition AMReX_SmallMatrix.H:407
static constexpr Order ordering
Definition AMReX_SmallMatrix.H:40
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:304
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 subtraction, lhs-rhs.
Definition AMReX_SmallMatrix.H:336
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:354
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:158
static constexpr int starting_index
Definition AMReX_SmallMatrix.H:41
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr SmallMatrix(Ts... vs)
Constructs column- or row-vector.
Definition AMReX_SmallMatrix.H:62
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE const T * end() const noexcept
Definition AMReX_SmallMatrix.H:193
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr SmallMatrix()=default
Default constructor.
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
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T * begin() noexcept
Definition AMReX_SmallMatrix.H:200
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:223
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T * end() noexcept
Definition AMReX_SmallMatrix.H:207
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T sum() const
Returns the sum of all elements in the matrix.
Definition AMReX_SmallMatrix.H:282
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
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:235
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:389
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:315
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:186
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T product() const
Returns the product of all elements in the matrix.
Definition AMReX_SmallMatrix.H:271
T & reference_type
Definition AMReX_SmallMatrix.H:37
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:365
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE SmallMatrix< T, NRows, NCols, ORDER, StartIndex > & transposeInPlace()
Transposes a square matrix in-place.
Definition AMReX_SmallMatrix.H:258
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T trace() const
Returns the trace of a square matrix.
Definition AMReX_SmallMatrix.H:294
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:325