v0.14.0
Loading...
Searching...
No Matches
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
7namespace 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
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 private:
257 template <int I>
258 Christof(const Christof<PackPtr<T *, I>, Tensor_Dim0, Tensor_Dim12> &) =
259 delete;
260 };
261
262 template <class T, int Tensor_Dim0, int Tensor_Dim12, int I>
263 class Christof<PackPtr<T *, I>, Tensor_Dim0, Tensor_Dim12>:
265 {
266
267 public:
268 template <class... U>
269 Christof(U *... d) : Christof<T *, Tensor_Dim0, Tensor_Dim12>(d...) {}
270
271 Christof() : Christof<T *, Tensor_Dim0, Tensor_Dim12>() {}
272
273 /* The ++ operator increments the pointer, not the number that the
274 pointer points to. This allows iterating over a grid. */
275
276 const Christof<PackPtr<T *, I>, Tensor_Dim0, Tensor_Dim12> &
277 operator++() const {
278 for (int i = 0; i < Tensor_Dim0; ++i)
279 for (int j = 0; j < (Tensor_Dim12 * (Tensor_Dim12 + 1)) / 2; ++j)
281 return *this;
282 }
283 };
284}
static Number< 2 > N2
static Number< 1 > N1
const Christof< PackPtr< T *, I >, Tensor_Dim0, Tensor_Dim12 > & operator++() 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
T & operator()(const int N1, const int N2, const int N3)
const Christof< T *, Tensor_Dim0, Tensor_Dim12 > & operator++() const
T * ptr(const int N1, const int N2, const int N3) const
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
Christof(const Christof< PackPtr< T *, I >, Tensor_Dim0, Tensor_Dim12 > &)=delete
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)
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 int N, const Index< j, Dim2 > index2) 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)
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
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)
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
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)
T operator()(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
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const double T
Tensors class implemented by Walter Landry.
Definition FTensor.hpp:51
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
constexpr IntegrationType I
const int N
Definition speed_test.cpp:3