v0.14.0
Loading...
Searching...
No Matches
Dg_value.hpp
Go to the documentation of this file.
1/* A general version, not for pointers. */
2
3#pragma once
4
5#include "../Tensor3.hpp"
6
7namespace FTensor
8{
9 template <class T, int Tensor_Dim01, int Tensor_Dim2> class Dg
10 {
11 T data[(Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2][Tensor_Dim2];
12
13 public:
14 template <class... U> Dg(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
21 Dg() {}
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_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
30 || N3 >= Tensor_Dim2 || N3 < 0)
31 {
32 std::stringstream s;
33 s << "Bad index in Dg<T," << Tensor_Dim01 << "," << Tensor_Dim2
34 << ">.operator(" << N1 << "," << N2 << "," << N3 << ")"
35 << std::endl;
36 throw std::out_of_range(s.str());
37 }
38#endif
39 return N1 > N2 ? data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2][N3]
40 : data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2][N3];
41 }
42
43 T operator()(const int N1, const int N2, const int N3) const
44 {
45#ifdef FTENSOR_DEBUG
46 if(N1 >= Tensor_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
47 || N3 >= Tensor_Dim2 || N3 < 0)
48 {
49 std::stringstream s;
50 s << "Bad index in Dg<T," << Tensor_Dim01 << "," << Tensor_Dim2
51 << ">.operator(" << N1 << "," << N2 << "," << N3 << ") const"
52 << std::endl;
53 throw std::out_of_range(s.str());
54 }
55#endif
56 return N1 > N2 ? data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2][N3]
57 : data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2][N3];
58 }
59
60 /* These operator()'s are the first part in constructing template
61 expressions. */
62
63 template <char i, char j, char k, int Dim01, int Dim2>
64 typename std::enable_if<
65 (Tensor_Dim01 >= Dim01 && Tensor_Dim2 >= Dim2),
66 Dg_Expr<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, Dim01, Dim2, i, j, k>>::type
67 operator()(const Index<i, Dim01>, const Index<j, Dim01>,
68 const Index<k, Dim2>)
69 {
70 return Dg_Expr<Dg<T, Tensor_Dim01, Tensor_Dim2>, T, Dim01, Dim2, i, j, k>(
71 *this);
72 }
73 // TODO Add support for multiple dimensions ^^^^^
74 template <char i, char j, char k, int Dim01, int Dim2>
75 typename std::enable_if<(Tensor_Dim01 >= Dim01 && Tensor_Dim2 >= Dim2),
77 Dim01, Dim2, i, j, k>>::type
78 operator()(const Index<i, Dim01>, const Index<j, Dim01>,
79 const Index<k, Dim2>) const
80 {
82 j, k>(*this);
83 }
84
85 /* These operators are for internal contractions. We can not do
86 something like A(i,j,j) where i and j have different dimensions,
87 because it becomes ambiguous. */
88
89 /* A(i,j,j) */
90
91 template <char i, char j, int Dim>
92 typename std::enable_if<
93 (Tensor_Dim01 >= Dim && Tensor_Dim2 >= Dim),
96 Dim, i>>::type
97 operator()(const Index<i, Dim>, const Index<j, Dim>,
98 const Index<j, Dim>) const
99 {
100 using TensorExpr
102 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
103 }
104
105 /* A(j,i,j) */
106
107 template <char i, char j, int Dim>
108 typename std::enable_if<
109 (Tensor_Dim01 >= Dim && Tensor_Dim2 >= Dim),
112 Dim, i>>::type
113 operator()(const Index<j, Dim>, const Index<i, Dim>,
114 const Index<j, Dim>) const
115 {
116 using TensorExpr
118 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
119 }
120
121 /* A(j,j,i) */
122
123 template <char i, char j, int Dim01, int Dim2>
124 typename std::enable_if<
125 (Tensor_Dim01 >= Dim01 && Tensor_Dim2 >= Dim2),
128 Dim2, i>>::type
129 operator()(const Index<j, Dim01>, const Index<j, Dim01>,
130 const Index<i, Dim2>) const
131 {
132 using TensorExpr
134 return Tensor1_Expr<TensorExpr, T, Dim2, i>(TensorExpr(*this));
135 }
136
137 /* This is for expressions where a number is used for one slot, and
138 indices for the others, yielding a Tensor2_Expr or
139 Tensor2_symmetric_Expr. The non-const versions don't actually
140 create a Dg_number_rhs_* object, while the const versions
141 do create a Dg_number_*. */
142
143 /* First slot. */
144
145 template <char i, char j, int N, int Dim1, int Dim2>
146 typename std::enable_if<
147 (Tensor_Dim01 > N && Tensor_Dim01 >= Dim1 && Tensor_Dim2 >= Dim2),
149 Dim1, Dim2, i, j>>::type
150 operator()(const Number<N>, const Index<i, Dim1>, const Index<j, Dim2>)
151 {
152 using TensorExpr
155 }
156
157 template <char i, char j, int N, int Dim1, int Dim2>
158 typename std::enable_if<
159 (Tensor_Dim01 > N && Tensor_Dim01 >= Dim1 && Tensor_Dim2 >= Dim2),
161 T, Dim1, Dim2, i, j>>::type
162 operator()(const Number<N>, const Index<i, Dim1>,
163 const Index<j, Dim2>) const
164 {
165 using TensorExpr
167 return Tensor2_Expr<TensorExpr, T, Dim1, Dim2, i, j>(TensorExpr(*this));
168 }
169
170 /* Second slot. */
171
172 template <char i, char j, int N, int Dim0, int Dim2>
173 typename std::enable_if<
174 (Tensor_Dim01 >= Dim0 && Tensor_Dim01 > N && Tensor_Dim2 >= Dim2),
176 Dim0, Dim2, i, j>>::type
177 operator()(const Index<i, Dim0>, const Number<N>, const Index<j, Dim2>)
178 {
179 using TensorExpr
182 }
183
184 template <char i, char j, int N, int Dim0, int Dim2>
185 typename std::enable_if<
186 (Tensor_Dim01 >= Dim0 && Tensor_Dim01 > N && Tensor_Dim2 >= Dim2),
188 T, Dim0, Dim2, i, j>>::type
189 operator()(const Index<i, Dim0>, const Number<N>,
190 const Index<j, Dim2>) const
191 {
192 using TensorExpr
194 return Tensor2_Expr<TensorExpr, T, Dim0, Dim2, i, j>(TensorExpr(*this));
195 }
196
197 /* Third slot. */
198
199 template <char i, char j, int N, int Dim>
200 typename std::enable_if<
201 (Tensor_Dim01 >= Dim && Tensor_Dim2 > N),
204 j>>::type
205 operator()(const Index<i, Dim>, const Index<j, Dim>, const Number<N>)
206 {
207 using TensorExpr
210 }
211
212 template <char i, char j, int N, int Dim>
213 typename std::enable_if<
214 (Tensor_Dim01 >= Dim && Tensor_Dim2 > N),
217 j>>::type
218 operator()(const Index<i, Dim>, const Index<j, Dim>, const Number<N>) const
219 {
220 using TensorExpr
223 TensorExpr(*this));
224 }
225 // TODO add multiple dimensions up here ^^^^^
226 /* This is for expressions where a number is used for two slots, and
227 an Index for the other, yielding a Tensor1_Expr. The non-const
228 versions don't actually create a Dg_number_rhs_* object,
229 while the const versions do create a Dg_number_*. */
230
231 /* Index in first slot. */
232
233 template <char i, int N1, int N2, int Dim>
234 typename std::enable_if<
235 (Tensor_Dim01 >= Dim && Tensor_Dim01 > N1 && Tensor_Dim2 > N2),
237 T, Dim, i>>::type
238 operator()(const Index<i, Dim> index, const Number<N1>, const Number<N2>)
239 {
240 using TensorExpr
243 }
244
245 template <char i, int N1, int N2, int Dim>
246 typename std::enable_if<
247 (Tensor_Dim01 >= Dim && Tensor_Dim01 > N1 && Tensor_Dim2 > N2),
250 Dim, i>>::type
251 operator()(const Index<i, Dim> index, const Number<N1>,
252 const Number<N2>) const
253 {
254 using TensorExpr
256 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
257 }
258
259 /* Index in second slot. I use the same structures as for the Index
260 in the first slot since the tensor is symmetric on the first two
261 indices. */
262
263 template <char i, int N0, int N2, int Dim>
264 typename std::enable_if<
265 (Tensor_Dim01 > N0 && Tensor_Dim01 >= Dim && Tensor_Dim2 > N2),
267 T, Dim, i>>::type
268 operator()(const Number<N0>, const Index<i, Dim>, const Number<N2>)
269 {
270 using TensorExpr
273 }
274
275 template <char i, int N0, int N2, int Dim>
276 typename std::enable_if<
277 (Tensor_Dim01 > N0 && Tensor_Dim01 >= Dim && Tensor_Dim2 > N2),
280 Dim, i>>::type
281 operator()(const Number<N0>, const Index<i, Dim> index,
282 const Number<N2>) const
283 {
284 using TensorExpr
286 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
287 }
288
289 /* Index in third slot. */
290
291 template <char i, int N0, int N1, int Dim>
292 typename std::enable_if<
293 (Tensor_Dim01 > N0 && Tensor_Dim01 > N1 && Tensor_Dim2 >= Dim),
295 T, Dim, i>>::type
296 operator()(const Number<N0>, const Number<N1>, const Index<i, Dim> index)
297 {
298 using TensorExpr
301 }
302
303 template <char i, int N0, int N1, int Dim>
304 typename std::enable_if<
305 (Tensor_Dim01 > N0 && Tensor_Dim01 > N1 && Tensor_Dim2 >= Dim),
308 Dim, i>>::type
309 operator()(const Number<N0>, const Number<N1>, const Index<i, Dim>) const
310 {
311 using TensorExpr
313 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
314 }
315
316 /* Specializations for using actual numbers instead of Number<>.
317 This is for expressions where an actual number is used for one
318 slot, and indices for the others, yielding a Tensor2_Expr or
319 Tensor2_symmetric_Expr. I only define the const versions. */
320
321 /* First slot. */
322
323 template <char i, char j, int Dim1, int Dim2>
324 typename std::enable_if<
325 (Tensor_Dim01 >= Dim1 && Tensor_Dim2 >= Dim2),
327 Dim1, Dim2, i, j>>::type
328 operator()(const int N, const Index<i, Dim1>, const Index<j, Dim2>) const
329 {
330 using TensorExpr
333 TensorExpr(*this, N));
334 }
335
336 /* Second slot. */
337
338 template <char i, char j, int Dim0, int Dim2>
339 typename std::enable_if<
340 (Tensor_Dim01 >= Dim0 && Tensor_Dim2 >= Dim2),
342 Dim0, Dim2, i, j>>::type
343 operator()(const Index<i, Dim0>, const int N, const Index<j, Dim2>) const
344 {
345 using TensorExpr
348 TensorExpr(*this, N));
349 }
350
351 /* Third slot. */
352
353 template <char i, char j, int Dim>
354 typename std::enable_if<
355 (Tensor_Dim01 >= Dim),
358 j>>::type
359 operator()(const Index<i, Dim>, const Index<j, Dim>, const int N) const
360 {
361 using TensorExpr
364 TensorExpr(*this, N));
365 }
366 // TODO Allow two different dimensions in 1 and 2
367 /* This is for expressions where a numeral is used for two slots, and
368 an Index for the other, yielding a Tensor1_Expr. */
369
370 /* Index in first slot. */
371
372 template <char i, int Dim>
373 typename std::enable_if<
374 (Tensor_Dim01 >= Dim),
376 Dim, i>>::type
377 operator()(const Index<i, Dim> index, const int N1, const int N2) const
378 {
379 using TensorExpr
381 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
382 }
383
384 /* Index in second slot. I use the same structures as for the Index
385 in the first slot since the tensor is symmetric on the first two
386 indices. */
387
388 template <char i, int Dim>
389 typename std::enable_if<
390 (Tensor_Dim01 >= Dim),
392 Dim, i>>::type
393 operator()(const int N1, const Index<i, Dim> index, const int N2) const
394 {
395 using TensorExpr
397 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
398 }
399
400 /* Index in third slot. */
401
402 template <char i, int Dim>
403 typename std::enable_if<
404 (Tensor_Dim2 >= Dim),
406 Dim, i>>::type
407 operator()(const int N1, const int N2, const Index<i, Dim>) const
408 {
409 using TensorExpr
411 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N1, N2));
412 }
413 };
414}
static Number< 2 > N2
static Number< 1 > N1
static Number< 0 > N0
T & operator()(const int N1, const int N2, const int N3)
Definition Dg_value.hpp:26
T data[(Tensor_Dim01 *(Tensor_Dim01+1))/2][Tensor_Dim2]
Definition Dg_value.hpp:11
T operator()(const int N1, const int N2, const int N3) const
Definition Dg_value.hpp:43
Dg(U... d)
Definition Dg_value.hpp:14
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