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>
54 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
55 static void run(Index rows, Index cols, Index depth,
56 const LhsScalar* _lhs, Index lhsStride,
57 const RhsScalar* _rhs, Index rhsStride,
58 ResScalar* res, Index resStride,
60 level3_blocking<LhsScalar,RhsScalar>& blocking,
61 GemmParallelInfo<Index>* info = 0)
63 const_blas_data_mapper<LhsScalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
64 const_blas_data_mapper<RhsScalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
66 typedef gebp_traits<LhsScalar,RhsScalar> Traits;
68 Index kc = blocking.kc();
69 Index mc = (std::min)(rows,blocking.mc());
72 gemm_pack_lhs<LhsScalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
73 gemm_pack_rhs<RhsScalar, Index, Traits::nr, RhsStorageOrder> pack_rhs;
74 gebp_kernel<LhsScalar, RhsScalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp;
76 #ifdef EIGEN_HAS_OPENMP
80 Index tid = omp_get_thread_num();
81 Index threads = omp_get_num_threads();
83 std::size_t sizeA = kc*mc;
84 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
85 ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, sizeA, 0);
86 ei_declare_aligned_stack_constructed_variable(RhsScalar, w, sizeW, 0);
88 RhsScalar* blockB = blocking.blockB();
89 eigen_internal_assert(blockB!=0);
92 for(Index k=0; k<depth; k+=kc)
94 const Index actual_kc = (std::min)(k+kc,depth)-k;
98 pack_lhs(blockA, &lhs(0,k), lhsStride, actual_kc, mc);
106 while(info[tid].users!=0) {}
107 info[tid].users += threads;
109 pack_rhs(blockB+info[tid].rhs_start*actual_kc, &rhs(k,info[tid].rhs_start), rhsStride, actual_kc, info[tid].rhs_length);
115 for(Index shift=0; shift<threads; ++shift)
117 Index j = (tid+shift)%threads;
123 while(info[j].sync!=k) {}
125 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);
129 for(Index i=mc; i<rows; i+=mc)
131 const Index actual_mc = (std::min)(i+mc,rows)-i;
134 pack_lhs(blockA, &lhs(i,k), lhsStride, actual_kc, actual_mc);
137 gebp(res+i, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha, -1,-1,0,0, w);
142 for(Index j=0; j<threads; ++j)
148 #endif // EIGEN_HAS_OPENMP
150 EIGEN_UNUSED_VARIABLE(info);
153 std::size_t sizeA = kc*mc;
154 std::size_t sizeB = kc*cols;
155 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
157 ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, sizeA, blocking.blockA());
158 ei_declare_aligned_stack_constructed_variable(RhsScalar, blockB, sizeB, blocking.blockB());
159 ei_declare_aligned_stack_constructed_variable(RhsScalar, blockW, sizeW, blocking.blockW());
163 for(Index k2=0; k2<depth; k2+=kc)
165 const Index actual_kc = (std::min)(k2+kc,depth)-k2;
171 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);
198 template<
typename Lhs,
typename Rhs>
199 struct traits<GeneralProduct<Lhs,Rhs,GemmProduct> >
200 : traits<ProductBase<GeneralProduct<Lhs,Rhs,GemmProduct>, Lhs, Rhs> >
203 template<
typename Scalar,
typename Index,
typename Gemm,
typename Lhs,
typename Rhs,
typename Dest,
typename BlockingType>
206 gemm_functor(
const Lhs& lhs,
const Rhs& rhs, Dest& dest,
const Scalar& actualAlpha,
207 BlockingType& blocking)
208 : m_lhs(lhs), m_rhs(rhs), m_dest(dest), m_actualAlpha(actualAlpha), m_blocking(blocking)
211 void initParallelSession()
const
213 m_blocking.allocateB();
216 void operator() (Index row, Index rows, Index col=0, Index cols=-1, GemmParallelInfo<Index>* info=0)
const
221 Gemm::run(rows, cols, m_lhs.cols(),
222 &m_lhs.coeffRef(row,0), m_lhs.outerStride(),
223 &m_rhs.coeffRef(0,col), m_rhs.outerStride(),
224 (Scalar*)&(m_dest.coeffRef(row,col)), m_dest.outerStride(),
225 m_actualAlpha, m_blocking, info);
232 Scalar m_actualAlpha;
233 BlockingType& m_blocking;
236 template<
int StorageOrder,
typename LhsScalar,
typename RhsScalar,
int MaxRows,
int MaxCols,
int MaxDepth,
int KcFactor=1,
237 bool FiniteAtCompileTime = MaxRows!=Dynamic && MaxCols!=Dynamic && MaxDepth != Dynamic>
class gemm_blocking_space;
239 template<
typename _LhsScalar,
typename _RhsScalar>
240 class level3_blocking
242 typedef _LhsScalar LhsScalar;
243 typedef _RhsScalar RhsScalar;
257 : m_blockA(0), m_blockB(0), m_blockW(0), m_mc(0), m_nc(0), m_kc(0)
260 inline DenseIndex mc()
const {
return m_mc; }
261 inline DenseIndex nc()
const {
return m_nc; }
262 inline DenseIndex kc()
const {
return m_kc; }
264 inline LhsScalar* blockA() {
return m_blockA; }
265 inline RhsScalar* blockB() {
return m_blockB; }
266 inline RhsScalar* blockW() {
return m_blockW; }
269 template<
int StorageOrder,
typename _LhsScalar,
typename _RhsScalar,
int MaxRows,
int MaxCols,
int MaxDepth,
int KcFactor>
270 class gemm_blocking_space<StorageOrder,_LhsScalar,_RhsScalar,MaxRows, MaxCols, MaxDepth, KcFactor, true>
271 :
public level3_blocking<
272 typename conditional<StorageOrder==RowMajor,_RhsScalar,_LhsScalar>::type,
273 typename conditional<StorageOrder==RowMajor,_LhsScalar,_RhsScalar>::type>
277 ActualRows = Transpose ? MaxCols : MaxRows,
278 ActualCols = Transpose ? MaxRows : MaxCols
280 typedef typename conditional<Transpose,_RhsScalar,_LhsScalar>::type LhsScalar;
281 typedef typename conditional<Transpose,_LhsScalar,_RhsScalar>::type RhsScalar;
282 typedef gebp_traits<LhsScalar,RhsScalar> Traits;
284 SizeA = ActualRows * MaxDepth,
285 SizeB = ActualCols * MaxDepth,
286 SizeW = MaxDepth * Traits::WorkSpaceFactor
289 EIGEN_ALIGN16 LhsScalar m_staticA[SizeA];
290 EIGEN_ALIGN16 RhsScalar m_staticB[SizeB];
291 EIGEN_ALIGN16 RhsScalar m_staticW[SizeW];
295 gemm_blocking_space(DenseIndex , DenseIndex , DenseIndex )
297 this->m_mc = ActualRows;
298 this->m_nc = ActualCols;
299 this->m_kc = MaxDepth;
300 this->m_blockA = m_staticA;
301 this->m_blockB = m_staticB;
302 this->m_blockW = m_staticW;
305 inline void allocateA() {}
306 inline void allocateB() {}
307 inline void allocateW() {}
308 inline void allocateAll() {}
311 template<
int StorageOrder,
typename _LhsScalar,
typename _RhsScalar,
int MaxRows,
int MaxCols,
int MaxDepth,
int KcFactor>
312 class gemm_blocking_space<StorageOrder,_LhsScalar,_RhsScalar,MaxRows, MaxCols, MaxDepth, KcFactor, false>
313 :
public level3_blocking<
314 typename conditional<StorageOrder==RowMajor,_RhsScalar,_LhsScalar>::type,
315 typename conditional<StorageOrder==RowMajor,_LhsScalar,_RhsScalar>::type>
320 typedef typename conditional<Transpose,_RhsScalar,_LhsScalar>::type LhsScalar;
321 typedef typename conditional<Transpose,_LhsScalar,_RhsScalar>::type RhsScalar;
322 typedef gebp_traits<LhsScalar,RhsScalar> Traits;
330 gemm_blocking_space(DenseIndex rows, DenseIndex cols, DenseIndex depth)
332 this->m_mc = Transpose ? cols : rows;
333 this->m_nc = Transpose ? rows : cols;
336 computeProductBlockingSizes<LhsScalar,RhsScalar,KcFactor>(this->m_kc, this->m_mc, this->m_nc);
337 m_sizeA = this->m_mc * this->m_kc;
338 m_sizeB = this->m_kc * this->m_nc;
339 m_sizeW = this->m_kc*Traits::WorkSpaceFactor;
344 if(this->m_blockA==0)
345 this->m_blockA = aligned_new<LhsScalar>(m_sizeA);
350 if(this->m_blockB==0)
351 this->m_blockB = aligned_new<RhsScalar>(m_sizeB);
356 if(this->m_blockW==0)
357 this->m_blockW = aligned_new<RhsScalar>(m_sizeW);
367 ~gemm_blocking_space()
369 aligned_delete(this->m_blockA, m_sizeA);
370 aligned_delete(this->m_blockB, m_sizeB);
371 aligned_delete(this->m_blockW, m_sizeW);
377 template<
typename Lhs,
typename Rhs>
378 class GeneralProduct<Lhs, Rhs, GemmProduct>
379 :
public ProductBase<GeneralProduct<Lhs,Rhs,GemmProduct>, Lhs, Rhs>
382 MaxDepthAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(Lhs::MaxColsAtCompileTime,Rhs::MaxRowsAtCompileTime)
385 EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
387 typedef typename Lhs::Scalar LhsScalar;
388 typedef typename Rhs::Scalar RhsScalar;
389 typedef Scalar ResScalar;
391 GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
393 typedef internal::scalar_product_op<LhsScalar,RhsScalar> BinOp;
394 EIGEN_CHECK_BINARY_COMPATIBILIY(BinOp,LhsScalar,RhsScalar);
397 template<
typename Dest>
void scaleAndAddTo(Dest& dst,
const Scalar& alpha)
const
399 eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
401 typename internal::add_const_on_value_type<ActualLhsType>::type lhs = LhsBlasTraits::extract(m_lhs);
402 typename internal::add_const_on_value_type<ActualRhsType>::type rhs = RhsBlasTraits::extract(m_rhs);
404 Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
405 * RhsBlasTraits::extractScalarFactor(m_rhs);
408 Dest::MaxRowsAtCompileTime,Dest::MaxColsAtCompileTime,MaxDepthAtCompileTime> BlockingType;
410 typedef internal::gemm_functor<
412 internal::general_matrix_matrix_product<
414 LhsScalar, (_ActualLhsType::Flags&
RowMajorBit) ?
RowMajor : ColMajor,
bool(LhsBlasTraits::NeedToConjugate),
415 RhsScalar, (_ActualRhsType::Flags&
RowMajorBit) ?
RowMajor : ColMajor,
bool(RhsBlasTraits::NeedToConjugate),
417 _ActualLhsType, _ActualRhsType, Dest, BlockingType> GemmFunctor;
419 BlockingType blocking(dst.rows(), dst.cols(), lhs.cols());
421 internal::parallelize_gemm<(Dest::MaxRowsAtCompileTime>32 || Dest::MaxRowsAtCompileTime==Dynamic)>(GemmFunctor(lhs, rhs, dst, actualAlpha, blocking), this->rows(), this->cols(), Dest::Flags&
RowMajorBit);
427 #endif // EIGEN_GENERAL_MATRIX_MATRIX_H
Definition: Eigen_Colamd.h:54
const unsigned int RowMajorBit
Definition: Constants.h:53
Definition: Constants.h:266
Definition: Constants.h:264