ublas/tensor/algorithms.hpp

0001 //
0002 //  Copyright (c) 2018-2019, Cem Bassoy, cem.bassoy@gmail.com
0003 //
0004 //  Distributed under the Boost Software License, Version 1.0. (See
0005 //  accompanying file LICENSE_1_0.txt or copy at
0006 //  http://www.boost.org/LICENSE_1_0.txt)
0007 //
0008 //  The authors gratefully acknowledge the support of
0009 //  Fraunhofer IOSB, Ettlingen, Germany
0010 //
0011
0012
0013 #ifndef _BOOST_UBLAS_TENSOR_ALGORITHMS_HPP
0014 #define _BOOST_UBLAS_TENSOR_ALGORITHMS_HPP
0015
0016
0017 #include <stdexcept>
0018 #include <complex>
0019 #include <functional>
0020
0021 namespace boost {
0022 namespace numeric {
0023 namespace ublas {
0024
0025
0026
0027 /** @brief Copies a tensor to another tensor with different layouts
0028  *
0029  * Implements C[i1,i2,...,ip] = A[i1,i2,...,ip]
0030  *
0031  * @param[in]  p rank of input and output tensor
0032  * @param[in]  n pointer to the extents of input or output tensor of length p
0033  * @param[in] pi pointer to a one-based permutation tuple of length p
0034  * @param[out] c pointer to the output tensor
0035  * @param[in] wc pointer to the strides of output tensor c
0036  * @param[in]  a pointer to the input tensor
0037  * @param[in] wa pointer to the strides of input tensor a
0038 */
0039 template <class PointerOut, class PointerIn, class SizeType>
0040 void copy(const SizeType p, SizeType const*const n,
0041           PointerOut c, SizeType const*const wc,
0042           PointerIn a,  SizeType const*const wa)
0043 {
0044     static_assert( std::is_pointer<PointerOut>::value & std::is_pointer<PointerIn>::value,
0045                    "Static error in boost::numeric::ublas::copy: Argument types for pointers are not pointer types.");
0046     if( p == 0 )
0047         return;
0048
0049     if(c == nullptr || a == nullptr)
0050         throw std::length_error("Error in boost::numeric::ublas::copy: Pointers shall not be null pointers.");
0051
0052     if(wc == nullptr || wa == nullptr)
0053         throw std::length_error("Error in boost::numeric::ublas::copy: Pointers shall not be null pointers.");
0054
0055     if(n == nullptr)
0056         throw std::length_error("Error in boost::numeric::ublas::copy: Pointers shall not be null pointers.");
0057
0058
0059     std::function<void(SizeType r, PointerOut c, PointerIn  a)> lambda;
0060
0061     lambda = [&lambda, n, wc, wa](SizeType r, PointerOut c, PointerIn a)
0062     {
0063         if(r > 0)
0064             for(auto d = 0u; d < n[r]; c += wc[r], a += wa[r], ++d)
0065                 lambda(r-1, c, a );
0066         else
0067             for(auto d = 0u; d < n[0]; c += wc[0], a += wa[0], ++d)
0068                 *c = *a;
0069     };
0070
0071     lambda( p-1, c, a );
0072 }
0073
0074
0075
0076 /** @brief Copies a tensor to another tensor with different layouts applying a unary operation
0077  *
0078  * Implements C[i1,i2,...,ip] = op ( A[i1,i2,...,ip] )
0079  *
0080  * @param[in]  p rank of input and output tensor
0081  * @param[in]  n pointer to the extents of input or output tensor of length p
0082  * @param[in] pi pointer to a one-based permutation tuple of length p
0083  * @param[out] c pointer to the output tensor
0084  * @param[in] wc pointer to the strides of output tensor c
0085  * @param[in]  a pointer to the input tensor
0086  * @param[in] wa pointer to the strides of input tensor a
0087  * @param[in] op unary operation
0088 */
0089 template <class PointerOut, class PointerIn, class SizeType, class UnaryOp>
0090 void transform(const SizeType p,
0091                SizeType const*const n,
0092                PointerOut c, SizeType const*const wc,
0093                PointerIn a,  SizeType const*const wa,
0094                UnaryOp op)
0095 {
0096     static_assert( std::is_pointer<PointerOut>::value & std::is_pointer<PointerIn>::value,
0097                    "Static error in boost::numeric::ublas::transform: Argument types for pointers are not pointer types.");
0098     if( p == 0 )
0099         return;
0100
0101     if(c == nullptr || a == nullptr)
0102         throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
0103
0104     if(wc == nullptr || wa == nullptr)
0105         throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
0106
0107     if(n == nullptr)
0108         throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
0109
0110
0111     std::function<void(SizeType r, PointerOut c, PointerIn a)> lambda;
0112
0113     lambda = [&lambda, n, wc, wa, op](SizeType r, PointerOut c, PointerIn a)
0114     {
0115         if(r > 0)
0116             for(auto d = 0u; d < n[r]; c += wc[r], a += wa[r], ++d)
0117                 lambda(r-1, c, a);
0118         else
0119             for(auto d = 0u; d < n[0]; c += wc[0], a += wa[0], ++d)
0120                 *c = op(*a);
0121     };
0122
0123     lambda( p-1, c, a );
0124
0125 }
0126
0127
0128 /** @brief Performs a reduce operation with all elements of the tensor and an initial value
0129  *
0130  * Implements k = sum_{i1,..,ip} A[i1,i2,...,ip]
0131  *
0132  * @param[in] r  zero-based recursion level starting with r=p-1
0133  * @param[in] n  pointer to the extents of input or output tensor
0134  * @param[in] a  pointer to the first input tensor
0135  * @param[in] w  pointer to the strides of input tensor a
0136  * @param[in] k  accumulated value
0137 */
0138 template <class PointerIn, class ValueType, class SizeType>
0139 ValueType accumulate(SizeType const p, SizeType const*const n,
0140                      PointerIn a, SizeType const*const w,
0141                      ValueType k)
0142 {
0143     static_assert(std::is_pointer<PointerIn>::value,
0144                   "Static error in boost::numeric::ublas::transform: Argument types for pointers are not pointer types.");
0145
0146     if( p == 0 )
0147         return k;
0148
0149     if(a == nullptr)
0150         throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
0151
0152     if(w == nullptr)
0153         throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
0154
0155     if(n == nullptr)
0156         throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
0157
0158
0159     std::function<ValueType(SizeType r, PointerIn a, ValueType k)> lambda;
0160
0161     lambda = [&lambda, n, w](SizeType r, PointerIn a, ValueType k)
0162     {
0163         if(r > 0u)
0164             for(auto d = 0u; d < n[r]; a += w[r], ++d)
0165                 k = lambda(r-1, a, k);
0166         else
0167             for(auto d = 0u; d < n[0]; a += w[0], ++d)
0168                 k += *a;
0169         return k;
0170     };
0171
0172     return lambda( p-1, a,  k );
0173 }
0174
0175 /** @brief Performs a reduce operation with all elements of the tensor and an initial value
0176  *
0177  * Implements k = op ( k , A[i1,i2,...,ip] ), for all ir
0178  *
0179  * @param[in] r  zero-based recursion level starting with r=p-1
0180  * @param[in] n  pointer to the extents of input or output tensor
0181  * @param[in] a  pointer to the first input tensor
0182  * @param[in] w  pointer to the strides of input tensor a
0183  * @param[in] k  accumulated value
0184  * @param[in] op binary operation
0185 */
0186
0187 template <class PointerIn, class ValueType, class SizeType, class BinaryOp>
0188 ValueType accumulate(SizeType const p, SizeType const*const n,
0189                      PointerIn a, SizeType const*const w,
0190                      ValueType k, BinaryOp op)
0191 {
0192     static_assert(std::is_pointer<PointerIn>::value,
0193                   "Static error in boost::numeric::ublas::transform: Argument types for pointers are not pointer types.");
0194
0195
0196     if( p == 0 )
0197         return k;
0198
0199     if(a == nullptr)
0200         throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
0201
0202     if(w == nullptr)
0203         throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
0204
0205     if(n == nullptr)
0206         throw std::length_error("Error in boost::numeric::ublas::transform: Pointers shall not be null pointers.");
0207
0208
0209     std::function<ValueType(SizeType r, PointerIn a, ValueType k)> lambda;
0210
0211     lambda = [&lambda, n, w, op](SizeType r, PointerIn a, ValueType k)
0212     {
0213         if(r > 0u)
0214             for(auto d = 0u; d < n[r]; a += w[r], ++d)
0215                 k = lambda(r-1, a, k);
0216         else
0217             for(auto d = 0u; d < n[0]; a += w[0], ++d)
0218                 k = op ( k, *a );
0219         return k;
0220     };
0221
0222     return lambda( p-1, a,  k );
0223 }
0224
0225 /** @brief Transposes a tensor
0226  *
0227  * Implements C[tau[i1],tau[i2],...,tau[ip]] = A[i1,i2,...,ip]
0228  *
0229  * @note is used in function trans
0230  *
0231  * @param[in]  p rank of input and output tensor
0232  * @param[in] na pointer to the extents of the input tensor a of length p
0233  * @param[in] pi pointer to a one-based permutation tuple of length p
0234  * @param[out] c pointer to the output tensor
0235  * @param[in] wc pointer to the strides of output tensor c
0236  * @param[in]  a pointer to the input tensor
0237  * @param[in] wa pointer to the strides of input tensor a
0238 */
0239
0240 template <class PointerOut, class PointerIn, class SizeType>
0241 void trans( SizeType const p,  SizeType const*const na, SizeType const*const pi,
0242             PointerOut c,      SizeType const*const wc,
0243             PointerIn a,       SizeType const*const wa)
0244 {
0245
0246     static_assert( std::is_pointer<PointerOut>::value & std::is_pointer<PointerIn>::value,
0247                    "Static error in boost::numeric::ublas::trans: Argument types for pointers are not pointer types.");
0248
0249     if( p < 2)
0250         return;
0251
0252     if(c == nullptr || a == nullptr)
0253         throw std::runtime_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
0254
0255     if(na == nullptr)
0256         throw std::runtime_error("Error in boost::numeric::ublas::trans: Pointers shall not be null.");
0257
0258     if(wc == nullptr || wa == nullptr)
0259         throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
0260
0261     if(na == nullptr)
0262         throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
0263
0264     if(pi == nullptr)
0265         throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
0266
0267
0268     std::function<void(SizeType r, PointerOut c, PointerIn a)> lambda;
0269
0270     lambda = [&lambda, na, wc, wa, pi](SizeType r, PointerOut c, PointerIn a)
0271     {
0272         if(r > 0)
0273             for(auto d = 0u; d < na[r]; c += wc[pi[r]-1], a += wa[r], ++d)
0274                 lambda(r-1, c, a);
0275         else
0276             for(auto d = 0u; d < na[0]; c += wc[pi[0]-1], a += wa[0], ++d)
0277                 *c = *a;
0278     };
0279
0280     lambda( p-1, c, a );
0281 }
0282
0283
0284 /** @brief Transposes a tensor
0285  *
0286  * Implements C[tau[i1],tau[i2],...,tau[ip]] = A[i1,i2,...,ip]
0287  *
0288  * @note is used in function trans
0289  *
0290  * @param[in]  p rank of input and output tensor
0291  * @param[in] na pointer to the extents of the input tensor a of length p
0292  * @param[in] pi pointer to a one-based permutation tuple of length p
0293  * @param[out] c pointer to the output tensor
0294  * @param[in] wc pointer to the strides of output tensor c
0295  * @param[in]  a pointer to the input tensor
0296  * @param[in] wa pointer to the strides of input tensor a
0297 */
0298
0299 template <class ValueType, class SizeType>
0300 void trans( SizeType const p,
0301             SizeType const*const na,
0302             SizeType const*const pi,
0303             std::complex<ValueType>* c,  SizeType const*const wc,
0304             std::complex<ValueType>* a,  SizeType const*const wa)
0305 {
0306     if( p < 2)
0307         return;
0308
0309     if(c == nullptr || a == nullptr)
0310         throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
0311
0312     if(wc == nullptr || wa == nullptr)
0313         throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
0314
0315     if(na == nullptr)
0316         throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
0317
0318     if(pi == nullptr)
0319         throw std::length_error("Error in boost::numeric::ublas::trans: Pointers shall not be null pointers.");
0320
0321
0322     std::function<void(SizeType r, std::complex<ValueType>* c, std::complex<ValueType>* a)> lambda;
0323
0324     lambda = [&lambda, na, wc, wa, pi](SizeType r, std::complex<ValueType>* c, std::complex<ValueType>* a)
0325     {
0326         if(r > 0)
0327             for(auto d = 0u; d < na[r]; c += wc[pi[r]-1], a += wa[r], ++d)
0328                 lambda(r-1, c, a);
0329         else
0330             for(auto d = 0u; d < na[0]; c += wc[pi[0]-1], a += wa[0], ++d)
0331                 *c = std::conj(*a);
0332     };
0333
0334     lambda( p-1, c, a );
0335
0336 }
0337
0338
0339
0340
0341 }
0342 }
0343 }
0344
0345 #endif