10 #ifndef EIGEN_GENERAL_MATRIX_MATRIX_H
11 #define EIGEN_GENERAL_MATRIX_MATRIX_H
17 template<
typename _LhsScalar,
typename _RhsScalar>
class level3_blocking;
22 typename LhsScalar,
int LhsStorageOrder,
bool ConjugateLhs,
23 typename RhsScalar,
int RhsStorageOrder,
bool ConjugateRhs>
24 struct general_matrix_matrix_product<Index,LhsScalar,LhsStorageOrder,ConjugateLhs,RhsScalar,RhsStorageOrder,ConjugateRhs,
RowMajor>
26 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
27 static EIGEN_STRONG_INLINE
void run(
28 Index rows, Index cols, Index depth,
29 const LhsScalar* lhs, Index lhsStride,
30 const RhsScalar* rhs, Index rhsStride,
31 ResScalar* res, Index resStride,
33 level3_blocking<RhsScalar,LhsScalar>& blocking,
34 GemmParallelInfo<Index>* info = 0)
37 general_matrix_matrix_product<Index,
41 ::run(cols,rows,depth,rhs,rhsStride,lhs,lhsStride,res,resStride,alpha,blocking,info);
49 typename LhsScalar,
int LhsStorageOrder,
bool ConjugateLhs,
50 typename RhsScalar,
int RhsStorageOrder,
bool ConjugateRhs>
51 struct general_matrix_matrix_product<Index,LhsScalar,LhsStorageOrder,ConjugateLhs,RhsScalar,RhsStorageOrder,ConjugateRhs,
ColMajor>
53 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
54 static void run(Index rows, Index cols, Index depth,
55 const LhsScalar* _lhs, Index lhsStride,
56 const RhsScalar* _rhs, Index rhsStride,
57 ResScalar* res, Index resStride,
59 level3_blocking<LhsScalar,RhsScalar>& blocking,
60 GemmParallelInfo<Index>* info = 0)
62 const_blas_data_mapper<LhsScalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
63 const_blas_data_mapper<RhsScalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
65 typedef gebp_traits<LhsScalar,RhsScalar> Traits;
67 Index kc = blocking.kc();
68 Index mc = (std::min)(rows,blocking.mc());
71 gemm_pack_lhs<LhsScalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
72 gemm_pack_rhs<RhsScalar, Index, Traits::nr, RhsStorageOrder> pack_rhs;
73 gebp_kernel<LhsScalar, RhsScalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp;
75 #ifdef EIGEN_HAS_OPENMP
79 Index tid = omp_get_thread_num();
80 Index threads = omp_get_num_threads();
82 std::size_t sizeA = kc*mc;
83 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
84 ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, sizeA, 0);
85 ei_declare_aligned_stack_constructed_variable(RhsScalar, w, sizeW, 0);
87 RhsScalar* blockB = blocking.blockB();
88 eigen_internal_assert(blockB!=0);
91 for(Index k=0; k<depth; k+=kc)
93 const Index actual_kc = (std::min)(k+kc,depth)-k;
97 pack_lhs(blockA, &lhs(0,k), lhsStride, actual_kc, mc);
105 while(info[tid].users!=0) {}
106 info[tid].users += threads;
108 pack_rhs(blockB+info[tid].rhs_start*actual_kc, &rhs(k,info[tid].rhs_start), rhsStride, actual_kc, info[tid].rhs_length);
114 for(Index shift=0; shift<threads; ++shift)
116 Index j = (tid+shift)%threads;
122 while(info[j].sync!=k) {}
124 gebp(res+info[j].rhs_start*resStride, resStride, blockA, blockB+info[j].rhs_start*actual_kc, mc, actual_kc, info[j].rhs_length, alpha, -1,-1,0,0, w);
128 for(Index i=mc; i<rows; i+=mc)
130 const Index actual_mc = (std::min)(i+mc,rows)-i;
133 pack_lhs(blockA, &lhs(i,k), lhsStride, actual_kc, actual_mc);
136 gebp(res+i, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha, -1,-1,0,0, w);
141 for(Index j=0; j<threads; ++j)
147 #endif // EIGEN_HAS_OPENMP
149 EIGEN_UNUSED_VARIABLE(info);
152 std::size_t sizeA = kc*mc;
153 std::size_t sizeB = kc*cols;
154 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
156 ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, sizeA, blocking.blockA());
157 ei_declare_aligned_stack_constructed_variable(RhsScalar, blockB, sizeB, blocking.blockB());
158 ei_declare_aligned_stack_constructed_variable(RhsScalar, blockW, sizeW, blocking.blockW());
162 for(Index k2=0; k2<depth; k2+=kc)
164 const Index actual_kc = (std::min)(k2+kc,depth)-k2;
170 pack_rhs(blockB, &rhs(k2,0), rhsStride, actual_kc, cols);
175 for(Index i2=0; i2<rows; i2+=mc)
177 const Index actual_mc = (std::min)(i2+mc,rows)-i2;
182 pack_lhs(blockA, &lhs(i2,k2), lhsStride, actual_kc, actual_mc);
185 gebp(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha, -1, -1, 0, 0, blockW);
199 template<
typename Lhs,
typename Rhs>
200 struct traits<GeneralProduct<Lhs,Rhs,GemmProduct> >
201 : traits<ProductBase<GeneralProduct<Lhs,Rhs,GemmProduct>, Lhs, Rhs> >
204 template<
typename Scalar,
typename Index,
typename Gemm,
typename Lhs,
typename Rhs,
typename Dest,
typename BlockingType>
207 gemm_functor(
const Lhs& lhs,
const Rhs& rhs, Dest& dest, Scalar actualAlpha,
208 BlockingType& blocking)
209 : m_lhs(lhs), m_rhs(rhs), m_dest(dest), m_actualAlpha(actualAlpha), m_blocking(blocking)
212 void initParallelSession()
const
214 m_blocking.allocateB();
217 void operator() (Index row, Index rows, Index col=0, Index cols=-1, GemmParallelInfo<Index>* info=0)
const
222 Gemm::run(rows, cols, m_lhs.cols(),
223 &m_lhs.coeffRef(row,0), m_lhs.outerStride(),
224 &m_rhs.coeffRef(0,col), m_rhs.outerStride(),
225 (Scalar*)&(m_dest.coeffRef(row,col)), m_dest.outerStride(),
226 m_actualAlpha, m_blocking, info);
233 Scalar m_actualAlpha;
234 BlockingType& m_blocking;
237 template<
int StorageOrder,
typename LhsScalar,
typename RhsScalar,
int MaxRows,
int MaxCols,
int MaxDepth,
int KcFactor=1,
238 bool FiniteAtCompileTime = MaxRows!=
Dynamic && MaxCols!=
Dynamic && MaxDepth !=
Dynamic>
class gemm_blocking_space;
240 template<
typename _LhsScalar,
typename _RhsScalar>
241 class level3_blocking
243 typedef _LhsScalar LhsScalar;
244 typedef _RhsScalar RhsScalar;
258 : m_blockA(0), m_blockB(0), m_blockW(0), m_mc(0), m_nc(0), m_kc(0)
261 inline DenseIndex mc()
const {
return m_mc; }
262 inline DenseIndex nc()
const {
return m_nc; }
263 inline DenseIndex kc()
const {
return m_kc; }
265 inline LhsScalar* blockA() {
return m_blockA; }
266 inline RhsScalar* blockB() {
return m_blockB; }
267 inline RhsScalar* blockW() {
return m_blockW; }
270 template<
int StorageOrder,
typename _LhsScalar,
typename _RhsScalar,
int MaxRows,
int MaxCols,
int MaxDepth,
int KcFactor>
271 class gemm_blocking_space<StorageOrder,_LhsScalar,_RhsScalar,MaxRows, MaxCols, MaxDepth, KcFactor, true>
272 :
public level3_blocking<
273 typename conditional<StorageOrder==RowMajor,_RhsScalar,_LhsScalar>::type,
274 typename conditional<StorageOrder==RowMajor,_LhsScalar,_RhsScalar>::type>
278 ActualRows = Transpose ? MaxCols : MaxRows,
279 ActualCols = Transpose ? MaxRows : MaxCols
281 typedef typename conditional<Transpose,_RhsScalar,_LhsScalar>::type LhsScalar;
282 typedef typename conditional<Transpose,_LhsScalar,_RhsScalar>::type RhsScalar;
283 typedef gebp_traits<LhsScalar,RhsScalar> Traits;
285 SizeA = ActualRows * MaxDepth,
286 SizeB = ActualCols * MaxDepth,
287 SizeW = MaxDepth * Traits::WorkSpaceFactor
290 EIGEN_ALIGN16 LhsScalar m_staticA[SizeA];
291 EIGEN_ALIGN16 RhsScalar m_staticB[SizeB];
292 EIGEN_ALIGN16 RhsScalar m_staticW[SizeW];
296 gemm_blocking_space(DenseIndex , DenseIndex , DenseIndex )
298 this->m_mc = ActualRows;
299 this->m_nc = ActualCols;
300 this->m_kc = MaxDepth;
301 this->m_blockA = m_staticA;
302 this->m_blockB = m_staticB;
303 this->m_blockW = m_staticW;
306 inline void allocateA() {}
307 inline void allocateB() {}
308 inline void allocateW() {}
309 inline void allocateAll() {}
312 template<
int StorageOrder,
typename _LhsScalar,
typename _RhsScalar,
int MaxRows,
int MaxCols,
int MaxDepth,
int KcFactor>
313 class gemm_blocking_space<StorageOrder,_LhsScalar,_RhsScalar,MaxRows, MaxCols, MaxDepth, KcFactor, false>
314 :
public level3_blocking<
315 typename conditional<StorageOrder==RowMajor,_RhsScalar,_LhsScalar>::type,
316 typename conditional<StorageOrder==RowMajor,_LhsScalar,_RhsScalar>::type>
321 typedef typename conditional<Transpose,_RhsScalar,_LhsScalar>::type LhsScalar;
322 typedef typename conditional<Transpose,_LhsScalar,_RhsScalar>::type RhsScalar;
323 typedef gebp_traits<LhsScalar,RhsScalar> Traits;
331 gemm_blocking_space(DenseIndex rows, DenseIndex cols, DenseIndex depth)
333 this->m_mc = Transpose ? cols : rows;
334 this->m_nc = Transpose ? rows : cols;
337 computeProductBlockingSizes<LhsScalar,RhsScalar,KcFactor>(this->m_kc, this->m_mc, this->m_nc);
338 m_sizeA = this->m_mc * this->m_kc;
339 m_sizeB = this->m_kc * this->m_nc;
340 m_sizeW = this->m_kc*Traits::WorkSpaceFactor;
345 if(this->m_blockA==0)
346 this->m_blockA = aligned_new<LhsScalar>(m_sizeA);
351 if(this->m_blockB==0)
352 this->m_blockB = aligned_new<RhsScalar>(m_sizeB);
357 if(this->m_blockW==0)
358 this->m_blockW = aligned_new<RhsScalar>(m_sizeW);
368 ~gemm_blocking_space()
370 aligned_delete(this->m_blockA, m_sizeA);
371 aligned_delete(this->m_blockB, m_sizeB);
372 aligned_delete(this->m_blockW, m_sizeW);
378 template<
typename Lhs,
typename Rhs>
379 class GeneralProduct<Lhs, Rhs, GemmProduct>
380 :
public ProductBase<GeneralProduct<Lhs,Rhs,GemmProduct>, Lhs, Rhs>
383 MaxDepthAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(Lhs::MaxColsAtCompileTime,Rhs::MaxRowsAtCompileTime)
386 EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
388 typedef typename Lhs::Scalar LhsScalar;
389 typedef typename Rhs::Scalar RhsScalar;
390 typedef Scalar ResScalar;
392 GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
394 typedef internal::scalar_product_op<LhsScalar,RhsScalar> BinOp;
395 EIGEN_CHECK_BINARY_COMPATIBILIY(BinOp,LhsScalar,RhsScalar);
398 template<
typename Dest>
void scaleAndAddTo(Dest& dst, Scalar alpha)
const
400 eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
402 typename internal::add_const_on_value_type<ActualLhsType>::type lhs = LhsBlasTraits::extract(m_lhs);
403 typename internal::add_const_on_value_type<ActualRhsType>::type rhs = RhsBlasTraits::extract(m_rhs);
405 Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
406 * RhsBlasTraits::extractScalarFactor(m_rhs);
409 Dest::MaxRowsAtCompileTime,Dest::MaxColsAtCompileTime,MaxDepthAtCompileTime> BlockingType;
411 typedef internal::gemm_functor<
413 internal::general_matrix_matrix_product<
415 LhsScalar, (_ActualLhsType::Flags&
RowMajorBit) ?
RowMajor : ColMajor,
bool(LhsBlasTraits::NeedToConjugate),
416 RhsScalar, (_ActualRhsType::Flags&
RowMajorBit) ?
RowMajor : ColMajor,
bool(RhsBlasTraits::NeedToConjugate),
418 _ActualLhsType, _ActualRhsType, Dest, BlockingType> GemmFunctor;
420 BlockingType blocking(dst.rows(), dst.cols(), lhs.cols());
422 internal::parallelize_gemm<(Dest::MaxRowsAtCompileTime>32 || Dest::MaxRowsAtCompileTime==
Dynamic)>(GemmFunctor(lhs, rhs, dst, actualAlpha, blocking), this->rows(), this->cols(), Dest::Flags&
RowMajorBit);
428 #endif // EIGEN_GENERAL_MATRIX_MATRIX_H