Line data Source code
1 : #ifndef SET_MATRIX4_MAJOR_H
2 : #define SET_MATRIX4_MAJOR_H
3 :
4 : #include "Base.H"
5 : #include "Matrix4.H"
6 : #include "Matrix4_MajorMinor.H"
7 : #include "Matrix4_Isotropic.H"
8 : #include "Matrix4_Diagonal.H"
9 :
10 : namespace Set
11 : {
12 :
13 : template <>
14 : class Matrix4<2, Sym::Major>
15 : {
16 : Scalar data[10] = {NAN, NAN, NAN, NAN, NAN,
17 : NAN, NAN, NAN, NAN, NAN};
18 :
19 : public:
20 5646988 : AMREX_GPU_HOST_DEVICE Matrix4(){};
21 :
22 : AMREX_FORCE_INLINE
23 : const Scalar &operator()(const int i, const int j, const int k, const int l) const
24 : {
25 15504336 : int uid = i + 2 * j + 4 * k + 8 * l;
26 15504336 : if (uid==0 ) return data[0]; // [0000]
27 14535315 : else if (uid==8 || uid==2 ) return data[1]; // [0001] [0100]
28 12597273 : else if (uid==4 || uid==1 ) return data[2]; // [0010] [1000]
29 10659231 : else if (uid==12 || uid==3 ) return data[3]; // [0011] [1100]
30 8721189 : else if (uid==10 ) return data[4]; // [0101]
31 7752168 : else if (uid==6 || uid==9 ) return data[5]; // [0110] [1001]
32 5814126 : else if (uid==14 || uid==11 ) return data[6]; // [0111] [1101]
33 3876084 : else if (uid==5 ) return data[7]; // [1010]
34 2907063 : else if (uid==13 || uid==7 ) return data[8]; // [1011] [1110]
35 969021 : else if (uid==15 ) return data[9]; // [1111]
36 : else
37 0 : Util::Abort(INFO, "Index out of range");
38 0 : return Set::Garbage;
39 : }
40 :
41 : AMREX_FORCE_INLINE
42 : Scalar &operator()(const int i, const int j, const int k, const int l)
43 : {
44 1146944 : int uid = i + 2 * j + 4 * k + 8 * l;
45 1144464 : if (uid==0 ) return data[0]; // [0000]
46 1075260 : else if (uid==8 || uid==2 ) return data[1]; // [0001] [0100]
47 931892 : else if (uid==4 || uid==1 ) return data[2]; // [0010] [1000]
48 788524 : else if (uid==12 || uid==3 ) return data[3]; // [0011] [1100]
49 645156 : else if (uid==10 ) return data[4]; // [0101]
50 573472 : else if (uid==6 || uid==9 ) return data[5]; // [0110] [1001]
51 430104 : else if (uid==14 || uid==11 ) return data[6]; // [0111] [1101]
52 286736 : else if (uid==5 ) return data[7]; // [1010]
53 215052 : else if (uid==13 || uid==7 ) return data[8]; // [1011] [1110]
54 71684 : else if (uid==15 ) return data[9]; // [1111]
55 : else
56 0 : Util::Abort(INFO, "Index out of range");
57 0 : return Set::Garbage;
58 : }
59 :
60 :
61 : #if AMREX_SPACEDIM == 2
62 40 : Matrix4(const Matrix4<2,Sym::Isotropic> &in)
63 40 : {
64 : // TODO: very inefficient, need to optimize
65 120 : for (int i = 0; i < 2; i++)
66 240 : for (int j = 0; j < 2; j++)
67 480 : for (int k = 0; k < 2; k++)
68 960 : for (int l = 0; l < 2; l++)
69 1280 : (*this)(i,j,k,l) = in(i,j,k,l);
70 40 : }
71 20 : Matrix4(const Matrix4<2,Sym::Diagonal> &in)
72 20 : {
73 : // TODO: very inefficient, need to optimize
74 60 : for (int i = 0; i < 2; i++)
75 120 : for (int j = 0; j < 2; j++)
76 240 : for (int k = 0; k < 2; k++)
77 480 : for (int l = 0; l < 2; l++)
78 640 : (*this)(i,j,k,l) = in(i,j,k,l);
79 20 : }
80 100 : Matrix4(const Matrix4<2,Sym::MajorMinor> &in)
81 100 : {
82 : // TODO: very inefficient, need to optimize
83 300 : for (int i = 0; i < 2; i++)
84 600 : for (int j = 0; j < 2; j++)
85 1200 : for (int k = 0; k < 2; k++)
86 2400 : for (int l = 0; l < 2; l++)
87 3200 : (*this)(i,j,k,l) = in(i,j,k,l);
88 100 : }
89 : #endif
90 :
91 0 : void Print(std::ostream &os)
92 : {
93 0 : os.precision(4);
94 0 : for (int k = 0; k < 2; k++)
95 : {
96 0 : for (int i = -1; i < 3; i++)
97 : {
98 0 : for (int l = 0; l < 2; l++)
99 : {
100 0 : if (i == -1)
101 0 : os << "┌ ┐ ";
102 0 : else if (i == 2)
103 0 : os << "└ ┘ ";
104 : else
105 : {
106 0 : os << "│ ";
107 0 : for (int j = 0; j < 2; j++)
108 : {
109 0 : const Set::Scalar &val = (*this)(i, j, k, l);
110 0 : os << std::scientific << std::setw(11) << val; //(fabs(val)>1E-10 ? val : 0);
111 0 : os << " ";
112 : }
113 0 : os << "│ ";
114 : }
115 : }
116 0 : os << std::endl;
117 : }
118 : }
119 0 : }
120 : //AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
121 : //void operator=(const Matrix4<2, Sym::Major> &a) { for (int i = 0; i < 10; i++) data[i] = a.data[i]; }
122 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
123 0 : void operator+=(const Matrix4<2, Sym::Major> &a) { for (int i = 0; i < 10; i++) data[i] += a.data[i]; }
124 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
125 : void operator-=(const Matrix4<2, Sym::Major> &a) { for (int i = 0; i < 10; i++) data[i] -= a.data[i]; }
126 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
127 : void operator*=(const Matrix4<2, Sym::Major> &a) { for (int i = 0; i < 10; i++) data[i] *= a.data[i]; }
128 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
129 : void operator/=(const Matrix4<2, Sym::Major> &a) { for (int i = 0; i < 10; i++) data[i] /= a.data[i]; }
130 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
131 : void operator*=(const Set::Scalar &alpha) { for (int i = 0; i < 10; i++) data[i] *= alpha; }
132 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
133 : void operator/=(const Set::Scalar &alpha) { for (int i = 0; i < 10; i++) data[i] /= alpha; }
134 :
135 : static Matrix4<2, Sym::Major> Increment()
136 : {
137 : Matrix4<2, Sym::Major> ret;
138 : for (int i = 0; i < 10; i++) ret.data[i] = (Set::Scalar)i;
139 : return ret;
140 : }
141 1 : void Randomize()
142 : {
143 11 : for (int i = 0; i < 10; i++) data[i] = (Util::Random());
144 1 : }
145 1 : static Matrix4<2, Sym::Major> Random()
146 : {
147 1 : Matrix4<2, Sym::Major> ret;
148 1 : ret.Randomize();
149 1 : return ret;
150 : }
151 183 : static Matrix4<2, Sym::Major> Zero()
152 : {
153 183 : Matrix4<2, Sym::Major> ret;
154 2013 : for (int i = 0; i < 10; i++) ret.data[i] = 0.0;
155 183 : return ret;
156 : }
157 240 : Set::Scalar Norm()
158 : {
159 240 : Set::Scalar normsq = 0.0;
160 2640 : for (int i = 0; i < 10; i++) normsq += data[i]*data[i];
161 240 : return std::sqrt(normsq);
162 : }
163 : bool contains_nan() const
164 : {
165 : if (std::isnan(data[ 0])) return true;
166 : if (std::isnan(data[ 1])) return true;
167 : if (std::isnan(data[ 2])) return true;
168 : if (std::isnan(data[ 3])) return true;
169 : if (std::isnan(data[ 4])) return true;
170 : if (std::isnan(data[ 5])) return true;
171 : if (std::isnan(data[ 6])) return true;
172 : if (std::isnan(data[ 7])) return true;
173 : if (std::isnan(data[ 8])) return true;
174 : if (std::isnan(data[ 9])) return true;
175 : return false;
176 : }
177 :
178 : friend Matrix4<2, Sym::Major> operator+(const Matrix4<2, Sym::Major> &a, const Matrix4<2, Sym::Major> &b);
179 : friend Matrix4<2, Sym::Major> operator-(const Matrix4<2, Sym::Major> &a, const Matrix4<2, Sym::Major> &b);
180 : friend Matrix4<2, Sym::Major> operator*(const Matrix4<2, Sym::Major> &a, const Set::Scalar &b);
181 : friend Matrix4<2, Sym::Major> operator*(const Set::Scalar &b,const Matrix4<2, Sym::Major> &a);
182 : friend Matrix4<2, Sym::Major> operator/(const Matrix4<2, Sym::Major> &a, const Set::Scalar &b);
183 : friend Set::Matrix operator*(const Matrix4<2, Sym::Major> &a, const Set::Matrix &b);
184 : };
185 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
186 : Matrix4<2, Sym::Major> operator+(const Matrix4<2, Sym::Major> &a, const Matrix4<2, Sym::Major> &b)
187 : {
188 7760 : Matrix4<2,Sym::Major> ret;
189 85360 : for (int i = 0; i < 10; i++) ret.data[i] = a.data[i] + b.data[i];
190 7760 : return ret;
191 : }
192 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
193 : Matrix4<2, Sym::Major> operator-(const Matrix4<2, Sym::Major> &a, const Matrix4<2, Sym::Major> &b)
194 : {
195 1766118 : Matrix4<2,Sym::Major> ret;
196 19427298 : for (int i = 0; i < 10; i++) ret.data[i] = a.data[i] - b.data[i];
197 1766118 : return ret;
198 : }
199 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
200 : Matrix4<2, Sym::Major> operator*(const Matrix4<2, Sym::Major> &a, const Set::Scalar &b)
201 : {
202 1772395 : Matrix4<2,Sym::Major> ret;
203 19496345 : for (int i = 0; i < 10; i++) ret.data[i] = a.data[i] * b;
204 1772395 : return ret;
205 : }
206 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
207 : Matrix4<2, Sym::Major> operator*(const Set::Scalar &b, const Matrix4<2, Sym::Major> &a)
208 : {
209 0 : return a*b;
210 : }
211 :
212 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
213 : Matrix4<2, Sym::Major> operator/(const Matrix4<2, Sym::Major> &a, const Set::Scalar &b)
214 : {
215 1766747 : Matrix4<2,Sym::Major> ret;
216 19434217 : for (int i = 0; i < 10; i++) ret.data[i] = a.data[i] / b;
217 1766747 : return ret;
218 : }
219 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
220 : Set::Matrix operator*(const Matrix4<2, Sym::Major> &a, const Set::Matrix &b)
221 : {
222 4652526 : Set::Matrix ret;
223 4652526 : ret(0,0) = a.data[0]*b(0,0) + a.data[1]*b(0,1) + a.data[2]*b(1,0) + a.data[3]*b(1,1);
224 4652526 : ret(0,1) = a.data[1]*b(0,0) + a.data[4]*b(0,1) + a.data[5]*b(1,0) + a.data[6]*b(1,1);
225 4652526 : ret(1,0) = a.data[2]*b(0,0) + a.data[5]*b(0,1) + a.data[7]*b(1,0) + a.data[8]*b(1,1);
226 4652526 : ret(1,1) = a.data[3]*b(0,0) + a.data[6]*b(0,1) + a.data[8]*b(1,0) + a.data[9]*b(1,1);
227 4652526 : return ret;
228 : }
229 :
230 :
231 :
232 : template <>
233 : class Matrix4<3, Sym::Major>
234 : {
235 : Scalar data[45] = {NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN,
236 : NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN,
237 : NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN,
238 : NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN,
239 : NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN};
240 :
241 : public:
242 71265 : AMREX_GPU_HOST_DEVICE Matrix4(){};
243 :
244 : #if AMREX_SPACEDIM == 3
245 40 : Matrix4(const Matrix4<3,Sym::Isotropic> &in)
246 40 : {
247 : // TODO: very inefficient, need to optimize
248 160 : for (int i = 0; i < 3; i++)
249 480 : for (int j = 0; j < 3; j++)
250 1440 : for (int k = 0; k < 3; k++)
251 4320 : for (int l = 0; l < 3; l++)
252 3240 : (*this)(i,j,k,l) = in(i,j,k,l);
253 40 : }
254 20 : Matrix4(const Matrix4<3,Sym::Diagonal> &in)
255 20 : {
256 : // TODO: very inefficient, need to optimize
257 80 : for (int i = 0; i < 3; i++)
258 240 : for (int j = 0; j < 3; j++)
259 720 : for (int k = 0; k < 3; k++)
260 2160 : for (int l = 0; l < 3; l++)
261 1620 : (*this)(i,j,k,l) = in(i,j,k,l);
262 20 : }
263 100 : Matrix4(const Matrix4<3,Sym::MajorMinor> &in)
264 100 : {
265 : // TODO: very inefficient, need to optimize
266 400 : for (int i = 0; i < 3; i++)
267 1200 : for (int j = 0; j < 3; j++)
268 3600 : for (int k = 0; k < 3; k++)
269 10800 : for (int l = 0; l < 3; l++)
270 8100 : (*this)(i,j,k,l) = in(i,j,k,l);
271 100 : }
272 : #endif
273 :
274 : AMREX_FORCE_INLINE
275 : Scalar &operator()(const int i, const int j, const int k, const int l)
276 : {
277 12734282 : int uid = i + 3 * j + 9 * k + 27 * l;
278 :
279 12734102 : if (uid == 0) return data[0]; // [0000]
280 12520050 : else if (uid == 27 || uid == 3) return data[1]; // [0001] [0100]
281 12091586 : else if (uid == 54 || uid == 6) return data[2]; // [0002] [0200]
282 11805230 : else if (uid == 9 || uid == 1) return data[3]; // [0010] [1000]
283 11376766 : else if (uid == 36 || uid == 4) return data[4]; // [0011] [1100]
284 10948302 : else if (uid == 63 || uid == 7) return data[5]; // [0012] [1200]
285 10661946 : else if (uid == 18 || uid == 2) return data[6]; // [0020] [2000]
286 10375590 : else if (uid == 45 || uid == 5) return data[7]; // [0021] [2100]
287 10089234 : else if (uid == 72 || uid == 8) return data[8]; // [0022] [2200]
288 9802878 : else if (uid == 30) return data[9]; // [0101]
289 9588646 : else if (uid == 57 || uid == 33) return data[10]; // [0102] [0201]
290 9302290 : else if (uid == 12 || uid == 28) return data[11]; // [0110] [1001]
291 8873826 : else if (uid == 39 || uid == 31) return data[12]; // [0111] [1101]
292 8445362 : else if (uid == 66 || uid == 34) return data[13]; // [0112] [1201]
293 8159006 : else if (uid == 21 || uid == 29) return data[14]; // [0120] [2001]
294 7872650 : else if (uid == 48 || uid == 32) return data[15]; // [0121] [2101]
295 7586294 : else if (uid == 75 || uid == 35) return data[16]; // [0122] [2201]
296 7299938 : else if (uid == 60) return data[17]; // [0202]
297 7156760 : else if (uid == 15 || uid == 55) return data[18]; // [0210] [1002]
298 6870404 : else if (uid == 42 || uid == 58) return data[19]; // [0211] [1102]
299 6584048 : else if (uid == 69 || uid == 61) return data[20]; // [0212] [1202]
300 6297692 : else if (uid == 24 || uid == 56) return data[21]; // [0220] [2002]
301 6011336 : else if (uid == 51 || uid == 59) return data[22]; // [0221] [2102]
302 5724980 : else if (uid == 78 || uid == 62) return data[23]; // [0222] [2202]
303 5438624 : else if (uid == 10) return data[24]; // [1010]
304 5224392 : else if (uid == 37 || uid == 13) return data[25]; // [1011] [1110]
305 4795928 : else if (uid == 64 || uid == 16) return data[26]; // [1012] [1210]
306 4509572 : else if (uid == 19 || uid == 11) return data[27]; // [1020] [2010]
307 4223216 : else if (uid == 46 || uid == 14) return data[28]; // [1021] [2110]
308 3936860 : else if (uid == 73 || uid == 17) return data[29]; // [1022] [2210]
309 3650504 : else if (uid == 40) return data[30]; // [1111]
310 3436272 : else if (uid == 67 || uid == 43) return data[31]; // [1112] [1211]
311 3149916 : else if (uid == 22 || uid == 38) return data[32]; // [1120] [2011]
312 2863560 : else if (uid == 49 || uid == 41) return data[33]; // [1121] [2111]
313 2577204 : else if (uid == 76 || uid == 44) return data[34]; // [1122] [2211]
314 2290848 : else if (uid == 70) return data[35]; // [1212]
315 2147670 : else if (uid == 25 || uid == 65) return data[36]; // [1220] [2012]
316 1861314 : else if (uid == 52 || uid == 68) return data[37]; // [1221] [2112]
317 1574958 : else if (uid == 79 || uid == 71) return data[38]; // [1222] [2212]
318 1288602 : else if (uid == 20) return data[39]; // [2020]
319 1145424 : else if (uid == 47 || uid == 23) return data[40]; // [2021] [2120]
320 859068 : else if (uid == 74 || uid == 26) return data[41]; // [2022] [2220]
321 572712 : else if (uid == 50) return data[42]; // [2121]
322 429534 : else if (uid == 77 || uid == 53) return data[43]; // [2122] [2221]
323 143178 : else if (uid == 80) return data[44]; // [2222]
324 : else
325 0 : Util::Abort(INFO, "Index out of range");
326 0 : return Set::Garbage;
327 : }
328 :
329 :
330 : AMREX_FORCE_INLINE
331 : const Scalar &operator()(const int i, const int j, const int k, const int l) const
332 : {
333 0 : int uid = i + 3 * j + 9 * k + 27 * l;
334 :
335 0 : if (uid == 0) return data[0]; // [0000]
336 0 : else if (uid == 27 || uid == 3) return data[1]; // [0001] [0100]
337 0 : else if (uid == 54 || uid == 6) return data[2]; // [0002] [0200]
338 0 : else if (uid == 9 || uid == 1) return data[3]; // [0010] [1000]
339 0 : else if (uid == 36 || uid == 4) return data[4]; // [0011] [1100]
340 0 : else if (uid == 63 || uid == 7) return data[5]; // [0012] [1200]
341 0 : else if (uid == 18 || uid == 2) return data[6]; // [0020] [2000]
342 0 : else if (uid == 45 || uid == 5) return data[7]; // [0021] [2100]
343 0 : else if (uid == 72 || uid == 8) return data[8]; // [0022] [2200]
344 0 : else if (uid == 30) return data[9]; // [0101]
345 0 : else if (uid == 57 || uid == 33) return data[10]; // [0102] [0201]
346 0 : else if (uid == 12 || uid == 28) return data[11]; // [0110] [1001]
347 0 : else if (uid == 39 || uid == 31) return data[12]; // [0111] [1101]
348 0 : else if (uid == 66 || uid == 34) return data[13]; // [0112] [1201]
349 0 : else if (uid == 21 || uid == 29) return data[14]; // [0120] [2001]
350 0 : else if (uid == 48 || uid == 32) return data[15]; // [0121] [2101]
351 0 : else if (uid == 75 || uid == 35) return data[16]; // [0122] [2201]
352 0 : else if (uid == 60) return data[17]; // [0202]
353 0 : else if (uid == 15 || uid == 55) return data[18]; // [0210] [1002]
354 0 : else if (uid == 42 || uid == 58) return data[19]; // [0211] [1102]
355 0 : else if (uid == 69 || uid == 61) return data[20]; // [0212] [1202]
356 0 : else if (uid == 24 || uid == 56) return data[21]; // [0220] [2002]
357 0 : else if (uid == 51 || uid == 59) return data[22]; // [0221] [2102]
358 0 : else if (uid == 78 || uid == 62) return data[23]; // [0222] [2202]
359 0 : else if (uid == 10) return data[24]; // [1010]
360 0 : else if (uid == 37 || uid == 13) return data[25]; // [1011] [1110]
361 0 : else if (uid == 64 || uid == 16) return data[26]; // [1012] [1210]
362 0 : else if (uid == 19 || uid == 11) return data[27]; // [1020] [2010]
363 0 : else if (uid == 46 || uid == 14) return data[28]; // [1021] [2110]
364 0 : else if (uid == 73 || uid == 17) return data[29]; // [1022] [2210]
365 0 : else if (uid == 40) return data[30]; // [1111]
366 0 : else if (uid == 67 || uid == 43) return data[31]; // [1112] [1211]
367 0 : else if (uid == 22 || uid == 38) return data[32]; // [1120] [2011]
368 0 : else if (uid == 49 || uid == 41) return data[33]; // [1121] [2111]
369 0 : else if (uid == 76 || uid == 44) return data[34]; // [1122] [2211]
370 0 : else if (uid == 70) return data[35]; // [1212]
371 0 : else if (uid == 25 || uid == 65) return data[36]; // [1220] [2012]
372 0 : else if (uid == 52 || uid == 68) return data[37]; // [1221] [2112]
373 0 : else if (uid == 79 || uid == 71) return data[38]; // [1222] [2212]
374 0 : else if (uid == 20) return data[39]; // [2020]
375 0 : else if (uid == 47 || uid == 23) return data[40]; // [2021] [2120]
376 0 : else if (uid == 74 || uid == 26) return data[41]; // [2022] [2220]
377 0 : else if (uid == 50) return data[42]; // [2121]
378 0 : else if (uid == 77 || uid == 53) return data[43]; // [2122] [2221]
379 0 : else if (uid == 80) return data[44]; // [2222]
380 : else
381 0 : Util::Abort(INFO, "Index out of range");
382 0 : return Set::Garbage;
383 : }
384 :
385 240 : Set::Scalar Norm()
386 : {
387 240 : Set::Scalar retsq = 0;
388 11040 : for (int i = 0; i < 45; i++) retsq += data[i]*data[i];
389 240 : return std::sqrt(retsq);
390 : }
391 :
392 :
393 0 : void Print(std::ostream &os)
394 : {
395 0 : for (int i = 0; i < 45; i++)
396 0 : os << "i = " << i << " " << data[i] << std::endl;
397 :
398 0 : os.precision(4);
399 0 : for (int k = 0; k < 3; k++)
400 : {
401 0 : for (int i = -1; i < 4; i++)
402 : {
403 0 : for (int l = 0; l < 3; l++)
404 : {
405 0 : if (i == -1)
406 0 : os << "┌ ┐ ";
407 0 : else if (i == 3)
408 0 : os << "└ ┘ ";
409 : else
410 : {
411 0 : os << "│ ";
412 0 : for (int j = 0; j < 3; j++)
413 : {
414 0 : const Set::Scalar &val = (*this)(i, j, k, l);
415 0 : os << std::scientific << std::setw(11) << val; //(fabs(val)>1E-10 ? val : 0);
416 0 : os << " ";
417 : }
418 0 : os << "│ ";
419 : }
420 : }
421 0 : os << std::endl;
422 : }
423 : }
424 0 : }
425 : //AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
426 : //void operator=(const Matrix4<3, Sym::Major> &a) { for (int i = 0; i < 45; i++) data[i] = a.data[i]; }
427 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
428 0 : void operator+=(const Matrix4<3, Sym::Major> &a) { for (int i = 0; i < 45; i++) data[i] += a.data[i]; }
429 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
430 : void operator-=(const Matrix4<3, Sym::Major> &a) { for (int i = 0; i < 45; i++) data[i] -= a.data[i]; }
431 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
432 : void operator*=(const Matrix4<3, Sym::Major> &a) { for (int i = 0; i < 45; i++) data[i] *= a.data[i]; }
433 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
434 : void operator/=(const Matrix4<3, Sym::Major> &a) { for (int i = 0; i < 45; i++) data[i] /= a.data[i]; }
435 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
436 : void operator*=(const Set::Scalar &alpha) { for (int i = 0; i < 45; i++) data[i] *= alpha; }
437 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
438 : void operator/=(const Set::Scalar &alpha) { for (int i = 0; i < 45; i++) data[i] /= alpha; }
439 :
440 : static Matrix4<3, Sym::Major> Increment()
441 : {
442 : Matrix4<3, Sym::Major> ret;
443 : for (int i = 0; i < 45; i++) ret.data[i] = (Set::Scalar)i;
444 : return ret;
445 : }
446 1 : void Randomize()
447 : {
448 46 : for (int i = 0; i < 45; i++) data[i] = (Util::Random());
449 1 : }
450 1 : static Matrix4<3, Sym::Major> Random()
451 : {
452 1 : Matrix4<3, Sym::Major> ret;
453 1 : ret.Randomize();
454 1 : return ret;
455 : }
456 50 : static Matrix4<3, Sym::Major> Zero()
457 : {
458 50 : Matrix4<3, Sym::Major> ret;
459 2300 : for (int i = 0; i < 45; i++) ret.data[i] = 0.0;
460 50 : return ret;
461 : }
462 : Set::Scalar norm()
463 : {
464 : Set::Scalar normsq = 0.0;
465 : for (int i = 0; i < 45; i++) normsq += data[i]*data[i];
466 : return std::sqrt(normsq);
467 : }
468 : bool contains_nan() const
469 : {
470 : for (int i = 0; i < 45; i++)
471 : if (std::isnan(data[i])) return true;
472 : return false;
473 : }
474 :
475 : friend Matrix4<3, Sym::Major> operator-(const Matrix4<3, Sym::Major> &a, const Matrix4<3, Sym::Major> &b);
476 : friend Matrix4<3, Sym::Major> operator+(const Matrix4<3, Sym::Major> &a, const Matrix4<3, Sym::Major> &b);
477 : friend Set::Matrix operator*(const Matrix4<3, Sym::Major> &a, const Set::Matrix &b);
478 : friend Matrix4<3, Sym::Major> operator*(const Matrix4<3, Sym::Major> &a, const Set::Scalar &b);
479 : friend Matrix4<3, Sym::Major> operator*(const Set::Scalar &b,const Matrix4<3, Sym::Major> &a);
480 : friend Matrix4<3, Sym::Major> operator/(const Matrix4<3, Sym::Major> &a, const Set::Scalar &b);
481 : };
482 :
483 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
484 : Matrix4<3, Sym::Major> operator+(const Matrix4<3, Sym::Major> &a, const Matrix4<3, Sym::Major> &b)
485 : {
486 0 : Matrix4<3,Sym::Major> ret;
487 0 : for (int i = 0; i < 45; i++) ret.data[i] = a.data[i] + b.data[i];
488 0 : return ret;
489 : }
490 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
491 : Matrix4<3, Sym::Major> operator-(const Matrix4<3, Sym::Major> &a, const Matrix4<3, Sym::Major> &b)
492 : {
493 120 : Matrix4<3,Sym::Major> ret;
494 5520 : for (int i = 0; i < 45; i++) ret.data[i] = a.data[i] - b.data[i];
495 120 : return ret;
496 : }
497 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
498 : Matrix4<3, Sym::Major> operator*(const Matrix4<3, Sym::Major> &a, const Set::Scalar &b)
499 : {
500 0 : Matrix4<3,Sym::Major> ret;
501 0 : for (int i = 0; i < 45; i++) ret.data[i] = a.data[i] * b;
502 0 : return ret;
503 : }
504 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
505 : Matrix4<3, Sym::Major> operator*(const Set::Scalar &b, const Matrix4<3, Sym::Major> &a)
506 : {
507 0 : return a*b;
508 : }
509 :
510 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
511 : Matrix4<3, Sym::Major> operator/(const Matrix4<3, Sym::Major> &a, const Set::Scalar &b)
512 : {
513 0 : Matrix4<3,Sym::Major> ret;
514 0 : for (int i = 0; i < 45; i++) ret.data[i] = a.data[i] / b;
515 0 : return ret;
516 : }
517 :
518 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
519 : Set::Matrix operator*(const Matrix4<3, Sym::Major> &a, const Set::Matrix &b)
520 : {
521 0 : Set::Matrix ret;
522 0 : ret(0,0) = a.data[ 0]*b(0,0) + a.data[ 1]*b(0,1) + a.data[ 2]*b(0,2) + a.data[ 3]*b(1,0) + a.data[ 4]*b(1,1) + a.data[ 5]*b(1,2) + a.data[ 6]*b(2,0) + a.data[ 7]*b(2,1) + a.data[ 8]*b(2,2);
523 0 : ret(0,1) = a.data[ 1]*b(0,0) + a.data[ 9]*b(0,1) + a.data[10]*b(0,2) + a.data[11]*b(1,0) + a.data[12]*b(1,1) + a.data[13]*b(1,2) + a.data[14]*b(2,0) + a.data[15]*b(2,1) + a.data[16]*b(2,2);
524 0 : ret(0,2) = a.data[ 2]*b(0,0) + a.data[10]*b(0,1) + a.data[17]*b(0,2) + a.data[18]*b(1,0) + a.data[19]*b(1,1) + a.data[20]*b(1,2) + a.data[21]*b(2,0) + a.data[22]*b(2,1) + a.data[23]*b(2,2);
525 0 : ret(1,0) = a.data[ 3]*b(0,0) + a.data[11]*b(0,1) + a.data[18]*b(0,2) + a.data[24]*b(1,0) + a.data[25]*b(1,1) + a.data[26]*b(1,2) + a.data[27]*b(2,0) + a.data[28]*b(2,1) + a.data[29]*b(2,2);
526 0 : ret(1,1) = a.data[ 4]*b(0,0) + a.data[12]*b(0,1) + a.data[19]*b(0,2) + a.data[25]*b(1,0) + a.data[30]*b(1,1) + a.data[31]*b(1,2) + a.data[32]*b(2,0) + a.data[33]*b(2,1) + a.data[34]*b(2,2);
527 0 : ret(1,2) = a.data[ 5]*b(0,0) + a.data[13]*b(0,1) + a.data[20]*b(0,2) + a.data[26]*b(1,0) + a.data[31]*b(1,1) + a.data[35]*b(1,2) + a.data[36]*b(2,0) + a.data[37]*b(2,1) + a.data[38]*b(2,2);
528 0 : ret(2,0) = a.data[ 6]*b(0,0) + a.data[14]*b(0,1) + a.data[21]*b(0,2) + a.data[27]*b(1,0) + a.data[32]*b(1,1) + a.data[36]*b(1,2) + a.data[39]*b(2,0) + a.data[40]*b(2,1) + a.data[41]*b(2,2);
529 0 : ret(2,1) = a.data[ 7]*b(0,0) + a.data[15]*b(0,1) + a.data[22]*b(0,2) + a.data[28]*b(1,0) + a.data[33]*b(1,1) + a.data[37]*b(1,2) + a.data[40]*b(2,0) + a.data[42]*b(2,1) + a.data[43]*b(2,2);
530 0 : ret(2,2) = a.data[ 8]*b(0,0) + a.data[16]*b(0,1) + a.data[23]*b(0,2) + a.data[29]*b(1,0) + a.data[34]*b(1,1) + a.data[38]*b(1,2) + a.data[41]*b(2,0) + a.data[43]*b(2,1) + a.data[44]*b(2,2);
531 0 : return ret;
532 : }
533 :
534 : AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE
535 : Set::Vector operator * (const Matrix4<AMREX_SPACEDIM,Sym::Major> &a, const Set::Matrix3 &b)
536 : {
537 : // TODO: improve efficiency of this method
538 969021 : Set::Vector ret = Set::Vector::Zero();
539 2907063 : for (int i = 0; i < AMREX_SPACEDIM; i++)
540 : {
541 5814126 : for (int J = 0; J < AMREX_SPACEDIM; J++)
542 11628252 : for (int k = 0; k < AMREX_SPACEDIM; k++)
543 23256504 : for (int L = 0; L < AMREX_SPACEDIM; L++)
544 31008672 : ret(i) += a(i,J,k,L) * b(k,L,J);
545 : }
546 969021 : return ret;
547 : }
548 :
549 : }
550 : #endif
|
