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