1 /* ----------------------------------------------------------------------
2 SPARTA - Stochastic PArallel Rarefied-gas Time-accurate Analyzer
3 http://sparta.sandia.gov
4 Steve Plimpton, sjplimp@sandia.gov, Michael Gallis, magalli@sandia.gov
5 Sandia National Laboratories
6
7 Copyright (2014) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level SPARTA directory.
13 ------------------------------------------------------------------------- */
14
15 /* ----------------------------------------------------------------------
16 Contributing author: Mike Brown (SNL)
17 ------------------------------------------------------------------------- */
18
19 #ifndef SPARTA_MATH_EXTRA_KOKKOS_H
20 #define SPARTA_MATH_EXTRA_KOKKOS_H
21
22 #include "spatype.h"
23 #include "math.h"
24 #include "stdio.h"
25 #include "string.h"
26
27 namespace MathExtraKokkos {
28
29 // 3 vector operations
30
31 KOKKOS_INLINE_FUNCTION void norm3(double *v);
32 KOKKOS_INLINE_FUNCTION void normalize3(const double *v, double *ans);
33 KOKKOS_INLINE_FUNCTION void snorm3(const double, double *v);
34 KOKKOS_INLINE_FUNCTION void snormalize3(const double, const double *v, double *ans);
35 KOKKOS_INLINE_FUNCTION void negate3(double *v);
36 KOKKOS_INLINE_FUNCTION void scale3(double s, double *v);
37 KOKKOS_INLINE_FUNCTION void scale3(double s, const double *v, double *ans);
38 KOKKOS_INLINE_FUNCTION void axpy3(double alpha, const double *x, double *y);
39 KOKKOS_INLINE_FUNCTION void axpy3(double alpha, const double *x, const double *y,
40 double *ynew);
41 KOKKOS_INLINE_FUNCTION void add3(const double *v1, const double *v2, double *ans);
42 KOKKOS_INLINE_FUNCTION void sub3(const double *v1, const double *v2, double *ans);
43 KOKKOS_INLINE_FUNCTION double len3(const double *v);
44 KOKKOS_INLINE_FUNCTION double lensq3(const double *v);
45 KOKKOS_INLINE_FUNCTION double dot3(const double *v1, const double *v2);
46 KOKKOS_INLINE_FUNCTION void cross3(const double *v1, const double *v2, double *ans);
47 KOKKOS_INLINE_FUNCTION void reflect3(double *v, const double *norm);
48
49 // 3x3 matrix operations
50
51 KOKKOS_INLINE_FUNCTION double det3(const double mat[3][3]);
52 KOKKOS_INLINE_FUNCTION void diag_times3(const double *diagonal, const double mat[3][3],
53 double ans[3][3]);
54 KOKKOS_INLINE_FUNCTION void plus3(const double m[3][3], const double m2[3][3],
55 double ans[3][3]);
56 KOKKOS_INLINE_FUNCTION void times3(const double m[3][3], const double m2[3][3],
57 double ans[3][3]);
58 KOKKOS_INLINE_FUNCTION void transpose_times3(const double mat1[3][3],
59 const double mat2[3][3],
60 double ans[3][3]);
61 KOKKOS_INLINE_FUNCTION void times3_transpose(const double mat1[3][3],
62 const double mat2[3][3],
63 double ans[3][3]);
64 KOKKOS_INLINE_FUNCTION void invert3(const double mat[3][3], double ans[3][3]);
65 KOKKOS_INLINE_FUNCTION void matvec(const double mat[3][3], const double*vec, double *ans);
66 KOKKOS_INLINE_FUNCTION void matvec(const double *ex, const double *ey, const double *ez,
67 const double *vec, double *ans);
68 KOKKOS_INLINE_FUNCTION void transpose_matvec(const double mat[3][3], const double*vec,
69 double *ans);
70 KOKKOS_INLINE_FUNCTION void transpose_matvec(const double *ex, const double *ey,
71 const double *ez, const double *v,
72 double *ans);
73 KOKKOS_INLINE_FUNCTION void transpose_diag3(const double mat[3][3], const double*vec,
74 double ans[3][3]);
75 KOKKOS_INLINE_FUNCTION void vecmat(const double *v, const double m[3][3], double *ans);
76 KOKKOS_INLINE_FUNCTION void scalar_times3(const double f, double m[3][3]);
77
78 // quaternion operations
79
80 KOKKOS_INLINE_FUNCTION void axisangle_to_quat(const double *v, const double angle,
81 double *quat);
82 KOKKOS_INLINE_FUNCTION void quat_to_mat(const double *quat, double mat[3][3]);
83 }
84
85 /* ----------------------------------------------------------------------
86 normalize a vector in place
87 ------------------------------------------------------------------------- */
88 KOKKOS_INLINE_FUNCTION
norm3(double * v)89 void MathExtraKokkos::norm3(double *v)
90 {
91 double scale = 1.0/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
92 v[0] *= scale;
93 v[1] *= scale;
94 v[2] *= scale;
95 }
96
97 /* ----------------------------------------------------------------------
98 normalize a vector, return in ans
99 ------------------------------------------------------------------------- */
100
101 KOKKOS_INLINE_FUNCTION
normalize3(const double * v,double * ans)102 void MathExtraKokkos::normalize3(const double *v, double *ans)
103 {
104 double scale = 1.0/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
105 ans[0] = v[0]*scale;
106 ans[1] = v[1]*scale;
107 ans[2] = v[2]*scale;
108 }
109
110 /* ----------------------------------------------------------------------
111 scale a vector to length in place
112 ------------------------------------------------------------------------- */
113
114 KOKKOS_INLINE_FUNCTION
snorm3(const double length,double * v)115 void MathExtraKokkos::snorm3(const double length, double *v)
116 {
117 double scale = length/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
118 v[0] *= scale;
119 v[1] *= scale;
120 v[2] *= scale;
121 }
122
123 /* ----------------------------------------------------------------------
124 scale a vector to length
125 ------------------------------------------------------------------------- */
126
127 KOKKOS_INLINE_FUNCTION
snormalize3(const double length,const double * v,double * ans)128 void MathExtraKokkos::snormalize3(const double length, const double *v, double *ans)
129 {
130 double scale = length/sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
131 ans[0] = v[0]*scale;
132 ans[1] = v[1]*scale;
133 ans[2] = v[2]*scale;
134 }
135
136 /* ----------------------------------------------------------------------
137 negate vector v in place
138 ------------------------------------------------------------------------- */
139
140 KOKKOS_INLINE_FUNCTION
negate3(double * v)141 void MathExtraKokkos::negate3(double *v)
142 {
143 v[0] = -v[0];
144 v[1] = -v[1];
145 v[2] = -v[2];
146 }
147
148 /* ----------------------------------------------------------------------
149 scale vector v by s in place
150 ------------------------------------------------------------------------- */
151
152 KOKKOS_INLINE_FUNCTION
scale3(double s,double * v)153 void MathExtraKokkos::scale3(double s, double *v)
154 {
155 v[0] *= s;
156 v[1] *= s;
157 v[2] *= s;
158 }
159
160 /* ----------------------------------------------------------------------
161 scale vector v by s, return in ans
162 ------------------------------------------------------------------------- */
163
164 KOKKOS_INLINE_FUNCTION
scale3(double s,const double * v,double * ans)165 void MathExtraKokkos::scale3(double s, const double *v, double *ans)
166 {
167 ans[0] = s*v[0];
168 ans[1] = s*v[1];
169 ans[2] = s*v[2];
170 }
171
172 /* ----------------------------------------------------------------------
173 axpy: y = alpha*x + y
174 y is replaced by result
175 ------------------------------------------------------------------------- */
176
177 KOKKOS_INLINE_FUNCTION
axpy3(double alpha,const double * x,double * y)178 void MathExtraKokkos::axpy3(double alpha, const double *x, double *y)
179 {
180 y[0] += alpha*x[0];
181 y[1] += alpha*x[1];
182 y[2] += alpha*x[2];
183 }
184
185 /* ----------------------------------------------------------------------
186 axpy: ynew = alpha*x + y
187 ------------------------------------------------------------------------- */
188
189 KOKKOS_INLINE_FUNCTION
axpy3(double alpha,const double * x,const double * y,double * ynew)190 void MathExtraKokkos::axpy3(double alpha, const double *x, const double *y,
191 double *ynew)
192 {
193 ynew[0] += alpha*x[0] + y[0];
194 ynew[1] += alpha*x[1] + y[1];
195 ynew[2] += alpha*x[2] + y[2];
196 }
197
198 /* ----------------------------------------------------------------------
199 ans = v1 + v2
200 ------------------------------------------------------------------------- */
201
202 KOKKOS_INLINE_FUNCTION
add3(const double * v1,const double * v2,double * ans)203 void MathExtraKokkos::add3(const double *v1, const double *v2, double *ans)
204 {
205 ans[0] = v1[0] + v2[0];
206 ans[1] = v1[1] + v2[1];
207 ans[2] = v1[2] + v2[2];
208 }
209
210 /* ----------------------------------------------------------------------
211 ans = v1 - v2
212 ------------------------------------------------------------------------- */
213
214 KOKKOS_INLINE_FUNCTION
sub3(const double * v1,const double * v2,double * ans)215 void MathExtraKokkos::sub3(const double *v1, const double *v2, double *ans)
216 {
217 ans[0] = v1[0] - v2[0];
218 ans[1] = v1[1] - v2[1];
219 ans[2] = v1[2] - v2[2];
220 }
221
222 /* ----------------------------------------------------------------------
223 length of vector v
224 ------------------------------------------------------------------------- */
225
226 KOKKOS_INLINE_FUNCTION
len3(const double * v)227 double MathExtraKokkos::len3(const double *v)
228 {
229 return sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
230 }
231
232 /* ----------------------------------------------------------------------
233 squared length of vector v, or dot product of v with itself
234 ------------------------------------------------------------------------- */
235
236 KOKKOS_INLINE_FUNCTION
lensq3(const double * v)237 double MathExtraKokkos::lensq3(const double *v)
238 {
239 return v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
240 }
241
242 /* ----------------------------------------------------------------------
243 dot product of 2 vectors
244 ------------------------------------------------------------------------- */
245
246 KOKKOS_INLINE_FUNCTION
dot3(const double * v1,const double * v2)247 double MathExtraKokkos::dot3(const double *v1, const double *v2)
248 {
249 return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
250 }
251
252 /* ----------------------------------------------------------------------
253 cross product of 2 vectors
254 ------------------------------------------------------------------------- */
255
256 KOKKOS_INLINE_FUNCTION
cross3(const double * v1,const double * v2,double * ans)257 void MathExtraKokkos::cross3(const double *v1, const double *v2, double *ans)
258 {
259 ans[0] = v1[1]*v2[2] - v1[2]*v2[1];
260 ans[1] = v1[2]*v2[0] - v1[0]*v2[2];
261 ans[2] = v1[0]*v2[1] - v1[1]*v2[0];
262 }
263
264 /* ----------------------------------------------------------------------
265 reflect vector v around unit normal n
266 return updated v of same length = v - 2(v dot n)n
267 ------------------------------------------------------------------------- */
268
269 KOKKOS_INLINE_FUNCTION
reflect3(double * v,const double * n)270 void MathExtraKokkos::reflect3(double *v, const double *n)
271 {
272 double dot = dot3(v,n);
273 v[0] -= 2.0*dot*n[0];
274 v[1] -= 2.0*dot*n[1];
275 v[2] -= 2.0*dot*n[2];
276 }
277
278 /* ----------------------------------------------------------------------
279 determinant of a matrix
280 ------------------------------------------------------------------------- */
281
282 KOKKOS_INLINE_FUNCTION
det3(const double m[3][3])283 double MathExtraKokkos::det3(const double m[3][3])
284 {
285 double ans = m[0][0]*m[1][1]*m[2][2] - m[0][0]*m[1][2]*m[2][1] -
286 m[1][0]*m[0][1]*m[2][2] + m[1][0]*m[0][2]*m[2][1] +
287 m[2][0]*m[0][1]*m[1][2] - m[2][0]*m[0][2]*m[1][1];
288 return ans;
289 }
290
291 /* ----------------------------------------------------------------------
292 diagonal matrix times a full matrix
293 ------------------------------------------------------------------------- */
294
295 KOKKOS_INLINE_FUNCTION
diag_times3(const double * d,const double m[3][3],double ans[3][3])296 void MathExtraKokkos::diag_times3(const double *d, const double m[3][3],
297 double ans[3][3])
298 {
299 ans[0][0] = d[0]*m[0][0];
300 ans[0][1] = d[0]*m[0][1];
301 ans[0][2] = d[0]*m[0][2];
302 ans[1][0] = d[1]*m[1][0];
303 ans[1][1] = d[1]*m[1][1];
304 ans[1][2] = d[1]*m[1][2];
305 ans[2][0] = d[2]*m[2][0];
306 ans[2][1] = d[2]*m[2][1];
307 ans[2][2] = d[2]*m[2][2];
308 }
309
310 /* ----------------------------------------------------------------------
311 add two matrices
312 ------------------------------------------------------------------------- */
313
314 KOKKOS_INLINE_FUNCTION
plus3(const double m[3][3],const double m2[3][3],double ans[3][3])315 void MathExtraKokkos::plus3(const double m[3][3], const double m2[3][3],
316 double ans[3][3])
317 {
318 ans[0][0] = m[0][0]+m2[0][0];
319 ans[0][1] = m[0][1]+m2[0][1];
320 ans[0][2] = m[0][2]+m2[0][2];
321 ans[1][0] = m[1][0]+m2[1][0];
322 ans[1][1] = m[1][1]+m2[1][1];
323 ans[1][2] = m[1][2]+m2[1][2];
324 ans[2][0] = m[2][0]+m2[2][0];
325 ans[2][1] = m[2][1]+m2[2][1];
326 ans[2][2] = m[2][2]+m2[2][2];
327 }
328
329 /* ----------------------------------------------------------------------
330 multiply mat1 times mat2
331 ------------------------------------------------------------------------- */
332
333 KOKKOS_INLINE_FUNCTION
times3(const double m[3][3],const double m2[3][3],double ans[3][3])334 void MathExtraKokkos::times3(const double m[3][3], const double m2[3][3],
335 double ans[3][3])
336 {
337 ans[0][0] = m[0][0]*m2[0][0] + m[0][1]*m2[1][0] + m[0][2]*m2[2][0];
338 ans[0][1] = m[0][0]*m2[0][1] + m[0][1]*m2[1][1] + m[0][2]*m2[2][1];
339 ans[0][2] = m[0][0]*m2[0][2] + m[0][1]*m2[1][2] + m[0][2]*m2[2][2];
340 ans[1][0] = m[1][0]*m2[0][0] + m[1][1]*m2[1][0] + m[1][2]*m2[2][0];
341 ans[1][1] = m[1][0]*m2[0][1] + m[1][1]*m2[1][1] + m[1][2]*m2[2][1];
342 ans[1][2] = m[1][0]*m2[0][2] + m[1][1]*m2[1][2] + m[1][2]*m2[2][2];
343 ans[2][0] = m[2][0]*m2[0][0] + m[2][1]*m2[1][0] + m[2][2]*m2[2][0];
344 ans[2][1] = m[2][0]*m2[0][1] + m[2][1]*m2[1][1] + m[2][2]*m2[2][1];
345 ans[2][2] = m[2][0]*m2[0][2] + m[2][1]*m2[1][2] + m[2][2]*m2[2][2];
346 }
347
348 /* ----------------------------------------------------------------------
349 multiply the transpose of mat1 times mat2
350 ------------------------------------------------------------------------- */
351
352 KOKKOS_INLINE_FUNCTION
transpose_times3(const double m[3][3],const double m2[3][3],double ans[3][3])353 void MathExtraKokkos::transpose_times3(const double m[3][3], const double m2[3][3],
354 double ans[3][3])
355 {
356 ans[0][0] = m[0][0]*m2[0][0] + m[1][0]*m2[1][0] + m[2][0]*m2[2][0];
357 ans[0][1] = m[0][0]*m2[0][1] + m[1][0]*m2[1][1] + m[2][0]*m2[2][1];
358 ans[0][2] = m[0][0]*m2[0][2] + m[1][0]*m2[1][2] + m[2][0]*m2[2][2];
359 ans[1][0] = m[0][1]*m2[0][0] + m[1][1]*m2[1][0] + m[2][1]*m2[2][0];
360 ans[1][1] = m[0][1]*m2[0][1] + m[1][1]*m2[1][1] + m[2][1]*m2[2][1];
361 ans[1][2] = m[0][1]*m2[0][2] + m[1][1]*m2[1][2] + m[2][1]*m2[2][2];
362 ans[2][0] = m[0][2]*m2[0][0] + m[1][2]*m2[1][0] + m[2][2]*m2[2][0];
363 ans[2][1] = m[0][2]*m2[0][1] + m[1][2]*m2[1][1] + m[2][2]*m2[2][1];
364 ans[2][2] = m[0][2]*m2[0][2] + m[1][2]*m2[1][2] + m[2][2]*m2[2][2];
365 }
366
367 /* ----------------------------------------------------------------------
368 multiply mat1 times transpose of mat2
369 ------------------------------------------------------------------------- */
370
371 KOKKOS_INLINE_FUNCTION
times3_transpose(const double m[3][3],const double m2[3][3],double ans[3][3])372 void MathExtraKokkos::times3_transpose(const double m[3][3], const double m2[3][3],
373 double ans[3][3])
374 {
375 ans[0][0] = m[0][0]*m2[0][0] + m[0][1]*m2[0][1] + m[0][2]*m2[0][2];
376 ans[0][1] = m[0][0]*m2[1][0] + m[0][1]*m2[1][1] + m[0][2]*m2[1][2];
377 ans[0][2] = m[0][0]*m2[2][0] + m[0][1]*m2[2][1] + m[0][2]*m2[2][2];
378 ans[1][0] = m[1][0]*m2[0][0] + m[1][1]*m2[0][1] + m[1][2]*m2[0][2];
379 ans[1][1] = m[1][0]*m2[1][0] + m[1][1]*m2[1][1] + m[1][2]*m2[1][2];
380 ans[1][2] = m[1][0]*m2[2][0] + m[1][1]*m2[2][1] + m[1][2]*m2[2][2];
381 ans[2][0] = m[2][0]*m2[0][0] + m[2][1]*m2[0][1] + m[2][2]*m2[0][2];
382 ans[2][1] = m[2][0]*m2[1][0] + m[2][1]*m2[1][1] + m[2][2]*m2[1][2];
383 ans[2][2] = m[2][0]*m2[2][0] + m[2][1]*m2[2][1] + m[2][2]*m2[2][2];
384 }
385
386 /* ----------------------------------------------------------------------
387 invert a matrix
388 does NOT checks for singular or badly scaled matrix
389 ------------------------------------------------------------------------- */
390
391 KOKKOS_INLINE_FUNCTION
invert3(const double m[3][3],double ans[3][3])392 void MathExtraKokkos::invert3(const double m[3][3], double ans[3][3])
393 {
394 double den = m[0][0]*m[1][1]*m[2][2]-m[0][0]*m[1][2]*m[2][1];
395 den += -m[1][0]*m[0][1]*m[2][2]+m[1][0]*m[0][2]*m[2][1];
396 den += m[2][0]*m[0][1]*m[1][2]-m[2][0]*m[0][2]*m[1][1];
397
398 ans[0][0] = (m[1][1]*m[2][2]-m[1][2]*m[2][1]) / den;
399 ans[0][1] = -(m[0][1]*m[2][2]-m[0][2]*m[2][1]) / den;
400 ans[0][2] = (m[0][1]*m[1][2]-m[0][2]*m[1][1]) / den;
401 ans[1][0] = -(m[1][0]*m[2][2]-m[1][2]*m[2][0]) / den;
402 ans[1][1] = (m[0][0]*m[2][2]-m[0][2]*m[2][0]) / den;
403 ans[1][2] = -(m[0][0]*m[1][2]-m[0][2]*m[1][0]) / den;
404 ans[2][0] = (m[1][0]*m[2][1]-m[1][1]*m[2][0]) / den;
405 ans[2][1] = -(m[0][0]*m[2][1]-m[0][1]*m[2][0]) / den;
406 ans[2][2] = (m[0][0]*m[1][1]-m[0][1]*m[1][0]) / den;
407 }
408
409 /* ----------------------------------------------------------------------
410 matrix times vector
411 ------------------------------------------------------------------------- */
412
413 KOKKOS_INLINE_FUNCTION
matvec(const double m[3][3],const double * v,double * ans)414 void MathExtraKokkos::matvec(const double m[3][3], const double *v, double *ans)
415 {
416 ans[0] = m[0][0]*v[0] + m[0][1]*v[1] + m[0][2]*v[2];
417 ans[1] = m[1][0]*v[0] + m[1][1]*v[1] + m[1][2]*v[2];
418 ans[2] = m[2][0]*v[0] + m[2][1]*v[1] + m[2][2]*v[2];
419 }
420
421 /* ----------------------------------------------------------------------
422 matrix times vector
423 ------------------------------------------------------------------------- */
424
425 KOKKOS_INLINE_FUNCTION
matvec(const double * ex,const double * ey,const double * ez,const double * v,double * ans)426 void MathExtraKokkos::matvec(const double *ex, const double *ey, const double *ez,
427 const double *v, double *ans)
428 {
429 ans[0] = ex[0]*v[0] + ey[0]*v[1] + ez[0]*v[2];
430 ans[1] = ex[1]*v[0] + ey[1]*v[1] + ez[1]*v[2];
431 ans[2] = ex[2]*v[0] + ey[2]*v[1] + ez[2]*v[2];
432 }
433
434 /* ----------------------------------------------------------------------
435 transposed matrix times vector
436 ------------------------------------------------------------------------- */
437
438 KOKKOS_INLINE_FUNCTION
transpose_matvec(const double m[3][3],const double * v,double * ans)439 void MathExtraKokkos::transpose_matvec(const double m[3][3], const double *v,
440 double *ans)
441 {
442 ans[0] = m[0][0]*v[0] + m[1][0]*v[1] + m[2][0]*v[2];
443 ans[1] = m[0][1]*v[0] + m[1][1]*v[1] + m[2][1]*v[2];
444 ans[2] = m[0][2]*v[0] + m[1][2]*v[1] + m[2][2]*v[2];
445 }
446
447 /* ----------------------------------------------------------------------
448 transposed matrix times vector
449 ------------------------------------------------------------------------- */
450
451 KOKKOS_INLINE_FUNCTION
transpose_matvec(const double * ex,const double * ey,const double * ez,const double * v,double * ans)452 void MathExtraKokkos::transpose_matvec(const double *ex, const double *ey,
453 const double *ez, const double *v,
454 double *ans)
455 {
456 ans[0] = ex[0]*v[0] + ex[1]*v[1] + ex[2]*v[2];
457 ans[1] = ey[0]*v[0] + ey[1]*v[1] + ey[2]*v[2];
458 ans[2] = ez[0]*v[0] + ez[1]*v[1] + ez[2]*v[2];
459 }
460
461 /* ----------------------------------------------------------------------
462 transposed matrix times diagonal matrix
463 ------------------------------------------------------------------------- */
464
465 KOKKOS_INLINE_FUNCTION
transpose_diag3(const double m[3][3],const double * d,double ans[3][3])466 void MathExtraKokkos::transpose_diag3(const double m[3][3], const double *d,
467 double ans[3][3])
468 {
469 ans[0][0] = m[0][0]*d[0];
470 ans[0][1] = m[1][0]*d[1];
471 ans[0][2] = m[2][0]*d[2];
472 ans[1][0] = m[0][1]*d[0];
473 ans[1][1] = m[1][1]*d[1];
474 ans[1][2] = m[2][1]*d[2];
475 ans[2][0] = m[0][2]*d[0];
476 ans[2][1] = m[1][2]*d[1];
477 ans[2][2] = m[2][2]*d[2];
478 }
479
480 /* ----------------------------------------------------------------------
481 row vector times matrix
482 ------------------------------------------------------------------------- */
483
484 KOKKOS_INLINE_FUNCTION
vecmat(const double * v,const double m[3][3],double * ans)485 void MathExtraKokkos::vecmat(const double *v, const double m[3][3], double *ans)
486 {
487 ans[0] = v[0]*m[0][0] + v[1]*m[1][0] + v[2]*m[2][0];
488 ans[1] = v[0]*m[0][1] + v[1]*m[1][1] + v[2]*m[2][1];
489 ans[2] = v[0]*m[0][2] + v[1]*m[1][2] + v[2]*m[2][2];
490 }
491
492 /* ----------------------------------------------------------------------
493 matrix times scalar, in place
494 ------------------------------------------------------------------------- */
495
496 KOKKOS_INLINE_FUNCTION
scalar_times3(const double f,double m[3][3])497 void MathExtraKokkos::scalar_times3(const double f, double m[3][3])
498 {
499 m[0][0] *= f; m[0][1] *= f; m[0][2] *= f;
500 m[1][0] *= f; m[1][1] *= f; m[1][2] *= f;
501 m[2][0] *= f; m[2][1] *= f; m[2][2] *= f;
502 }
503
504 /* ----------------------------------------------------------------------
505 compute quaternion from axis-angle rotation
506 v MUST be a unit vector
507 ------------------------------------------------------------------------- */
508
509 KOKKOS_INLINE_FUNCTION
axisangle_to_quat(const double * v,const double angle,double * quat)510 void MathExtraKokkos::axisangle_to_quat(const double *v, const double angle,
511 double *quat)
512 {
513 double halfa = 0.5*angle;
514 double sina = sin(halfa);
515 quat[0] = cos(halfa);
516 quat[1] = v[0]*sina;
517 quat[2] = v[1]*sina;
518 quat[3] = v[2]*sina;
519 }
520
521 /* ----------------------------------------------------------------------
522 compute rotation matrix from quaternion
523 quat = [w i j k]
524 ------------------------------------------------------------------------- */
525
526 KOKKOS_INLINE_FUNCTION
quat_to_mat(const double * quat,double mat[3][3])527 void MathExtraKokkos::quat_to_mat(const double *quat, double mat[3][3])
528 {
529 double w2 = quat[0]*quat[0];
530 double i2 = quat[1]*quat[1];
531 double j2 = quat[2]*quat[2];
532 double k2 = quat[3]*quat[3];
533 double twoij = 2.0*quat[1]*quat[2];
534 double twoik = 2.0*quat[1]*quat[3];
535 double twojk = 2.0*quat[2]*quat[3];
536 double twoiw = 2.0*quat[1]*quat[0];
537 double twojw = 2.0*quat[2]*quat[0];
538 double twokw = 2.0*quat[3]*quat[0];
539
540 mat[0][0] = w2+i2-j2-k2;
541 mat[0][1] = twoij-twokw;
542 mat[0][2] = twojw+twoik;
543
544 mat[1][0] = twoij+twokw;
545 mat[1][1] = w2-i2+j2-k2;
546 mat[1][2] = twojk-twoiw;
547
548 mat[2][0] = twoik-twojw;
549 mat[2][1] = twojk+twoiw;
550 mat[2][2] = w2-i2-j2+k2;
551 }
552
553 #endif
554