v0.14.0
Loading...
Searching...
No Matches
Ddg_pointer.hpp
Go to the documentation of this file.
1/* A version for pointers. */
2
3#pragma once
4
5namespace FTensor
6{
7 template <class T, int Tensor_Dim01, int Tensor_Dim23>
8 class Ddg<T *, Tensor_Dim01, Tensor_Dim23>
9 {
10 const int inc;
11
12 protected:
13
14 mutable T *restrict data[(Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2]
15 [(Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2];
16
17 public:
18 Ddg(T *d0000, T *d0010, T *d0011, T *d1000, T *d1010, T *d1011, T *d1100,
19 T *d1110, T *d1111, const int i = 1)
20 : inc(i)
21 {
22 ptr(0, 0, 0, 0) = d0000;
23 ptr(0, 0, 1, 0) = d0010;
24 ptr(0, 0, 1, 1) = d0011;
25 ptr(1, 0, 0, 0) = d1000;
26 ptr(1, 0, 1, 0) = d1010;
27 ptr(1, 0, 1, 1) = d1011;
28 ptr(1, 1, 0, 0) = d1100;
29 ptr(1, 1, 1, 0) = d1110;
30 ptr(1, 1, 1, 1) = d1111;
31 }
32
33 Ddg(T *d0000, T *d0001, T *d0002, T *d0011, T *d0012, T *d0022, T *d0100,
34 T *d0101, T *d0102, T *d0111, T *d0112, T *d0122, T *d0200, T *d0201,
35 T *d0202, T *d0211, T *d0212, T *d0222, T *d1100, T *d1101, T *d1102,
36 T *d1111, T *d1112, T *d1122, T *d1200, T *d1201, T *d1202, T *d1211,
37 T *d1212, T *d1222, T *d2200, T *d2201, T *d2202, T *d2211, T *d2212,
38 T *d2222, const int i = 1)
39 : inc(i)
40 {
41 ptr(0, 0, 0, 0) = d0000;
42 ptr(0, 0, 0, 1) = d0001;
43 ptr(0, 0, 0, 2) = d0002;
44 ptr(0, 0, 1, 1) = d0011;
45 ptr(0, 0, 1, 2) = d0012;
46 ptr(0, 0, 2, 2) = d0022;
47 ptr(0, 1, 0, 0) = d0100;
48 ptr(0, 1, 0, 1) = d0101;
49 ptr(0, 1, 0, 2) = d0102;
50 ptr(0, 1, 1, 1) = d0111;
51 ptr(0, 1, 1, 2) = d0112;
52 ptr(0, 1, 2, 2) = d0122;
53 ptr(0, 2, 0, 0) = d0200;
54 ptr(0, 2, 0, 1) = d0201;
55 ptr(0, 2, 0, 2) = d0202;
56 ptr(0, 2, 1, 1) = d0211;
57 ptr(0, 2, 1, 2) = d0212;
58 ptr(0, 2, 2, 2) = d0222;
59 ptr(1, 1, 0, 0) = d1100;
60 ptr(1, 1, 0, 1) = d1101;
61 ptr(1, 1, 0, 2) = d1102;
62 ptr(1, 1, 1, 1) = d1111;
63 ptr(1, 1, 1, 2) = d1112;
64 ptr(1, 1, 2, 2) = d1122;
65 ptr(1, 2, 0, 0) = d1200;
66 ptr(1, 2, 0, 1) = d1201;
67 ptr(1, 2, 0, 2) = d1202;
68 ptr(1, 2, 1, 1) = d1211;
69 ptr(1, 2, 1, 2) = d1212;
70 ptr(1, 2, 2, 2) = d1222;
71 ptr(2, 2, 0, 0) = d2200;
72 ptr(2, 2, 0, 1) = d2201;
73 ptr(2, 2, 0, 2) = d2202;
74 ptr(2, 2, 1, 1) = d2211;
75 ptr(2, 2, 1, 2) = d2212;
76 ptr(2, 2, 2, 2) = d2222;
77 }
78
79 /* Initializations for varying numbers of elements. */
80 template <class... U> Ddg(U *... d) : data{d...}, inc(1) {}
81
82 /* There are two operator(int,int,int,int)'s, one for non-consts
83 that lets you change the value, and one for consts that
84 doesn't. */
85
86 T &operator()(const int N1, const int N2, const int N3, const int N4)
87 {
88#ifdef FTENSOR_DEBUG
89 if(N1 >= Tensor_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
90 || N3 >= Tensor_Dim23 || N3 < 0 || N4 >= Tensor_Dim23 || N4 < 0)
91 {
92 std::stringstream s;
93 s << "Bad index in Dg<T*," << Tensor_Dim01 << "," << Tensor_Dim23
94 << ">.operator(" << N1 << "," << N2 << "," << N3 << "," << N4
95 << ")" << std::endl;
96 throw std::out_of_range(s.str());
97 }
98#endif
99 return N1 > N2
100 ? (N3 > N4 ? *data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2]
101 [N3 + (N4 * (2 * Tensor_Dim23 - N4 - 1)) / 2]
102 : *data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2]
103 [N4 + (N3 * (2 * Tensor_Dim23 - N3 - 1)) / 2])
104 : (N3 > N4
105 ? *data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2]
106 [N3 + (N4 * (2 * Tensor_Dim23 - N4 - 1)) / 2]
107 : *data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2]
108 [N4 + (N3 * (2 * Tensor_Dim23 - N3 - 1)) / 2]);
109 }
110
111 T operator()(const int N1, const int N2, const int N3, const int N4) const
112 {
113#ifdef FTENSOR_DEBUG
114 if(N1 >= Tensor_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
115 || N3 >= Tensor_Dim23 || N3 < 0 || N4 >= Tensor_Dim23 || N4 < 0)
116 {
117 std::stringstream s;
118 s << "Bad index in Dg<T*," << Tensor_Dim01 << "," << Tensor_Dim23
119 << ">.operator(" << N1 << "," << N2 << "," << N3 << "," << N4
120 << ") const" << std::endl;
121 throw std::out_of_range(s.str());
122 }
123#endif
124 return N1 > N2
125 ? (N3 > N4 ? *data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2]
126 [N3 + (N4 * (2 * Tensor_Dim23 - N4 - 1)) / 2]
127 : *data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2]
128 [N4 + (N3 * (2 * Tensor_Dim23 - N3 - 1)) / 2])
129 : (N3 > N4
130 ? *data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2]
131 [N3 + (N4 * (2 * Tensor_Dim23 - N4 - 1)) / 2]
132 : *data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2]
133 [N4 + (N3 * (2 * Tensor_Dim23 - N3 - 1)) / 2]);
134 }
135
136 T *restrict &ptr(const int N1, const int N2, const int N3, const int N4) const
137 {
138#ifdef FTENSOR_DEBUG
139 if(N1 >= Tensor_Dim01 || N1 < 0 || N2 >= Tensor_Dim01 || N2 < 0
140 || N3 >= Tensor_Dim23 || N3 < 0 || N4 >= Tensor_Dim23 || N4 < 0)
141 {
142 std::stringstream s;
143 s << "Bad index in Dg<T," << Tensor_Dim01 << "," << Tensor_Dim23
144 << ">.ptr(" << N1 << "," << N2 << "," << N3 << "," << N4 << ")"
145 << std::endl;
146 throw std::out_of_range(s.str());
147 }
148#endif
149 return N1 > N2
150 ? (N3 > N4 ? data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2]
151 [N3 + (N4 * (2 * Tensor_Dim23 - N4 - 1)) / 2]
152 : data[N1 + (N2 * (2 * Tensor_Dim01 - N2 - 1)) / 2]
153 [N4 + (N3 * (2 * Tensor_Dim23 - N3 - 1)) / 2])
154 : (N3 > N4 ? data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2]
155 [N3 + (N4 * (2 * Tensor_Dim23 - N4 - 1)) / 2]
156 : data[N2 + (N1 * (2 * Tensor_Dim01 - N1 - 1)) / 2]
157 [N4 + (N3 * (2 * Tensor_Dim23 - N3 - 1)) / 2]);
158 }
159
160 /* These operator()'s are the first part in constructing template
161 expressions. They can be used to slice off lower dimensional
162 parts. They are not entirely safe, since you can accidently use a
163 higher dimension than what is really allowed (like Dim=5). */
164
165 template <char i, char j, char k, char l, int Dim01, int Dim23>
167 operator()(const Index<i, Dim01> index1, const Index<j, Dim01> index2,
168 const Index<k, Dim23> index3, const Index<l, Dim23> index4)
169 {
171 j, k, l>(*this);
172 }
173
174 template <char i, char j, char k, char l, int Dim01, int Dim23>
176 k, l>
177 operator()(const Index<i, Dim01> index1, const Index<j, Dim01> index2,
178 const Index<k, Dim23> index3,
179 const Index<l, Dim23> index4) const
180 {
182 Dim23, i, j, k, l>(*this);
183 }
184
185 /* This is for expressions where a number is used for two slots, and
186 an index for the other two, yielding a Tensor2_symmetric_Expr. */
187
188 template <char i, char j, int N0, int N1, int Dim>
191 Dim, i, j>
192 operator()(const Number<N0> n1, const Number<N1> n2,
193 const Index<i, Dim> index1, const Index<j, Dim> index2) const
194 {
195 using TensorExpr
198 TensorExpr(*this));
199 }
200
201 template <char i, char j, int N0, int N1, int Dim>
204 Dim, i, j>
205 operator()(const Number<N0> n1, const Number<N1> n2,
206 const Index<i, Dim> index1, const Index<j, Dim> index2)
207 {
208 using TensorExpr
211 }
212
213 /* This is for expressions where a number is used for one slot, and
214 an index for the other three, yielding a Dg_Expr. */
215
216 template <char i, char j, char k, int N0, int Dim1, int Dim23>
218 Dim23, Dim1, i, j, k>
219 operator()(const Number<N0> n1, const Index<k, Dim1> index3,
220 const Index<i, Dim23> index1,
221 const Index<j, Dim23> index2) const
222 {
223 using TensorExpr
225 return Dg_Expr<TensorExpr, T, Dim23, Dim1, i, j, k>(TensorExpr(*this));
226 }
227
228 template <char i, char j, char k, int N0, int Dim1, int Dim23>
230 Dim23, Dim1, i, j, k>
231 operator()(const Number<N0> n1, const Index<k, Dim1> index3,
232 const Index<i, Dim23> index1, const Index<j, Dim23> index2)
233 {
234 using TensorExpr
237 }
238
239 /* This is for expressions where an int (not a Number) is used for
240 two slots, and an index for the other two, yielding a
241 Tensor2_symmetric_Expr. */
242
243 template <char i, char j, int Dim>
246 j>
247 operator()(const int N0, const int N1, const Index<i, Dim> index1,
248 const Index<j, Dim> index2) const
249 {
250 using TensorExpr
253 TensorExpr(*this, N0, N1));
254 }
255
256 template <char i, char j, int Dim>
259 operator()(const Index<i, Dim> index1, const Index<j, Dim> index2,
260 const int N2, const int N3) const
261 {
262 using TensorExpr
265 TensorExpr(*this, N2, N3));
266 }
267
268 /* int in three slots, Index in the other yielding a Tensor1_Expr. */
269
270 template <char i, int Dim>
272 Dim, i>
273 operator()(const Index<i, Dim> index1, const int N1, const int N2,
274 const int N3)
275 {
276 using TensorExpr
279 TensorExpr(*this, N1, N2, N3));
280 }
281
282 template <char i, int Dim>
284 Dim, i>
285 operator()(const int N1, const Index<i, Dim> index1, const int N2,
286 const int N3)
287 {
288 using TensorExpr
291 TensorExpr(*this, N1, N2, N3));
292 }
293
294 /* This is for expressions where an int (not a Number) is used for
295 one slot, and an index for the other three, yielding a
296 Dg_Expr. */
297
298 template <char i, char j, char k, int Dim1, int Dim23>
300 Dim23, Dim1, i, j, k>
301 operator()(const int N0, const Index<k, Dim1> index3,
302 const Index<i, Dim23> index1,
303 const Index<j, Dim23> index2) const
304 {
305 using TensorExpr
308 TensorExpr(*this, N0));
309 }
310
311 /* The ++ operator increments the pointer, not the number that the
312 pointer points to. This allows iterating over a grid. */
313
314 const Ddg &operator++() const {
315 for(int i = 0; i < (Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2; ++i)
316 for(int j = 0; j < (Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2; ++j)
317 data[i][j] += inc;
318 return *this;
319 }
320
321 private:
322 template <int I>
323 Ddg(const Ddg<PackPtr<T *, I>, Tensor_Dim01, Tensor_Dim23> &) = delete;
324 };
325
326 template <class T, int Tensor_Dim01, int Tensor_Dim23, int I>
327 class Ddg<PackPtr<T *, I>, Tensor_Dim01, Tensor_Dim23>
329
330 public:
331 /* Initializations for varying numbers of elements. */
332 template <class... U>
333 Ddg(U *... d) : Ddg<T *, Tensor_Dim01, Tensor_Dim23>(d...) {}
334
335 /* The ++ operator increments the pointer, not the number that the
336 pointer points to. This allows iterating over a grid. */
337
338 const Ddg &operator++() const {
339 for (int i = 0; i < (Tensor_Dim01 * (Tensor_Dim01 + 1)) / 2; ++i)
340 for (int j = 0; j < (Tensor_Dim23 * (Tensor_Dim23 + 1)) / 2; ++j)
342 return *this;
343 }
344 };
345}
static Number< 2 > N2
static Number< 1 > N1
static Number< 0 > N0
Ddg_Expr< Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T, Dim01, Dim23, i, j, k, l > operator()(const Index< i, Dim01 > index1, const Index< j, Dim01 > index2, const Index< k, Dim23 > index3, const Index< l, Dim23 > index4)
T *restrict & ptr(const int N1, const int N2, const int N3, const int N4) const
Tensor1_Expr< Ddg_numeral_123< Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T >, T, Dim, i > operator()(const Index< i, Dim > index1, const int N1, const int N2, const int N3)
Tensor2_symmetric_Expr< Ddg_number_01< const Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T, N0, N1 >, T, Dim, i, j > operator()(const Number< N0 > n1, const Number< N1 > n2, const Index< i, Dim > index1, const Index< j, Dim > index2) const
Tensor1_Expr< Ddg_numeral_123< Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T >, T, Dim, i > operator()(const int N1, const Index< i, Dim > index1, const int N2, const int N3)
Ddg(T *d0000, T *d0001, T *d0002, T *d0011, T *d0012, T *d0022, T *d0100, T *d0101, T *d0102, T *d0111, T *d0112, T *d0122, T *d0200, T *d0201, T *d0202, T *d0211, T *d0212, T *d0222, T *d1100, T *d1101, T *d1102, T *d1111, T *d1112, T *d1122, T *d1200, T *d1201, T *d1202, T *d1211, T *d1212, T *d1222, T *d2200, T *d2201, T *d2202, T *d2211, T *d2212, T *d2222, const int i=1)
Definition: Ddg_pointer.hpp:33
Ddg(const Ddg< PackPtr< T *, I >, Tensor_Dim01, Tensor_Dim23 > &)=delete
T operator()(const int N1, const int N2, const int N3, const int N4) const
T & operator()(const int N1, const int N2, const int N3, const int N4)
Definition: Ddg_pointer.hpp:86
Dg_Expr< Ddg_numeral_0< const Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T >, T, Dim23, Dim1, i, j, k > operator()(const int N0, const Index< k, Dim1 > index3, const Index< i, Dim23 > index1, const Index< j, Dim23 > index2) const
Ddg(T *d0000, T *d0010, T *d0011, T *d1000, T *d1010, T *d1011, T *d1100, T *d1110, T *d1111, const int i=1)
Definition: Ddg_pointer.hpp:18
Tensor2_symmetric_Expr< Ddg_numeral_23< Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T >, T, Dim, i, j > operator()(const Index< i, Dim > index1, const Index< j, Dim > index2, const int N2, const int N3) const
Tensor2_symmetric_Expr< Ddg_numeral_01< const Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T >, T, Dim, i, j > operator()(const int N0, const int N1, const Index< i, Dim > index1, const Index< j, Dim > index2) const
Tensor2_symmetric_Expr< Ddg_number_rhs_01< Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T, N0, N1 >, T, Dim, i, j > operator()(const Number< N0 > n1, const Number< N1 > n2, const Index< i, Dim > index1, const Index< j, Dim > index2)
Dg_Expr< Ddg_number_rhs_0< Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T, N0 >, T, Dim23, Dim1, i, j, k > operator()(const Number< N0 > n1, const Index< k, Dim1 > index3, const Index< i, Dim23 > index1, const Index< j, Dim23 > index2)
Dg_Expr< Ddg_number_0< const Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T, N0 >, T, Dim23, Dim1, i, j, k > operator()(const Number< N0 > n1, const Index< k, Dim1 > index3, const Index< i, Dim23 > index1, const Index< j, Dim23 > index2) const
Ddg_Expr< const Ddg< T *, Tensor_Dim01, Tensor_Dim23 >, T, Dim01, Dim23, i, j, k, l > operator()(const Index< i, Dim01 > index1, const Index< j, Dim01 > index2, const Index< k, Dim23 > index3, const Index< l, Dim23 > index4) 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
constexpr IntegrationType I