123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491 |
- // 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()
|