Peano 4
Loading...
Searching...
No Matches
KernelUtils.h
Go to the documentation of this file.
1
14#pragma once
15
16#include "tarch/Assertions.h"
17#include "tarch/la/Vector.h"
19
20#include <sstream>
21
22namespace kernels {
23
29inline double pow1(const double& a) { return a; }
30inline double pow2(const double& a) { return a * a; }
31inline double pow3(const double& a) { return a * a * a; }
32inline double pow4(const double& a) { const double a2=a*a; return a2 * a2; }
33inline double pow5(const double& a) { const double a2=a*a; return a2 * a2 * a; }
34inline double pow6(const double& a) { const double a2=a*a; return a2 * a2 * a2; }
35inline double pow7(const double& a) { const double a2=a*a; return a2 * a2 * a2 * a;}
36inline double pow8(const double& a) { const double a2=a*a; const double a4=a2*a2; return a4 * a4; }
37inline double pow9(const double& a) { const double a3=a*a*a; return a3 * a3 * a3; }
38inline double pow10(const double& a) { const double a2=a*a; const double a4=a2*a2; return a4 * a4 * a2; }
39inline double pow11(const double& a) { const double a2=a*a; const double a4=a2*a2; return a4*a4*a2*a; }
40inline double pow12(const double& a) { const double a2=a*a; const double a4=a2*a2; return a4*a4*a4; }
41inline double pow13(const double& a) { const double a2=a*a; const double a4=a2*a2; return a4*a4*a4*a; }
42inline double pow14(const double& a) { const double a2=a*a; const double a4=a2*a2; return a4*a4*a4*a2; }
43inline double pow15(const double& a) { const double a2=a*a; const double a5=a2*a2*a; return a5*a5*a5; }
44inline double pow16(const double& a) { const double a2=a*a; const double a4=a2*a2; const double a8=a4*a4; return a8*a8; }
45
52typedef double (* UnivariateFunction) (const double& s);
53
62struct index {
63 // add further indices if you want to support more than maximal 6 indices.
64 const int i0, i1, i2, i3, i4, i5; // Lenghts
65 const int b0, b1, b2, b3, b4, b5; // Basis
66 const int size;
67
68 index(int j0=1, int j1=1, int j2=1, int j3=1, int j4=1, int j5=1) :
69 i0(j0), i1(j1), i2(j2), i3(j3), i4(j4), i5(j5),
70 b0(i1 * i2 * i3 * i4 * i5),
71 b1(i2 * i3 * i4 * i5),
72 b2(i3 * i4 * i5),
73 b3(i4 * i5),
74 b4(i5),
75 b5(1),
76 size(j0*j1*j2*j3*j4*j5) {}
77
78 index(const index&) = default; // implicit trivial copy constructor
79
83 int get(int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) const {
84 assertion2(j0 < i0, j1, i1);
85 assertion2(j1 < i1, j2, i2);
86 assertion2(j2 < i2, j3, i3);
87 assertion2(j3 < i3, j4, i4);
88 assertion2(j4 < i4, j5, i5);
89
90 return j0 * b0 + j1 * b1 + j2 * b2 + j3 * b3 + j4 * b4 + j5 * b5;
91 }
92
108 int rowMajor(const tarch::la::Vector<2,int>& j, int j2=0, int j3=0, int j4=0, int j5=0) const {
109 return get(j(1), j(0), j2, j3, j4, j5);
110 }
111 int colMajor(const tarch::la::Vector<2,int>& j, int j2=0, int j3=0, int j4=0, int j5=0) const {
112 return get(j(0), j(1), j2, j3, j4, j5);
113 }
114 int rowMajor(const tarch::la::Vector<3,int>& j, int j3=0, int j4=0, int j5=0) const {
115 return get(j(2), j(1), j(0), j3, j4, j5);
116 }
117 int colMajor(const tarch::la::Vector<3,int>& j, int j3=0, int j4=0, int j5=0) const {
118 return get(j(0), j(1), j(2), j3, j4, j5);
119 }
120
126 void rev(int pos, int& j0, int& j1, int& j2, int& j3, int& j4, int& j5) const {
127 // The algorithm exploits the short notation for int a,b that
128 // a/b == std::floor(a/b). The c's are the contributions at position
129 // i while the j's are the digits for position i. If you want to
130 // increase the overall amount of positions, extend logically.
131 j0 = (pos) / b0;
132 j1 = (pos-j0*b0) / b1;
133 j2 = (pos-j0*b0-j1*b1) / b2;
134 j3 = (pos-j0*b0-j1*b1-j2*b2) / b3;
135 j4 = (pos-j0*b0-j1*b1-j2*b2-j3*b3) / b4;
136 j5 = (pos-j0*b0-j1*b1-j2*b2-j3*b3-j4*b4) / b5;
137 }
138
139 /* As there are no reference default values, ie. we cannot write
140 * void rev(int pos, int& j0=NULL) or
141 * void rev(int pos, int& j1=int(), ...)
142 * we do it with pointers. Now you can write instead
143 * int i, j, k;
144 * myidx.rev(position, &i, &j, &k);
145 * which is much more convenient than tracking the non-used indices
146 * for your own.
147 */
148 void rev(int pos, int* const j0=NULL, int* const j1=NULL, int* const j2=NULL, int* const j3=NULL, int* const j4=NULL, int* const j5=NULL) const {
149 int a0=0, a1=0, a2=0, a3=0, a4=0, a5=0; // default storages
150 rev(pos, j0?*j0:a0, j1?*j1:a1, j2?*j2:a2, j3?*j3:a3, j4?*j4:a4, j5?*j5:a5);
151 }
152
153 // syntactic sugar:
154 int operator()(int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) const {
155 return get(j0, j1, j2, j3, j4, j5);
156 }
157
162 bool check(int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) const {
163 return (j0 < i0) && (j1 < i1) && (j2 < i2) && (j3 < i3) && (j4 < i4) && (j5 < i5);
164 }
165
169 static std::string strIndex(int min=0, int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) {
170 std::stringstream s;
171 s << "(";
172 s << j0;
173 if(j1 > min) s << ", " << j1;
174 if(j2 > min) s << ", " << j2;
175 if(j3 > min) s << ", " << j3;
176 if(j4 > min) s << ", " << j4;
177 if(j5 > min) s << ", " << j5;
178 s << ")";
179 return s.str();
180 }
181
183 std::string getStr(int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) const {
184 return strIndex(/*min*/ -1, j0,j1,j2,j3,j4,j5);
185 }
186 std::string revStr(int pos) const {
187 int j0,j1,j2,j3,j4,j5;
188 rev(pos, j0,j1,j2,j3,j4,j5);
189 return strIndex(/*min*/ -1, j0,j1,j2,j3,j4,j5);
190 }
191
193 std::string toString() const {
194 return strIndex(/*min*/ +1, i0,i1,i2,i3,i4,i5);
195 }
196};
197
203struct dindex : public index {
204 dindex(int max, int jN1=1, int jN2=1, int jN3=1 ) :
205 #if Dimensions == 2
206 kernels::index(max, max, jN1, jN2, jN3)
207 #elif Dimensions == 3
208 kernels::index(max, max, max, jN1, jN2, jN3)
209 #else
210 #error "ExaHyPE doesnt support chosen dimension. Only 2d and 3d are available"
211 #endif
212 {}
213};
214
226// In order to have a compile time sized class, we should have
227// template<typename T, std::size_t size> here or even better,
228// template<typename T, std::size j0, std::size j1, ...>
229// However, this will blow up the compile time/resulting binary, not that great.
230template<typename T>
231struct array {
234
235 array(int j0=1, int j1=1, int j2=1, int j3=1, int j4=1, int j5=1) : idx(j0,j1,j2,j3,j4,j5) {
236 data = new T[idx.size];
237 }
238
239 T& operator()(int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) const {
240 return data[idx(j0, j1, j2, j3, j4, j5)];
241 }
242
243 ~array() { delete data; }
244};
245
256template<typename T>
257struct shadow {
260
261 shadow(T* _data, int j0=1, int j1=1, int j2=1, int j3=1, int j4=1, int j5=1) : data(_data), idx(j0,j1,j2,j3,j4,j5) {}
262
263 T& operator()(int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) const {
264 return data[idx(j0, j1, j2, j3, j4, j5)];
265 }
266};
267
270
271
272struct idx2 {
273 idx2(int I, int J, int line = -1) : I_(I), J_(J), size(I*J), line_(line) {}
274
275 int operator()(int i, int j) const {
276 assertion3(i < I_, i, I_, line_);
277 assertion3(j < J_, j, J_, line_);
278 return i * J_ + j;
279 }
280
281 void rev(int pos, int& i, int& j) const {
282 assertion(pos < size);
283 i = pos % J_;
284 j = pos - i * J_;
285 }
286
287 const int I_, J_, size, line_;
288};
289
290struct idx3 {
291 idx3(int I, int J, int K, int line = -1) : I_(I), J_(J), K_(K), size(I*J*K), line_(line) {}
292
293 int operator()(int i, int j, int k) const {
294 assertion3(i < I_, i, I_, line_);
295 assertion3(j < J_, j, J_, line_);
296 assertion3(k < K_, k, K_, line_);
297 return i * (J_ * K_) + j * K_ + k;
298 }
299
300 const int I_, J_, K_, size, line_;
301};
302
303struct idx4 {
304 idx4(int I, int J, int K, int L, int line = -1)
305 : I_(I), J_(J), K_(K), L_(L), size(I*J*K*L), line_(line) {}
306
307 int operator()(int i, int j, int k, int l) const {
308 assertion3(i < I_, i, I_, line_);
309 assertion3(j < J_, j, J_, line_);
310 assertion3(k < K_, k, K_, line_);
311 assertion3(l < L_, l, L_, line_);
312 return i * (J_ * K_ * L_) + j * (K_ * L_) + k * L_ + l;
313 }
314
315 const int I_, J_, K_, L_, size, line_;
316};
317
318struct idx5 {
319 idx5(int I, int J, int K, int L, int M, int line = -1)
320 : I_(I), J_(J), K_(K), L_(L), M_(M), size(I*J*K*L*M), line_(line) {}
321
322 int operator()(int i, int j, int k, int l, int m) const {
323 assertion3(i < I_, i, I_, line_);
324 assertion3(j < J_, j, J_, line_);
325 assertion3(k < K_, k, K_, line_);
326 assertion3(l < L_, l, L_, line_);
327 assertion3(m < M_, m, M_, line_);
328 return i * (J_ * K_ * L_ * M_) + j * (K_ * L_ * M_) + k * (L_ * M_) +
329 l * M_ + m;
330 }
331
332 const int I_, J_, K_, L_, M_, size, line_;
333};
334
335struct idx6 {
336 idx6(int I, int J, int K, int L, int M, int N, int line = -1)
337 : I_(I), J_(J), K_(K), L_(L), M_(M), N_(N), size(I*J*K*L*M*N), line_(line) {}
338
339 int operator()(int i, int j, int k, int l, int m, int n) const {
340 assertion3(i < I_, i, I_, line_);
341 assertion3(j < J_, j, J_, line_);
342 assertion3(k < K_, k, K_, line_);
343 assertion3(l < L_, l, L_, line_);
344 assertion3(m < M_, m, M_, line_);
345 assertion3(n < N_, n, N_, line_);
346 return i * (J_ * K_ * L_ * M_ * N_) + j * (K_ * L_ * M_ * N_) +
347 k * (L_ * M_ * N_) + l * (M_ * N_) + m * N_ + n;
348 }
349
350 const int I_, J_, K_, L_, M_, N_, size, line_;
351};
352
353} // namespace kernels
354
#define assertion2(expr, param0, param1)
#define assertion3(expr, param0, param1, param2)
#define assertion(expr)
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
examples::exahype2::elastic::VariableShortcuts s
Definition loh.cpp:10
This file is part of the ExaHyPE project.
double pow7(const double &a)
Definition KernelUtils.h:35
double(* UnivariateFunction)(const double &s)
Univariate function.
Definition KernelUtils.h:52
double pow11(const double &a)
Definition KernelUtils.h:39
double pow8(const double &a)
Definition KernelUtils.h:36
shadow< double > dshadow
double pow15(const double &a)
Definition KernelUtils.h:43
double pow4(const double &a)
Definition KernelUtils.h:32
array< double > darray
double pow14(const double &a)
Definition KernelUtils.h:42
double pow10(const double &a)
Definition KernelUtils.h:38
double pow6(const double &a)
Definition KernelUtils.h:34
double pow1(const double &a)
Small integer powers with minimum number of multiplications.
Definition KernelUtils.h:29
double pow9(const double &a)
Definition KernelUtils.h:37
double pow13(const double &a)
Definition KernelUtils.h:41
double pow12(const double &a)
Definition KernelUtils.h:40
double pow16(const double &a)
Definition KernelUtils.h:44
double pow2(const double &a)
Definition KernelUtils.h:30
double pow5(const double &a)
Definition KernelUtils.h:33
double pow3(const double &a)
Definition KernelUtils.h:31
A tiny class for storing a continous n-dimensional array together with its sizes.
T & operator()(int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) const
array(int j0=1, int j1=1, int j2=1, int j3=1, int j4=1, int j5=1)
Returns conveniently an index with DIMENSIONS entries, each the argument max.
dindex(int max, int jN1=1, int jN2=1, int jN3=1)
idx2(int I, int J, int line=-1)
int operator()(int i, int j) const
const int J_
void rev(int pos, int &i, int &j) const
const int size
const int line_
const int I_
int operator()(int i, int j, int k) const
const int size
const int I_
idx3(int I, int J, int K, int line=-1)
const int J_
const int K_
const int line_
const int L_
int operator()(int i, int j, int k, int l) const
const int line_
const int I_
const int size
idx4(int I, int J, int K, int L, int line=-1)
const int K_
const int J_
const int I_
idx5(int I, int J, int K, int L, int M, int line=-1)
const int K_
const int M_
const int J_
const int L_
const int line_
int operator()(int i, int j, int k, int l, int m) const
const int size
const int J_
const int N_
const int line_
const int K_
idx6(int I, int J, int K, int L, int M, int N, int line=-1)
int operator()(int i, int j, int k, int l, int m, int n) const
const int I_
const int L_
const int M_
const int size
This is a single successor class for the idx2, idx3, idx4, idx5, idx6 classes.
Definition KernelUtils.h:62
std::string getStr(int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) const
Some index to string.
std::string revStr(int pos) const
const int size
Definition KernelUtils.h:66
bool check(int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) const
Checks if the given indices are in the valid range.
int rowMajor(const tarch::la::Vector< 2, int > &j, int j2=0, int j3=0, int j4=0, int j5=0) const
A nice feature for the Peano likers.
const int i2
Definition KernelUtils.h:64
const int b0
Definition KernelUtils.h:65
std::string toString() const
Object to string.
const int i4
Definition KernelUtils.h:64
void rev(int pos, int &j0, int &j1, int &j2, int &j3, int &j4, int &j5) const
Inverse index: Get the index tuple for a given index.
int get(int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) const
Compute a single index ("superindex", global index, ...) from the tuples.
Definition KernelUtils.h:83
int colMajor(const tarch::la::Vector< 3, int > &j, int j3=0, int j4=0, int j5=0) const
index(int j0=1, int j1=1, int j2=1, int j3=1, int j4=1, int j5=1)
Definition KernelUtils.h:68
static std::string strIndex(int min=0, int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0)
Allows to print index tuples.
const int b2
Definition KernelUtils.h:65
const int b1
Definition KernelUtils.h:65
void rev(int pos, int *const j0=NULL, int *const j1=NULL, int *const j2=NULL, int *const j3=NULL, int *const j4=NULL, int *const j5=NULL) const
const int i5
Definition KernelUtils.h:64
index(const index &)=default
const int i3
Definition KernelUtils.h:64
const int i0
Definition KernelUtils.h:64
const int b3
Definition KernelUtils.h:65
int colMajor(const tarch::la::Vector< 2, int > &j, int j2=0, int j3=0, int j4=0, int j5=0) const
const int b5
Definition KernelUtils.h:65
const int b4
Definition KernelUtils.h:65
int operator()(int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) const
const int i1
Definition KernelUtils.h:64
int rowMajor(const tarch::la::Vector< 3, int > &j, int j3=0, int j4=0, int j5=0) const
The same as array, just without local storage but instead shadowing data which is not owned.
shadow(T *_data, int j0=1, int j1=1, int j2=1, int j3=1, int j4=1, int j5=1)
T & operator()(int j0=0, int j1=0, int j2=0, int j3=0, int j4=0, int j5=0) const
Simple vector class.
Definition Vector.h:134