v0.9.0
Christof_pointer.hpp
Go to the documentation of this file.
1 /* A version for pointers */
2 
3 #pragma once
4 
5 #include "../Tensor3.hpp"
6 
7 namespace FTensor
8 {
9  template <class T, int Tensor_Dim0, int Tensor_Dim12>
10  class Christof<T *, Tensor_Dim0, Tensor_Dim12>
11  {
12  const int inc;
13 
14  protected:
15 
16  mutable T *restrict
17  data[Tensor_Dim0][(Tensor_Dim12 * (Tensor_Dim12 + 1)) / 2];
18 
19  public:
20  template <class... U> explicit Christof(U *... d) : inc(1), data{d...} {
21  static_assert(sizeof...(d) == sizeof(data) / sizeof(T),
22  "Incorrect number of Arguments. Constructor should "
23  "initialize the entire Tensor");
24  }
25 
26  template <class... U>
27  explicit Christof(const int i, U *... d) : data{d...}, inc(i) {
28  static_assert(sizeof...(d) == sizeof(data) / sizeof(T),
29  "Incorrect number of Arguments. Constructor should "
30  "initialize the entire Tensor");
31  }
32 
33  Christof() {}
34 
35  /* There are two operator(int,int,int)'s, one for non-consts that lets you
36  change the value, and one for consts that doesn't. */
37 
38  T &operator()(const int N1, const int N2, const int N3)
39  {
40 #ifdef FTENSOR_DEBUG
41  if(N1 >= Tensor_Dim0 || N1 < 0 || N2 >= Tensor_Dim12 || N2 < 0
42  || N3 >= Tensor_Dim12 || N3 < 0)
43  {
44  std::stringstream s;
45  s << "Bad index in Christof<T*," << Tensor_Dim0 << ","
46  << Tensor_Dim12 << ">.operator(" << N1 << "," << N2 << "," << N3
47  << ")" << std::endl;
48  throw std::out_of_range(s.str());
49  }
50 #endif
51  return N2 > N3 ? *data[N1][N2 + (N3 * (2 * Tensor_Dim12 - N3 - 1)) / 2]
52  : *data[N1][N3 + (N2 * (2 * Tensor_Dim12 - N2 - 1)) / 2];
53  }
54 
55  T operator()(const int N1, const int N2, const int N3) const
56  {
57 #ifdef FTENSOR_DEBUG
58  if(N1 >= Tensor_Dim0 || N1 < 0 || N2 >= Tensor_Dim12 || N2 < 0
59  || N3 >= Tensor_Dim12 || N3 < 0)
60  {
61  std::stringstream s;
62  s << "Bad index in Christof<T*," << Tensor_Dim0 << ","
63  << Tensor_Dim12 << ">.operator(" << N1 << "," << N2 << "," << N3
64  << ") const" << std::endl;
65  throw std::out_of_range(s.str());
66  }
67 #endif
68  return N2 > N3 ? *data[N1][N2 + (N3 * (2 * Tensor_Dim12 - N3 - 1)) / 2]
69  : *data[N1][N3 + (N2 * (2 * Tensor_Dim12 - N2 - 1)) / 2];
70  }
71 
72  T *ptr(const int N1, const int N2, const int N3) const
73  {
74 #ifdef FTENSOR_DEBUG
75  if(N1 >= Tensor_Dim0 || N1 < 0 || N2 >= Tensor_Dim12 || N2 < 0
76  || N3 >= Tensor_Dim12 || N3 < 0)
77  {
78  std::stringstream s;
79  s << "Bad index in Christof<T*," << Tensor_Dim0 << ","
80  << Tensor_Dim12 << ">.ptr(" << N1 << "," << N2 << "," << N3 << ")"
81  << std::endl;
82  throw std::out_of_range(s.str());
83  }
84 #endif
85  return N2 > N3 ? data[N1][N2 + (N3 * (2 * Tensor_Dim12 - N3 - 1)) / 2]
86  : data[N1][N3 + (N2 * (2 * Tensor_Dim12 - N2 - 1)) / 2];
87  }
88 
89  /* These operator()'s are the first part in constructing template
90  expressions. I mix up the indices here so that it behaves like a
91  Dg. That way I don't have to have a separate wrapper
92  class Christof_Expr, which simplifies things. */
93 
94  template <char i, char j, char k, int Dim0, int Dim12>
96  operator()(const Index<k, Dim0> index1, const Index<i, Dim12> index2,
97  const Index<j, Dim12> index3)
98  {
100  i, j, k>(*this);
101  }
102 
103  template <char i, char j, char k, int Dim0, int Dim12>
105  j, k>
106  operator()(const Index<k, Dim0> index1, const Index<i, Dim12> index2,
107  const Index<j, Dim12> index3) const
108  {
110  Dim0, i, j, k>(*this);
111  }
112 
113  /* These operators are for internal contractions. */
114 
115  /* const versions */
116 
117  template <char i, char j, int Dim0, int Dim12>
120  T, Dim0, i>
121  operator()(const Index<i, Dim0> index1, const Index<j, Dim12> index2,
122  const Index<j, Dim12> index3) const
123  {
124  using TensorExpr
126  T, Dim12>;
127  return Tensor1_Expr<TensorExpr, T, Dim0, i>(TensorExpr(*this));
128  }
129 
130  template <char i, char j, int Dim02, int Dim1>
133  T, Dim1, i>
134  operator()(const Index<j, Dim02> index1, const Index<i, Dim1> index2,
135  const Index<j, Dim02> index3) const
136  {
137  using TensorExpr
139  T, Dim02>;
140  return Tensor1_Expr<TensorExpr, T, Dim1, i>(TensorExpr(*this));
141  }
142 
143  template <char i, char j, int Dim01, int Dim2>
146  T, Dim2, i>
147  operator()(const Index<j, Dim01> index1, const Index<j, Dim01> index2,
148  const Index<i, Dim2> index3) const
149  {
150  using TensorExpr
152  T, Dim01>;
153  return Tensor1_Expr<TensorExpr, T, Dim2, i>(TensorExpr(*this));
154  }
155 
156  /* non-const versions */
157 
158  template <char i, char j, int Dim0, int Dim12>
161  T, Dim0, i>
162  operator()(const Index<i, Dim0> index1, const Index<j, Dim12> index2,
163  const Index<j, Dim12> index3)
164  {
165  using TensorExpr
167  T, Dim12>;
168  return Tensor1_Expr<TensorExpr, T, Dim0, i>(TensorExpr(*this));
169  }
170 
171  template <char i, char j, int Dim02, int Dim1>
174  T, Dim1, i>
175  operator()(const Index<j, Dim02> index1, const Index<i, Dim1> index2,
176  const Index<j, Dim02> index3)
177  {
178  using TensorExpr
180  T, Dim02>;
181  return Tensor1_Expr<TensorExpr, T, Dim1, i>(TensorExpr(*this));
182  }
183 
184  template <char i, char j, int Dim01, int Dim2>
187  T, Dim2, i>
188  operator()(const Index<j, Dim01> index1, const Index<j, Dim01> index2,
189  const Index<i, Dim2> index3)
190  {
191  using TensorExpr
193  T, Dim01>;
194  return Tensor1_Expr<TensorExpr, T, Dim2, i>(TensorExpr(*this));
195  }
196 
197  /* This is for expressions where a number is used for one slot, and
198  an index for the others, yielding a Tensor2_symmetric_Expr. */
199 
200  template <char i, char j, int N, int Dim12>
203  T, N>,
204  T, Dim12, i, j>
205  operator()(const Number<N> n1, const Index<i, Dim12> index1,
206  const Index<j, Dim12> index2) const
207  {
208  using TensorExpr
210  N>;
212  TensorExpr(*this));
213  }
214 
215  /* An int in one spot, Index for the others, yielding a Tensor2. I
216  can use the same structure for both, since Christof is
217  symmetric on the last two indices. */
218 
219  template <char i, char j, int Dim0, int Dim2>
222  T, Dim0, Dim2, i, j>
223  operator()(const Index<i, Dim0> index1, const int N,
224  const Index<j, Dim2> index2) const
225  {
226  using TensorExpr
229  TensorExpr(*this, N));
230  }
231 
232  template <char i, char j, int Dim0, int Dim2>
235  T, Dim0, Dim2, i, j>
236  operator()(const Index<i, Dim0> index1, const Index<j, Dim2> index2,
237  const int N) const
238  {
239  using TensorExpr
242  TensorExpr(*this, N));
243  }
244 
245  /* The ++ operator increments the pointer, not the number that the
246  pointer points to. This allows iterating over a grid. */
247 
249  {
250  for(int i = 0; i < Tensor_Dim0; ++i)
251  for(int j = 0; j < (Tensor_Dim12 * (Tensor_Dim12 + 1)) / 2; ++j)
252  data[i][j] += inc;
253  return *this;
254  }
255  };
256 
257  template <class T, int Tensor_Dim0, int Tensor_Dim12, int I>
258  class Christof<PackPtr<T *, I>, Tensor_Dim0, Tensor_Dim12>:
260  {
261 
262  public:
263  template <class... U>
264  Christof(U *... d) : Christof<T *, Tensor_Dim0, Tensor_Dim12>(d...) {}
265 
266  Christof() : Christof<T *, Tensor_Dim0, Tensor_Dim12>() {}
267 
268  /* The ++ operator increments the pointer, not the number that the
269  pointer points to. This allows iterating over a grid. */
270 
271  const Christof<PackPtr<T *, I>, Tensor_Dim0, Tensor_Dim12> &
272  operator++() const {
273  for (int i = 0; i < Tensor_Dim0; ++i)
274  for (int j = 0; j < (Tensor_Dim12 * (Tensor_Dim12 + 1)) / 2; ++j)
276  return *this;
277  }
278  };
279 }
const Christof< T *, Tensor_Dim0, Tensor_Dim12 > & operator++() const
T & operator()(const int N1, const int N2, const int N3)
Tensor1_Expr< const Tensor3_contracted_01< const Christof< T *, Tensor_Dim0, Tensor_Dim12 >, T, Dim01 >, T, Dim2, i > operator()(const Index< j, Dim01 > index1, const Index< j, Dim01 > index2, const Index< i, Dim2 > index3)
T operator()(const int N1, const int N2, const int N3) const
Tensor1_Expr< const Tensor3_contracted_02< const Christof< T *, Tensor_Dim0, Tensor_Dim12 >, T, Dim02 >, T, Dim1, i > operator()(const Index< j, Dim02 > index1, const Index< i, Dim1 > index2, const Index< j, Dim02 > index3)
Fully Antisymmetric Levi-Civita Tensor.
Tensor1_Expr< const Tensor3_contracted_01< const Christof< T *, Tensor_Dim0, Tensor_Dim12 >, T, Dim01 >, T, Dim2, i > operator()(const Index< j, Dim01 > index1, const Index< j, Dim01 > index2, const Index< i, Dim2 > index3) const
Tensor2_Expr< const Christof_numeral_1< const Christof< T *, Tensor_Dim0, Tensor_Dim12 >, T >, T, Dim0, Dim2, i, j > operator()(const Index< i, Dim0 > index1, const Index< j, Dim2 > index2, const int N) 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
T * ptr(const int N1, const int N2, const int N3) const
Tensor1_Expr< const Tensor3_contracted_12< const Christof< T *, Tensor_Dim0, Tensor_Dim12 >, T, Dim12 >, T, Dim0, i > operator()(const Index< i, Dim0 > index1, const Index< j, Dim12 > index2, const Index< j, Dim12 > index3) const
Tensor2_symmetric_Expr< const Christof_number_0< const Christof< T *, Tensor_Dim0, Tensor_Dim12 >, T, N >, T, Dim12, i, j > operator()(const Number< N > n1, const Index< i, Dim12 > index1, const Index< j, Dim12 > index2) const
const Christof< PackPtr< T *, I >, Tensor_Dim0, Tensor_Dim12 > & operator++() const
Tensor2_Expr< const Christof_numeral_1< const Christof< T *, Tensor_Dim0, Tensor_Dim12 >, T >, T, Dim0, Dim2, i, j > operator()(const Index< i, Dim0 > index1, const int N, const Index< j, Dim2 > index2) const
Tensor1_Expr< const Tensor3_contracted_02< const Christof< T *, Tensor_Dim0, Tensor_Dim12 >, T, Dim02 >, T, Dim1, i > operator()(const Index< j, Dim02 > index1, const Index< i, Dim1 > index2, const Index< j, Dim02 > index3) const
Dg_Expr< Christof< T *, Tensor_Dim0, Tensor_Dim12 >, T, Dim12, Dim0, i, j, k > operator()(const Index< k, Dim0 > index1, const Index< i, Dim12 > index2, const Index< j, Dim12 > index3)
const int N
Definition: speed_test.cpp:3
Dg_Expr< const Christof< T *, Tensor_Dim0, Tensor_Dim12 >, T, Dim12, Dim0, i, j, k > operator()(const Index< k, Dim0 > index1, const Index< i, Dim12 > index2, const Index< j, Dim12 > index3) const
Tensor1_Expr< const Tensor3_contracted_12< const Christof< T *, Tensor_Dim0, Tensor_Dim12 >, T, Dim12 >, T, Dim0, i > operator()(const Index< i, Dim0 > index1, const Index< j, Dim12 > index2, const Index< j, Dim12 > index3)