00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_GENERAL_MATRIX_MATRIX_H
00026 #define EIGEN_GENERAL_MATRIX_MATRIX_H
00027
00028 namespace Eigen {
00029
00030 namespace internal {
00031
00032 template<typename _LhsScalar, typename _RhsScalar> class level3_blocking;
00033
00034
00035 template<
00036 typename Index,
00037 typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs,
00038 typename RhsScalar, int RhsStorageOrder, bool ConjugateRhs>
00039 struct general_matrix_matrix_product<Index,LhsScalar,LhsStorageOrder,ConjugateLhs,RhsScalar,RhsStorageOrder,ConjugateRhs,RowMajor>
00040 {
00041 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
00042 static EIGEN_STRONG_INLINE void run(
00043 Index rows, Index cols, Index depth,
00044 const LhsScalar* lhs, Index lhsStride,
00045 const RhsScalar* rhs, Index rhsStride,
00046 ResScalar* res, Index resStride,
00047 ResScalar alpha,
00048 level3_blocking<RhsScalar,LhsScalar>& blocking,
00049 GemmParallelInfo<Index>* info = 0)
00050 {
00051
00052 general_matrix_matrix_product<Index,
00053 RhsScalar, RhsStorageOrder==RowMajor ? ColMajor : RowMajor, ConjugateRhs,
00054 LhsScalar, LhsStorageOrder==RowMajor ? ColMajor : RowMajor, ConjugateLhs,
00055 ColMajor>
00056 ::run(cols,rows,depth,rhs,rhsStride,lhs,lhsStride,res,resStride,alpha,blocking,info);
00057 }
00058 };
00059
00060
00061
00062 template<
00063 typename Index,
00064 typename LhsScalar, int LhsStorageOrder, bool ConjugateLhs,
00065 typename RhsScalar, int RhsStorageOrder, bool ConjugateRhs>
00066 struct general_matrix_matrix_product<Index,LhsScalar,LhsStorageOrder,ConjugateLhs,RhsScalar,RhsStorageOrder,ConjugateRhs,ColMajor>
00067 {
00068 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
00069 static void run(Index rows, Index cols, Index depth,
00070 const LhsScalar* _lhs, Index lhsStride,
00071 const RhsScalar* _rhs, Index rhsStride,
00072 ResScalar* res, Index resStride,
00073 ResScalar alpha,
00074 level3_blocking<LhsScalar,RhsScalar>& blocking,
00075 GemmParallelInfo<Index>* info = 0)
00076 {
00077 const_blas_data_mapper<LhsScalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
00078 const_blas_data_mapper<RhsScalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
00079
00080 typedef gebp_traits<LhsScalar,RhsScalar> Traits;
00081
00082 Index kc = blocking.kc();
00083 Index mc = (std::min)(rows,blocking.mc());
00084
00085
00086 gemm_pack_lhs<LhsScalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
00087 gemm_pack_rhs<RhsScalar, Index, Traits::nr, RhsStorageOrder> pack_rhs;
00088 gebp_kernel<LhsScalar, RhsScalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp;
00089
00090 #ifdef EIGEN_HAS_OPENMP
00091 if(info)
00092 {
00093
00094 Index tid = omp_get_thread_num();
00095 Index threads = omp_get_num_threads();
00096
00097 std::size_t sizeA = kc*mc;
00098 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
00099 ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, sizeA, 0);
00100 ei_declare_aligned_stack_constructed_variable(RhsScalar, w, sizeW, 0);
00101
00102 RhsScalar* blockB = blocking.blockB();
00103 eigen_internal_assert(blockB!=0);
00104
00105
00106 for(Index k=0; k<depth; k+=kc)
00107 {
00108 const Index actual_kc = (std::min)(k+kc,depth)-k;
00109
00110
00111
00112 pack_lhs(blockA, &lhs(0,k), lhsStride, actual_kc, mc);
00113
00114
00115
00116
00117
00118
00119
00120 while(info[tid].users!=0) {}
00121 info[tid].users += threads;
00122
00123 pack_rhs(blockB+info[tid].rhs_start*actual_kc, &rhs(k,info[tid].rhs_start), rhsStride, actual_kc, info[tid].rhs_length);
00124
00125
00126 info[tid].sync = k;
00127
00128
00129 for(Index shift=0; shift<threads; ++shift)
00130 {
00131 Index j = (tid+shift)%threads;
00132
00133
00134
00135
00136 if(shift>0)
00137 while(info[j].sync!=k) {}
00138
00139 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);
00140 }
00141
00142
00143 for(Index i=mc; i<rows; i+=mc)
00144 {
00145 const Index actual_mc = (std::min)(i+mc,rows)-i;
00146
00147
00148 pack_lhs(blockA, &lhs(i,k), lhsStride, actual_kc, actual_mc);
00149
00150
00151 gebp(res+i, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha, -1,-1,0,0, w);
00152 }
00153
00154
00155
00156 for(Index j=0; j<threads; ++j)
00157 #pragma omp atomic
00158 --(info[j].users);
00159 }
00160 }
00161 else
00162 #endif // EIGEN_HAS_OPENMP
00163 {
00164 EIGEN_UNUSED_VARIABLE(info);
00165
00166
00167 std::size_t sizeA = kc*mc;
00168 std::size_t sizeB = kc*cols;
00169 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
00170
00171 ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, sizeA, blocking.blockA());
00172 ei_declare_aligned_stack_constructed_variable(RhsScalar, blockB, sizeB, blocking.blockB());
00173 ei_declare_aligned_stack_constructed_variable(RhsScalar, blockW, sizeW, blocking.blockW());
00174
00175
00176
00177 for(Index k2=0; k2<depth; k2+=kc)
00178 {
00179 const Index actual_kc = (std::min)(k2+kc,depth)-k2;
00180
00181
00182
00183
00184
00185 pack_rhs(blockB, &rhs(k2,0), rhsStride, actual_kc, cols);
00186
00187
00188
00189
00190 for(Index i2=0; i2<rows; i2+=mc)
00191 {
00192 const Index actual_mc = (std::min)(i2+mc,rows)-i2;
00193
00194
00195
00196
00197 pack_lhs(blockA, &lhs(i2,k2), lhsStride, actual_kc, actual_mc);
00198
00199
00200 gebp(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha, -1, -1, 0, 0, blockW);
00201
00202 }
00203 }
00204 }
00205 }
00206
00207 };
00208
00209
00210
00211
00212
00213
00214 template<typename Lhs, typename Rhs>
00215 struct traits<GeneralProduct<Lhs,Rhs,GemmProduct> >
00216 : traits<ProductBase<GeneralProduct<Lhs,Rhs,GemmProduct>, Lhs, Rhs> >
00217 {};
00218
00219 template<typename Scalar, typename Index, typename Gemm, typename Lhs, typename Rhs, typename Dest, typename BlockingType>
00220 struct gemm_functor
00221 {
00222 gemm_functor(const Lhs& lhs, const Rhs& rhs, Dest& dest, Scalar actualAlpha,
00223 BlockingType& blocking)
00224 : m_lhs(lhs), m_rhs(rhs), m_dest(dest), m_actualAlpha(actualAlpha), m_blocking(blocking)
00225 {}
00226
00227 void initParallelSession() const
00228 {
00229 m_blocking.allocateB();
00230 }
00231
00232 void operator() (Index row, Index rows, Index col=0, Index cols=-1, GemmParallelInfo<Index>* info=0) const
00233 {
00234 if(cols==-1)
00235 cols = m_rhs.cols();
00236
00237 Gemm::run(rows, cols, m_lhs.cols(),
00238 &m_lhs.coeffRef(row,0), m_lhs.outerStride(),
00239 &m_rhs.coeffRef(0,col), m_rhs.outerStride(),
00240 (Scalar*)&(m_dest.coeffRef(row,col)), m_dest.outerStride(),
00241 m_actualAlpha, m_blocking, info);
00242 }
00243
00244 protected:
00245 const Lhs& m_lhs;
00246 const Rhs& m_rhs;
00247 Dest& m_dest;
00248 Scalar m_actualAlpha;
00249 BlockingType& m_blocking;
00250 };
00251
00252 template<int StorageOrder, typename LhsScalar, typename RhsScalar, int MaxRows, int MaxCols, int MaxDepth,
00253 bool FiniteAtCompileTime = MaxRows!=Dynamic && MaxCols!=Dynamic && MaxDepth != Dynamic> class gemm_blocking_space;
00254
00255 template<typename _LhsScalar, typename _RhsScalar>
00256 class level3_blocking
00257 {
00258 typedef _LhsScalar LhsScalar;
00259 typedef _RhsScalar RhsScalar;
00260
00261 protected:
00262 LhsScalar* m_blockA;
00263 RhsScalar* m_blockB;
00264 RhsScalar* m_blockW;
00265
00266 DenseIndex m_mc;
00267 DenseIndex m_nc;
00268 DenseIndex m_kc;
00269
00270 public:
00271
00272 level3_blocking()
00273 : m_blockA(0), m_blockB(0), m_blockW(0), m_mc(0), m_nc(0), m_kc(0)
00274 {}
00275
00276 inline DenseIndex mc() const { return m_mc; }
00277 inline DenseIndex nc() const { return m_nc; }
00278 inline DenseIndex kc() const { return m_kc; }
00279
00280 inline LhsScalar* blockA() { return m_blockA; }
00281 inline RhsScalar* blockB() { return m_blockB; }
00282 inline RhsScalar* blockW() { return m_blockW; }
00283 };
00284
00285 template<int StorageOrder, typename _LhsScalar, typename _RhsScalar, int MaxRows, int MaxCols, int MaxDepth>
00286 class gemm_blocking_space<StorageOrder,_LhsScalar,_RhsScalar,MaxRows, MaxCols, MaxDepth, true>
00287 : public level3_blocking<
00288 typename conditional<StorageOrder==RowMajor,_RhsScalar,_LhsScalar>::type,
00289 typename conditional<StorageOrder==RowMajor,_LhsScalar,_RhsScalar>::type>
00290 {
00291 enum {
00292 Transpose = StorageOrder==RowMajor,
00293 ActualRows = Transpose ? MaxCols : MaxRows,
00294 ActualCols = Transpose ? MaxRows : MaxCols
00295 };
00296 typedef typename conditional<Transpose,_RhsScalar,_LhsScalar>::type LhsScalar;
00297 typedef typename conditional<Transpose,_LhsScalar,_RhsScalar>::type RhsScalar;
00298 typedef gebp_traits<LhsScalar,RhsScalar> Traits;
00299 enum {
00300 SizeA = ActualRows * MaxDepth,
00301 SizeB = ActualCols * MaxDepth,
00302 SizeW = MaxDepth * Traits::WorkSpaceFactor
00303 };
00304
00305 EIGEN_ALIGN16 LhsScalar m_staticA[SizeA];
00306 EIGEN_ALIGN16 RhsScalar m_staticB[SizeB];
00307 EIGEN_ALIGN16 RhsScalar m_staticW[SizeW];
00308
00309 public:
00310
00311 gemm_blocking_space(DenseIndex , DenseIndex , DenseIndex )
00312 {
00313 this->m_mc = ActualRows;
00314 this->m_nc = ActualCols;
00315 this->m_kc = MaxDepth;
00316 this->m_blockA = m_staticA;
00317 this->m_blockB = m_staticB;
00318 this->m_blockW = m_staticW;
00319 }
00320
00321 inline void allocateA() {}
00322 inline void allocateB() {}
00323 inline void allocateW() {}
00324 inline void allocateAll() {}
00325 };
00326
00327 template<int StorageOrder, typename _LhsScalar, typename _RhsScalar, int MaxRows, int MaxCols, int MaxDepth>
00328 class gemm_blocking_space<StorageOrder,_LhsScalar,_RhsScalar,MaxRows, MaxCols, MaxDepth, false>
00329 : public level3_blocking<
00330 typename conditional<StorageOrder==RowMajor,_RhsScalar,_LhsScalar>::type,
00331 typename conditional<StorageOrder==RowMajor,_LhsScalar,_RhsScalar>::type>
00332 {
00333 enum {
00334 Transpose = StorageOrder==RowMajor
00335 };
00336 typedef typename conditional<Transpose,_RhsScalar,_LhsScalar>::type LhsScalar;
00337 typedef typename conditional<Transpose,_LhsScalar,_RhsScalar>::type RhsScalar;
00338 typedef gebp_traits<LhsScalar,RhsScalar> Traits;
00339
00340 DenseIndex m_sizeA;
00341 DenseIndex m_sizeB;
00342 DenseIndex m_sizeW;
00343
00344 public:
00345
00346 gemm_blocking_space(DenseIndex rows, DenseIndex cols, DenseIndex depth)
00347 {
00348 this->m_mc = Transpose ? cols : rows;
00349 this->m_nc = Transpose ? rows : cols;
00350 this->m_kc = depth;
00351
00352 computeProductBlockingSizes<LhsScalar,RhsScalar>(this->m_kc, this->m_mc, this->m_nc);
00353 m_sizeA = this->m_mc * this->m_kc;
00354 m_sizeB = this->m_kc * this->m_nc;
00355 m_sizeW = this->m_kc*Traits::WorkSpaceFactor;
00356 }
00357
00358 void allocateA()
00359 {
00360 if(this->m_blockA==0)
00361 this->m_blockA = aligned_new<LhsScalar>(m_sizeA);
00362 }
00363
00364 void allocateB()
00365 {
00366 if(this->m_blockB==0)
00367 this->m_blockB = aligned_new<RhsScalar>(m_sizeB);
00368 }
00369
00370 void allocateW()
00371 {
00372 if(this->m_blockW==0)
00373 this->m_blockW = aligned_new<RhsScalar>(m_sizeW);
00374 }
00375
00376 void allocateAll()
00377 {
00378 allocateA();
00379 allocateB();
00380 allocateW();
00381 }
00382
00383 ~gemm_blocking_space()
00384 {
00385 aligned_delete(this->m_blockA, m_sizeA);
00386 aligned_delete(this->m_blockB, m_sizeB);
00387 aligned_delete(this->m_blockW, m_sizeW);
00388 }
00389 };
00390
00391 }
00392
00393 template<typename Lhs, typename Rhs>
00394 class GeneralProduct<Lhs, Rhs, GemmProduct>
00395 : public ProductBase<GeneralProduct<Lhs,Rhs,GemmProduct>, Lhs, Rhs>
00396 {
00397 enum {
00398 MaxDepthAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(Lhs::MaxColsAtCompileTime,Rhs::MaxRowsAtCompileTime)
00399 };
00400 public:
00401 EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct)
00402
00403 typedef typename Lhs::Scalar LhsScalar;
00404 typedef typename Rhs::Scalar RhsScalar;
00405 typedef Scalar ResScalar;
00406
00407 GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs)
00408 {
00409 typedef internal::scalar_product_op<LhsScalar,RhsScalar> BinOp;
00410 EIGEN_CHECK_BINARY_COMPATIBILIY(BinOp,LhsScalar,RhsScalar);
00411 }
00412
00413 template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
00414 {
00415 eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
00416
00417 typename internal::add_const_on_value_type<ActualLhsType>::type lhs = LhsBlasTraits::extract(m_lhs);
00418 typename internal::add_const_on_value_type<ActualRhsType>::type rhs = RhsBlasTraits::extract(m_rhs);
00419
00420 Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs)
00421 * RhsBlasTraits::extractScalarFactor(m_rhs);
00422
00423 typedef internal::gemm_blocking_space<(Dest::Flags&RowMajorBit) ? RowMajor : ColMajor,LhsScalar,RhsScalar,
00424 Dest::MaxRowsAtCompileTime,Dest::MaxColsAtCompileTime,MaxDepthAtCompileTime> BlockingType;
00425
00426 typedef internal::gemm_functor<
00427 Scalar, Index,
00428 internal::general_matrix_matrix_product<
00429 Index,
00430 LhsScalar, (_ActualLhsType::Flags&RowMajorBit) ? RowMajor : ColMajor, bool(LhsBlasTraits::NeedToConjugate),
00431 RhsScalar, (_ActualRhsType::Flags&RowMajorBit) ? RowMajor : ColMajor, bool(RhsBlasTraits::NeedToConjugate),
00432 (Dest::Flags&RowMajorBit) ? RowMajor : ColMajor>,
00433 _ActualLhsType, _ActualRhsType, Dest, BlockingType> GemmFunctor;
00434
00435 BlockingType blocking(dst.rows(), dst.cols(), lhs.cols());
00436
00437 internal::parallelize_gemm<(Dest::MaxRowsAtCompileTime>32 || Dest::MaxRowsAtCompileTime==Dynamic)>(GemmFunctor(lhs, rhs, dst, actualAlpha, blocking), this->rows(), this->cols(), Dest::Flags&RowMajorBit);
00438 }
00439 };
00440
00441 }
00442
00443 #endif // EIGEN_GENERAL_MATRIX_MATRIX_H