492 lines
13 KiB
C++
492 lines
13 KiB
C++
// Copyright (c) 2018-2019 Cem Bassoy
|
|
//
|
|
// Distributed under the Boost Software License, Version 1.0. (See
|
|
// accompanying file LICENSE_1_0.txt or copy at
|
|
// http://www.boost.org/LICENSE_1_0.txt)
|
|
//
|
|
// The authors gratefully acknowledge the support of
|
|
// Fraunhofer and Google in producing this work
|
|
// which started as a Google Summer of Code project.
|
|
//
|
|
|
|
|
|
#include <iostream>
|
|
#include <algorithm>
|
|
#include <vector>
|
|
|
|
#include <boost/numeric/ublas/tensor/multiplication.hpp>
|
|
#include <boost/numeric/ublas/tensor/extents.hpp>
|
|
#include <boost/numeric/ublas/tensor/strides.hpp>
|
|
#include "utility.hpp"
|
|
|
|
#include <boost/test/unit_test.hpp>
|
|
|
|
|
|
BOOST_AUTO_TEST_SUITE (test_tensor_contraction)
|
|
|
|
|
|
using test_types = zip<int,long,float,double,std::complex<float>>::with_t<boost::numeric::ublas::first_order, boost::numeric::ublas::last_order>;
|
|
|
|
//using test_types = zip<int>::with_t<boost::numeric::ublas::first_order>;
|
|
|
|
|
|
struct fixture
|
|
{
|
|
using extents_type = boost::numeric::ublas::shape;
|
|
fixture()
|
|
: extents {
|
|
extents_type{1,1}, // 1
|
|
extents_type{1,2}, // 2
|
|
extents_type{2,1}, // 3
|
|
extents_type{2,3}, // 4
|
|
extents_type{5,4}, // 5
|
|
extents_type{2,3,1}, // 6
|
|
extents_type{4,1,3}, // 7
|
|
extents_type{1,2,3}, // 8
|
|
extents_type{4,2,3}, // 9
|
|
extents_type{4,2,3,5}} // 10
|
|
{
|
|
}
|
|
std::vector<extents_type> extents;
|
|
};
|
|
|
|
|
|
|
|
BOOST_FIXTURE_TEST_CASE_TEMPLATE(test_tensor_mtv, value, test_types, fixture )
|
|
{
|
|
using namespace boost::numeric;
|
|
using value_type = typename value::first_type;
|
|
using layout_type = typename value::second_type;
|
|
using strides_type = ublas::strides<layout_type>;
|
|
using vector_type = std::vector<value_type>;
|
|
using extents_type = ublas::shape;
|
|
using extents_type_base = typename extents_type::base_type;
|
|
using size_type = typename extents_type_base::value_type;
|
|
|
|
|
|
for(auto const& na : extents) {
|
|
|
|
if(na.size() > 2)
|
|
continue;
|
|
|
|
auto a = vector_type(na.product(), value_type{2});
|
|
auto wa = strides_type(na);
|
|
for(auto m = 0u; m < na.size(); ++m){
|
|
auto nb = extents_type {na[m],1};
|
|
auto wb = strides_type (nb);
|
|
auto b = vector_type (nb.product(), value_type{1} );
|
|
|
|
auto nc_base = extents_type_base(std::max(na.size()-1, size_type{2}), 1);
|
|
|
|
for(auto i = 0u, j = 0u; i < na.size(); ++i)
|
|
if(i != m)
|
|
nc_base[j++] = na[i];
|
|
|
|
auto nc = extents_type (nc_base);
|
|
auto wc = strides_type (nc);
|
|
auto c = vector_type (nc.product(), value_type{0});
|
|
|
|
ublas::detail::recursive::mtv(
|
|
size_type(m),
|
|
c.data(), nc.data(), wc.data(),
|
|
a.data(), na.data(), wa.data(),
|
|
b.data());
|
|
|
|
|
|
for(auto i = 0u; i < c.size(); ++i)
|
|
BOOST_CHECK_EQUAL( c[i] , value_type(na[m]) * a[i] );
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_mtm, value, test_types, fixture )
|
|
{
|
|
using namespace boost::numeric;
|
|
using value_type = typename value::first_type;
|
|
using layout_type = typename value::second_type;
|
|
using strides_type = ublas::strides<layout_type>;
|
|
using vector_type = std::vector<value_type>;
|
|
using extents_type = ublas::shape;
|
|
// using extents_type_base = typename extents_type::base_type;
|
|
|
|
|
|
for(auto const& na : extents) {
|
|
|
|
if(na.size() != 2)
|
|
continue;
|
|
|
|
auto a = vector_type (na.product(), value_type{2});
|
|
auto wa = strides_type (na);
|
|
|
|
auto nb = extents_type {na[1],na[0]};
|
|
auto wb = strides_type (nb);
|
|
auto b = vector_type (nb.product(), value_type{1} );
|
|
|
|
auto nc = extents_type {na[0],nb[1]};
|
|
auto wc = strides_type (nc);
|
|
auto c = vector_type (nc.product());
|
|
|
|
|
|
ublas::detail::recursive::mtm(
|
|
c.data(), nc.data(), wc.data(),
|
|
a.data(), na.data(), wa.data(),
|
|
b.data(), nb.data(), wb.data());
|
|
|
|
|
|
for(auto i = 0u; i < c.size(); ++i)
|
|
BOOST_CHECK_EQUAL( c[i] , value_type(na[1]) * a[0] );
|
|
|
|
|
|
}
|
|
}
|
|
|
|
|
|
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttv, value, test_types, fixture )
|
|
{
|
|
using namespace boost::numeric;
|
|
using value_type = typename value::first_type;
|
|
using layout_type = typename value::second_type;
|
|
using strides_type = ublas::strides<layout_type>;
|
|
using vector_type = std::vector<value_type>;
|
|
using extents_type = ublas::shape;
|
|
using extents_type_base = typename extents_type::base_type;
|
|
using size_type = typename extents_type_base::value_type;
|
|
|
|
|
|
for(auto const& na : extents) {
|
|
|
|
auto a = vector_type(na.product(), value_type{2});
|
|
auto wa = strides_type(na);
|
|
for(auto m = 0u; m < na.size(); ++m){
|
|
auto b = vector_type (na[m], value_type{1} );
|
|
auto nb = extents_type {na[m],1};
|
|
auto wb = strides_type (nb);
|
|
|
|
auto nc_base = extents_type_base(std::max(na.size()-1, size_type(2)),1);
|
|
|
|
for(auto i = 0ul, j = 0ul; i < na.size(); ++i)
|
|
if(i != m)
|
|
nc_base[j++] = na[i];
|
|
|
|
auto nc = extents_type (nc_base);
|
|
auto wc = strides_type (nc);
|
|
auto c = vector_type (nc.product(), value_type{0});
|
|
|
|
ublas::ttv(size_type(m+1), na.size(),
|
|
c.data(), nc.data(), wc.data(),
|
|
a.data(), na.data(), wa.data(),
|
|
b.data(), nb.data(), wb.data());
|
|
|
|
|
|
for(auto i = 0u; i < c.size(); ++i)
|
|
BOOST_CHECK_EQUAL( c[i] , value_type(na[m]) * a[i] );
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttm, value, test_types, fixture )
|
|
{
|
|
using namespace boost::numeric;
|
|
using value_type = typename value::first_type;
|
|
using layout_type = typename value::second_type;
|
|
using strides_type = ublas::strides<layout_type>;
|
|
using vector_type = std::vector<value_type>;
|
|
using extents_type = ublas::shape;
|
|
using size_type = typename extents_type::value_type;
|
|
|
|
|
|
for(auto const& na : extents) {
|
|
|
|
auto a = vector_type(na.product(), value_type{2});
|
|
auto wa = strides_type(na);
|
|
for(auto m = 0u; m < na.size(); ++m){
|
|
auto nb = extents_type {na[m], na[m] };
|
|
auto b = vector_type (nb.product(), value_type{1} );
|
|
auto wb = strides_type (nb);
|
|
|
|
|
|
auto nc = na;
|
|
auto wc = strides_type (nc);
|
|
auto c = vector_type (nc.product(), value_type{0});
|
|
|
|
ublas::ttm(size_type(m+1), na.size(),
|
|
c.data(), nc.data(), wc.data(),
|
|
a.data(), na.data(), wa.data(),
|
|
b.data(), nb.data(), wb.data());
|
|
|
|
for(auto i = 0u; i < c.size(); ++i)
|
|
BOOST_CHECK_EQUAL( c[i] , value_type(na[m]) * a[i] );
|
|
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttt_permutation, value, test_types, fixture )
|
|
{
|
|
using namespace boost::numeric;
|
|
using value_type = typename value::first_type;
|
|
using layout_type = typename value::second_type;
|
|
using strides_type = ublas::strides<layout_type>;
|
|
using vector_type = std::vector<value_type>;
|
|
using extents_type = ublas::shape;
|
|
using size_type = typename strides_type::value_type;
|
|
|
|
|
|
auto compute_factorial = [](auto const& p){
|
|
auto f = 1ul;
|
|
for(auto i = 1u; i <= p; ++i)
|
|
f *= i;
|
|
return f;
|
|
};
|
|
|
|
|
|
auto compute_inverse_permutation = [](auto const& pi){
|
|
auto pi_inv = pi;
|
|
for(auto j = 0u; j < pi.size(); ++j)
|
|
pi_inv[pi[j]-1] = j+1;
|
|
return pi_inv;
|
|
};
|
|
|
|
auto permute_extents = [](auto const& pi, auto const& na){
|
|
auto nb = na;
|
|
assert(pi.size() == na.size());
|
|
for(auto j = 0u; j < pi.size(); ++j)
|
|
nb[j] = na[pi[j]-1];
|
|
return nb;
|
|
};
|
|
|
|
|
|
// left-hand and right-hand side have the
|
|
// the same number of elements
|
|
|
|
// computing the inner product with
|
|
// different permutation tuples for
|
|
// right-hand side
|
|
|
|
for(auto const& na : extents) {
|
|
|
|
auto wa = strides_type(na);
|
|
auto a = vector_type(na.product(), value_type{2});
|
|
auto pa = na.size();
|
|
auto pia = std::vector<size_type>(pa);
|
|
std::iota( pia.begin(), pia.end(), 1 );
|
|
|
|
auto pib = pia;
|
|
auto pib_inv = compute_inverse_permutation(pib);
|
|
|
|
auto f = compute_factorial(pa);
|
|
|
|
// for the number of possible permutations
|
|
// only permutation tuple pib is changed.
|
|
for(auto i = 0u; i < f; ++i) {
|
|
|
|
auto nb = permute_extents( pib, na );
|
|
auto wb = strides_type(nb);
|
|
auto b = vector_type(nb.product(), value_type{3});
|
|
auto pb = nb.size();
|
|
|
|
// the number of contractions is changed.
|
|
for( auto q = size_type(0); q <= pa; ++q) {
|
|
|
|
auto r = pa - q;
|
|
auto s = pb - q;
|
|
|
|
auto pc = r+s > 0 ? std::max(r+s,size_type(2)) : size_type(2);
|
|
|
|
auto nc_base = std::vector<size_type>( pc , 1 );
|
|
|
|
for(auto i = 0u; i < r; ++i)
|
|
nc_base[ i ] = na[ pia[i]-1 ];
|
|
|
|
for(auto i = 0u; i < s; ++i)
|
|
nc_base[ r + i ] = nb[ pib_inv[i]-1 ];
|
|
|
|
auto nc = extents_type ( nc_base );
|
|
auto wc = strides_type ( nc );
|
|
auto c = vector_type ( nc.product(), value_type(0) );
|
|
|
|
ublas::ttt(pa,pb,q,
|
|
pia.data(), pib_inv.data(),
|
|
c.data(), nc.data(), wc.data(),
|
|
a.data(), na.data(), wa.data(),
|
|
b.data(), nb.data(), wb.data());
|
|
|
|
|
|
auto acc = value_type(1);
|
|
for(auto i = r; i < pa; ++i)
|
|
acc *= value_type(na[pia[i]-1]);
|
|
|
|
for(auto i = 0ul; i < c.size(); ++i)
|
|
BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
|
|
|
|
}
|
|
|
|
std::next_permutation(pib.begin(), pib.end());
|
|
pib_inv = compute_inverse_permutation(pib);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
|
|
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_ttt, value, test_types, fixture )
|
|
{
|
|
using namespace boost::numeric;
|
|
using value_type = typename value::first_type;
|
|
using layout_type = typename value::second_type;
|
|
using strides_type = ublas::strides<layout_type>;
|
|
using vector_type = std::vector<value_type>;
|
|
using extents_type = ublas::shape;
|
|
using size_type = typename strides_type::value_type;
|
|
|
|
// left-hand and right-hand side have the
|
|
// the same number of elements
|
|
|
|
// computing the inner product with
|
|
// different permutation tuples for
|
|
// right-hand side
|
|
|
|
for(auto const& na : extents) {
|
|
|
|
auto wa = strides_type(na);
|
|
auto a = vector_type(na.product(), value_type{2});
|
|
auto pa = na.size();
|
|
|
|
auto nb = na;
|
|
auto wb = strides_type(nb);
|
|
auto b = vector_type(nb.product(), value_type{3});
|
|
auto pb = nb.size();
|
|
|
|
// std::cout << "na = ";
|
|
// std::copy(na.begin(), na.end(), std::ostream_iterator<size_type>(std::cout, " "));
|
|
// std::cout << std::endl;
|
|
|
|
// std::cout << "nb = ";
|
|
// std::copy(nb.begin(), nb.end(), std::ostream_iterator<size_type>(std::cout, " "));
|
|
// std::cout << std::endl;
|
|
|
|
|
|
// the number of contractions is changed.
|
|
for( auto q = size_type(0); q <= pa; ++q) { // pa
|
|
|
|
auto r = pa - q;
|
|
auto s = pb - q;
|
|
|
|
auto pc = r+s > 0 ? std::max(r+s, size_type(2)) : size_type(2);
|
|
|
|
auto nc_base = std::vector<size_type>( pc , 1 );
|
|
|
|
for(auto i = 0u; i < r; ++i)
|
|
nc_base[ i ] = na[ i ];
|
|
|
|
for(auto i = 0u; i < s; ++i)
|
|
nc_base[ r + i ] = nb[ i ];
|
|
|
|
auto nc = extents_type ( nc_base );
|
|
auto wc = strides_type ( nc );
|
|
auto c = vector_type ( nc.product(), value_type{0} );
|
|
|
|
// std::cout << "nc = ";
|
|
// std::copy(nc.begin(), nc.end(), std::ostream_iterator<size_type>(std::cout, " "));
|
|
// std::cout << std::endl;
|
|
|
|
ublas::ttt(pa,pb,q,
|
|
c.data(), nc.data(), wc.data(),
|
|
a.data(), na.data(), wa.data(),
|
|
b.data(), nb.data(), wb.data());
|
|
|
|
|
|
auto acc = value_type(1);
|
|
for(auto i = r; i < pa; ++i)
|
|
acc *= value_type(na[i]);
|
|
|
|
for(auto i = 0u; i < c.size(); ++i)
|
|
BOOST_CHECK_EQUAL( c[i] , acc * a[0] * b[0] );
|
|
|
|
}
|
|
|
|
}
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_inner, value, test_types, fixture )
|
|
{
|
|
using namespace boost::numeric;
|
|
using value_type = typename value::first_type;
|
|
using layout_type = typename value::second_type;
|
|
using strides_type = ublas::strides<layout_type>;
|
|
using vector_type = std::vector<value_type>;
|
|
|
|
|
|
for(auto const& n : extents) {
|
|
|
|
auto a = vector_type(n.product(), value_type{2});
|
|
auto b = vector_type(n.product(), value_type{3});
|
|
auto w = strides_type(n);
|
|
|
|
auto c = ublas::inner(n.size(), n.data(), a.data(), w.data(), b.data(), w.data(), value_type(0));
|
|
auto cref = std::inner_product(a.begin(), a.end(), b.begin(), value_type(0));
|
|
|
|
|
|
BOOST_CHECK_EQUAL( c , cref );
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_tensor_outer, value, test_types, fixture )
|
|
{
|
|
using namespace boost::numeric;
|
|
using value_type = typename value::first_type;
|
|
using layout_type = typename value::second_type;
|
|
using extents_type = ublas::shape;
|
|
using strides_type = ublas::strides<layout_type>;
|
|
using vector_type = std::vector<value_type>;
|
|
|
|
|
|
for(auto const& na : extents) {
|
|
|
|
auto a = vector_type(na.product(), value_type{2});
|
|
auto wa = strides_type(na);
|
|
|
|
for(auto const& nb : extents) {
|
|
|
|
auto b = vector_type(nb.product(), value_type{3});
|
|
auto wb = strides_type(nb);
|
|
|
|
auto c = vector_type(nb.product()*na.product());
|
|
auto nc = typename extents_type::base_type(na.size()+nb.size());
|
|
|
|
for(auto i = 0u; i < na.size(); ++i)
|
|
nc[i] = na[i];
|
|
for(auto i = 0u; i < nb.size(); ++i)
|
|
nc[i+na.size()] = nb[i];
|
|
|
|
auto wc = strides_type(extents_type(nc));
|
|
|
|
ublas::outer(c.data(), nc.size(), nc.data(), wc.data(),
|
|
a.data(), na.size(), na.data(), wa.data(),
|
|
b.data(), nb.size(), nb.data(), wb.data());
|
|
|
|
for(auto const& cc : c)
|
|
BOOST_CHECK_EQUAL( cc , a[0]*b[0] );
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
BOOST_AUTO_TEST_SUITE_END()
|
|
|