From 2d4a091beb9e55664c1475137af7166d524cbc1d Mon Sep 17 00:00:00 2001 From: Mehdi Goli Date: Wed, 14 Dec 2016 15:30:37 +0000 Subject: Adding tensor contraction operation backend for Sycl; adding test for contractionOp sycl backend; adding temporary solution to prevent memory leak in buffer; cleaning up cxx11_tensor_buildins_sycl.h --- .../Eigen/CXX11/src/Tensor/TensorContractionSycl.h | 355 +++++++++++++++++++++ 1 file changed, 355 insertions(+) create mode 100644 unsupported/Eigen/CXX11/src/Tensor/TensorContractionSycl.h (limited to 'unsupported/Eigen/CXX11/src/Tensor/TensorContractionSycl.h') diff --git a/unsupported/Eigen/CXX11/src/Tensor/TensorContractionSycl.h b/unsupported/Eigen/CXX11/src/Tensor/TensorContractionSycl.h new file mode 100644 index 000000000..7e3c73caf --- /dev/null +++ b/unsupported/Eigen/CXX11/src/Tensor/TensorContractionSycl.h @@ -0,0 +1,355 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Mehdi Goli Codeplay Software Ltd. +// Ralph Potter Codeplay Software Ltd. +// Luke Iwanski Codeplay Software Ltd. +// Contact: +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +/***************************************************************** + * TensorSyclConvertToDeviceExpression.h + * + * \brief: + * TensorContractionsycl + * +*****************************************************************/ + +#ifndef EIGEN_CXX11_TENSOR_TENSOR_CONTRACTION_SYCL_H +#define EIGEN_CXX11_TENSOR_TENSOR_CONTRACTION_SYCL_H +namespace Eigen { + +template struct LaunchSyclKernels; +template +struct TensorEvaluator, const Eigen::SyclDevice> : + public TensorContractionEvaluatorBase, const Eigen::SyclDevice> > { + + typedef const Eigen::SyclDevice Device; + + typedef TensorEvaluator, Device> Self; + typedef TensorContractionEvaluatorBase Base; + typedef TensorContractionOp XprType; + typedef typename internal::remove_const::type Scalar; + typedef typename XprType::Index Index; + typedef typename XprType::CoeffReturnType CoeffReturnType; + typedef typename PacketType::type PacketReturnType; + + enum { + Layout = TensorEvaluator::Layout, + }; + + // Most of the code is assuming that both input tensors are ColMajor. If the + // inputs are RowMajor, we will "cheat" by swapping the LHS and RHS: + // If we want to compute A * B = C, where A is LHS and B is RHS, the code + // will pretend B is LHS and A is RHS. + typedef typename internal::conditional< + static_cast(Layout) == static_cast(ColMajor), LeftArgType, RightArgType>::type EvalLeftArgType; + typedef typename internal::conditional< + static_cast(Layout) == static_cast(ColMajor), RightArgType, LeftArgType>::type EvalRightArgType; + + static const int LDims = + internal::array_size::Dimensions>::value; + static const int RDims = + internal::array_size::Dimensions>::value; + static const int ContractDims = internal::array_size::value; + + typedef array left_dim_mapper_t; + typedef array right_dim_mapper_t; + + typedef array contract_t; + typedef array left_nocontract_t; + typedef array right_nocontract_t; + + static const int NumDims = LDims + RDims - 2 * ContractDims; + + typedef DSizes Dimensions; + + // typedefs needed in evalTo + typedef typename internal::remove_const::type LhsScalar; + typedef typename internal::remove_const::type RhsScalar; + + typedef TensorEvaluator LeftEvaluator; + typedef TensorEvaluator RightEvaluator; + + typedef typename LeftEvaluator::Dimensions LeftDimensions; + typedef typename RightEvaluator::Dimensions RightDimensions; + + EIGEN_DEVICE_FUNC TensorEvaluator(const XprType& op, const Device& device) : + Base(op, device) {} + + // We need to redefine this method to make nvcc happy + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool evalSubExprsIfNeeded(Scalar* data) { + this->m_leftImpl.evalSubExprsIfNeeded(NULL); + this->m_rightImpl.evalSubExprsIfNeeded(NULL); + if (data) { + evalTo(data); + return false; + } else { + this->m_result = static_cast(this->m_device.allocate(this->dimensions().TotalSize() * sizeof(Scalar))); + evalTo(this->m_result); + return true; + } + } + const Eigen::SyclDevice& device() const {return this->m_device;} + void evalTo(Scalar* buffer) const { + // Here is the result + if (this->m_lhs_inner_dim_contiguous) { + if (this->m_rhs_inner_dim_contiguous) { + if (this->m_rhs_inner_dim_reordered) { + evalTyped(buffer); + } + else { + evalTyped(buffer); + } + } + else { + if (this->m_rhs_inner_dim_reordered) { + evalTyped(buffer); + } + else { + evalTyped(buffer); + } + } + } + else { + if (this->m_rhs_inner_dim_contiguous) { + if (this->m_rhs_inner_dim_reordered) { + evalTyped(buffer); + } + else { + evalTyped(buffer); + } + } + else { + if (this->m_rhs_inner_dim_reordered) { + evalTyped(buffer); + } + else { + evalTyped(buffer); + } + } + } + } + + template + void evalTyped(Scalar* buffer) const { + // columns in left side, rows in right side + const Index k = this->m_k_size; + EIGEN_UNUSED_VARIABLE(k) + // rows in left side + const Index m = this->m_i_size; + // columns in right side + const Index n = this->m_j_size; + + // zero out the result buffer (which must be of size at least m * n * sizeof(Scalar) + this->m_device.memset(buffer, 0, m * n * sizeof(Scalar)); + LaunchSyclKernels::Run(*this, buffer, m, n, k, + this->m_k_strides, this->m_left_contracting_strides, this->m_right_contracting_strides, + this->m_i_strides, this->m_j_strides, this->m_left_nocontract_strides, this->m_right_nocontract_strides); + } + // required by sycl to construct the expr on the device. Returns original left_impl + const TensorEvaluator& left_impl() const { + return choose(Cond(Layout) == static_cast(ColMajor)>(), this->m_leftImpl, this->m_rightImpl); + } + // required by sycl to construct the expr on the device. Returns original right_impl + const TensorEvaluator& right_impl() const { + return choose(Cond(Layout) == static_cast(ColMajor)>(), this->m_rightImpl, this->m_leftImpl); + } + // required by sycl to construct the expr on the device + const Indices& indices() const {return this->m_expr_indices;} +}; + +/// Dummy container on the device. This is used to avoid calling the constructor of TensorEvaluator for TensorContractionOp. This makes the code much faster. +template struct TensorEvaluatorContainer; +template +struct TensorEvaluatorContainer>{ + typedef Eigen::DefaultDevice Device; + typedef TensorContractionOp XprType; + typedef typename internal::remove_const::type Scalar; + typedef typename XprType::Index Index; + typedef typename XprType::CoeffReturnType CoeffReturnType; + typedef typename PacketType::type PacketReturnType; + enum { + Layout = TensorEvaluator::Layout, + }; + + typedef typename internal::conditional(Layout) == static_cast(ColMajor), LeftArgType, RightArgType>::type EvalLeftArgType; + typedef typename internal::conditional(Layout) == static_cast(ColMajor), RightArgType, LeftArgType>::type EvalRightArgType; + typedef TensorEvaluator LeftEvaluator; + typedef TensorEvaluator RightEvaluator; + + TensorEvaluatorContainer(const XprType& op, const Eigen::DefaultDevice& device) + : m_leftImpl(choose(Cond(Layout) == static_cast(ColMajor)>(), + op.lhsExpression(), op.rhsExpression()), device), + m_rightImpl(choose(Cond(Layout) == static_cast(ColMajor)>(), + op.rhsExpression(), op.lhsExpression()), device){} +LeftEvaluator m_leftImpl; +RightEvaluator m_rightImpl; +}; + +#define TileSizeDimM 32 // Tile size for dimension M +#define TileSizeDimN 32 // Tile size for dimension N +#define TileSizeDimK 16 // Tile size for dimension K +#define WorkLoadPerThreadM 4 // Work load per thread in dimension M +#define WorkLoadPerThreadN 4 // work load per thread in dimension N +#define LocalThreadSizeM (TileSizeDimM/WorkLoadPerThreadM) // Local thread size for the first dimension (M here) +#define LocalThreadSizeN (TileSizeDimN/WorkLoadPerThreadN) // Local thread size for the second dimension (N here) +#define LoadPerThreadLhs ((TileSizeDimK*WorkLoadPerThreadM*WorkLoadPerThreadN)/(TileSizeDimN)) // workload per thread for Lhs expression +#define LoadPerThreadRhs ((TileSizeDimK*WorkLoadPerThreadM*WorkLoadPerThreadN)/(TileSizeDimM)) // workload per thread for Rhs expression +#define RoundUp(x,y) ((((x) + (y) - 1) / (y))*(y)) // RoundUp function to make sure that the global threadId is dividabe by local threadId + +template struct KernelNameConstructor; +template struct LaunchSyclKernels { +template< typename Self, typename Output, typename Index, typename ContractT, typename LeftNocontractT, typename RightNocontractT> + static void Run(const Self& self, Output* buffer, Index M, Index N, Index K, + ContractT m_k_strides, ContractT m_left_contracting_strides, ContractT m_right_contracting_strides, + LeftNocontractT m_i_strides, RightNocontractT m_j_strides, LeftNocontractT m_left_nocontract_strides, RightNocontractT m_right_nocontract_strides){ + // create a tuple of accessors from Evaluator + typedef typename Eigen::TensorSycl::internal::createPlaceHolderExpression::Type PlaceHolderExpr; + typedef KernelNameConstructor KernelName; + auto functors = Eigen::TensorSycl::internal::extractFunctors(self); + Index roundUpK = RoundUp(K, TileSizeDimK); + Index roundUpM = RoundUp(M, TileSizeDimM); + Index roundUpN = RoundUp(N, TileSizeDimN); + self.device().sycl_queue().submit([&](cl::sycl::handler &cgh) { + auto tuple_of_accessors = Eigen::TensorSycl::internal::createTupleOfAccessors(cgh, self); + // Local memory for elements of Lhs + cl::sycl::accessor localLhs(cl::sycl::range<1>(2* TileSizeDimM * TileSizeDimK), cgh); + // Local memory for elements of Rhs + cl::sycl::accessor localRhs(cl::sycl::range<1>(2* TileSizeDimK * TileSizeDimN), cgh); + //Output memory + auto out_privateRes= self.device(). template get_sycl_accessor(cgh, buffer); + // sycl parallel for + cgh.parallel_for( cl::sycl::nd_range<2>(cl::sycl::range<2>(roundUpM/WorkLoadPerThreadM, roundUpN/WorkLoadPerThreadN), cl::sycl::range<2>(LocalThreadSizeM, LocalThreadSizeN)), [=](cl::sycl::nd_item<2> itemID) { + typedef typename Eigen::TensorSycl::internal::ConvertToDeviceExpression::Type DevExpr; + auto device_expr =Eigen::TensorSycl::internal::createDeviceExpression(functors, tuple_of_accessors); + auto device_evaluator = TensorEvaluatorContainer(device_expr.expr, Eigen::DefaultDevice()); + typedef TensorEvaluatorContainer DevEvaluator; + typedef internal::TensorContractionInputMapper LhsMapper; + + typedef internal::TensorContractionInputMapper RhsMapper; + // initialize data mappers must happen inside the kernel for device eval + LhsMapper lhs(device_evaluator.m_leftImpl, m_left_nocontract_strides, m_i_strides, m_left_contracting_strides, m_k_strides); + RhsMapper rhs(device_evaluator.m_rightImpl, m_right_nocontract_strides, m_j_strides, m_right_contracting_strides, m_k_strides); + auto out_ptr = ConvertToActualTypeSycl(Output, out_privateRes); + // Matmul Kernel + // Thread identifiers + const int mLocalThreadId = itemID.get_local(0); // Local ID row + const int nLocalThreadId = itemID.get_local(1); // Local ID col + const int mGroupId = itemID.get_group(0); // Work-group ID row + const int nGroupId = itemID.get_group(1); // Work-group ID localCol + const int linearLocalThreadId = nLocalThreadId*LocalThreadSizeM + mLocalThreadId; // linear local thread ID + // Allocate register space + float privateLhs; + float privateRhs[WorkLoadPerThreadN]; + float privateRes[WorkLoadPerThreadM][WorkLoadPerThreadN]; + // Initialise the privateResumulation registers + for (int wLPTM=0; wLPTM(0); + } + // Tile Rhs + for (int lPTR=0; lPTR(0); + + } + // Loop over all tiles + const int numTiles = roundUpK/TileSizeDimK; + int firstHalf=0; + do { + // Synchronise + itemID.barrier(cl::sycl::access::fence_space::local_space); + // Load the next tile of Lhs and Rhs into local memory + int nextHalf = firstHalf + 1; + if (nextHalf < numTiles) { + // Tile A + for (int lPTL=0; lPTL(0); + } + // Tile B + for (int lPTR=0; lPTR(0); + } + } + // Loop over the values of a single tile + for (int k=0; k