CMS 3D CMS Logo

LinearGridInterpolator3D.cc
Go to the documentation of this file.
2 #include "Grid1D.h"
3 #include "Grid3D.h"
5 
9 }
10 
12  /*
13  int i = grida.index(a);
14  int j = gridb.index(b);
15  int k = gridc.index(c);
16 
17  if (i==-1 || j==-1 || k==-1) {
18  // point outside of grid validity!
19  throwGridInterpolator3DException();
20  }
21 
22 
23  Scalar s = (a - grida.node(i)) / grida.step();
24  Scalar t = (b - gridb.node(j)) / gridb.step();
25  Scalar u = (c - gridc.node(k)) / gridc.step();
26 
27  */
28 
29  Scalar s, t, u;
30  int i = grida.index(a, s);
31  int j = gridb.index(b, t);
32  int k = gridc.index(c, u);
33 
34  // test range??
35 
36  grida.normalize(i, s);
37  gridb.normalize(j, t);
38  gridc.normalize(k, u);
39 
40 #ifdef DEBUG_LinearGridInterpolator3D
42  using std::cout;
43  using std::endl;
44  cout << "LinearGridInterpolator3D called with a,b,c " << a << "," << b << "," << c << endl;
45  cout << " i,j,k = " << i << "," << j << "," << k << " s,t,u = " << s << "," << t << "," << u << endl;
46  cout << "Node positions for " << i << "," << j << "," << k << " : " << grida.node(i) << "," << gridb.node(j) << ","
47  << gridc.node(k) << endl;
48  cout << "Node positions for " << i + 1 << "," << j + 1 << "," << k + 1 << " : " << grida.node(i + 1) << ","
49  << gridb.node(j + 1) << "," << gridc.node(k + 1) << endl;
50  cout << "Grid(" << i << "," << j << "," << k << ") = " << grid(i, j, k) << " ";
51  cout << "Grid(" << i << "," << j << "," << k + 1 << ") = " << grid(i, j, k + 1) << endl;
52  cout << "Grid(" << i << "," << j + 1 << "," << k << ") = " << grid(i, j + 1, k) << " ";
53  cout << "Grid(" << i << "," << j + 1 << "," << k + 1 << ") = " << grid(i, j + 1, k + 1) << endl;
54  cout << "Grid(" << i + 1 << "," << j << "," << k << ") = " << grid(i + 1, j, k) << " ";
55  cout << "Grid(" << i + 1 << "," << j << "," << k + 1 << ") = " << grid(i + 1, j, k + 1) << endl;
56  cout << "Grid(" << i + 1 << "," << j + 1 << "," << k << ") = " << grid(i + 1, j + 1, k) << " ";
57  cout << "Grid(" << i + 1 << "," << j + 1 << "," << k + 1 << ") = " << grid(i + 1, j + 1, k + 1) << endl;
58  cout << "number of nodes: " << grida.nodes() << "," << gridb.nodes() << "," << gridc.nodes() << endl;
59 
60  // cout << (1-s)*(1-t)*(1-u)*grid(i, j, k) << " " << (1-s)*(1-t)*u*grid(i, j, k+1) << endl
61  // << (1-s)* t *(1-u)*grid(i, j+1,k) << " " << (1-s)* t *u*grid(i, j+1,k+1) << endl
62  // << s *(1-t)*(1-u)*grid(i+1,j, k) << " " << s *(1-t)*u*grid(i+1,j, k+1) << endl
63  // << s * t *(1-u)*grid(i+1,j+1,k) << " " << s * t *u*grid(i+1,j+1,k+1) << endl;
64  }
65 
66 #endif
67 
68  int ind = grid.index(i, j, k);
69  int s1 = grid.stride1();
70  int s2 = grid.stride2();
71  int s3 = grid.stride3();
72  //chances are this is more numerically precise this way
73 
74  // this code for test to check properly inline of wrapped math...
75 #if defined(CMS_TEST_RAWSSE)
76 
77  __m128 resultSIMD =
78  _mm_mul_ps(_mm_set1_ps((1.f - s) * (1.f - t) * u), _mm_sub_ps(grid(ind + s3).v.vec, grid(ind).v.vec));
79  resultSIMD = _mm_add_ps(
80  resultSIMD,
81  _mm_mul_ps(_mm_set1_ps((1.f - s) * t * u), _mm_sub_ps(grid(ind + s2 + s3).v.vec, grid(ind + s2).v.vec)));
82  resultSIMD = _mm_add_ps(
83  resultSIMD,
84  _mm_mul_ps(_mm_set1_ps(s * (1.f - t) * u), _mm_sub_ps(grid(ind + s1 + s3).v.vec, grid(ind + s1).v.vec)));
85  resultSIMD = _mm_add_ps(
86  resultSIMD,
87  _mm_mul_ps(_mm_set1_ps(s * t * u), _mm_sub_ps(grid(ind + s1 + s2 + s3).v.vec, grid(ind + s1 + s2).v.vec)));
88  resultSIMD =
89  _mm_add_ps(resultSIMD, _mm_mul_ps(_mm_set1_ps((1.f - s) * t), _mm_sub_ps(grid(ind + s2).v.vec, grid(ind).v.vec)));
90  resultSIMD = _mm_add_ps(resultSIMD,
91  _mm_mul_ps(_mm_set1_ps(s * t), _mm_sub_ps(grid(ind + s1 + s2).v.vec, grid(ind + s1).v.vec)));
92  resultSIMD = _mm_add_ps(resultSIMD, _mm_mul_ps(_mm_set1_ps(s), _mm_sub_ps(grid(ind + s1).v.vec, grid(ind).v.vec)));
93  resultSIMD = _mm_add_ps(resultSIMD, grid(ind).v.vec);
94 
96  result.v = resultSIMD;
97 
98 #else
99 
100  ValueType result = ((1.f - s) * (1.f - t) * u) * (grid(ind + s3) - grid(ind));
101  result = result + ((1.f - s) * t * u) * (grid(ind + s2 + s3) - grid(ind + s2));
102  result = result + (s * (1.f - t) * u) * (grid(ind + s1 + s3) - grid(ind + s1));
103  result = result + (s * t * u) * (grid(ind + s1 + s2 + s3) - grid(ind + s1 + s2));
104  result = result + ((1.f - s) * t) * (grid(ind + s2) - grid(ind));
105  result = result + (s * t) * (grid(ind + s1 + s2) - grid(ind + s1));
106  result = result + (s) * (grid(ind + s1) - grid(ind));
107  result = result + grid(ind);
108 
109 #endif
110 
111  // (1-s)*(1-t)*(1-u)*grid(i, j, k) + (1-s)*(1-t)*u*grid(i, j, k+1) +
112  // (1-s)* t *(1-u)*grid(i, j+1,k) + (1-s)* t *u*grid(i, j+1,k+1) +
113  // s *(1-t)*(1-u)*grid(i+1,j, k) + s *(1-t)*u*grid(i+1,j, k+1) +
114  // s * t *(1-u)*grid(i+1,j+1,k) + s * t *u*grid(i+1,j+1,k+1);
115 
116  return result;
117 }
Scalar upper() const
Definition: Grid1D.h:20
void normalize(int &ind, Scalar &f) const
Definition: Grid1D.h:36
Scalar lower() const
Definition: Grid1D.h:19
int stride1() const
Definition: Grid3D.h:52
int index(Scalar a, Scalar &f) const
Definition: Grid1D.h:29
ReturnType interpolate(Scalar a, Scalar b, Scalar c)
Scalar node(int i) const
Definition: Grid1D.h:24
double f[11][100]
int stride2() const
Definition: Grid3D.h:53
static const bool debug
int stride3() const
Definition: Grid3D.h:54
int nodes() const
Definition: Grid1D.h:21
int index(int i, int j, int k) const
Definition: Grid3D.h:51
double b
Definition: hdecay.h:120
double a
Definition: hdecay.h:121