GMatTensor 0.10.6
Loading...
Searching...
No Matches
detail.hpp
Go to the documentation of this file.
1
9#ifndef GMATTENSOR_DETAIL_HPP
10#define GMATTENSOR_DETAIL_HPP
11
12#include "config.h"
13
14namespace GMatTensor {
15namespace detail {
16
17// --------------------------
18// Array of 2nd-order tensors
19// --------------------------
20
21template <class T, size_t nd>
22struct impl_A2 {
23 using value_type = typename T::value_type;
24 using shape_type = typename T::shape_type;
25 static_assert(xt::has_fixed_rank_t<T>::value, "Only fixed rank allowed.");
26 static_assert(xt::get_rank<T>::value >= 2, "Rank too low.");
27 constexpr static size_t rank = xt::get_rank<T>::value - 2; // rank of the 'pure' matrix
28 constexpr static size_t stride0 = 1;
29 constexpr static size_t stride2 = nd * nd;
30 constexpr static size_t stride4 = nd * nd * nd * nd;
31
32 // Get the size of the underlying array (strips the last to dimensions)
33
34 template <class S>
35 static size_t toSizeT0(const S& shape)
36 {
37 using ST = typename S::value_type;
38 return std::accumulate(shape.cbegin(), shape.cbegin() + rank, ST(1), std::multiplies<ST>());
39 }
40
41 // Extract shape of the underlying array (strips the last to dimensions)
42 // and optionally add the dimensions
43
44 template <class S>
45 static std::array<size_t, rank> toShapeT0(const S& shape)
46 {
47 std::array<size_t, rank> ret;
48 std::copy(shape.cbegin(), shape.cbegin() + rank, ret.begin());
49 return ret;
50 }
51
52 template <class S>
53 static std::array<size_t, rank + 2> toShapeT2(const S& shape)
54 {
55 std::array<size_t, rank + 2> ret;
56 std::copy(shape.cbegin(), shape.cbegin() + rank, ret.begin());
57 ret[rank + 0] = nd;
58 ret[rank + 1] = nd;
59 return ret;
60 }
61
62 template <class S>
63 static std::array<size_t, rank + 4> toShapeT4(const S& shape)
64 {
65 std::array<size_t, rank + 4> ret;
66 std::copy(shape.cbegin(), shape.cbegin() + rank, ret.begin());
67 ret[rank + 0] = nd;
68 ret[rank + 1] = nd;
69 ret[rank + 2] = nd;
70 ret[rank + 3] = nd;
71 return ret;
72 }
73
78 template <class R, typename F>
79 static void ret0(const T& A, R& ret, F func)
80 {
81 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT2(A.shape())));
82 GMATTENSOR_ASSERT(xt::has_shape(ret, toShapeT0(A.shape())));
83#pragma omp parallel for
84 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
85 ret.flat(i) = func(&A.flat(i * stride2));
86 }
87 }
88
93 template <typename F>
94 static auto ret0(const T& A, F func) -> typename allocate<rank, T>::type
95 {
96 using return_type = typename allocate<rank, T>::type;
97 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT2(A.shape())));
98 return_type ret = return_type::from_shape(toShapeT0(A.shape()));
99#pragma omp parallel for
100 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
101 ret.flat(i) = func(&A.flat(i * stride2));
102 }
103 return ret;
104 }
105
106 // Argument: A2, B2
107 // Return: scalar
108
109 template <class R, typename F>
110 static void B2_ret0(const T& A, const T& B, R& ret, F func)
111 {
112 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT2(A.shape())));
113 GMATTENSOR_ASSERT(xt::has_shape(B, toShapeT2(B.shape())));
114 GMATTENSOR_ASSERT(xt::has_shape(ret, toShapeT0(A.shape())));
115#pragma omp parallel for
116 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
117 ret.flat(i) = func(&A.flat(i * stride2), &B.flat(i * stride2));
118 }
119 }
120
121 template <typename F>
122 static auto B2_ret0(const T& A, const T& B, F func) -> typename allocate<rank, T>::type
123 {
124 using return_type = typename allocate<rank, T>::type;
125 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT2(A.shape())));
126 GMATTENSOR_ASSERT(xt::has_shape(B, toShapeT2(A.shape())));
127 return_type ret = return_type::from_shape(toShapeT0(A.shape()));
128#pragma omp parallel for
129 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
130 ret.flat(i) = func(&A.flat(i * stride2), &B.flat(i * stride2));
131 }
132 return ret;
133 }
134
135 // Argument: A2, B2
136 // Return: 2nd-order tensor
137
138 template <class R, typename F>
139 static void B2_ret2(const T& A, const T& B, R& ret, F func)
140 {
141 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT2(A.shape())));
142 GMATTENSOR_ASSERT(xt::has_shape(B, toShapeT2(B.shape())));
143 GMATTENSOR_ASSERT(xt::has_shape(ret, toShapeT2(A.shape())));
144#pragma omp parallel for
145 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
146 func(&A.flat(i * stride2), &B.flat(i * stride2), &ret.flat(i * stride2));
147 }
148 }
149
150 template <typename F>
151 static auto B2_ret2(const T& A, const T& B, F func) -> typename allocate<rank + 2, T>::type
152 {
153 using return_type = typename allocate<rank + 2, T>::type;
154 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT2(A.shape())));
155 GMATTENSOR_ASSERT(xt::has_shape(B, toShapeT2(A.shape())));
156 return_type ret = return_type::from_shape(toShapeT2(A.shape()));
157#pragma omp parallel for
158 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
159 func(&A.flat(i * stride2), &B.flat(i * stride2), &ret.flat(i * stride2));
160 }
161 return ret;
162 }
163
164 // Argument: A2, B2
165 // Return: 4th-order tensor
166
167 template <class R, typename F>
168 static void B2_ret4(const T& A, const T& B, R& ret, F func)
169 {
170 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT2(A.shape())));
171 GMATTENSOR_ASSERT(xt::has_shape(B, toShapeT2(B.shape())));
172 GMATTENSOR_ASSERT(xt::has_shape(ret, toShapeT4(A.shape())));
173#pragma omp parallel for
174 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
175 func(&A.flat(i * stride2), &B.flat(i * stride2), &ret.flat(i * stride4));
176 }
177 }
178
179 template <typename F>
180 static auto B2_ret4(const T& A, const T& B, F func) -> typename allocate<rank + 4, T>::type
181 {
182 using return_type = typename allocate<rank + 4, T>::type;
183 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT2(A.shape())));
184 GMATTENSOR_ASSERT(xt::has_shape(B, toShapeT2(A.shape())));
185 return_type ret = return_type::from_shape(toShapeT4(A.shape()));
186#pragma omp parallel for
187 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
188 func(&A.flat(i * stride2), &B.flat(i * stride2), &ret.flat(i * stride4));
189 }
190 return ret;
191 }
192
193 // Argument: A2
194 // Return: 2nd-order tensor
195
196 template <class R, typename F>
197 static void ret2(const T& A, R& ret, F func)
198 {
199 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT2(A.shape())));
200 GMATTENSOR_ASSERT(xt::has_shape(ret, toShapeT2(A.shape())));
201#pragma omp parallel for
202 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
203 func(&A.flat(i * stride2), &ret.flat(i * stride2));
204 }
205 }
206
207 template <typename F>
208 static auto ret2(const T& A, F func) -> typename allocate<rank + 2, T>::type
209 {
210 using return_type = typename allocate<rank + 2, T>::type;
211 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT2(A.shape())));
212 return_type ret = return_type::from_shape(toShapeT2(A.shape()));
213#pragma omp parallel for
214 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
215 func(&A.flat(i * stride2), &ret.flat(i * stride2));
216 }
217 return ret;
218 }
219};
220
221// --------------------------
222// Array of 4th-order tensors
223// --------------------------
224
225template <class T, size_t nd>
226struct impl_A4 {
227 using value_type = typename T::value_type;
228 using shape_type = typename T::shape_type;
229 static_assert(xt::has_fixed_rank_t<T>::value, "Only fixed rank allowed.");
230 static_assert(xt::get_rank<T>::value >= 4, "Rank too low.");
231 constexpr static size_t rank = xt::get_rank<T>::value - 4;
232 constexpr static size_t stride0 = 1;
233 constexpr static size_t stride2 = nd * nd;
234 constexpr static size_t stride4 = nd * nd * nd * nd;
235
236 // Get the size of the underlying array (strips the last to dimensions)
237
238 template <class S>
239 static size_t toSizeT0(const S& shape)
240 {
241 using ST = typename S::value_type;
242 return std::accumulate(shape.cbegin(), shape.cbegin() + rank, ST(1), std::multiplies<ST>());
243 }
244
245 // Extract shape of the underlying array (strips the last to dimensions)
246 // and optionally add the dimensions
247
248 template <class S>
249 static std::array<size_t, rank> toShapeT0(const S& shape)
250 {
251 std::array<size_t, rank> ret;
252 std::copy(shape.cbegin(), shape.cbegin() + rank, ret.begin());
253 return ret;
254 }
255
256 template <class S>
257 static std::array<size_t, rank + 2> toShapeT2(const S& shape)
258 {
259 std::array<size_t, rank + 2> ret;
260 std::copy(shape.cbegin(), shape.cbegin() + rank, ret.begin());
261 ret[rank + 0] = nd;
262 ret[rank + 1] = nd;
263 return ret;
264 }
265
266 template <class S>
267 static std::array<size_t, rank + 4> toShapeT4(const S& shape)
268 {
269 std::array<size_t, rank + 4> ret;
270 std::copy(shape.cbegin(), shape.cbegin() + rank, ret.begin());
271 ret[rank + 0] = nd;
272 ret[rank + 1] = nd;
273 ret[rank + 2] = nd;
274 ret[rank + 3] = nd;
275 return ret;
276 }
277
278 // Argument: A4, B2
279 // Return: 2nd-order tensor
280
281 template <class U, class R, typename F>
282 static void B2_ret2(const T& A, const U& B, R& ret, F func)
283 {
284 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT4(A.shape())));
285 GMATTENSOR_ASSERT(xt::has_shape(B, toShapeT2(A.shape())));
286 GMATTENSOR_ASSERT(xt::has_shape(ret, toShapeT2(A.shape())));
287#pragma omp parallel for
288 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
289 func(&A.flat(i * stride4), &B.flat(i * stride2), &ret.flat(i * stride2));
290 }
291 }
292
293 template <class U, typename F>
294 static auto B2_ret2(const T& A, const U& B, F func) -> typename allocate<rank + 2, T>::type
295 {
296 using return_type = typename allocate<rank + 2, T>::type;
297 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT4(A.shape())));
298 GMATTENSOR_ASSERT(xt::has_shape(B, toShapeT2(A.shape())));
299 return_type ret = return_type::from_shape(toShapeT2(A.shape()));
300#pragma omp parallel for
301 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
302 func(&A.flat(i * stride4), &B.flat(i * stride2), &ret.flat(i * stride2));
303 }
304 return ret;
305 }
306
307 // Argument: A4, B2
308 // Return: 4th-order tensor
309
310 template <class U, class R, typename F>
311 static void B2_ret4(const T& A, const U& B, R& ret, F func)
312 {
313 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT4(A.shape())));
314 GMATTENSOR_ASSERT(xt::has_shape(B, toShapeT2(A.shape())));
315 GMATTENSOR_ASSERT(xt::has_shape(ret, toShapeT4(A.shape())));
316#pragma omp parallel for
317 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
318 func(&A.flat(i * stride4), &B.flat(i * stride2), &ret.flat(i * stride4));
319 }
320 }
321
322 template <class U, typename F>
323 static auto B2_ret4(const T& A, const U& B, F func) -> typename allocate<rank + 4, T>::type
324 {
325 using return_type = typename allocate<rank + 4, T>::type;
326 GMATTENSOR_ASSERT(xt::has_shape(A, toShapeT4(A.shape())));
327 GMATTENSOR_ASSERT(xt::has_shape(B, toShapeT2(A.shape())));
328 return_type ret = return_type::from_shape(toShapeT4(A.shape()));
329#pragma omp parallel for
330 for (size_t i = 0; i < toSizeT0(A.shape()); ++i) {
331 func(&A.flat(i * stride4), &B.flat(i * stride2), &ret.flat(i * stride4));
332 }
333 return ret;
334 }
335};
336
337} // namespace detail
338} // namespace GMatTensor
339
340#endif
config.h
Macros used in the library.
GMATTENSOR_ASSERT
#define GMATTENSOR_ASSERT(expr)
All assertions are implementation as:
Definition: config.h:59
GMatTensor
Tensor products / operations.
Definition: Cartesian2d.h:20