[PATCH 2/3] boost::ublas: gcc support for optimal matrix multiplication
Imre Palik <[email protected]>
| Newsgroups | gmane.comp.lib.boost.ublas |
|---|---|
| Message-ID | <[email protected]> |
This patch contains an optimised, vectorised kernel, using gcc's SIMD vectors. This reaches matrix multiplication speeds comparable to hand-crafted assembly kernels on x86 Haswell microarchitecture. The patch also contains optimised kerenl parameters for float, double, and long double. Signed-off-by: Imre Palik <[email protected]> Cc: Michael Lehn <[email protected]> --- include/boost/numeric/ublas/detail/block_sizes.hpp | 78 +++++++ include/boost/numeric/ublas/detail/gemm.hpp | 223 +++++++++++---------- include/boost/numeric/ublas/detail/vector.hpp | 27 +++ include/boost/numeric/ublas/matrix_expression.hpp | 1 + include/boost/numeric/ublas/operation.hpp | 1 + 5 files changed, 223 insertions(+), 107 deletions(-) create mode 100644 include/boost/numeric/ublas/detail/block_sizes.hpp create mode 100644 include/boost/numeric/ublas/detail/vector.hpp diff --git a/include/boost/numeric/ublas/detail/block_sizes.hpp b/include/boost/numeric/ublas/detail/block_sizes.hpp new file mode 100644 index 0000000..a146f10 --- /dev/null +++ b/include/boost/numeric/ublas/detail/block_sizes.hpp @@ -0,0 +1,78 @@ +// +// Copyright (c) 2016 +// Michael Lehn, Imre Palik +// +// 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) + +#ifndef _BOOST_UBLAS_BLOCK_SIZES_ +#define _BOOST_UBLAS_BLOCK_SIZES_ + +#include <boost/numeric/ublas/detail/vector.hpp> + +namespace boost { namespace numeric { namespace ublas { namespace detail { + + template <typename T> + struct prod_block_size { + static const unsigned vector_length = + _BOOST_UBLAS_VECTOR_SIZE > sizeof(T)? _BOOST_UBLAS_VECTOR_SIZE/sizeof(T) : 1; // Number of elements in a vector register + static const unsigned mc = 256; + static const unsigned kc = 512; // stripe length + static const unsigned nc = (4096/(3 * vector_length)) * (3 * vector_length); + static const unsigned mr = 4; // stripe width for lhs + static const unsigned nr = 3 * vector_length; // stripe width for rhs + static const unsigned align = 64; // align temporary arrays to this boundary + static const unsigned limit = 27; // Use gemm from this size + BOOST_STATIC_ASSERT_MSG(mc>0 && kc>0 && nc>0 && mr>0 && nr>0, "Invalid block size."); + BOOST_STATIC_ASSERT_MSG(mc % mr == 0, "MC must be a multiple of MR."); + BOOST_STATIC_ASSERT_MSG(nc % nr == 0, "NC must be a multiple of NR."); + }; + + template <> + struct prod_block_size<float> { + static const unsigned mc = 256; + static const unsigned kc = 512; // stripe length + static const unsigned nc = 4096; + static const unsigned mr = 4; // stripe width for lhs + static const unsigned nr = 16; // stripe width for rhs + static const unsigned align = 64; // align temporary arrays to this boundary + static const unsigned limit = 27; // Use gemm from this size + static const unsigned vector_length = _BOOST_UBLAS_VECTOR_SIZE/sizeof(float); // Number of elements in a vector register + BOOST_STATIC_ASSERT_MSG(mc>0 && kc>0 && nc>0 && mr>0 && nr>0, "Invalid block size."); + BOOST_STATIC_ASSERT_MSG(mc % mr == 0, "MC must be a multiple of MR."); + BOOST_STATIC_ASSERT_MSG(nc % nr == 0, "NC must be a multiple of NR."); + }; + + template <> + struct prod_block_size<long double> { + static const unsigned mc = 256; + static const unsigned kc = 512; // stripe length + static const unsigned nc = 4096; + static const unsigned mr = 1; // stripe width for lhs + static const unsigned nr = 4; // stripe width for rhs + static const unsigned align = 64; // align temporary arrays to this boundary + static const unsigned limit = 71; // Use gemm from this size + static const unsigned vector_length = 1; // Number of elements in a vector register + BOOST_STATIC_ASSERT_MSG(mc>0 && kc>0 && nc>0 && mr>0 && nr>0, "Invalid block size."); + BOOST_STATIC_ASSERT_MSG(mc % mr == 0, "MC must be a multiple of MR."); + BOOST_STATIC_ASSERT_MSG(nc % nr == 0, "NC must be a multiple of NR."); + }; + + template<typename T> + struct is_blocksize { + struct fallback { static const int nr = 0; }; + struct derived : T, fallback {}; + template<int C1> + struct nonr { + static const bool value = false; + typedef false_type type; + }; + + template<typename C> static char (&f(typename nonr<C::nr>::type*))[1]; + template<typename C> static char (&f(...))[2]; + + static bool const value = sizeof(f<derived>(0)) == 2; + }; +}}}} +#endif diff --git a/include/boost/numeric/ublas/detail/gemm.hpp b/include/boost/numeric/ublas/detail/gemm.hpp index f66e80f..7fd0ab9 100644 --- a/include/boost/numeric/ublas/detail/gemm.hpp +++ b/include/boost/numeric/ublas/detail/gemm.hpp @@ -10,80 +10,26 @@ #define _BOOST_UBLAS_GEMM_ #include <boost/type_traits/common_type.hpp> +#include <boost/type_traits/aligned_storage.hpp> +#include <boost/type_traits/is_arithmetic.hpp> #include <boost/align/aligned_alloc.hpp> #include <boost/align/assume_aligned.hpp> #include <boost/static_assert.hpp> #include <boost/numeric/ublas/matrix_proxy.hpp> +#include <boost/numeric/ublas/detail/vector.hpp> +#include <boost/predef/compiler.h> namespace boost { namespace numeric { namespace ublas { namespace detail { - template <typename T> - struct prod_block_size { - static const unsigned mc = 256; - static const unsigned kc = 512; // stripe length - static const unsigned nc = 4092; - static const unsigned mr = 4; // stripe width for lhs - static const unsigned nr = 12; // stripe width for rhs - static const unsigned align = 64; // align temporary arrays to this boundary - static const unsigned limit = 26; // Use gemm from this size - BOOST_STATIC_ASSERT_MSG(mc>0 && kc>0 && nc>0 && mr>0 && nr>0, "Invalid block size."); - BOOST_STATIC_ASSERT_MSG(mc % mr == 0, "MC must be a multiple of MR."); - BOOST_STATIC_ASSERT_MSG(nc % nr == 0, "NC must be a multiple of NR."); - }; - - template <> - struct prod_block_size<float> { - static const unsigned mc = 256; - static const unsigned kc = 512; // stripe length - static const unsigned nc = 4080; - static const unsigned mr = 4; // stripe width for lhs - static const unsigned nr = 24; // stripe width for rhs - static const unsigned align = 64; // align temporary arrays to this boundary - static const unsigned limit = 26; // Use gemm from this size - BOOST_STATIC_ASSERT_MSG(mc>0 && kc>0 && nc>0 && mr>0 && nr>0, "Invalid block size."); - BOOST_STATIC_ASSERT_MSG(mc % mr == 0, "MC must be a multiple of MR."); - BOOST_STATIC_ASSERT_MSG(nc % nr == 0, "NC must be a multiple of NR."); - }; - - template <> - struct prod_block_size<long double> { - static const unsigned mc = 256; - static const unsigned kc = 512; // stripe length - static const unsigned nc = 4096; - static const unsigned mr = 2; // stripe width for lhs - static const unsigned nr = 2; // stripe width for rhs - static const unsigned align = 64; // align temporary arrays to this boundary - static const unsigned limit = 26; // Use gemm from this size - BOOST_STATIC_ASSERT_MSG(mc>0 && kc>0 && nc>0 && mr>0 && nr>0, "Invalid block size."); - BOOST_STATIC_ASSERT_MSG(mc % mr == 0, "MC must be a multiple of MR."); - BOOST_STATIC_ASSERT_MSG(nc % nr == 0, "NC must be a multiple of NR."); - }; - - template<typename T> - struct is_blocksize { - struct fallback { static const int nr = 0; }; - struct derived : T, fallback {}; - template<int C1> - struct nonr { - static const bool value = false; - typedef false_type type; - }; - - template<typename C> static char (&f(typename nonr<C::nr>::type*))[1]; - template<typename C> static char (&f(...))[2]; - - static bool const value = sizeof(f<derived>(0)) == 2; - }; - template <typename E> void gescal(const typename E::value_type &alpha, matrix_expression<E> &X) { typedef typename E::size_type size_type; - for (size_type i=0; i<X().size1(); ++i) { + for (size_type i=0; i<X().size1(); ++i) { for (size_type j=0; j<X().size2(); ++j) { - X()(i,j) *= alpha; + X()(i,j) *= alpha; } } } @@ -114,9 +60,73 @@ namespace boost { namespace numeric { namespace ublas { namespace detail { } } - //-- Micro Kernel -------------------------------------------------------------- - template <typename Index, typename T, typename TC, typename BlockSize> + //-- Micro Kernel ---------------------------------------------------------- +#if defined(BOOST_COMP_GNUC_DETECTION) \ + && BOOST_COMP_GNUC_DETECTION >= BOOST_VERSION_NUMBER(4,8,0) + template <typename Index, typename T, typename TC, + typename BlockSize> + typename enable_if_c<is_arithmetic<T>::value + && (1 < BlockSize::vector_length), + void>::type + ugemm(Index kc, TC alpha, const T *A, const T *B, + TC beta, TC *C, Index incRowC, Index incColC) + { + BOOST_ALIGN_ASSUME_ALIGNED (A, BlockSize::align); + BOOST_ALIGN_ASSUME_ALIGNED (B, BlockSize::align); + static const unsigned vector_length = BlockSize::vector_length; + static const Index MR = BlockSize::mr; + static const Index NR = BlockSize::nr/vector_length; + + typedef T vx __attribute__((vector_size (_BOOST_UBLAS_VECTOR_SIZE))); + + vx P[MR*NR] = {}; + + const vx *b = (const vx *)B; + for (Index l=0; l<kc; ++l) { + for (Index i=0; i<MR; ++i) { + for (Index j=0; j<NR; ++j) { + P[i*NR+j] += A[i]*b[j]; + } + } + A += MR; + b += NR; + } + + if (alpha!=TC(1)) { + for (Index i=0; i<MR; ++i) { + for (Index j=0; j<NR; ++j) { + P[i*NR+j] *= alpha; + } + } + } + + const T *p = (const T *) P; + if (beta == TC(0)) { + for (Index i=0; i<MR; ++i) { + for (Index j=0; j<NR * vector_length; ++j) { + C[i*incRowC+j*incColC] = p[i*NR * vector_length +j]; + } + } + } else { + for (Index i=0; i<MR; ++i) { + for (Index j=0; j<NR * vector_length; ++j) { + C[i*incRowC+j*incColC] *= beta; + C[i*incRowC+j*incColC] += p[i*NR * vector_length+j]; + } + } + } + } + + template <typename Index, typename T, typename TC, + typename BlockSize> + typename enable_if_c<!is_arithmetic<T>::value + || (BlockSize::vector_length <= 1), + void>::type +#else + template <typename Index, typename T, typename TC, + typename BlockSize> void +#endif ugemm(Index kc, TC alpha, const T *A, const T *B, TC beta, TC *C, Index incRowC, Index incColC) { @@ -124,47 +134,47 @@ namespace boost { namespace numeric { namespace ublas { namespace detail { BOOST_ALIGN_ASSUME_ALIGNED(B, BlockSize::align); static const Index MR = BlockSize::mr; static const Index NR = BlockSize::nr; - typename boost::aligned_storage<sizeof(T[MR*NR]),BlockSize::align>::type Pa; - T *P = reinterpret_cast<T*>(Pa.address()); - for (unsigned c = 0; c < MR * NR; c++) - P[c] = 0; + typename aligned_storage<sizeof(T[MR*NR]),BlockSize::align>::type Pa; + T *P = reinterpret_cast<T*>(Pa.address()); + for (unsigned c = 0; c < MR * NR; c++) + P[c] = 0; for (Index l=0; l<kc; ++l) { - for (Index i=0; i<MR; ++i) { + for (Index i=0; i<MR; ++i) { for (Index j=0; j<NR; ++j) { P[i* NR+j] += A[i]*B[j]; } } - A += MR; - B += NR; + A += MR; + B += NR; + } + + if (alpha!=TC(1)) { + for (Index i=0; i<MR; ++i) { + for (Index j=0; j<NR; ++j) { + P[i*NR+j] *= alpha; + } + } } - if (alpha!=TC(1)) { - for (Index i=0; i<MR; ++i) { - for (Index j=0; j<NR; ++j) { - P[i*NR+j] *= alpha; - } - } - } - - if (beta == TC(0)) { - for (Index i=0; i<MR; ++i) { - for (Index j=0; j<NR; ++j) { + if (beta == TC(0)) { + for (Index i=0; i<MR; ++i) { + for (Index j=0; j<NR; ++j) { C[i*incRowC+j*incColC] = P[i*NR+j]; - } - } - } else { - for (Index i=0; i<MR; ++i) { - for (Index j=0; j<NR; ++j) { + } + } + } else { + for (Index i=0; i<MR; ++i) { + for (Index j=0; j<NR; ++j) { C[i*incRowC+j*incColC] *= beta; C[i*incRowC+j*incColC] += P[i*NR+j]; - } - } - } + } + } + } } - //-- Macro Kernel -------------------------------------------------------------- - template <typename Index, typename T, typename TC, typename BlockSize> + //-- Macro Kernel ---------------------------------------------------------- + template <typename Index, typename T, typename TC, typename BlockSize> void mgemm(Index mc, Index nc, Index kc, TC alpha, const T *A, const T *B, TC beta, @@ -181,9 +191,6 @@ namespace boost { namespace numeric { namespace ublas { namespace detail { // #pragma omp parallel for // #endif for (Index j=0; j<np; ++j) { - // __builtin_prefetch(B + j * kc * NR, 0); - // __builtin_prefetch(B + j * kc * NR + 8, 0); - // __builtin_prefetch(B + j * kc * NR + 16, 0); const Index nr = (j!=np-1 || nr_==0) ? NR : nr_; T C_[MR*NR]; @@ -191,7 +198,7 @@ namespace boost { namespace numeric { namespace ublas { namespace detail { const Index mr = (i!=mp-1 || mr_==0) ? MR : mr_; if (mr==MR && nr==NR) { - ugemm<Index, T, TC, BlockSize>(kc, alpha, + ugemm<Index, T, TC, BlockSize>(kc, alpha, &A[i*kc*MR], &B[j*kc*NR], beta, &C[i*MR*incRowC+j*NR*incColC], @@ -214,11 +221,12 @@ namespace boost { namespace numeric { namespace ublas { namespace detail { } //-- Packing blocks ------------------------------------------------------------ - template <typename E, typename T, typename BlockSize> + template <typename E, typename T, typename BlockSize> void pack_A(const matrix_expression<E> &A, T *p) { typedef typename E::size_type size_type; + BOOST_ALIGN_ASSUME_ALIGNED(p, BlockSize::align); const size_type mc = A ().size1(); const size_type kc = A ().size2(); @@ -236,11 +244,12 @@ namespace boost { namespace numeric { namespace ublas { namespace detail { } } - template <typename E, typename T, typename BlockSize> + template <typename E, typename T, typename BlockSize> void pack_B(const matrix_expression<E> &B, T *p) { typedef typename E::size_type size_type; + BOOST_ALIGN_ASSUME_ALIGNED(p, BlockSize::align); const size_type kc = B ().size1(); const size_type nc = B ().size2(); @@ -262,9 +271,9 @@ namespace boost { namespace numeric { namespace ublas { namespace detail { template <typename E1, typename E2, typename E3, typename BlockSize> void gemm(typename E3::value_type alpha, const matrix_expression<E1> &e1, - const matrix_expression<E2> &e2, + const matrix_expression<E2> &e2, typename E3::value_type beta, matrix_expression<E3> &e3, - BlockSize) + BlockSize) { typedef typename E3::size_type size_type; typedef typename E1::value_type value_type1; @@ -276,7 +285,7 @@ namespace boost { namespace numeric { namespace ublas { namespace detail { static const size_type NC = BlockSize::nc; const size_type m = BOOST_UBLAS_SAME (e3 ().size1 (), e1 ().size1 ()); - const size_type n = BOOST_UBLAS_SAME (e3 ().size2 (), e2 ().size2 ()); + const size_type n = BOOST_UBLAS_SAME (e3 ().size2 (), e2 ().size2 ()); const size_type k = BOOST_UBLAS_SAME (e1 ().size2 (), e2 ().size1 ()); static const size_type KC = BlockSize::kc; @@ -291,11 +300,11 @@ namespace boost { namespace numeric { namespace ublas { namespace detail { const size_type incRowC = &e3()(1,0) - &e3()(0,0); const size_type incColC = &e3()(0,1) - &e3()(0,0); value_type_i *A = - static_cast<value_type_i *>(boost::alignment::aligned_alloc(BlockSize::align, - sizeof(value_type_i) * MC*KC)); + static_cast<value_type_i *>(::boost::alignment::aligned_alloc(BlockSize::align, + sizeof(value_type_i) * MC*KC)); value_type_i *B = - static_cast<value_type_i *>(boost::alignment::aligned_alloc(BlockSize::align, - sizeof(value_type_i) * KC*NC)); + static_cast<value_type_i *>(::boost::alignment::aligned_alloc(BlockSize::align, + sizeof(value_type_i) * KC*NC)); if (alpha==value_type3(0) || k==0) { gescal(beta, e3); @@ -319,13 +328,13 @@ namespace boost { namespace numeric { namespace ublas { namespace detail { pack_A<matrix_range<const E1>, value_type_i, BlockSize>(As, A); mgemm<size_type, value_type_i, value_type3, BlockSize>(mc, nc, kc, alpha, A, B, beta_, - &C_[i*MC*incRowC+j*NC*incColC], - incRowC, incColC); + &C_[i*MC*incRowC+j*NC*incColC], + incRowC, incColC); } } } - boost::alignment::aligned_free(A); - boost::alignment::aligned_free(B); + ::boost::alignment::aligned_free(A); + ::boost::alignment::aligned_free(B); } }}}} #endif diff --git a/include/boost/numeric/ublas/detail/vector.hpp b/include/boost/numeric/ublas/detail/vector.hpp new file mode 100644 index 0000000..1681eda --- /dev/null +++ b/include/boost/numeric/ublas/detail/vector.hpp @@ -0,0 +1,27 @@ +// +// Copyright (c) 2016 +// Imre Palik +// +// 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) + +#ifndef _BOOST_UBLAS_DETAIL_VECTOR_ +#define _BOOST_UBLAS_DETAIL_VECTOR_ +#include <boost/predef/hardware/simd.h> + +#if defined(BOOST_HW_SIMD) && (BOOST_HW_SIMD == BOOST_HW_SIMD_X86 || BOOST_HW_SIMD == BOOST_HW_SIMD_X86_AMD) +#if BOOST_HW_SIMD >= BOOST_HW_SIMD_X86_AVX_VERSION +#define _BOOST_UBLAS_VECTOR_SIZE 32 +#elif BOOST_HW_SIMD >= BOOST_HW_SIMD_X86_SSE_VERSION +#define _BOOST_UBLAS_VECTOR_SIZE 16 +#else +#define _BOOST_UBLAS_VECTOR_SIZE 8 +#endif +#endif + +#ifndef _BOOST_UBLAS_VECTOR_SIZE +#define _BOOST_UBLAS_VECTOR_SIZE 1 +#endif + +#endif diff --git a/include/boost/numeric/ublas/matrix_expression.hpp b/include/boost/numeric/ublas/matrix_expression.hpp index 78daf2b..c4e0ac8 100644 --- a/include/boost/numeric/ublas/matrix_expression.hpp +++ b/include/boost/numeric/ublas/matrix_expression.hpp @@ -15,6 +15,7 @@ #include <boost/numeric/ublas/vector_expression.hpp> #include <boost/numeric/ublas/detail/gemm.hpp> +#include <boost/numeric/ublas/detail/block_sizes.hpp> // Expression templates based on ideas of Todd Veldhuizen and Geoffrey Furnish // Iterators based on ideas of Jeremy Siek diff --git a/include/boost/numeric/ublas/operation.hpp b/include/boost/numeric/ublas/operation.hpp index 393a9bd..8fd3b63 100644 --- a/include/boost/numeric/ublas/operation.hpp +++ b/include/boost/numeric/ublas/operation.hpp @@ -15,6 +15,7 @@ #include <boost/numeric/ublas/matrix_proxy.hpp> #include <boost/numeric/ublas/detail/gemm.hpp> +#include <boost/numeric/ublas/detail/block_sizes.hpp> /** \file operation.hpp * \brief This file contains some specialized products. */ -- 1.9.1