v0.9.0
Dg_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_Dim2>
8  class Dg<T *, Tensor_Dim01, Tensor_Dim2>
9  {
10 
11  protected:
12 
13  mutable T *restrict
14  data[(Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2][Tensor_Dim2];
15 
16  public:
17 
18  template <class... U> Dg(U *... d) : data{d...}
19  {
20  static_assert(sizeof...(d) == sizeof(data) / sizeof(T),
21  "Incorrect number of Arguments. Constructor should "
22  "initialize the entire Tensor");
23  }
24 
25  Dg() {}
26 
27  /* There are two operator(int,int,int)'s, one for non-consts that lets you
28  change the value, and one for consts that doesn't. */
29 
30  T &operator()(const int N1, const int N2, const int N3)
31  {
32 #ifdef FTENSOR_DEBUG
33  if(N1 >= Tensor_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
34  || N3 >= Tensor_Dim2 || N3 < 0)
35  {
36  std::stringstream s;
37  s << "Bad index in Dg<T*," << Tensor_Dim01 << "," << Tensor_Dim2
38  << ">.operator(" << N1 << "," << N2 << "," << N3 << ")"
39  << std::endl;
40  throw std::out_of_range(s.str());
41  }
42 #endif
43  return N1 > N2 ? *data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2][N3]
44  : *data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2][N3];
45  }
46 
47  T operator()(const int N1, const int N2, const int N3) const
48  {
49 #ifdef FTENSOR_DEBUG
50  if(N1 >= Tensor_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
51  || N3 >= Tensor_Dim2 || N3 < 0)
52  {
53  std::stringstream s;
54  s << "Bad index in Dg<T*," << Tensor_Dim01 << "," << Tensor_Dim2
55  << ">.operator(" << N1 << "," << N2 << "," << N3 << ") const"
56  << std::endl;
57  throw std::out_of_range(s.str());
58  }
59 #endif
60  return N1 > N2 ? *data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2][N3]
61  : *data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2][N3];
62  }
63 
64  T *ptr(const int N1, const int N2, const int N3) const
65  {
66 #ifdef FTENSOR_DEBUG
67  if(N1 >= Tensor_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
68  || N3 >= Tensor_Dim2 || N3 < 0)
69  {
70  std::stringstream s;
71  s << "Bad index in Dg<T," << Tensor_Dim01 << "," << Tensor_Dim2
72  << ">.ptr(" << N1 << "," << N2 << "," << N3 << ")" << std::endl;
73  throw std::out_of_range(s.str());
74  }
75 #endif
76  return N1 > N2 ? data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2][N3]
77  : data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2][N3];
78  }
79 
80  /* These operator()'s are the first part in constructing template
81  expressions. */
82 
83  template <char i, char j, char k, int Dim01, int Dim2>
84  Dg_Expr<Dg<T *, Tensor_Dim01, Tensor_Dim2>, T, Dim01, Dim2, i, j, k>
85  operator()(const Index<i, Dim01> index1, const Index<j, Dim01> index2,
86  const Index<k, Dim2> index3)
87  {
88  return Dg_Expr<Dg<T *, Tensor_Dim01, Tensor_Dim2>, T, Dim01, Dim2, i, j,
89  k>(*this);
90  }
91 
92  template <char i, char j, char k, int Dim01, int Dim2>
94  operator()(const Index<i, Dim01> index1, const Index<j, Dim01> index2,
95  const Index<k, Dim2> index3) const
96  {
98  i, j, k>(*this);
99  }
100 
101  /* These operators are for internal contractions. The commented out
102  versions are more general, but are ambiguous. This means,
103  unfortunately, that you can't do something like A(i,j,j) where i
104  and j have different dimensions. */
105 
106  template <char i, char j, int Dim>
107  Tensor1_Expr<
109  T, Dim, i>
110  operator()(const Index<i, Dim> index1, const Index<j, Dim> index2,
111  const Index<j, Dim> index3) const
112  {
113  using TensorExpr
115  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
116  }
117 
118  template <char i, char j, int Dim>
119  Tensor1_Expr<
121  T, Dim, i>
122  operator()(const Index<j, Dim> index1, const Index<i, Dim> index2,
123  const Index<j, Dim> index3) const
124  {
125  using TensorExpr
127  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
128  }
129 
130  // template<char i, char j, int Dim0, int Dim12>
131  // Tensor1_Expr<const Tensor3_contracted_12<Dg
132  // <T*,Tensor_Dim01,Tensor_Dim2>,T,Dim12>,T,Dim0,i>
133  // operator()(const Index<i,Dim0> index1, const Index<j,Dim12> index2,
134  // const Index<j,Dim12> index3) const
135  // {
136  // typedef Tensor3_contracted_12<Dg<T*,Tensor_Dim01,Tensor_Dim2>,
137  // T,Dim12> TensorExpr;
138  // return Tensor1_Expr<TensorExpr,T,Dim0,i>(TensorExpr(*this));
139  // }
140 
141  // template<char i, char j, int Dim02, int Dim1>
142  // Tensor1_Expr<const Tensor3_contracted_02<Dg
143  // <T*,Tensor_Dim01,Tensor_Dim2>,T,Dim02>,T,Dim1,i>
144  // operator()(const Index<j,Dim02> index1, const Index<i,Dim1> index2,
145  // const Index<j,Dim02> index3) const
146  // {
147  // typedef Tensor3_contracted_02<Dg<T*,Tensor_Dim01,Tensor_Dim2>,
148  // T,Dim02> TensorExpr;
149  // return Tensor1_Expr<TensorExpr,T,Dim1,i>(TensorExpr(*this));
150  // }
151 
152  template <char i, char j, int Dim01, int Dim2>
153  Tensor1_Expr<
155  T, Dim2, i>
156  operator()(const Index<j, Dim01> index1, const Index<j, Dim01> index2,
157  const Index<i, Dim2> index3) const
158  {
159  using TensorExpr
161  return Tensor1_Expr<TensorExpr, T, Dim2, i>(TensorExpr(*this));
162  }
163 
164  /* This is for expressions where a number is used for one slot, and
165  indices for the others, yielding a Tensor2_Expr or
166  Tensor2_symmetric_Expr. The non-const versions don't actually
167  create a Dg_number_rhs_* object, while the const versions
168  do create a Dg_number_*. */
169 
170  /* First slot. */
171 
172  template <char i, char j, int N, int Dim1, int Dim2>
174  Dim1, Dim2, i, j>
175  operator()(const Number<N> n1, const Index<i, Dim1> index1,
176  const Index<j, Dim2> index2)
177  {
178  using TensorExpr
181  }
182 
183  template <char i, char j, int N, int Dim1, int Dim2>
184  const Tensor2_Expr<
186  Dim1, Dim2, i, j>
187  operator()(const Number<N> n1, const Index<i, Dim1> index1,
188  const Index<j, Dim2> index2) const
189  {
190  using TensorExpr
192  return Tensor2_Expr<TensorExpr, T, Dim1, Dim2, i, j>(TensorExpr(*this));
193  }
194 
195  /* Second slot. */
196 
197  template <char i, char j, int N, int Dim0, int Dim2>
199  Dim0, Dim2, i, j>
200  operator()(const Index<i, Dim0> index1, const Number<N> n1,
201  const Index<j, Dim2> index2)
202  {
203  using TensorExpr
206  }
207 
208  template <char i, char j, int N, int Dim0, int Dim2>
209  const Tensor2_Expr<
211  Dim0, Dim2, i, j>
212  operator()(const Index<i, Dim0> index1, const Number<N> n1,
213  const Index<j, Dim2> index2) const
214  {
215  using TensorExpr
217  return Tensor2_Expr<TensorExpr, T, Dim0, Dim2, i, j>(TensorExpr(*this));
218  }
219 
220  /* Third slot. */
221 
222  template <char i, char j, int N, int Dim>
225  operator()(const Index<i, Dim> index1, const Index<j, Dim> index2,
226  const Number<N> n1)
227  {
228  using TensorExpr
231  }
232 
233  template <char i, char j, int N, int Dim>
236  Dim, i, j>
237  operator()(const Index<i, Dim> index1, const Index<j, Dim> index2,
238  const Number<N> n1) const
239  {
240  using TensorExpr
243  TensorExpr(*this));
244  }
245 
246  /* This is for expressions where a number is used for two slots, and
247  an Index for the other, yielding a Tensor1_Expr. The non-const
248  versions don't actually create a Dg_number_rhs_* object,
249  while the const versions do create a Dg_number_*. */
250 
251  /* Index in first slot. */
252 
253  template <char i, int N1, int N2, int Dim>
255  T, Dim, i>
256  operator()(const Index<i, Dim> index, const Number<N1> n1,
257  const Number<N2> n2)
258  {
259  using TensorExpr
262  }
263 
264  template <char i, int N1, int N2, int Dim>
265  const Tensor1_Expr<
267  T, Dim, i>
268  operator()(const Index<i, Dim> index, const Number<N1> n1,
269  const Number<N2> n2) const
270  {
271  using TensorExpr
273  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
274  }
275 
276  /* Index in second slot. I use the same structures as for the Index
277  in the first slot since the tensor is symmetric on the first two
278  indices. */
279 
280  template <char i, int N1, int N2, int Dim>
282  T, Dim, i>
283  operator()(const Number<N1> n1, const Index<i, Dim> index,
284  const Number<N2> n2)
285  {
286  using TensorExpr
289  }
290 
291  template <char i, int N1, int N2, int Dim>
292  const Tensor1_Expr<
294  T, Dim, i>
295  operator()(const Number<N1> n1, const Index<i, Dim> index,
296  const Number<N2> n2) const
297  {
298  using TensorExpr
300  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
301  }
302 
303  /* Index in third slot. */
304 
305  template <char i, int N1, int N2, int Dim>
307  T, Dim, i>
308  operator()(const Number<N1> n1, const Number<N2> n2,
309  const Index<i, Dim> index)
310  {
311  using TensorExpr
314  }
315 
316  template <char i, int N1, int N2, int Dim>
317  const Tensor1_Expr<
319  T, Dim, i>
320  operator()(const Number<N1> n1, const Number<N2> n2,
321  const Index<i, Dim> index) const
322  {
323  using TensorExpr
325  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
326  }
327 
328  /* Specializations for using actual numbers instead of Number<>.
329  This is for expressions where an actual number is used for one
330  slot, and indices for the others, yielding a Tensor2_Expr or
331  Tensor2_symmetric_Expr. I only define the const versions. */
332 
333  /* First slot. */
334 
335  template <char i, char j, int Dim1, int Dim2>
336  const Tensor2_Expr<
338  Dim2, i, j>
339  operator()(const int N, const Index<i, Dim1> index1,
340  const Index<j, Dim2> index2) const
341  {
342  using TensorExpr
345  TensorExpr(*this, N));
346  }
347 
348  /* Second slot. */
349 
350  template <char i, char j, int Dim0, int Dim2>
351  const Tensor2_Expr<
353  Dim2, i, j>
354  operator()(const Index<i, Dim0> index1, const int N,
355  const Index<j, Dim2> index2) const
356  {
357  using TensorExpr
360  TensorExpr(*this, N));
361  }
362 
363  /* Third slot. */
364 
365  template <char i, char j, int Dim>
368  i, j>
369  operator()(const Index<i, Dim> index1, const Index<j, Dim> index2,
370  const int N) const
371  {
372  using TensorExpr
375  TensorExpr(*this, N));
376  }
377 
378  /* This is for expressions where a numeral is used for two slots, and
379  an Index for the other, yielding a Tensor1_Expr. */
380 
381  /* Index in first slot. */
382 
383  template <char i, int Dim>
384  const Tensor1_Expr<
386  i>
387  operator()(const Index<i, Dim> index, const int N1, const int N2) const
388  {
389  using TensorExpr
391  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
392  }
393 
394  /* Index in second slot. I use the same structures as for the Index
395  in the first slot since the tensor is symmetric on the first two
396  indices. */
397 
398  template <char i, int Dim>
399  const Tensor1_Expr<
401  i>
402  operator()(const int N1, const Index<i, Dim> index, const int N2) const
403  {
404  using TensorExpr
406  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
407  }
408 
409  /* Index in third slot. */
410 
411  template <char i, int Dim>
412  const Tensor1_Expr<
414  i>
415  operator()(const int N1, const int N2, const Index<i, Dim> index) const
416  {
417  using TensorExpr
419  return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
420  }
421 
422  /* The ++ operator increments the pointer, not the number that the
423  pointer points to. This allows iterating over a grid. */
424 
426  {
427  for(int i = 0; i < (Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2; ++i)
428  for(int j = 0; j < Tensor_Dim2; ++j)
429  ++data[i][j];
430  return *this;
431  }
432  };
433 
434  template <class T, int Tensor_Dim01, int Tensor_Dim2, int I>
435  class Dg<PackPtr<T *, I>, Tensor_Dim01, Tensor_Dim2>
437 
438  public:
439  template <class... U>
440  Dg(U *... d) : Dg<T *, Tensor_Dim01, Tensor_Dim2>(d...) {}
441 
442  Dg() : Dg<T *, Tensor_Dim01, Tensor_Dim2>() {}
443 
444  /* The ++ operator increments the pointer, not the number that the
445  pointer points to. This allows iterating over a grid. */
446 
447  const Dg &operator++() const {
448  for (int i = 0; i < (Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2; ++i)
449  for (int j = 0; j < Tensor_Dim2; ++j)
451  return *this;
452  }
453  };
454  } // namespace FTensor
const Tensor1_Expr< const Dg_number_12< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, N1, N2 >, T, Dim, i > operator()(const Index< i, Dim > index, const Number< N1 > n1, const Number< N2 > n2) const
Definition: Dg_pointer.hpp:268
const Tensor2_Expr< const Dg_numeral_0< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T >, T, Dim0, Dim2, i, j > operator()(const Index< i, Dim0 > index1, const int N, const Index< j, Dim2 > index2) const
Definition: Dg_pointer.hpp:354
const Tensor2_symmetric_Expr< const Dg_number_2< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, N >, T, Dim, i, j > operator()(const Index< i, Dim > index1, const Index< j, Dim > index2, const Number< N > n1) const
Definition: Dg_pointer.hpp:237
const Dg< T *, Tensor_Dim01, Tensor_Dim2 > & operator++() const
Definition: Dg_pointer.hpp:425
Tensor2_Expr< Dg_number_rhs_0< Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, N >, T, Dim1, Dim2, i, j > operator()(const Number< N > n1, const Index< i, Dim1 > index1, const Index< j, Dim2 > index2)
Definition: Dg_pointer.hpp:175
Tensor1_Expr< const Tensor3_contracted_01< Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, Dim01 >, T, Dim2, i > operator()(const Index< j, Dim01 > index1, const Index< j, Dim01 > index2, const Index< i, Dim2 > index3) const
Definition: Dg_pointer.hpp:156
Tensor1_Expr< Dg_number_rhs_12< Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, N1, N2 >, T, Dim, i > operator()(const Number< N1 > n1, const Index< i, Dim > index, const Number< N2 > n2)
Definition: Dg_pointer.hpp:283
Dg_Expr< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, Dim01, Dim2, i, j, k > operator()(const Index< i, Dim01 > index1, const Index< j, Dim01 > index2, const Index< k, Dim2 > index3) const
Definition: Dg_pointer.hpp:94
Tensor2_Expr< Dg_number_rhs_0< Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, N >, T, Dim0, Dim2, i, j > operator()(const Index< i, Dim0 > index1, const Number< N > n1, const Index< j, Dim2 > index2)
Definition: Dg_pointer.hpp:200
const Tensor2_Expr< const Dg_numeral_0< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T >, T, Dim1, Dim2, i, j > operator()(const int N, const Index< i, Dim1 > index1, const Index< j, Dim2 > index2) const
Definition: Dg_pointer.hpp:339
T operator()(const int N1, const int N2, const int N3) const
Definition: Dg_pointer.hpp:47
Tensor1_Expr< const Tensor3_contracted_12< Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, Dim >, T, Dim, i > operator()(const Index< i, Dim > index1, const Index< j, Dim > index2, const Index< j, Dim > index3) const
Definition: Dg_pointer.hpp:110
Fully Antisymmetric Levi-Civita Tensor.
const Tensor1_Expr< const Dg_number_01< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, N1, N2 >, T, Dim, i > operator()(const Number< N1 > n1, const Number< N2 > n2, const Index< i, Dim > index) const
Definition: Dg_pointer.hpp:320
Tensor1_Expr< Dg_number_rhs_12< Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, N1, N2 >, T, Dim, i > operator()(const Index< i, Dim > index, const Number< N1 > n1, const Number< N2 > n2)
Definition: Dg_pointer.hpp:256
Tensor1_Expr< const Tensor3_contracted_02< Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, Dim >, T, Dim, i > operator()(const Index< j, Dim > index1, const Index< i, Dim > index2, const Index< j, Dim > index3) const
Definition: Dg_pointer.hpp:122
Tensor1_Expr< Dg_number_rhs_01< Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, N1, N2 >, T, Dim, i > operator()(const Number< N1 > n1, const Number< N2 > n2, const Index< i, Dim > index)
Definition: Dg_pointer.hpp:308
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
T * ptr(const int N1, const int N2, const int N3) const
Definition: Dg_pointer.hpp:64
const Tensor2_symmetric_Expr< const Dg_numeral_2< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T >, T, Dim, i, j > operator()(const Index< i, Dim > index1, const Index< j, Dim > index2, const int N) const
Definition: Dg_pointer.hpp:369
const Tensor2_Expr< const Dg_number_0< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, N >, T, Dim1, Dim2, i, j > operator()(const Number< N > n1, const Index< i, Dim1 > index1, const Index< j, Dim2 > index2) const
Definition: Dg_pointer.hpp:187
Dg_Expr< Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, Dim01, Dim2, i, j, k > operator()(const Index< i, Dim01 > index1, const Index< j, Dim01 > index2, const Index< k, Dim2 > index3)
Definition: Dg_pointer.hpp:85
T & operator()(const int N1, const int N2, const int N3)
Definition: Dg_pointer.hpp:30
const Tensor1_Expr< const Dg_number_12< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, N1, N2 >, T, Dim, i > operator()(const Number< N1 > n1, const Index< i, Dim > index, const Number< N2 > n2) const
Definition: Dg_pointer.hpp:295
const Tensor1_Expr< const Dg_numeral_01< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T >, T, Dim, i > operator()(const int N1, const int N2, const Index< i, Dim > index) const
Definition: Dg_pointer.hpp:415
Tensor2_symmetric_Expr< Dg_number_rhs_2< Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, N >, T, Dim, i, j > operator()(const Index< i, Dim > index1, const Index< j, Dim > index2, const Number< N > n1)
Definition: Dg_pointer.hpp:225
const int N
Definition: speed_test.cpp:3
const Tensor1_Expr< const Dg_numeral_12< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T >, T, Dim, i > operator()(const int N1, const Index< i, Dim > index, const int N2) const
Definition: Dg_pointer.hpp:402
const Tensor1_Expr< const Dg_numeral_12< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T >, T, Dim, i > operator()(const Index< i, Dim > index, const int N1, const int N2) const
Definition: Dg_pointer.hpp:387
const Tensor2_Expr< const Dg_number_0< const Dg< T *, Tensor_Dim01, Tensor_Dim2 >, T, N >, T, Dim0, Dim2, i, j > operator()(const Index< i, Dim0 > index1, const Number< N > n1, const Index< j, Dim2 > index2) const
Definition: Dg_pointer.hpp:212