1 //
2 // Copyright (c) 2018-2019, Cem Bassoy, [email protected]
3 //
4 // Distributed under the Boost Software License, Version 1.0. (See
5 // accompanying file LICENSE_1_0.txt or copy at
6 // http://www.boost.org/LICENSE_1_0.txt)
7 //
8 // The authors gratefully acknowledge the support of
9 // Fraunhofer IOSB, Ettlingen, Germany
10 //
11
12
13 #ifndef _BOOST_UBLAS_TENSOR_ALGORITHMS_HPP
14 #define _BOOST_UBLAS_TENSOR_ALGORITHMS_HPP
15
16
17 #include <stdexcept>
18 #include <complex>
19 #include <functional>
20
21 namespace boost {
22 namespace numeric {
23 namespace ublas {
24
25
26
27 /** @brief Copies a tensor to another tensor with different layouts
28 *
29 * Implements C[i1,i2,...,ip] = A[i1,i2,...,ip]
30 *
31 * @param[in] p rank of input and output tensor
32 * @param[in] n pointer to the extents of input or output tensor of length p
33 * @param[in] pi pointer to a one-based permutation tuple of length p
34 * @param[out] c pointer to the output tensor
35 * @param[in] wc pointer to the strides of output tensor c
36 * @param[in] a pointer to the input tensor
37 * @param[in] wa pointer to the strides of input tensor a
38 */
39 template <class PointerOut, class PointerIn, class SizeType>
copy(const SizeType p,SizeType const * const n,PointerOut c,SizeType const * const wc,PointerIn a,SizeType const * const wa)40 void copy(const SizeType p, SizeType const*const n,
41 PointerOut c, SizeType const*const wc,
42 PointerIn a, SizeType const*const wa)
43 {
44 static_assert( std::is_pointer<PointerOut>::value & std::is_pointer<PointerIn>::value,
45 "Static error in boost::numeric::ublas::copy: Argument types for pointers are not pointer types.");
46 if( p == 0 )
47 return;
48
49 if(c == nullptr || a == nullptr)
50 throw std::length_error("Error in boost::numeric::ublas::copy: Pointers shall not be null pointers.");
51
52 if(wc == nullptr || wa == nullptr)
53 throw std::length_error("Error in boost::numeric::ublas::copy: Pointers shall not be null pointers.");
54
55 if(n == nullptr)
56 throw std::length_error("Error in boost::numeric::ublas::copy: Pointers shall not be null pointers.");
57
58
59 std::function<void(SizeType r, PointerOut c, PointerIn a)> lambda;
60
61 lambda = [&lambda, n, wc, wa](SizeType r, PointerOut c, PointerIn a)
62 {
63 if(r > 0)
64 for(auto d = 0u; d < n[r]; c += wc[r], a += wa[r], ++d)
65 lambda(r-1, c, a );
66 else
67 for(auto d = 0u; d < n[0]; c += wc[0], a += wa[0], ++d)
68 *c = *a;
69 };
70
71 lambda( p-1, c, a );
72 }
73
74
75
76 /** @brief Copies a tensor to another tensor with different layouts applying a unary operation
77 *
78 * Implements C[i1,i2,...,ip] = op ( A[i1,i2,...,ip] )
79 *
80 * @param[in] p rank of input and output tensor
81 * @param[in] n pointer to the extents of input or output tensor of length p
82 * @param[in] pi pointer to a one-based permutation tuple of length p
83 * @param[out] c pointer to the output tensor
84 * @param[in] wc pointer to the strides of output tensor c
85 * @param[in] a pointer to the input tensor
86 * @param[in] wa pointer to the strides of input tensor a
87 * @param[in] op unary operation
88 */
89 template <class PointerOut, class PointerIn, class SizeType, class UnaryOp>
transform(const SizeType p,SizeType const * const n,PointerOut c,SizeType const * const wc,PointerIn a,SizeType const * const wa,UnaryOp op)90 void transform(const SizeType p,
91 SizeType const*const n,
92 PointerOut c, SizeType const*const wc,
93 PointerIn a, SizeType const*const wa,
94 UnaryOp op)
95 {
96 static_assert( std::is_pointer<PointerOut>::value & std::is_pointer<PointerIn>::value,
97 "Static error in boost::numeric::ublas::transform: Argument types for pointers are not pointer types.");
98 if( p == 0 )
99 return;
100
101 if(c == nullptr || a == nullptr)
102 throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
103
104 if(wc == nullptr || wa == nullptr)
105 throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
106
107 if(n == nullptr)
108 throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
109
110
111 std::function<void(SizeType r, PointerOut c, PointerIn a)> lambda;
112
113 lambda = [&lambda, n, wc, wa, op](SizeType r, PointerOut c, PointerIn a)
114 {
115 if(r > 0)
116 for(auto d = 0u; d < n[r]; c += wc[r], a += wa[r], ++d)
117 lambda(r-1, c, a);
118 else
119 for(auto d = 0u; d < n[0]; c += wc[0], a += wa[0], ++d)
120 *c = op(*a);
121 };
122
123 lambda( p-1, c, a );
124
125 }
126
127
128 /** @brief Performs a reduce operation with all elements of the tensor and an initial value
129 *
130 * Implements k = sum_{i1,..,ip} A[i1,i2,...,ip]
131 *
132 * @param[in] r zero-based recursion level starting with r=p-1
133 * @param[in] n pointer to the extents of input or output tensor
134 * @param[in] a pointer to the first input tensor
135 * @param[in] w pointer to the strides of input tensor a
136 * @param[in] k accumulated value
137 */
138 template <class PointerIn, class ValueType, class SizeType>
accumulate(SizeType const p,SizeType const * const n,PointerIn a,SizeType const * const w,ValueType k)139 ValueType accumulate(SizeType const p, SizeType const*const n,
140 PointerIn a, SizeType const*const w,
141 ValueType k)
142 {
143 static_assert(std::is_pointer<PointerIn>::value,
144 "Static error in boost::numeric::ublas::transform: Argument types for pointers are not pointer types.");
145
146 if( p == 0 )
147 return k;
148
149 if(a == nullptr)
150 throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
151
152 if(w == nullptr)
153 throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
154
155 if(n == nullptr)
156 throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
157
158
159 std::function<ValueType(SizeType r, PointerIn a, ValueType k)> lambda;
160
161 lambda = [&lambda, n, w](SizeType r, PointerIn a, ValueType k)
162 {
163 if(r > 0u)
164 for(auto d = 0u; d < n[r]; a += w[r], ++d)
165 k = lambda(r-1, a, k);
166 else
167 for(auto d = 0u; d < n[0]; a += w[0], ++d)
168 k += *a;
169 return k;
170 };
171
172 return lambda( p-1, a, k );
173 }
174
175 /** @brief Performs a reduce operation with all elements of the tensor and an initial value
176 *
177 * Implements k = op ( k , A[i1,i2,...,ip] ), for all ir
178 *
179 * @param[in] r zero-based recursion level starting with r=p-1
180 * @param[in] n pointer to the extents of input or output tensor
181 * @param[in] a pointer to the first input tensor
182 * @param[in] w pointer to the strides of input tensor a
183 * @param[in] k accumulated value
184 * @param[in] op binary operation
185 */
186
187 template <class PointerIn, class ValueType, class SizeType, class BinaryOp>
accumulate(SizeType const p,SizeType const * const n,PointerIn a,SizeType const * const w,ValueType k,BinaryOp op)188 ValueType accumulate(SizeType const p, SizeType const*const n,
189 PointerIn a, SizeType const*const w,
190 ValueType k, BinaryOp op)
191 {
192 static_assert(std::is_pointer<PointerIn>::value,
193 "Static error in boost::numeric::ublas::transform: Argument types for pointers are not pointer types.");
194
195
196 if( p == 0 )
197 return k;
198
199 if(a == nullptr)
200 throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
201
202 if(w == nullptr)
203 throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
204
205 if(n == nullptr)
206 throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
207
208
209 std::function<ValueType(SizeType r, PointerIn a, ValueType k)> lambda;
210
211 lambda = [&lambda, n, w, op](SizeType r, PointerIn a, ValueType k)
212 {
213 if(r > 0u)
214 for(auto d = 0u; d < n[r]; a += w[r], ++d)
215 k = lambda(r-1, a, k);
216 else
217 for(auto d = 0u; d < n[0]; a += w[0], ++d)
218 k = op ( k, *a );
219 return k;
220 };
221
222 return lambda( p-1, a, k );
223 }
224
225 /** @brief Transposes a tensor
226 *
227 * Implements C[tau[i1],tau[i2],...,tau[ip]] = A[i1,i2,...,ip]
228 *
229 * @note is used in function trans
230 *
231 * @param[in] p rank of input and output tensor
232 * @param[in] na pointer to the extents of the input tensor a of length p
233 * @param[in] pi pointer to a one-based permutation tuple of length p
234 * @param[out] c pointer to the output tensor
235 * @param[in] wc pointer to the strides of output tensor c
236 * @param[in] a pointer to the input tensor
237 * @param[in] wa pointer to the strides of input tensor a
238 */
239
240 template <class PointerOut, class PointerIn, class SizeType>
trans(SizeType const p,SizeType const * const na,SizeType const * const pi,PointerOut c,SizeType const * const wc,PointerIn a,SizeType const * const wa)241 void trans( SizeType const p, SizeType const*const na, SizeType const*const pi,
242 PointerOut c, SizeType const*const wc,
243 PointerIn a, SizeType const*const wa)
244 {
245
246 static_assert( std::is_pointer<PointerOut>::value & std::is_pointer<PointerIn>::value,
247 "Static error in boost::numeric::ublas::trans: Argument types for pointers are not pointer types.");
248
249 if( p < 2)
250 return;
251
252 if(c == nullptr || a == nullptr)
253 throw std::runtime_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
254
255 if(na == nullptr)
256 throw std::runtime_error("Error in boost::numeric::ublas::trans: Pointers shall not be null.");
257
258 if(wc == nullptr || wa == nullptr)
259 throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
260
261 if(na == nullptr)
262 throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
263
264 if(pi == nullptr)
265 throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
266
267
268 std::function<void(SizeType r, PointerOut c, PointerIn a)> lambda;
269
270 lambda = [&lambda, na, wc, wa, pi](SizeType r, PointerOut c, PointerIn a)
271 {
272 if(r > 0)
273 for(auto d = 0u; d < na[r]; c += wc[pi[r]-1], a += wa[r], ++d)
274 lambda(r-1, c, a);
275 else
276 for(auto d = 0u; d < na[0]; c += wc[pi[0]-1], a += wa[0], ++d)
277 *c = *a;
278 };
279
280 lambda( p-1, c, a );
281 }
282
283
284 /** @brief Transposes a tensor
285 *
286 * Implements C[tau[i1],tau[i2],...,tau[ip]] = A[i1,i2,...,ip]
287 *
288 * @note is used in function trans
289 *
290 * @param[in] p rank of input and output tensor
291 * @param[in] na pointer to the extents of the input tensor a of length p
292 * @param[in] pi pointer to a one-based permutation tuple of length p
293 * @param[out] c pointer to the output tensor
294 * @param[in] wc pointer to the strides of output tensor c
295 * @param[in] a pointer to the input tensor
296 * @param[in] wa pointer to the strides of input tensor a
297 */
298
299 template <class ValueType, class SizeType>
trans(SizeType const p,SizeType const * const na,SizeType const * const pi,std::complex<ValueType> * c,SizeType const * const wc,std::complex<ValueType> * a,SizeType const * const wa)300 void trans( SizeType const p,
301 SizeType const*const na,
302 SizeType const*const pi,
303 std::complex<ValueType>* c, SizeType const*const wc,
304 std::complex<ValueType>* a, SizeType const*const wa)
305 {
306 if( p < 2)
307 return;
308
309 if(c == nullptr || a == nullptr)
310 throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
311
312 if(wc == nullptr || wa == nullptr)
313 throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
314
315 if(na == nullptr)
316 throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
317
318 if(pi == nullptr)
319 throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
320
321
322 std::function<void(SizeType r, std::complex<ValueType>* c, std::complex<ValueType>* a)> lambda;
323
324 lambda = [&lambda, na, wc, wa, pi](SizeType r, std::complex<ValueType>* c, std::complex<ValueType>* a)
325 {
326 if(r > 0)
327 for(auto d = 0u; d < na[r]; c += wc[pi[r]-1], a += wa[r], ++d)
328 lambda(r-1, c, a);
329 else
330 for(auto d = 0u; d < na[0]; c += wc[pi[0]-1], a += wa[0], ++d)
331 *c = std::conj(*a);
332 };
333
334 lambda( p-1, c, a );
335
336 }
337
338
339
340
341 }
342 }
343 }
344
345 #endif
346