v0.14.0
Loading...
Searching...
No Matches
Tensor2_symmetric_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_Dim> class Tensor2_symmetric<T *, Tensor_Dim>
8 {
9 const int inc;
10
11 protected:
12
13 mutable T *restrict data[(Tensor_Dim * (Tensor_Dim + 1)) / 2];
14
15 public:
16 template <class... U>
17 explicit Tensor2_symmetric(U *... d) : data{d...}, inc(1) {
18 static_assert(sizeof...(d) == sizeof(data) / sizeof(T),
19 "Incorrect number of Arguments. Constructor should "
20 "initialize the entire Tensor");
21 }
22
23 template <class... U>
24 explicit Tensor2_symmetric(const int i, U *... d) : data{d...}, inc(i) {
25 static_assert(sizeof...(d) == sizeof(data) / sizeof(T),
26 "Incorrect number of Arguments. Constructor should "
27 "initialize the entire Tensor");
28 }
29
31
32 /* Tensor_Dim=2 */
33 Tensor2_symmetric(T *d00, T *d01, T *d11, const int i = 1) : inc(i)
34 {
36 d11);
37 }
38
39 /* Tensor_Dim=3 */
40 Tensor2_symmetric(T *d00, T *d01, T *d02, T *d11, T *d12, T *d22,
41 const int i = 1)
42 : inc(i) {
44 data, d00, d01, d02, d11, d12, d22);
45 }
46
47 /* There are two operator(int,int)'s, one for non-consts that lets you
48 change the value, and one for consts that doesn't. */
49
50 T &operator()(const int N1, const int N2)
51 {
52#ifdef FTENSOR_DEBUG
53 if(N1 >= Tensor_Dim || N1 < 0 || N2 >= Tensor_Dim || N2 < 0)
54 {
55 std::stringstream s;
56 s << "Bad index in Tensor2_symmetric<T*," << Tensor_Dim
57 << ">.operator(" << N1 << "," << N2 << ")" << std::endl;
58 throw std::out_of_range(s.str());
59 }
60#endif
61 return N1 > N2 ? *data[N1 + (N2 * (2 * Tensor_Dim - N2 - 1)) / 2]
62 : *data[N2 + (N1 * (2 * Tensor_Dim - N1 - 1)) / 2];
63 }
64
65 T operator()(const int N1, const int N2) const
66 {
67#ifdef FTENSOR_DEBUG
68 if(N1 >= Tensor_Dim || N1 < 0 || N2 >= Tensor_Dim || N2 < 0)
69 {
70 std::stringstream s;
71 s << "Bad index in Tensor2_symmetric<T*," << Tensor_Dim
72 << ">.operator(" << N1 << "," << N2 << ") const" << std::endl;
73 throw std::out_of_range(s.str());
74 }
75#endif
76 return N1 > N2 ? *data[N1 + (N2 * (2 * Tensor_Dim - N2 - 1)) / 2]
77 : *data[N2 + (N1 * (2 * Tensor_Dim - N1 - 1)) / 2];
78 }
79
80 T *ptr(const int N1, const int N2) const
81 {
82#ifdef FTENSOR_DEBUG
83 if(N1 >= Tensor_Dim || N1 < 0 || N2 >= Tensor_Dim || N2 < 0)
84 {
85 std::stringstream s;
86 s << "Bad index in Tensor2_symmetric<T*," << Tensor_Dim << ">.ptr("
87 << N1 << "," << N2 << ")" << std::endl;
88 throw std::out_of_range(s.str());
89 }
90#endif
91 return N1 > N2 ? data[N1 + (N2 * (2 * Tensor_Dim - N2 - 1)) / 2]
92 : data[N2 + (N1 * (2 * Tensor_Dim - N1 - 1)) / 2];
93 }
94
95 /* These operator()'s are the first part in constructing template
96 expressions. They can be used to slice off lower dimensional
97 parts. They are not entirely safe, since you can accidentaly use a
98 higher dimension than what is really allowed (like Dim=5). */
99
100 /* This returns a Tensor2_Expr, since the indices are not really
101 symmetric anymore since they cover different dimensions. */
102
103 template <char i, char j, int Dim0, int Dim1>
105 operator()(const Index<i, Dim0> index1, const Index<j, Dim1> index2)
106 {
108 j>(*this);
109 }
110
111 template <char i, char j, int Dim0, int Dim1>
113 operator()(const Index<i, Dim0> index1, const Index<j, Dim1> index2) const
114 {
116 Dim1, i, j>(*this);
117 }
118
119 /* This returns a Tensor2_symmetric_Expr, since the indices are still
120 symmetric on the lower dimensions. */
121
122 template <char i, char j, int Dim>
124 operator()(const Index<i, Dim> index1, const Index<j, Dim> index2)
125 {
127 i, j>(*this);
128 }
129
130 template <char i, char j, int Dim>
132 j>
133 operator()(const Index<i, Dim> index1, const Index<j, Dim> index2) const
134 {
136 T, Dim, i, j>(*this);
137 }
138
139 /* This is for expressions where a number is used for one slot, and
140 an index for another, yielding a Tensor1_Expr. The non-const
141 versions don't actually create a Tensor2_number_rhs_[01] object.
142 They create a Tensor1_Expr directly, which provides the
143 appropriate indexing operators. The const versions do create a
144 Tensor2_number_[01]. */
145
146 template <char i, int N, int Dim>
148 T, Dim, i>
149 operator()(const Index<i, Dim> index1, const Number<N> &n1)
150 {
151 using TensorExpr
154 }
155
156 template <char i, int N, int Dim>
159 T, Dim, i>
160 operator()(const Index<i, Dim> index1, const Number<N> &n1) const
161 {
162 using TensorExpr
164 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
165 }
166
167 template <char i, int N, int Dim>
169 T, Dim, i>
170 operator()(const Number<N> &n1, const Index<i, Dim> index1)
171 {
172 using TensorExpr
175 }
176
177 template <char i, int N, int Dim>
180 T, Dim, i>
181 operator()(const Number<N> &n1, const Index<i, Dim> index1) const
182 {
183 using TensorExpr
185 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this));
186 }
187
188 /* Specializations for using actual numbers instead of Number<> */
189
190 template <char i, int Dim>
193 Dim, i>
194 operator()(const Index<i, Dim> index1, const int N) const
195 {
196 using TensorExpr
198 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N));
199 }
200
201 template <char i, int Dim>
204 Dim, i>
205 operator()(const int N, const Index<i, Dim> index1) const
206 {
207 using TensorExpr
209 return Tensor1_Expr<TensorExpr, T, Dim, i>(TensorExpr(*this, N));
210 }
211
212 /* The ++ operator increments the pointer, not the number that the
213 pointer points to. This allows iterating over a grid. */
214
216 {
217 for(int i = 0; i < (Tensor_Dim * (Tensor_Dim + 1)) / 2; ++i)
218 data[i] += inc;
219 return *this;
220 }
221
222 /* These two operator()'s return the Tensor2 with internal
223 contractions, yielding a T. I have to specify one for both
224 const and non-const because otherwise the compiler will use the
225 operator() which gives a Tensor2_Expr<>. */
226
227 template <char i, int Dim>
228 T operator()(const Index<i, Dim> index1, const Index<i, Dim> index2)
229 {
230 return internal_contract(Number<Dim>());
231 }
232
233 template <char i, int Dim>
234 T operator()(const Index<i, Dim> index1, const Index<i, Dim> index2) const
235 {
236 return internal_contract(Number<Dim>());
237 }
238
239 private:
240 template <int N> T internal_contract(Number<N>) const
241 {
242 return *data[N - 1 + ((N - 1) * (2 * Tensor_Dim - N)) / 2]
243 + internal_contract(Number<N - 1>());
244 }
245
246 T internal_contract(Number<1>) const { return *data[0]; }
247 };
248
249 template <class T, int Tensor_Dim, int I>
250 class Tensor2_symmetric<PackPtr<T *, I>, Tensor_Dim>
252
253 public:
254 template <class... U>
255 Tensor2_symmetric(U *... d) : Tensor2_symmetric<T *, Tensor_Dim>(d...) {}
256
257 Tensor2_symmetric(): Tensor2_symmetric<T *, Tensor_Dim>() {}
258
259 /* The ++ operator increments the pointer, not the number that the
260 pointer points to. This allows iterating over a grid. */
261
262 const Tensor2_symmetric<PackPtr<T *, I>, Tensor_Dim> &
264 {
265 for(int i = 0; i < (Tensor_Dim * (Tensor_Dim + 1)) / 2; ++i)
267 return *this;
268 }
269
270 };
271
272}
static Number< 2 > N2
static Number< 1 > N1
const Tensor2_symmetric< PackPtr< T *, I >, Tensor_Dim > & operator++() const
Tensor1_Expr< const Tensor2_number_1< const Tensor2_symmetric< T *, Tensor_Dim >, T, N >, T, Dim, i > operator()(const Index< i, Dim > index1, const Number< N > &n1) const
Tensor2_symmetric_Expr< const Tensor2_symmetric< T *, Tensor_Dim >, T, Dim, i, j > operator()(const Index< i, Dim > index1, const Index< j, Dim > index2) const
Tensor1_Expr< const Tensor2_number_0< const Tensor2_symmetric< T *, Tensor_Dim >, T, N >, T, Dim, i > operator()(const Number< N > &n1, const Index< i, Dim > index1) const
Tensor1_Expr< Tensor2_number_rhs_1< Tensor2_symmetric< T *, Tensor_Dim >, T, N >, T, Dim, i > operator()(const Index< i, Dim > index1, const Number< N > &n1)
Tensor1_Expr< Tensor2_number_rhs_0< Tensor2_symmetric< T *, Tensor_Dim >, T, N >, T, Dim, i > operator()(const Number< N > &n1, const Index< i, Dim > index1)
T operator()(const Index< i, Dim > index1, const Index< i, Dim > index2)
Tensor2_Expr< const Tensor2_symmetric< T *, Tensor_Dim >, T, Dim0, Dim1, i, j > operator()(const Index< i, Dim0 > index1, const Index< j, Dim1 > index2) const
T operator()(const Index< i, Dim > index1, const Index< i, Dim > index2) const
const Tensor2_symmetric< T *, Tensor_Dim > & operator++() const
Tensor2_symmetric(T *d00, T *d01, T *d02, T *d11, T *d12, T *d22, const int i=1)
Tensor2_Expr< Tensor2_symmetric< T *, Tensor_Dim >, T, Dim0, Dim1, i, j > operator()(const Index< i, Dim0 > index1, const Index< j, Dim1 > index2)
Tensor1_Expr< const Tensor2_numeral_0< const Tensor2_symmetric< T *, Tensor_Dim >, T >, T, Dim, i > operator()(const int N, const Index< i, Dim > index1) const
Tensor2_symmetric_Expr< Tensor2_symmetric< T *, Tensor_Dim >, T, Dim, i, j > operator()(const Index< i, Dim > index1, const Index< j, Dim > index2)
Tensor2_symmetric(T *d00, T *d01, T *d11, const int i=1)
Tensor1_Expr< const Tensor2_numeral_1< const Tensor2_symmetric< T *, Tensor_Dim >, T >, T, Dim, i > operator()(const Index< i, Dim > index1, const int N) const
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'j', 3 > j
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