v0.8.23
Ddg_pointer.hpp
Go to the documentation of this file.
1 /* A version for pointers. */
2 
3 #pragma once
4 
5 namespace FTensor
6 {
7  template <class T, int Tensor_Dim01, int Tensor_Dim23>
8  class Ddg<T *, Tensor_Dim01, Tensor_Dim23>
9  {
10  const int inc;
11 
12  protected:
13 
14  mutable T *restrict data[(Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2]
15  [(Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2];
16 
17  public:
18  Ddg(T *d0000, T *d0010, T *d0011, T *d1000, T *d1010, T *d1011, T *d1100,
19  T *d1110, T *d1111, const int i = 1)
20  : inc(i)
21  {
22  ptr(0, 0, 0, 0) = d0000;
23  ptr(0, 0, 1, 0) = d0010;
24  ptr(0, 0, 1, 1) = d0011;
25  ptr(1, 0, 0, 0) = d1000;
26  ptr(1, 0, 1, 0) = d1010;
27  ptr(1, 0, 1, 1) = d1011;
28  ptr(1, 1, 0, 0) = d1100;
29  ptr(1, 1, 1, 0) = d1110;
30  ptr(1, 1, 1, 1) = d1111;
31  }
32 
33  Ddg(T *d0000, T *d0001, T *d0002, T *d0011, T *d0012, T *d0022, T *d0100,
34  T *d0101, T *d0102, T *d0111, T *d0112, T *d0122, T *d0200, T *d0201,
35  T *d0202, T *d0211, T *d0212, T *d0222, T *d1100, T *d1101, T *d1102,
36  T *d1111, T *d1112, T *d1122, T *d1200, T *d1201, T *d1202, T *d1211,
37  T *d1212, T *d1222, T *d2200, T *d2201, T *d2202, T *d2211, T *d2212,
38  T *d2222, const int i = 1)
39  : inc(i)
40  {
41  ptr(0, 0, 0, 0) = d0000;
42  ptr(0, 0, 0, 1) = d0001;
43  ptr(0, 0, 0, 2) = d0002;
44  ptr(0, 0, 1, 1) = d0011;
45  ptr(0, 0, 1, 2) = d0012;
46  ptr(0, 0, 2, 2) = d0022;
47  ptr(0, 1, 0, 0) = d0100;
48  ptr(0, 1, 0, 1) = d0101;
49  ptr(0, 1, 0, 2) = d0102;
50  ptr(0, 1, 1, 1) = d0111;
51  ptr(0, 1, 1, 2) = d0112;
52  ptr(0, 1, 2, 2) = d0122;
53  ptr(0, 2, 0, 0) = d0200;
54  ptr(0, 2, 0, 1) = d0201;
55  ptr(0, 2, 0, 2) = d0202;
56  ptr(0, 2, 1, 1) = d0211;
57  ptr(0, 2, 1, 2) = d0212;
58  ptr(0, 2, 2, 2) = d0222;
59  ptr(1, 1, 0, 0) = d1100;
60  ptr(1, 1, 0, 1) = d1101;
61  ptr(1, 1, 0, 2) = d1102;
62  ptr(1, 1, 1, 1) = d1111;
63  ptr(1, 1, 1, 2) = d1112;
64  ptr(1, 1, 2, 2) = d1122;
65  ptr(1, 2, 0, 0) = d1200;
66  ptr(1, 2, 0, 1) = d1201;
67  ptr(1, 2, 0, 2) = d1202;
68  ptr(1, 2, 1, 1) = d1211;
69  ptr(1, 2, 1, 2) = d1212;
70  ptr(1, 2, 2, 2) = d1222;
71  ptr(2, 2, 0, 0) = d2200;
72  ptr(2, 2, 0, 1) = d2201;
73  ptr(2, 2, 0, 2) = d2202;
74  ptr(2, 2, 1, 1) = d2211;
75  ptr(2, 2, 1, 2) = d2212;
76  ptr(2, 2, 2, 2) = d2222;
77  }
78 
79  /* Initializations for varying numbers of elements. */
80  template <class... U> Ddg(U *... d) : data{d...}, inc(1) {}
81 
82  /* There are two operator(int,int,int,int)'s, one for non-consts
83  that lets you change the value, and one for consts that
84  doesn't. */
85 
86  T &operator()(const int N1, const int N2, const int N3, const int N4)
87  {
88 #ifdef FTENSOR_DEBUG
89  if(N1 >= Tensor_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
90  || N3 >= Tensor_Dim23 || N3 < 0 || N4 >= Tensor_Dim23 || N4 < 0)
91  {
92  std::stringstream s;
93  s << "Bad index in Dg<T*," << Tensor_Dim01 << "," << Tensor_Dim23
94  << ">.operator(" << N1 << "," << N2 << "," << N3 << "," << N4
95  << ")" << std::endl;
96  throw std::out_of_range(s.str());
97  }
98 #endif
99  return N1 > N2
100  ? (N3 > N4 ? *data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2]
101  [N3 + (N4 * (2 * Tensor_Dim23 - N4 - 1)) / 2]
102  : *data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2]
103  [N4 + (N3 * (2 * Tensor_Dim23 - N3 - 1)) / 2])
104  : (N3 > N4
105  ? *data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2]
106  [N3 + (N4 * (2 * Tensor_Dim23 - N4 - 1)) / 2]
107  : *data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2]
108  [N4 + (N3 * (2 * Tensor_Dim23 - N3 - 1)) / 2]);
109  }
110 
111  T operator()(const int N1, const int N2, const int N3, const int N4) const
112  {
113 #ifdef FTENSOR_DEBUG
114  if(N1 >= Tensor_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
115  || N3 >= Tensor_Dim23 || N3 < 0 || N4 >= Tensor_Dim23 || N4 < 0)
116  {
117  std::stringstream s;
118  s << "Bad index in Dg<T*," << Tensor_Dim01 << "," << Tensor_Dim23
119  << ">.operator(" << N1 << "," << N2 << "," << N3 << "," << N4
120  << ") const" << std::endl;
121  throw std::out_of_range(s.str());
122  }
123 #endif
124  return N1 > N2
125  ? (N3 > N4 ? *data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2]
126  [N3 + (N4 * (2 * Tensor_Dim23 - N4 - 1)) / 2]
127  : *data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2]
128  [N4 + (N3 * (2 * Tensor_Dim23 - N3 - 1)) / 2])
129  : (N3 > N4
130  ? *data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2]
131  [N3 + (N4 * (2 * Tensor_Dim23 - N4 - 1)) / 2]
132  : *data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2]
133  [N4 + (N3 * (2 * Tensor_Dim23 - N3 - 1)) / 2]);
134  }
135 
136  T *restrict &ptr(const int N1, const int N2, const int N3, const int N4) const
137  {
138 #ifdef FTENSOR_DEBUG
139  if(N1 >= Tensor_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
140  || N3 >= Tensor_Dim23 || N3 < 0 || N4 >= Tensor_Dim23 || N4 < 0)
141  {
142  std::stringstream s;
143  s << "Bad index in Dg<T," << Tensor_Dim01 << "," << Tensor_Dim23
144  << ">.ptr(" << N1 << "," << N2 << "," << N3 << "," << N4 << ")"
145  << std::endl;
146  throw std::out_of_range(s.str());
147  }
148 #endif
149  return N1 > N2
150  ? (N3 > N4 ? data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2]
151  [N3 + (N4 * (2 * Tensor_Dim23 - N4 - 1)) / 2]
152  : data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2]
153  [N4 + (N3 * (2 * Tensor_Dim23 - N3 - 1)) / 2])
154  : (N3 > N4 ? data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2]
155  [N3 + (N4 * (2 * Tensor_Dim23 - N4 - 1)) / 2]
156  : data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2]
157  [N4 + (N3 * (2 * Tensor_Dim23 - N3 - 1)) / 2]);
158  }
159 
160  /* These operator()'s are the first part in constructing template
161  expressions. They can be used to slice off lower dimensional
162  parts. They are not entirely safe, since you can accidently use a
163  higher dimension than what is really allowed (like Dim=5). */
164 
165  template <char i, char j, char k, char l, int Dim01, int Dim23>
166  Ddg_Expr<Ddg<T *, Tensor_Dim01, Tensor_Dim23>, T, Dim01, Dim23, i, j, k, l>
167  operator()(const Index<i, Dim01> index1, const Index<j, Dim01> index2,
168  const Index<k, Dim23> index3, const Index<l, Dim23> index4)
169  {
170  return Ddg_Expr<Ddg<T *, Tensor_Dim01, Tensor_Dim23>, T, Dim01, Dim23, i,
171  j, k, l>(*this);
172  }
173 
174  template <char i, char j, char k, char l, int Dim01, int Dim23>
176  k, l>
177  operator()(const Index<i, Dim01> index1, const Index<j, Dim01> index2,
178  const Index<k, Dim23> index3,
179  const Index<l, Dim23> index4) const
180  {
182  Dim23, i, j, k, l>(*this);
183  }
184 
185  /* This is for expressions where a number is used for two slots, and
186  an index for the other two, yielding a Tensor2_symmetric_Expr. */
187 
188  template <char i, char j, int N0, int N1, int Dim>
191  Dim, i, j>
192  operator()(const Number<N0> n1, const Number<N1> n2,
193  const Index<i, Dim> index1, const Index<j, Dim> index2) const
194  {
195  using TensorExpr
198  TensorExpr(*this));
199  }
200 
201  template <char i, char j, int N0, int N1, int Dim>
204  Dim, i, j>
205  operator()(const Number<N0> n1, const Number<N1> n2,
206  const Index<i, Dim> index1, const Index<j, Dim> index2)
207  {
208  using TensorExpr
211  }
212 
213  /* This is for expressions where a number is used for one slot, and
214  an index for the other three, yielding a Dg_Expr. */
215 
216  template <char i, char j, char k, int N0, int Dim1, int Dim23>
218  Dim23, Dim1, i, j, k>
219  operator()(const Number<N0> n1, const Index<k, Dim1> index3,
220  const Index<i, Dim23> index1,
221  const Index<j, Dim23> index2) const
222  {
223  using TensorExpr
225  return Dg_Expr<TensorExpr, T, Dim23, Dim1, i, j, k>(TensorExpr(*this));
226  }
227 
228  template <char i, char j, char k, int N0, int Dim1, int Dim23>
230  Dim23, Dim1, i, j, k>
231  operator()(const Number<N0> n1, const Index<k, Dim1> index3,
232  const Index<i, Dim23> index1, const Index<j, Dim23> index2)
233  {
234  using TensorExpr
237  }
238 
239  /* This is for expressions where an int (not a Number) is used for
240  two slots, and an index for the other two, yielding a
241  Tensor2_symmetric_Expr. */
242 
243  template <char i, char j, int Dim>
246  j>
247  operator()(const int N0, const int N1, const Index<i, Dim> index1,
248  const Index<j, Dim> index2) const
249  {
250  using TensorExpr
253  TensorExpr(*this, N0, N1));
254  }
255 
256  template <char i, char j, int Dim>
259  operator()(const Index<i, Dim> index1, const Index<j, Dim> index2,
260  const int N2, const int N3) const
261  {
262  using TensorExpr
265  TensorExpr(*this, N2, N3));
266  }
267 
268  /* int in three slots, Index in the other yielding a Tensor1_Expr. */
269 
270  template <char i, int Dim>
272  Dim, i>
273  operator()(const Index<i, Dim> index1, const int N1, const int N2,
274  const int N3)
275  {
276  using TensorExpr
279  TensorExpr(*this, N1, N2, N3));
280  }
281 
282  template <char i, int Dim>
284  Dim, i>
285  operator()(const int N1, const Index<i, Dim> index1, const int N2,
286  const int N3)
287  {
288  using TensorExpr
291  TensorExpr(*this, N1, N2, N3));
292  }
293 
294  /* This is for expressions where an int (not a Number) is used for
295  one slot, and an index for the other three, yielding a
296  Dg_Expr. */
297 
298  template <char i, char j, char k, int Dim1, int Dim23>
300  Dim23, Dim1, i, j, k>
301  operator()(const int N0, const Index<k, Dim1> index3,
302  const Index<i, Dim23> index1,
303  const Index<j, Dim23> index2) const
304  {
305  using TensorExpr
308  TensorExpr(*this, N0));
309  }
310 
311  /* The ++ operator increments the pointer, not the number that the
312  pointer points to. This allows iterating over a grid. */
313 
314  const Ddg &operator++() const {
315  for(int i = 0; i < (Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2; ++i)
316  for(int j = 0; j < (Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2; ++j)
317  data[i][j] += inc;
318  return *this;
319  }
320  };
321 
322  template <class T, int Tensor_Dim01, int Tensor_Dim23, int I>
323  class Ddg<PackPtr<T *, I>, Tensor_Dim01, Tensor_Dim23>
325 
326  public:
327  /* Initializations for varying numbers of elements. */
328  template <class... U>
329  Ddg(U *... d) : Ddg<T *, Tensor_Dim01, Tensor_Dim23>(d...) {}
330 
331  /* The ++ operator increments the pointer, not the number that the
332  pointer points to. This allows iterating over a grid. */
333 
334  const Ddg &operator++() const {
335  for (int i = 0; i < (Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2; ++i)
336  for (int j = 0; j < (Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2; ++j)
338  return *this;
339  }
340  };
341 }
Tensor1_Expr< Ddg_numeral_123< Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T >, T, Dim, i > operator()(const int N1, const Index< i, Dim > index1, const int N2, const int N3)
Dg_Expr< Ddg_number_rhs_0< Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T, N0 >, T, Dim23, Dim1, i, j, k > operator()(const Number< N0 > n1, const Index< k, Dim1 > index3, const Index< i, Dim23 > index1, const Index< j, Dim23 > index2)
Dg_Expr< Ddg_number_0< const Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T, N0 >, T, Dim23, Dim1, i, j, k > operator()(const Number< N0 > n1, const Index< k, Dim1 > index3, const Index< i, Dim23 > index1, const Index< j, Dim23 > index2) const
Ddg(T *d0000, T *d0001, T *d0002, T *d0011, T *d0012, T *d0022, T *d0100, T *d0101, T *d0102, T *d0111, T *d0112, T *d0122, T *d0200, T *d0201, T *d0202, T *d0211, T *d0212, T *d0222, T *d1100, T *d1101, T *d1102, T *d1111, T *d1112, T *d1122, T *d1200, T *d1201, T *d1202, T *d1211, T *d1212, T *d1222, T *d2200, T *d2201, T *d2202, T *d2211, T *d2212, T *d2222, const int i=1)
Definition: Ddg_pointer.hpp:33
T & operator()(const int N1, const int N2, const int N3, const int N4)
Definition: Ddg_pointer.hpp:86
Fully Antisymmetric Levi-Civita Tensor.
Dg_Expr< Ddg_numeral_0< const Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T >, T, Dim23, Dim1, i, j, k > operator()(const int N0, const Index< k, Dim1 > index3, const Index< i, Dim23 > index1, const Index< j, Dim23 > index2) const
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
Ddg_Expr< Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T, Dim01, Dim23, i, j, k, l > operator()(const Index< i, Dim01 > index1, const Index< j, Dim01 > index2, const Index< k, Dim23 > index3, const Index< l, Dim23 > index4)
Tensor2_symmetric_Expr< Ddg_numeral_23< Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T >, T, Dim, i, j > operator()(const Index< i, Dim > index1, const Index< j, Dim > index2, const int N2, const int N3) const
Ddg_Expr< const Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T, Dim01, Dim23, i, j, k, l > operator()(const Index< i, Dim01 > index1, const Index< j, Dim01 > index2, const Index< k, Dim23 > index3, const Index< l, Dim23 > index4) const
T operator()(const int N1, const int N2, const int N3, const int N4) const
Tensor2_symmetric_Expr< Ddg_number_rhs_01< Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T, N0, N1 >, T, Dim, i, j > operator()(const Number< N0 > n1, const Number< N1 > n2, const Index< i, Dim > index1, const Index< j, Dim > index2)
T *restrict & ptr(const int N1, const int N2, const int N3, const int N4) const
Tensor1_Expr< Ddg_numeral_123< Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T >, T, Dim, i > operator()(const Index< i, Dim > index1, const int N1, const int N2, const int N3)
Tensor2_symmetric_Expr< Ddg_numeral_01< const Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T >, T, Dim, i, j > operator()(const int N0, const int N1, const Index< i, Dim > index1, const Index< j, Dim > index2) const
Tensor2_symmetric_Expr< Ddg_number_01< const Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T, N0, N1 >, T, Dim, i, j > operator()(const Number< N0 > n1, const Number< N1 > n2, const Index< i, Dim > index1, const Index< j, Dim > index2) const
Ddg(T *d0000, T *d0010, T *d0011, T *d1000, T *d1010, T *d1011, T *d1100, T *d1110, T *d1111, const int i=1)
Definition: Ddg_pointer.hpp:18