Added support for RowMajor inputs to the contraction code.

This commit is contained in:
Benoit Steiner 2015-01-14 12:36:57 -08:00
parent 0a0ab6dd15
commit 71676eaddd
3 changed files with 219 additions and 93 deletions

View File

@ -320,6 +320,8 @@ class TensorContractionInputMapper
}; };
template<typename Scalar, typename Index, int side, template<typename Scalar, typename Index, int side,
typename Tensor, typename Tensor,
typename nocontract_t, typename contract_t, typename nocontract_t, typename contract_t,
@ -362,6 +364,14 @@ class TensorContractionInputMapper<Scalar, Index, side, Tensor, nocontract_t, co
}; };
template <size_t n> struct max_n_1 {
static const size_t size = n;
};
template <> struct max_n_1<0> {
static const size_t size = 1;
};
template<typename Dimensions, typename LhsXprType, typename RhsXprType> template<typename Dimensions, typename LhsXprType, typename RhsXprType>
struct traits<TensorContractionOp<Dimensions, LhsXprType, RhsXprType> > struct traits<TensorContractionOp<Dimensions, LhsXprType, RhsXprType> >
{ {
@ -378,6 +388,10 @@ struct traits<TensorContractionOp<Dimensions, LhsXprType, RhsXprType> >
typedef typename remove_reference<LhsNested>::type _LhsNested; typedef typename remove_reference<LhsNested>::type _LhsNested;
typedef typename remove_reference<RhsNested>::type _RhsNested; typedef typename remove_reference<RhsNested>::type _RhsNested;
// From NumDims below.
static const int NumDimensions = max_n_1<traits<RhsXprType>::NumDimensions + traits<RhsXprType>::NumDimensions - 2 * array_size<Dimensions>::value>::size;
static const int Layout = traits<LhsXprType>::Layout;
enum { enum {
Flags = 0, Flags = 0,
}; };
@ -401,19 +415,19 @@ struct traits<TensorEvaluator<const TensorContractionOp<Indices_, LeftArgType_,
typedef LeftArgType_ LeftArgType; typedef LeftArgType_ LeftArgType;
typedef RightArgType_ RightArgType; typedef RightArgType_ RightArgType;
typedef Device_ Device; typedef Device_ Device;
// From NumDims below.
static const int NumDimensions = max_n_1<traits<LeftArgType_>::NumDimensions + traits<RightArgType_>::NumDimensions - 2 * array_size<Indices_>::value>::size;
}; };
} // end namespace internal } // end namespace internal
template<typename Indices, typename LhsXprType, typename RhsXprType> template<typename Indices, typename LhsXprType, typename RhsXprType>
class TensorContractionOp : public TensorBase<TensorContractionOp<Indices, LhsXprType, RhsXprType> > class TensorContractionOp : public TensorBase<TensorContractionOp<Indices, LhsXprType, RhsXprType> >
{ {
public: public:
typedef typename Eigen::internal::traits<TensorContractionOp>::Scalar Scalar; typedef typename Eigen::internal::traits<TensorContractionOp>::Scalar Scalar;
typedef typename Eigen::internal::traits<TensorContractionOp>::Packet Packet; typedef typename Eigen::internal::traits<TensorContractionOp>::Packet Packet;
typedef typename Eigen::NumTraits<Scalar>::Real RealScalar;
typedef typename internal::promote_storage_type<typename LhsXprType::CoeffReturnType, typedef typename internal::promote_storage_type<typename LhsXprType::CoeffReturnType,
typename RhsXprType::CoeffReturnType>::ret CoeffReturnType; typename RhsXprType::CoeffReturnType>::ret CoeffReturnType;
typedef typename internal::promote_storage_type<typename LhsXprType::PacketReturnType, typedef typename internal::promote_storage_type<typename LhsXprType::PacketReturnType,
@ -422,20 +436,21 @@ class TensorContractionOp : public TensorBase<TensorContractionOp<Indices, LhsXp
typedef typename Eigen::internal::traits<TensorContractionOp>::StorageKind StorageKind; typedef typename Eigen::internal::traits<TensorContractionOp>::StorageKind StorageKind;
typedef typename Eigen::internal::traits<TensorContractionOp>::Index Index; typedef typename Eigen::internal::traits<TensorContractionOp>::Index Index;
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorContractionOp(const LhsXprType& lhs, const RhsXprType& rhs, const Indices& dims) EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorContractionOp(
const LhsXprType& lhs, const RhsXprType& rhs, const Indices& dims)
: m_lhs_xpr(lhs), m_rhs_xpr(rhs), m_indices(dims) {} : m_lhs_xpr(lhs), m_rhs_xpr(rhs), m_indices(dims) {}
EIGEN_DEVICE_FUNC EIGEN_DEVICE_FUNC
const Indices& indices() const { return m_indices; } const Indices& indices() const { return m_indices; }
/** \returns the nested expressions */ /** \returns the nested expressions */
EIGEN_DEVICE_FUNC EIGEN_DEVICE_FUNC
const typename internal::remove_all<typename LhsXprType::Nested>::type& const typename internal::remove_all<typename LhsXprType::Nested>::type&
lhsExpression() const { return m_lhs_xpr; } lhsExpression() const { return m_lhs_xpr; }
EIGEN_DEVICE_FUNC EIGEN_DEVICE_FUNC
const typename internal::remove_all<typename RhsXprType::Nested>::type& const typename internal::remove_all<typename RhsXprType::Nested>::type&
rhsExpression() const { return m_rhs_xpr; } rhsExpression() const { return m_rhs_xpr; }
protected: protected:
typename LhsXprType::Nested m_lhs_xpr; typename LhsXprType::Nested m_lhs_xpr;
@ -444,12 +459,17 @@ class TensorContractionOp : public TensorBase<TensorContractionOp<Indices, LhsXp
}; };
template <size_t n> struct max_n_1 { template<bool cond> struct Cond {};
static const size_t size = n;
}; template<typename T1, typename T2> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
template <> struct max_n_1<0> { const T1& choose(Cond<true>, const T1& first, const T2&) {
static const size_t size = 1; return first;
}; }
template<typename T1, typename T2> EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE
const T2& choose(Cond<false>, const T1&, const T2& second) {
return second;
}
template<typename Derived> template<typename Derived>
@ -467,37 +487,94 @@ struct TensorContractionEvaluatorBase
typedef typename XprType::CoeffReturnType CoeffReturnType; typedef typename XprType::CoeffReturnType CoeffReturnType;
typedef typename XprType::PacketReturnType PacketReturnType; typedef typename XprType::PacketReturnType PacketReturnType;
typedef array<Index, TensorEvaluator<LeftArgType, Device>::Dimensions::count> left_dim_mapper_t;
typedef array<Index, TensorEvaluator<RightArgType, Device>::Dimensions::count> right_dim_mapper_t;
typedef array<Index, internal::array_size<Indices>::value> contract_t;
typedef array<Index, max_n_1<TensorEvaluator<LeftArgType, Device>::Dimensions::count - internal::array_size<Indices>::value>::size> left_nocontract_t;
typedef array<Index, max_n_1<TensorEvaluator<RightArgType, Device>::Dimensions::count - internal::array_size<Indices>::value>::size> right_nocontract_t;
static const int NumDims = max_n_1<TensorEvaluator<LeftArgType, Device>::Dimensions::count + TensorEvaluator<RightArgType, Device>::Dimensions::count - 2 * internal::array_size<Indices>::value>::size;
typedef DSizes<Index, NumDims> Dimensions;
enum { enum {
IsAligned = true, IsAligned = true,
PacketAccess = (internal::packet_traits<Scalar>::size > 1), PacketAccess = (internal::packet_traits<Scalar>::size > 1),
Layout = TensorEvaluator<LeftArgType, Device>::Layout,
CoordAccess = false, // to be implemented
}; };
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TensorContractionEvaluatorBase(const XprType& op, const Device& device) // Most of the code is assuming that both input tensors are ColMajor. If the
: m_leftImpl(op.lhsExpression(), device), m_rightImpl(op.rhsExpression(), device), m_device(device), m_result(NULL) // 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<
Layout == ColMajor, LeftArgType, RightArgType>::type EvalLeftArgType;
typedef typename internal::conditional<
Layout == ColMajor, RightArgType, LeftArgType>::type EvalRightArgType;
static const int LDims =
internal::array_size<typename TensorEvaluator<EvalLeftArgType, Device>::Dimensions>::value;
static const int RDims =
internal::array_size<typename TensorEvaluator<EvalRightArgType, Device>::Dimensions>::value;
static const int ContractDims = internal::array_size<Indices>::value;
static const int NumDims = internal::max_n_1<LDims + RDims - 2 * ContractDims>::size;
typedef array<Index, LDims> left_dim_mapper_t;
typedef array<Index, RDims> right_dim_mapper_t;
typedef array<Index, ContractDims> contract_t;
typedef array<Index, internal::max_n_1<LDims - ContractDims>::size> left_nocontract_t;
typedef array<Index, internal::max_n_1<RDims - ContractDims>::size> right_nocontract_t;
typedef DSizes<Index, NumDims> Dimensions;
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
TensorContractionEvaluatorBase(const XprType& op, const Device& device)
: m_leftImpl(choose(Cond<Layout == ColMajor>(),
op.lhsExpression(), op.rhsExpression()), device),
m_rightImpl(choose(Cond<Layout == ColMajor>(),
op.rhsExpression(), op.lhsExpression()), device),
m_device(device),
m_result(NULL) {
EIGEN_STATIC_ASSERT((TensorEvaluator<LeftArgType, Device>::Layout ==
TensorEvaluator<RightArgType, Device>::Layout),
YOU_MADE_A_PROGRAMMING_MISTAKE);
eigen_assert((internal::array_size<contract_t>::value > 0) && "Must contract on some indices"); eigen_assert((internal::array_size<contract_t>::value > 0) && "Must contract on some indices");
array<Index, TensorEvaluator<LeftArgType, Device>::Dimensions::count> lhs_strides;
lhs_strides[0] = 1; DSizes<Index, LDims> eval_left_dims;
for (int i = 0; i < TensorEvaluator<LeftArgType, Device>::Dimensions::count-1; ++i) { DSizes<Index, RDims> eval_right_dims;
lhs_strides[i+1] = lhs_strides[i] * m_leftImpl.dimensions()[i]; array<IndexPair<Index>, ContractDims> eval_op_indices;
if (Layout == ColMajor) {
// For ColMajor, we keep using the existing dimensions
for (int i = 0; i < LDims; i++) {
eval_left_dims[i] = m_leftImpl.dimensions()[i];
}
for (int i = 0; i < RDims; i++) {
eval_right_dims[i] = m_rightImpl.dimensions()[i];
}
// We keep the pairs of contracting indices.
for (int i = 0; i < ContractDims; i++) {
eval_op_indices[i].first = op.indices()[i].first;
eval_op_indices[i].second = op.indices()[i].second;
}
} else {
// For RowMajor, we need to reverse the existing dimensions
for (int i = 0; i < LDims; i++) {
eval_left_dims[i] = m_leftImpl.dimensions()[LDims - i - 1];
}
for (int i = 0; i < RDims; i++) {
eval_right_dims[i] = m_rightImpl.dimensions()[RDims - i - 1];
}
// We need to flip all the pairs of contracting indices as well as
// reversing the dimensions.
for (int i = 0; i < ContractDims; i++) {
eval_op_indices[i].first = LDims - 1 - op.indices()[i].second;
eval_op_indices[i].second = RDims - 1 - op.indices()[i].first;
}
} }
array<Index, TensorEvaluator<RightArgType, Device>::Dimensions::count> rhs_strides; array<Index, LDims> lhs_strides;
lhs_strides[0] = 1;
for (int i = 0; i < LDims-1; ++i) {
lhs_strides[i+1] = lhs_strides[i] * eval_left_dims[i];
}
array<Index, RDims> rhs_strides;
rhs_strides[0] = 1; rhs_strides[0] = 1;
for (int i = 0; i < TensorEvaluator<RightArgType, Device>::Dimensions::count-1; ++i) { for (int i = 0; i < RDims-1; ++i) {
rhs_strides[i+1] = rhs_strides[i] * m_rightImpl.dimensions()[i]; rhs_strides[i+1] = rhs_strides[i] * eval_right_dims[i];
} }
m_i_strides[0] = 1; m_i_strides[0] = 1;
@ -515,27 +592,28 @@ struct TensorContractionEvaluatorBase
m_lhs_inner_dim_contiguous = true; m_lhs_inner_dim_contiguous = true;
int dim_idx = 0; int dim_idx = 0;
int nocontract_idx = 0; int nocontract_idx = 0;
const typename TensorEvaluator<LeftArgType, Device>::Dimensions& left_dims = m_leftImpl.dimensions();
for (int i = 0; i < TensorEvaluator<LeftArgType, Device>::Dimensions::count; i++) { for (int i = 0; i < LDims; i++) {
// find if we are contracting on index i of left tensor // find if we are contracting on index i of left tensor
bool contracting = false; bool contracting = false;
for (int j = 0; j < internal::array_size<Indices>::value; j++) { for (int j = 0; j < ContractDims; j++) {
if (op.indices()[j].first == i) { if (eval_op_indices[j].first == i) {
contracting = true; contracting = true;
break; break;
} }
} }
if (!contracting) { if (!contracting) {
// add dimension size to output dimensions // add dimension size to output dimensions
m_dimensions[dim_idx] = left_dims[i]; m_dimensions[dim_idx] = eval_left_dims[i];
m_left_nocontract_strides[nocontract_idx] = lhs_strides[i]; m_left_nocontract_strides[nocontract_idx] = lhs_strides[i];
if (dim_idx != i) { if (dim_idx != i) {
m_lhs_inner_dim_contiguous = false; m_lhs_inner_dim_contiguous = false;
} }
if (nocontract_idx+1 < internal::array_size<left_nocontract_t>::value) { if (nocontract_idx+1 < internal::array_size<left_nocontract_t>::value) {
m_i_strides[nocontract_idx+1] = m_i_strides[nocontract_idx] * left_dims[i]; m_i_strides[nocontract_idx+1] =
m_i_strides[nocontract_idx] * eval_left_dims[i];
} else { } else {
m_i_size = m_i_strides[nocontract_idx] * left_dims[i]; m_i_size = m_i_strides[nocontract_idx] * eval_left_dims[i];
} }
dim_idx++; dim_idx++;
nocontract_idx++; nocontract_idx++;
@ -543,22 +621,22 @@ struct TensorContractionEvaluatorBase
} }
nocontract_idx = 0; nocontract_idx = 0;
const typename TensorEvaluator<RightArgType, Device>::Dimensions& right_dims = m_rightImpl.dimensions(); for (int i = 0; i < RDims; i++) {
for (int i = 0; i < TensorEvaluator<RightArgType, Device>::Dimensions::count; i++) {
bool contracting = false; bool contracting = false;
// find if we are contracting on index i of right tensor // find if we are contracting on index i of right tensor
for (int j = 0; j < internal::array_size<Indices>::value; j++) { for (int j = 0; j < ContractDims; j++) {
if (op.indices()[j].second == i) { if (eval_op_indices[j].second == i) {
contracting = true; contracting = true;
break; break;
} }
} }
if (!contracting) { if (!contracting) {
m_dimensions[dim_idx] = right_dims[i]; m_dimensions[dim_idx] = eval_right_dims[i];
if (nocontract_idx+1 < internal::array_size<right_nocontract_t>::value) { if (nocontract_idx+1 < internal::array_size<right_nocontract_t>::value) {
m_j_strides[nocontract_idx+1] = m_j_strides[nocontract_idx] * right_dims[i]; m_j_strides[nocontract_idx+1] =
m_j_strides[nocontract_idx] * eval_right_dims[i];
} else { } else {
m_j_size = m_j_strides[nocontract_idx] * right_dims[i]; m_j_size = m_j_strides[nocontract_idx] * eval_right_dims[i];
} }
m_right_nocontract_strides[nocontract_idx] = rhs_strides[i]; m_right_nocontract_strides[nocontract_idx] = rhs_strides[i];
dim_idx++; dim_idx++;
@ -573,12 +651,13 @@ struct TensorContractionEvaluatorBase
// each tensor, we'll only look at the first tensor here. // each tensor, we'll only look at the first tensor here.
m_rhs_inner_dim_contiguous = true; m_rhs_inner_dim_contiguous = true;
m_rhs_inner_dim_reordered = false; m_rhs_inner_dim_reordered = false;
for (int i = 0; i < internal::array_size<Indices>::value; i++) { for (int i = 0; i < ContractDims; i++) {
Index left = op.indices()[i].first; Index left = eval_op_indices[i].first;
Index right = op.indices()[i].second; Index right = eval_op_indices[i].second;
Index size = left_dims[left]; Index size = eval_left_dims[left];
eigen_assert(size == right_dims[right] && "Contraction axes must be same size"); eigen_assert(size == eval_right_dims[right] &&
"Contraction axes must be same size");
if (i+1 < internal::array_size<contract_t>::value) { if (i+1 < internal::array_size<contract_t>::value) {
m_k_strides[i+1] = m_k_strides[i] * size; m_k_strides[i+1] = m_k_strides[i] * size;
@ -588,7 +667,7 @@ struct TensorContractionEvaluatorBase
m_left_contracting_strides[i] = lhs_strides[left]; m_left_contracting_strides[i] = lhs_strides[left];
m_right_contracting_strides[i] = rhs_strides[right]; m_right_contracting_strides[i] = rhs_strides[right];
if (i > 0 && right < op.indices()[i-1].second) { if (i > 0 && right < eval_op_indices[i-1].second) {
m_rhs_inner_dim_reordered = true; m_rhs_inner_dim_reordered = true;
} }
if (right != i) { if (right != i) {
@ -597,9 +676,16 @@ struct TensorContractionEvaluatorBase
} }
// Scalar case. We represent the result as a 1d tensor of size 1. // Scalar case. We represent the result as a 1d tensor of size 1.
if (TensorEvaluator<LeftArgType, Device>::Dimensions::count + TensorEvaluator<RightArgType, Device>::Dimensions::count == 2 * internal::array_size<Indices>::value) { if (LDims + RDims == 2 * ContractDims) {
m_dimensions[0] = 1; m_dimensions[0] = 1;
} }
// If the layout is RowMajor, we need to reverse the m_dimensions
if (Layout == RowMajor) {
for (int i = 0, j = NumDims - 1; i < j; i++, j--) {
std::swap(m_dimensions[i], m_dimensions[j]);
}
}
} }
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Dimensions& dimensions() const { return m_dimensions; } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Dimensions& dimensions() const { return m_dimensions; }
@ -661,10 +747,10 @@ struct TensorContractionEvaluatorBase
const Index rows = m_i_size; const Index rows = m_i_size;
const Index cols = m_k_size; const Index cols = m_k_size;
typedef typename internal::remove_const<typename LeftArgType::Scalar>::type LhsScalar; typedef typename internal::remove_const<typename EvalLeftArgType::Scalar>::type LhsScalar;
typedef typename internal::remove_const<typename RightArgType::Scalar>::type RhsScalar; typedef typename internal::remove_const<typename EvalRightArgType::Scalar>::type RhsScalar;
typedef TensorEvaluator<LeftArgType, Device> LeftEvaluator; typedef TensorEvaluator<EvalLeftArgType, Device> LeftEvaluator;
typedef TensorEvaluator<RightArgType, Device> RightEvaluator; typedef TensorEvaluator<EvalRightArgType, Device> RightEvaluator;
const int lhs_packet_size = internal::packet_traits<LhsScalar>::size; const int lhs_packet_size = internal::packet_traits<LhsScalar>::size;
const int rhs_packet_size = internal::packet_traits<RhsScalar>::size; const int rhs_packet_size = internal::packet_traits<RhsScalar>::size;
typedef internal::TensorContractionInputMapper<LhsScalar, Index, internal::Lhs, typedef internal::TensorContractionInputMapper<LhsScalar, Index, internal::Lhs,
@ -719,7 +805,6 @@ struct TensorContractionEvaluatorBase
protected: protected:
// Prevent assignment // Prevent assignment
TensorContractionEvaluatorBase& operator = (const TensorContractionEvaluatorBase&); TensorContractionEvaluatorBase& operator = (const TensorContractionEvaluatorBase&);
Dimensions m_dimensions; Dimensions m_dimensions;
contract_t m_k_strides; contract_t m_k_strides;
@ -739,16 +824,18 @@ struct TensorContractionEvaluatorBase
Index m_j_size; Index m_j_size;
Index m_k_size; Index m_k_size;
TensorEvaluator<LeftArgType, Device> m_leftImpl; TensorEvaluator<EvalLeftArgType, Device> m_leftImpl;
TensorEvaluator<RightArgType, Device> m_rightImpl; TensorEvaluator<EvalRightArgType, Device> m_rightImpl;
const Device& m_device; const Device& m_device;
Scalar* m_result; Scalar* m_result;
}; };
// evaluator for default device
template<typename Indices, typename LeftArgType, typename RightArgType, typename Device> template<typename Indices, typename LeftArgType, typename RightArgType, typename Device>
struct TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgType>, Device> : struct TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgType>, Device> :
public TensorContractionEvaluatorBase<TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgType>, Device> > { public TensorContractionEvaluatorBase<
TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgType>, Device> > {
typedef TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgType>, Device> Self; typedef TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgType>, Device> Self;
typedef TensorContractionEvaluatorBase<Self> Base; typedef TensorContractionEvaluatorBase<Self> Base;
@ -759,15 +846,35 @@ struct TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgT
typedef typename XprType::CoeffReturnType CoeffReturnType; typedef typename XprType::CoeffReturnType CoeffReturnType;
typedef typename XprType::PacketReturnType PacketReturnType; typedef typename XprType::PacketReturnType PacketReturnType;
typedef array<Index, TensorEvaluator<LeftArgType, Device>::Dimensions::count> left_dim_mapper_t; enum {
typedef array<Index, TensorEvaluator<RightArgType, Device>::Dimensions::count> right_dim_mapper_t; Layout = TensorEvaluator<LeftArgType, Device>::Layout,
};
typedef array<Index, internal::array_size<Indices>::value> contract_t; // Most of the code is assuming that both input tensors are ColMajor. If the
typedef array<Index, max_n_1<TensorEvaluator<LeftArgType, Device>::Dimensions::count - internal::array_size<Indices>::value>::size> left_nocontract_t; // inputs are RowMajor, we will "cheat" by swapping the LHS and RHS:
typedef array<Index, max_n_1<TensorEvaluator<RightArgType, Device>::Dimensions::count - internal::array_size<Indices>::value>::size> right_nocontract_t; // 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<
Layout == ColMajor, LeftArgType, RightArgType>::type EvalLeftArgType;
typedef typename internal::conditional<
Layout == ColMajor, RightArgType, LeftArgType>::type EvalRightArgType;
static const int NumDims = max_n_1<TensorEvaluator<LeftArgType, Device>::Dimensions::count + TensorEvaluator<RightArgType, Device>::Dimensions::count - 2 * internal::array_size<Indices>::value>::size; static const int LDims =
internal::array_size<typename TensorEvaluator<EvalLeftArgType, Device>::Dimensions>::value;
static const int RDims =
internal::array_size<typename TensorEvaluator<EvalRightArgType, Device>::Dimensions>::value;
static const int ContractDims = internal::array_size<Indices>::value;
typedef array<Index, LDims> left_dim_mapper_t;
typedef array<Index, RDims> right_dim_mapper_t;
typedef array<Index, ContractDims> contract_t;
typedef array<Index, internal::max_n_1<LDims - ContractDims>::size> left_nocontract_t;
typedef array<Index, internal::max_n_1<RDims - ContractDims>::size> right_nocontract_t;
static const int NumDims = internal::max_n_1<LDims + RDims - 2 * ContractDims>::size;
// Could we use NumDimensions here?
typedef DSizes<Index, NumDims> Dimensions; typedef DSizes<Index, NumDims> Dimensions;
@ -799,15 +906,15 @@ struct TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgT
this->m_device.memset(buffer, 0, m * n * sizeof(Scalar)); this->m_device.memset(buffer, 0, m * n * sizeof(Scalar));
// define mr, nr, and all of my data mapper types // define mr, nr, and all of my data mapper types
typedef typename internal::remove_const<typename LeftArgType::Scalar>::type LhsScalar; typedef typename internal::remove_const<typename EvalLeftArgType::Scalar>::type LhsScalar;
typedef typename internal::remove_const<typename RightArgType::Scalar>::type RhsScalar; typedef typename internal::remove_const<typename EvalRightArgType::Scalar>::type RhsScalar;
typedef typename internal::gebp_traits<LhsScalar, RhsScalar> Traits; typedef typename internal::gebp_traits<LhsScalar, RhsScalar> Traits;
const Index nr = Traits::nr; const Index nr = Traits::nr;
const Index mr = Traits::mr; const Index mr = Traits::mr;
typedef TensorEvaluator<LeftArgType, Device> LeftEvaluator; typedef TensorEvaluator<EvalLeftArgType, Device> LeftEvaluator;
typedef TensorEvaluator<RightArgType, Device> RightEvaluator; typedef TensorEvaluator<EvalRightArgType, Device> RightEvaluator;
const int lhs_packet_size = internal::packet_traits<LhsScalar>::size; const int lhs_packet_size = internal::packet_traits<LhsScalar>::size;
const int rhs_packet_size = internal::packet_traits<RhsScalar>::size; const int rhs_packet_size = internal::packet_traits<RhsScalar>::size;
@ -826,10 +933,10 @@ struct TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgT
typedef internal::blas_data_mapper<Scalar, Index, ColMajor> OutputMapper; typedef internal::blas_data_mapper<Scalar, Index, ColMajor> OutputMapper;
// Declare GEBP packing and kernel structs // Declare GEBP packing and kernel structs
internal::gemm_pack_lhs<LhsScalar, Index, typename LhsMapper::SubMapper, mr, Traits::LhsProgress, ColMajor> pack_lhs; internal::gemm_pack_lhs<LhsScalar, Index, typename LhsMapper::SubMapper, mr, Traits::LhsProgress, ColMajor> pack_lhs;
internal::gemm_pack_rhs<RhsScalar, Index, typename RhsMapper::SubMapper, nr, ColMajor> pack_rhs; internal::gemm_pack_rhs<RhsScalar, Index, typename RhsMapper::SubMapper, nr, ColMajor> pack_rhs;
internal::gebp_kernel<LhsScalar, RhsScalar, Index, OutputMapper, mr, nr, false, false> gebp; internal::gebp_kernel<LhsScalar, RhsScalar, Index, OutputMapper, mr, nr, false, false> gebp;
// initialize data mappers // initialize data mappers

View File

@ -1241,10 +1241,10 @@ struct TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgT
typedef array<Index, RDims> right_dim_mapper_t; typedef array<Index, RDims> right_dim_mapper_t;
typedef array<Index, ContractDims> contract_t; typedef array<Index, ContractDims> contract_t;
typedef array<Index, max_n_1<LDims - ContractDims>::size> left_nocontract_t; typedef array<Index, internal::max_n_1<LDims - ContractDims>::size> left_nocontract_t;
typedef array<Index, max_n_1<RDims - ContractDims>::size> right_nocontract_t; typedef array<Index, internal::max_n_1<RDims - ContractDims>::size> right_nocontract_t;
static const int NumDims = max_n_1<LDims + RDims - 2 * ContractDims>::size; static const int NumDims = internal::max_n_1<LDims + RDims - 2 * ContractDims>::size;
typedef DSizes<Index, NumDims> Dimensions; typedef DSizes<Index, NumDims> Dimensions;

View File

@ -70,24 +70,43 @@ struct TensorEvaluator<const TensorContractionOp<Indices, LeftArgType, RightArgT
typedef typename XprType::CoeffReturnType CoeffReturnType; typedef typename XprType::CoeffReturnType CoeffReturnType;
typedef typename XprType::PacketReturnType PacketReturnType; typedef typename XprType::PacketReturnType PacketReturnType;
typedef array<Index, TensorEvaluator<LeftArgType, Device>::Dimensions::count> left_dim_mapper_t; enum {
typedef array<Index, TensorEvaluator<RightArgType, Device>::Dimensions::count> right_dim_mapper_t; Layout = TensorEvaluator<LeftArgType, Device>::Layout,
};
typedef array<Index, internal::array_size<Indices>::value> contract_t; // Most of the code is assuming that both input tensors are ColMajor. If the
typedef array<Index, max_n_1<TensorEvaluator<LeftArgType, Device>::Dimensions::count - internal::array_size<Indices>::value>::size> left_nocontract_t; // inputs are RowMajor, we will "cheat" by swapping the LHS and RHS:
typedef array<Index, max_n_1<TensorEvaluator<RightArgType, Device>::Dimensions::count - internal::array_size<Indices>::value>::size> right_nocontract_t; // 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<
Layout == ColMajor, LeftArgType, RightArgType>::type EvalLeftArgType;
typedef typename internal::conditional<
Layout == ColMajor, RightArgType, LeftArgType>::type EvalRightArgType;
static const int NumDims = max_n_1<TensorEvaluator<LeftArgType, Device>::Dimensions::count + TensorEvaluator<RightArgType, Device>::Dimensions::count - 2 * internal::array_size<Indices>::value>::size; static const int LDims =
internal::array_size<typename TensorEvaluator<EvalLeftArgType, Device>::Dimensions>::value;
static const int RDims =
internal::array_size<typename TensorEvaluator<EvalRightArgType, Device>::Dimensions>::value;
static const int ContractDims = internal::array_size<Indices>::value;
typedef array<Index, LDims> left_dim_mapper_t;
typedef array<Index, RDims> right_dim_mapper_t;
typedef array<Index, ContractDims> contract_t;
typedef array<Index, max_n_1<LDims - ContractDims>::size> left_nocontract_t;
typedef array<Index, max_n_1<RDims - ContractDims>::size> right_nocontract_t;
static const int NumDims = max_n_1<LDims + RDims - 2 * ContractDims>::size;
typedef DSizes<Index, NumDims> Dimensions; typedef DSizes<Index, NumDims> Dimensions;
// typedefs needed in evalTo // typedefs needed in evalTo
typedef typename internal::remove_const<typename LeftArgType::Scalar>::type LhsScalar; typedef typename internal::remove_const<typename EvalLeftArgType::Scalar>::type LhsScalar;
typedef typename internal::remove_const<typename RightArgType::Scalar>::type RhsScalar; typedef typename internal::remove_const<typename EvalRightArgType::Scalar>::type RhsScalar;
typedef typename internal::gebp_traits<LhsScalar, RhsScalar> Traits; typedef typename internal::gebp_traits<LhsScalar, RhsScalar> Traits;
typedef TensorEvaluator<LeftArgType, Device> LeftEvaluator; typedef TensorEvaluator<EvalLeftArgType, Device> LeftEvaluator;
typedef TensorEvaluator<RightArgType, Device> RightEvaluator; typedef TensorEvaluator<EvalRightArgType, Device> RightEvaluator;
TensorEvaluator(const XprType& op, const Device& device) : TensorEvaluator(const XprType& op, const Device& device) :
Base(op, device) {} Base(op, device) {}