v0.14.0
Loading...
Searching...
No Matches
Christof_value.hpp
Go to the documentation of this file.
1/* A general version, not for pointers */
2
3#pragma once
4
5#include "../Dg.hpp"
6
7namespace FTensor
8{
9 template <class T, int Tensor_Dim0, int Tensor_Dim12> class Christof
10 {
11 T data[Tensor_Dim0][(Tensor_Dim12 * (Tensor_Dim12 + 1)) / 2];
12
13 public:
14 template <class... U> Christof(U... d) : data{d...}
15 {
16 static_assert(sizeof...(d) == sizeof(data) / sizeof(T),
17 "Incorrect number of Arguments. Constructor should "
18 "initialize the entire Tensor");
19 }
20
22
23 /* There are two operator(int,int,int)'s, one for non-consts that lets you
24 change the value, and one for consts that doesn't. */
25
26 T &operator()(const int N1, const int N2, const int N3)
27 {
28#ifdef FTENSOR_DEBUG
29 if(N1 >= Tensor_Dim0 || N1 < 0 || N2 >= Tensor_Dim12 || N2 < 0
30 || N3 >= Tensor_Dim12 || N3 < 0)
31 {
32 std::stringstream s;
33 s << "Bad index in Christof<T," << Tensor_Dim0 << "," << Tensor_Dim12
34 << ">.operator(" << N1 << "," << N2 << "," << N3 << ")"
35 << std::endl;
36 throw std::out_of_range(s.str());
37 }
38#endif
39 return N2 > N3 ? data[N1][N2 + (N3 * (2 * Tensor_Dim12 - N3 - 1)) / 2]
40 : data[N1][N3 + (N2 * (2 * Tensor_Dim12 - N2 - 1)) / 2];
41 }
42
43 T operator()(const int N1, const int N2, const int N3) const
44 {
45#ifdef FTENSOR_DEBUG
46 if(N1 >= Tensor_Dim0 || N1 < 0 || N2 >= Tensor_Dim12 || N2 < 0
47 || N3 >= Tensor_Dim12 || N3 < 0)
48 {
49 std::stringstream s;
50 s << "Bad index in Christof<T," << Tensor_Dim0 << "," << Tensor_Dim12
51 << ">.operator(" << N1 << "," << N2 << "," << N3 << ") const"
52 << std::endl;
53 throw std::out_of_range(s.str());
54 }
55#endif
56 return N2 > N3 ? data[N1][N2 + (N3 * (2 * Tensor_Dim12 - N3 - 1)) / 2]
57 : data[N1][N3 + (N2 * (2 * Tensor_Dim12 - N2 - 1)) / 2];
58 }
59
60 /* These operator()'s are the first part in constructing template
61 expressions. I mix up the indices here so that it behaves like a
62 Dg. That way I don't have to have a separate wrapper
63 class Christof_Expr, which simplifies things. */
64
65 template <char i, char j, char k, int Dim0, int Dim12>
66 typename std::enable_if<(Tensor_Dim0 >= Dim0 && Tensor_Dim12 >= Dim12),
68 Dim12, Dim0, i, j, k>>::type
69 operator()(const Index<k, Dim0>, const Index<i, Dim12>,
70 const Index<j, Dim12>)
71 {
73 j, k>(*this);
74 }
75
76 template <char i, char j, char k, int Dim0, int Dim12>
77 typename std::enable_if<(Tensor_Dim0 >= Dim0 && Tensor_Dim12 >= Dim12),
79 T, Dim12, Dim0, i, j, k>>::type
80 operator()(const Index<k, Dim0>, const Index<i, Dim12>,
81 const Index<j, Dim12>) const
82 {
84 Dim0, i, j, k>(*this);
85 }
86
87 /* This is for expressions where a number is used for two slots, and
88 an Index for the other, yielding a Tensor1_Expr. The non-const
89 versions don't actually create a Dg_number_rhs_* object,
90 while the const versions do create a Dg_number_*. */
91
92 /* Index in first slot. */
93
94 template <char i, int N1, int N2, int Dim>
95 typename std::enable_if<
96 (Tensor_Dim0 >= Dim && Tensor_Dim12 > N1 && Tensor_Dim12 > N2),
99 Dim, i>>::type
100 operator()(const Index<i, Dim>, const Number<N1>, const Number<N2>)
101 {
102 using TensorExpr
105 }
106
107 template <char i, int N1, int N2, int Dim>
108 typename std::enable_if<
109 (Tensor_Dim0 >= Dim && Tensor_Dim12 > N1 && Tensor_Dim12 > N2),
112 Dim, i>>::type
113 operator()(const Index<i, Dim>, const Number<N1>, const Number<N2>) const
114 {
115 using TensorExpr
117 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
118 }
119 // TODO index on second and third position
120
121 /* These operators are for internal contractions. Some of them are
122 less general, because, for example, A(j,i,j) with i and j having
123 different dimensions is ambiguous. */
124
125 /* const versions */
126
127 template <char i, char j, int Dim0, int Dim12>
128 typename std::enable_if<
129 (Tensor_Dim0 >= Dim0 && Tensor_Dim12 >= Dim12),
132 T, Dim0, i>>::type
133 operator()(const Index<i, Dim0>, const Index<j, Dim12>,
134 const Index<j, Dim12>) const
135 {
136 using TensorExpr
138 Dim12>;
139 return Tensor1_Expr<TensorExpr, T, Dim0, i>(TensorExpr(*this));
140 }
141
142 template <char i, char j, int Dim>
143 typename std::enable_if<
144 (Tensor_Dim0 >= Dim && Tensor_Dim12 >= Dim),
147 T, Dim, i>>::type
148 operator()(const Index<j, Dim>, const Index<i, Dim>,
149 const Index<j, Dim>) const
150 {
151 using TensorExpr
153 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
154 }
155 // TODO allow different dimensions here ^^^^ with the one below
156 // template<char i, char j, int Dim02, int Dim1>
157 // Tensor1_Expr<const Tensor3_contracted_02
158 // <const Christof<T,Tensor_Dim0,Tensor_Dim12>,T,Dim02>,T,Dim1,i>
159 // operator()(const Index<j,Dim02> index1, const Index<i,Dim1> index2,
160 // const Index<j,Dim02> index3) const
161 // {
162 // typedef Tensor3_contracted_02
163 // <const Christof<T,Tensor_Dim0,Tensor_Dim12>,T,Dim02>
164 // TensorExpr;
165 // return Tensor1_Expr<TensorExpr,T,Dim1,i>(TensorExpr(*this));
166 // }
167
168 template <char i, char j, int Dim>
169 typename std::enable_if<
170 (Tensor_Dim0 >= Dim && Tensor_Dim12 >= Dim),
173 T, Dim, i>>::type
174 operator()(const Index<j, Dim>, const Index<j, Dim>,
175 const Index<i, Dim>) const
176 {
177 using TensorExpr
179 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
180 }
181 // TODO allow different dimensions here ^^^^ with the one below
182 // template<char i, char j, int Dim01, int Dim2>
183 // Tensor1_Expr<const Tensor3_contracted_01
184 // <const Christof<T,Tensor_Dim0,Tensor_Dim12>,T,Dim01>,T,Dim2,i>
185 // operator()(const Index<j,Dim01> index1, const Index<j,Dim01> index2,
186 // const Index<i,Dim2> index3) const
187 // {
188 // typedef Tensor3_contracted_01
189 // <const Christof<T,Tensor_Dim0,Tensor_Dim12>,T,Dim01>
190 // TensorExpr;
191 // return Tensor1_Expr<TensorExpr,T,Dim2,i>(TensorExpr(*this));
192 // }
193
194 /* non-const versions, needed so that overload resolution doesn't
195 pick the more general indexing operator. */
196
197 template <char i, char j, int Dim0, int Dim12>
198 typename std::enable_if<
199 (Tensor_Dim0 >= Dim0 && Tensor_Dim12 >= Dim12),
202 T, Dim0, i>>::type
203 operator()(const Index<i, Dim0>, const Index<j, Dim12>,
204 const Index<j, Dim12>)
205 {
206 using TensorExpr
208 Dim12>;
209 return Tensor1_Expr<TensorExpr, T, Dim0, i>(TensorExpr(*this));
210 }
211
212 template <char i, char j, int Dim>
213 typename std::enable_if<
214 (Tensor_Dim0 >= Dim && Tensor_Dim12 >= Dim),
217 T, Dim, i>>::type
218 operator()(const Index<j, Dim>, const Index<i, Dim>, const Index<j, Dim>)
219 {
220 using TensorExpr
222 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
223 }
224 // TODO allow different dimensions here ^^^^ with the one below
225 // template<char i, char j, int Dim02, int Dim1>
226 // Tensor1_Expr<const Tensor3_contracted_02
227 // <const Christof<T,Tensor_Dim0,Tensor_Dim12>,T,Dim02>,T,Dim1,i>
228 // operator()(const Index<j,Dim02> index1, const Index<i,Dim1> index2,
229 // const Index<j,Dim02> index3)
230 // {
231 // typedef Tensor3_contracted_02
232 // <const Christof<T,Tensor_Dim0,Tensor_Dim12>,T,Dim02>
233 // TensorExpr;
234 // return Tensor1_Expr<TensorExpr,T,Dim1,i>(TensorExpr(*this));
235 // }
236
237 template <char i, char j, int Dim>
238 typename std::enable_if<
239 (Tensor_Dim0 >= Dim && Tensor_Dim12 >= Dim),
242 T, Dim, i>>::type
243 operator()(const Index<j, Dim>, const Index<j, Dim>, const Index<i, Dim>)
244 {
245 using TensorExpr
247 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
248 }
249
250 // template<char i, char j, int Dim01, int Dim2>
251 // Tensor1_Expr<const Tensor3_contracted_01
252 // <const Christof<T,Tensor_Dim0,Tensor_Dim12>,T,Dim01>,T,Dim2,i>
253 // operator()(const Index<j,Dim01> index1, const Index<j,Dim01> index2,
254 // const Index<i,Dim2> index3)
255 // {
256 // typedef Tensor3_contracted_01
257 // <const Christof<T,Tensor_Dim0,Tensor_Dim12>,T,Dim01>
258 // TensorExpr;
259 // return Tensor1_Expr<TensorExpr,T,Dim2,i>(TensorExpr(*this));
260 // }
261
262 /* This is for expressions where a Number<> is used for one slot, and
263 an index for the others, yielding a Tensor2_symmetric_Expr. */
264
265 template <char i, char j, int N, int Dim12>
266 typename std::enable_if<
267 (Tensor_Dim0 > N && Tensor_Dim12 >= Dim12),
270 Dim12, i, j>>::type
271 operator()(const Number<N>, const Index<i, Dim12>,
272 const Index<j, Dim12>) const
273 {
274 using TensorExpr
277 TensorExpr(*this));
278 }
279 // TODO Allow multiple dimensions here
280 /* This is for expressions where an int is used for one slot, and
281 an index for the others, yielding a Tensor2_symmetric_Expr. */
282
283 template <char i, char j, int Dim12>
284 typename std::enable_if<
285 (Tensor_Dim12 >= Dim12),
288 Dim12, i, j>>::type
289 operator()(const int N, const Index<i, Dim12>, const Index<j, Dim12>) const
290 {
291 using TensorExpr
294 TensorExpr(*this, N));
295 }
296 // TODO allow multiple dimensions here
297 /* An int in one spot, Index for the others, yielding a Tensor2. I
298 can use the same structure for both, since Christof is
299 symmetric on the last two indices. */
300
301 template <char i, char j, int Dim0, int Dim2>
302 typename std::enable_if<
303 (Tensor_Dim0 >= Dim0 && Tensor_Dim12 >= Dim2),
305 T, Dim0, Dim2, i, j>>::type
306 operator()(const Index<i, Dim0>, const int N, const Index<j, Dim2>) const
307 {
308 using TensorExpr
311 TensorExpr(*this, N));
312 }
313
314 template <char i, char j, int Dim0, int Dim2>
315 typename std::enable_if<
316 (Tensor_Dim0 >= Dim0 && Tensor_Dim12 >= Dim2),
318 T, Dim0, Dim2, i, j>>::type
319 operator()(const Index<i, Dim0>, const Index<j, Dim2>, const int N) const
320 {
321 using TensorExpr
324 TensorExpr(*this, N));
325 }
326 };
327}
static Number< 2 > N2
static Number< 1 > N1
T & operator()(const int N1, const int N2, const int N3)
T operator()(const int N1, const int N2, const int N3) const
T data[Tensor_Dim0][(Tensor_Dim12 *(Tensor_Dim12+1))/2]
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
const int N
Definition: speed_test.cpp:3