v0.13.1
Loading...
Searching...
No Matches
Tensor4_value.hpp
Go to the documentation of this file.
1/* A general version, not for pointers. */
2
3#pragma once
4
6
7#include <ostream>
8
9#ifdef FTENSOR_DEBUG
10#include <sstream>
11#include <stdexcept>
12#endif
13
14namespace FTensor
15{
16 template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2,
17 int Tensor_Dim3>
18 class Tensor4
19 {
20 T data[Tensor_Dim0][Tensor_Dim1][Tensor_Dim2][Tensor_Dim3];
21
22 public:
23 /* Initialization operators */
24 template <class... U> Tensor4(U... d) : data{d...}
25 {
26 static_assert(sizeof...(d) == sizeof(data) / sizeof(T),
27 "Incorrect number of Arguments. Constructor should "
28 "initialize the entire Tensor");
29 }
30
32 /* There are two operator(int,int,int,int)'s, one for non-consts
33 that lets you change the value, and one for consts that
34 doesn't. */
35 T &operator()(const int N1, const int N2, const int N3, const int N4)
36 {
37#ifdef FTENSOR_DEBUG
38 if(N1 >= Tensor_Dim0 || N1 < 0 || N2 >= Tensor_Dim1 || N2 < 0
39 || N3 >= Tensor_Dim2 || N3 < 0 || N4 >= Tensor_Dim3 || N4 < 0)
40 {
41 std::stringstream s;
42 s << "Bad index in Tensor4<T," << Tensor_Dim0 << "," << Tensor_Dim1
43 << "," << Tensor_Dim2 << "," << Tensor_Dim3 << ">.operator(" << N1
44 << "," << N2 << "," << N3 << "," << N4 << ") const" << std::endl;
45 throw std::out_of_range(s.str());
46 }
47#endif
48 return data[N1][N2][N3][N4];
49 }
50
51 T operator()(const int N1, const int N2, const int N3, const int N4) const
52 {
53#ifdef FTENSOR_DEBUG
54 if(N1 >= Tensor_Dim0 || N1 < 0 || N2 >= Tensor_Dim1 || N2 < 0
55 || N3 >= Tensor_Dim2 || N3 < 0 || N4 >= Tensor_Dim3 || N4 < 0)
56 {
57 std::stringstream s;
58 s << "Bad index in Tensor4<T," << Tensor_Dim0 << "," << Tensor_Dim1
59 << "," << Tensor_Dim2 << "," << Tensor_Dim3 << ">.operator(" << N1
60 << "," << N2 << "," << N3 << "," << N4 << ")" << std::endl;
61 throw std::out_of_range(s.str());
62 }
63#endif
64 return data[N1][N2][N3][N4];
65 }
66
67 /* These operator()'s are the first part in constructing template
68 expressions. They can be used to slice off lower dimensional
69 parts. */
70
71 template <char i, char j, char k, char l, int Dim0, int Dim1, int Dim2,
72 int Dim3>
73 typename std::enable_if<
74 (Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2
75 && Tensor_Dim3 >= Dim3),
77 T, Dim0, Dim1, Dim2, Dim3, i, j, k, l>>::type
78 operator()(const Index<i, Dim0>, const Index<j, Dim1>,
79 const Index<k, Dim2>, const Index<l, Dim3>)
80 {
81 return Tensor4_Expr<
83 Dim0, Dim1, Dim2, Dim3, i, j, k, l>(*this);
84 };
85
86 template <char i, char j, char k, char l, int Dim0, int Dim1, int Dim2,
87 int Dim3>
88 typename std::enable_if<
89 (Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2
90 && Tensor_Dim3 >= Dim3),
93 T, Dim0, Dim1, Dim2, Dim3, i, j, k, l>>::type
94 operator()(const Index<i, Dim0>, const Index<j, Dim1>,
95 const Index<k, Dim2>, const Index<l, Dim3>) const
96 {
97 return Tensor4_Expr<
99 T, Dim0, Dim1, Dim2, Dim3, i, j, k, l>(*this);
100 };
101
102 /* These operators are for internal contractions, resulting in a Tensor2.
103 * For example something like A(k,i,k,j) */
104
105 // (i,j,k,k)
106 template <char i, char j, char k, int Dim0, int Dim1, int Dim23>
108 const Index<k, Dim23>, const Index<k, Dim23>) const
109 {
110 static_assert(Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim1
111 && Tensor_Dim2 >= Dim23 && Tensor_Dim3 >= Dim23,
112 "Incompatible indices");
113
114 using TensorExpr = Tensor4_contracted_23<
116 Dim23>;
117 return Tensor2_Expr<TensorExpr, T, Dim0, Dim1, i, j>(TensorExpr(*this));
118 };
119
120 // (i,j,k,j)
121 template <char i, char j, char k, int Dim0, int Dim13, int Dim2>
123 const Index<k, Dim2>, const Index<j, Dim13>) const
124 {
125 static_assert(Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim13
126 && Tensor_Dim2 >= Dim2 && Tensor_Dim3 >= Dim13,
127 "Incompatible indices");
128
129 using TensorExpr = Tensor4_contracted_13<
131 Dim13>;
132 return Tensor2_Expr<TensorExpr, T, Dim0, Dim2, i, k>(TensorExpr(*this));
133 };
134
135 // (i,j,j,l)
136 template <char i, char j, char l, int Dim0, int Dim12, int Dim3>
138 const Index<j, Dim12>, const Index<l, Dim3>) const
139 {
140 static_assert(Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim12
141 && Tensor_Dim2 >= Dim12 && Tensor_Dim3 >= Dim3,
142 "Incompatible indices");
143
144 using TensorExpr = Tensor4_contracted_12<
146 Dim12>;
147 return Tensor2_Expr<TensorExpr, T, Dim0, Dim3, i, l>(TensorExpr(*this));
148 };
149
150 // (i,j,k,i)
151 template <char i, char j, char k, int Dim03, int Dim1, int Dim2>
153 const Index<k, Dim2>, const Index<i, Dim03>) const
154 {
155 static_assert(Tensor_Dim0 >= Dim03 && Tensor_Dim1 >= Dim1
156 && Tensor_Dim2 >= Dim2 && Tensor_Dim3 >= Dim03,
157 "Incompatible indices");
158
159 using TensorExpr = Tensor4_contracted_03<
161 Dim03>;
162 return Tensor2_Expr<TensorExpr, T, Dim1, Dim2, j, k>(TensorExpr(*this));
163 };
164
165 // (i,j,i,l)
166 template <char i, char j, char l, int Dim02, int Dim1, int Dim3>
168 const Index<i, Dim02>, const Index<l, Dim3>) const
169 {
170 static_assert(Tensor_Dim0 >= Dim02 && Tensor_Dim1 >= Dim1
171 && Tensor_Dim2 >= Dim02 && Tensor_Dim3 >= Dim3,
172 "Incompatible indices");
173
174 using TensorExpr = Tensor4_contracted_02<
176 Dim02>;
177 return Tensor2_Expr<TensorExpr, T, Dim1, Dim3, j, l>(TensorExpr(*this));
178 };
179
180 // (i,i,k,l)
181 template <char i, char k, char l, int Dim01, int Dim2, int Dim3>
183 const Index<k, Dim2>, const Index<l, Dim3>) const
184 {
185 static_assert(Tensor_Dim0 >= Dim01 && Tensor_Dim1 >= Dim01
186 && Tensor_Dim2 >= Dim2 && Tensor_Dim3 >= Dim3,
187 "Incompatible indices");
188
189 using TensorExpr = Tensor4_contracted_01<
191 Dim01>;
192 return Tensor2_Expr<TensorExpr, T, Dim2, Dim3, k, l>(TensorExpr(*this));
193 };
194
195 /* This is for expressions where a number is used for one slot, and
196 an index for another, yielding a Tensor3_Expr. The non-const
197 versions don't actually create a Tensor3_number_rhs_[0123] object.
198 They create a Tensor3_Expr directly, which provides the
199 appropriate indexing operators. The const versions do create a
200 Tensor3_number_[0123]. */
201
202 /* Third slot. */
203
204 template <char i, char j, char k, int N, int Dim0, int Dim1, int Dim3>
205 typename std::enable_if<
206 (Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim1 && Tensor_Dim2 > N &&
207 Tensor_Dim3 >= Dim3),
208 Tensor3_Expr<Tensor4_number_rhs_2<Tensor4<T, Tensor_Dim0, Tensor_Dim1,
209 Tensor_Dim2, Tensor_Dim3>,
210 T, N>,
211 T, Dim0, Dim1, Dim3, i, j, k>>::type
212 operator()(const Index<i, Dim0>, const Index<j, Dim1>, const Number<N>,
213 const Index<k, Dim3>) {
214 using TensorExpr = Tensor4_number_rhs_2<
217 }
218
219 template <char i, char j, char k, int N, int Dim0, int Dim1, int Dim3>
220 typename std::enable_if<
221 (Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim1 && Tensor_Dim2 > N &&
222 Tensor_Dim3 >= Dim3),
223 Tensor3_Expr<Tensor4_number_2<const Tensor4<T, Tensor_Dim0, Tensor_Dim1,
224 Tensor_Dim2, Tensor_Dim3>,
225 T, N>,
226 T, Dim0, Dim1, Dim3, i, j, k>>::type
227 operator()(const Index<i, Dim0>, const Index<j, Dim1>, const Number<N>,
228 const Index<k, Dim3>) const {
229 using TensorExpr = Tensor4_number_2<
231 T, N>;
233 TensorExpr(*this));
234 }
235
236 /* Forth slot. */
237
238 template <char i, char j, char k, int N, int Dim0, int Dim1, int Dim2>
239 typename std::enable_if<
240 (Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2 &&
241 Tensor_Dim3 > N),
242 Tensor3_Expr<Tensor4_number_rhs_3<Tensor4<T, Tensor_Dim0, Tensor_Dim1,
243 Tensor_Dim2, Tensor_Dim3>,
244 T, N>,
245 T, Dim0, Dim1, Dim2, i, j, k>>::type
246 operator()(const Index<i, Dim0>, const Index<j, Dim1>, const Index<k, Dim2>,
247 const Number<N>) {
248 using TensorExpr = Tensor4_number_rhs_3<
251 }
252
253 template <char i, char j, char k, int N, int Dim0, int Dim1, int Dim2>
254 typename std::enable_if<
255 (Tensor_Dim0 >= Dim0 && Tensor_Dim1 >= Dim1 && Tensor_Dim2 >= Dim2 &&
256 Tensor_Dim3 > N),
257 Tensor3_Expr<Tensor4_number_3<const Tensor4<T, Tensor_Dim0, Tensor_Dim1,
258 Tensor_Dim2, Tensor_Dim3>,
259 T, N>,
260 T, Dim0, Dim1, Dim2, i, j, k>>::type
261 operator()(const Index<i, Dim0>, const Index<j, Dim1>, const Index<k, Dim2>,
262 const Number<N>) const {
263 using TensorExpr = Tensor4_number_3<
265 T, N>;
267 TensorExpr(*this));
268 }
269 };
270}
271
272/// JSON compatible output
273
274namespace FTensor
275{
276 template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2,
277 int Tensor_Dim3>
278 std::ostream &Tensor4_0001(
279 std::ostream &os,
281 const int &iterator0, const int &iterator1, const int &iterator2)
282 {
283 os << '[';
284 for(int i = 0; i < Tensor_Dim3 - 1; ++i)
285 {
286 os << t(iterator0, iterator1, iterator2, i);
287 os << ',';
288 }
289 os << t(iterator0, iterator1, iterator2, Tensor_Dim3 - 1);
290 os << ']';
291
292 return os;
293 }
294
295 template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2,
296 int Tensor_Dim3>
297 std::ostream &Tensor4_0010(
298 std::ostream &os,
300 const int &iterator0, const int &iterator1)
301 {
302 os << '[';
303 for(int i = 0; i < Tensor_Dim2 - 1; ++i)
304 {
305 FTensor::Tensor4_0001(os, t, iterator0, iterator1, i);
306 os << ',';
307 }
308 FTensor::Tensor4_0001(os, t, iterator0, iterator1, Tensor_Dim2 - 1);
309 os << ']';
310
311 return os;
312 }
313
314 template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2,
315 int Tensor_Dim3>
316 std::ostream &Tensor4_0100(
317 std::ostream &os,
319 const int &iterator0)
320 {
321 os << '[';
322 for(int i = 0; i < Tensor_Dim1 - 1; ++i)
323 {
324 FTensor::Tensor4_0010(os, t, iterator0, i);
325 os << ',';
326 }
327 FTensor::Tensor4_0010(os, t, iterator0, Tensor_Dim1 - 1);
328 os << ']';
329
330 return os;
331 }
332}
333
334template <class T, int Tensor_Dim0, int Tensor_Dim1, int Tensor_Dim2,
335 int Tensor_Dim3>
336std::ostream &operator<<(std::ostream &os,
337 const FTensor::Tensor4<T, Tensor_Dim0, Tensor_Dim1,
338 Tensor_Dim2, Tensor_Dim3> &t)
339{
340 os << '[';
341 for(int i = 0; i < Tensor_Dim0 - 1; ++i)
342 {
344 os << ',';
345 }
346 FTensor::Tensor4_0100(os, t, Tensor_Dim0 - 1);
347 os << ']';
348
349 return os;
350}
static Number< 2 > N2
static Number< 1 > N1
T data[Tensor_Dim0][Tensor_Dim1][Tensor_Dim2][Tensor_Dim3]
T & operator()(const int N1, const int N2, const int N3, const int N4)
T operator()(const int N1, const int N2, const int N3, const int N4) const
auto operator()(const Index< i, Dim0 >, const Index< j, Dim13 >, const Index< k, Dim2 >, const Index< j, Dim13 >) const
auto operator()(const Index< i, Dim0 >, const Index< j, Dim12 >, const Index< j, Dim12 >, const Index< l, Dim3 >) const
auto operator()(const Index< i, Dim02 >, const Index< j, Dim1 >, const Index< i, Dim02 >, const Index< l, Dim3 >) const
auto operator()(const Index< i, Dim01 >, const Index< i, Dim01 >, const Index< k, Dim2 >, const Index< l, Dim3 >) const
auto operator()(const Index< i, Dim03 >, const Index< j, Dim1 >, const Index< k, Dim2 >, const Index< i, Dim03 >) const
auto operator()(const Index< i, Dim0 >, const Index< j, Dim1 >, const Index< k, Dim23 >, const Index< k, Dim23 >) const
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
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
std::ostream & Tensor4_0010(std::ostream &os, const Tensor4< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3 > &t, const int &iterator0, const int &iterator1)
std::ostream & Tensor4_0001(std::ostream &os, const Tensor4< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3 > &t, const int &iterator0, const int &iterator1, const int &iterator2)
std::ostream & Tensor4_0100(std::ostream &os, const Tensor4< T, Tensor_Dim0, Tensor_Dim1, Tensor_Dim2, Tensor_Dim3 > &t, const int &iterator0)
constexpr double t
plate stiffness
Definition: plate.cpp:59
const int N
Definition: speed_test.cpp:3