10 #ifndef EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H
11 #define EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H
25 template<
typename LhsScalar,
typename RhsScalar,
typename Index,
int mr,
int nr,
bool ConjLhs,
bool ConjRhs,
int UpLo>
29 template <
typename Index,
30 typename LhsScalar,
int LhsStorageOrder,
bool ConjugateLhs,
31 typename RhsScalar,
int RhsStorageOrder,
bool ConjugateRhs,
32 int ResStorageOrder,
int UpLo,
int Version = Specialized>
33 struct general_matrix_matrix_triangular_product;
36 template <
typename Index,
typename LhsScalar,
int LhsStorageOrder,
bool ConjugateLhs,
37 typename RhsScalar,
int RhsStorageOrder,
bool ConjugateRhs,
int UpLo,
int Version>
38 struct general_matrix_matrix_triangular_product<Index,LhsScalar,LhsStorageOrder,ConjugateLhs,RhsScalar,RhsStorageOrder,ConjugateRhs,
RowMajor,UpLo,Version>
40 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
41 static EIGEN_STRONG_INLINE
void run(Index size, Index depth,
const LhsScalar* lhs, Index lhsStride,
42 const RhsScalar* rhs, Index rhsStride, ResScalar* res, Index resStride, ResScalar alpha)
44 general_matrix_matrix_triangular_product<Index,
48 ::run(size,depth,rhs,rhsStride,lhs,lhsStride,res,resStride,alpha);
52 template <
typename Index,
typename LhsScalar,
int LhsStorageOrder,
bool ConjugateLhs,
53 typename RhsScalar,
int RhsStorageOrder,
bool ConjugateRhs,
int UpLo,
int Version>
54 struct general_matrix_matrix_triangular_product<Index,LhsScalar,LhsStorageOrder,ConjugateLhs,RhsScalar,RhsStorageOrder,ConjugateRhs,
ColMajor,UpLo,Version>
56 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
57 static EIGEN_STRONG_INLINE
void run(Index size, Index depth,
const LhsScalar* _lhs, Index lhsStride,
58 const RhsScalar* _rhs, Index rhsStride, ResScalar* res, Index resStride, ResScalar alpha)
60 const_blas_data_mapper<LhsScalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
61 const_blas_data_mapper<RhsScalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
63 typedef gebp_traits<LhsScalar,RhsScalar> Traits;
68 computeProductBlockingSizes<LhsScalar,RhsScalar>(kc, mc, nc);
71 mc = (mc/Traits::nr)*Traits::nr;
73 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
74 std::size_t sizeB = sizeW + kc*size;
75 ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, kc*mc, 0);
76 ei_declare_aligned_stack_constructed_variable(RhsScalar, allocatedBlockB, sizeB, 0);
77 RhsScalar* blockB = allocatedBlockB + sizeW;
79 gemm_pack_lhs<LhsScalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
80 gemm_pack_rhs<RhsScalar, Index, Traits::nr, RhsStorageOrder> pack_rhs;
81 gebp_kernel <LhsScalar, RhsScalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp;
82 tribb_kernel<LhsScalar, RhsScalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs, UpLo> sybb;
84 for(Index k2=0; k2<depth; k2+=kc)
86 const Index actual_kc = (std::min)(k2+kc,depth)-k2;
89 pack_rhs(blockB, &rhs(k2,0), rhsStride, actual_kc, size);
91 for(Index i2=0; i2<size; i2+=mc)
93 const Index actual_mc = (std::min)(i2+mc,size)-i2;
95 pack_lhs(blockA, &lhs(i2, k2), lhsStride, actual_kc, actual_mc);
102 gebp(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, (std::min)(size,i2), alpha,
103 -1, -1, 0, 0, allocatedBlockB);
105 sybb(res+resStride*i2 + i2, resStride, blockA, blockB + actual_kc*i2, actual_mc, actual_kc, alpha, allocatedBlockB);
109 Index j2 = i2+actual_mc;
110 gebp(res+resStride*j2+i2, resStride, blockA, blockB+actual_kc*j2, actual_mc, actual_kc, (std::max)(Index(0), size-j2), alpha,
111 -1, -1, 0, 0, allocatedBlockB);
127 template<
typename LhsScalar,
typename RhsScalar,
typename Index,
int mr,
int nr,
bool ConjLhs,
bool ConjRhs,
int UpLo>
130 typedef gebp_traits<LhsScalar,RhsScalar,ConjLhs,ConjRhs> Traits;
131 typedef typename Traits::ResScalar ResScalar;
134 BlockSize = EIGEN_PLAIN_ENUM_MAX(mr,nr)
136 void operator()(ResScalar* res, Index resStride,
const LhsScalar* blockA,
const RhsScalar* blockB, Index size, Index depth, ResScalar alpha, RhsScalar* workspace)
138 gebp_kernel<LhsScalar, RhsScalar, Index, mr, nr, ConjLhs, ConjRhs> gebp_kernel;
139 Matrix<ResScalar,BlockSize,BlockSize,ColMajor> buffer;
143 for (Index j=0; j<size; j+=BlockSize)
145 Index actualBlockSize = std::min<Index>(BlockSize,size - j);
146 const RhsScalar* actual_b = blockB+j*depth;
149 gebp_kernel(res+j*resStride, resStride, blockA, actual_b, j, depth, actualBlockSize, alpha,
150 -1, -1, 0, 0, workspace);
157 gebp_kernel(buffer.data(), BlockSize, blockA+depth*i, actual_b, actualBlockSize, depth, actualBlockSize, alpha,
158 -1, -1, 0, 0, workspace);
160 for(Index j1=0; j1<actualBlockSize; ++j1)
162 ResScalar* r = res + (j+j1)*resStride + i;
163 for(Index i1=UpLo==
Lower ? j1 : 0;
164 UpLo==
Lower ? i1<actualBlockSize : i1<=j1; ++i1)
165 r[i1] += buffer(i1,j1);
171 Index i = j+actualBlockSize;
172 gebp_kernel(res+j*resStride+i, resStride, blockA+depth*i, actual_b, size-i, depth, actualBlockSize, alpha,
173 -1, -1, 0, 0, workspace);
183 template<
typename MatrixType,
unsigned int UpLo>
184 template<
typename ProductDerived,
typename _Lhs,
typename _Rhs>
185 TriangularView<MatrixType,UpLo>& TriangularView<MatrixType,UpLo>::assignProduct(
const ProductBase<ProductDerived, _Lhs,_Rhs>& prod,
const Scalar& alpha)
187 typedef typename internal::remove_all<typename ProductDerived::LhsNested>::type Lhs;
188 typedef internal::blas_traits<Lhs> LhsBlasTraits;
189 typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhs;
190 typedef typename internal::remove_all<ActualLhs>::type _ActualLhs;
191 typename internal::add_const_on_value_type<ActualLhs>::type actualLhs = LhsBlasTraits::extract(prod.lhs());
193 typedef typename internal::remove_all<typename ProductDerived::RhsNested>::type Rhs;
194 typedef internal::blas_traits<Rhs> RhsBlasTraits;
195 typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhs;
196 typedef typename internal::remove_all<ActualRhs>::type _ActualRhs;
197 typename internal::add_const_on_value_type<ActualRhs>::type actualRhs = RhsBlasTraits::extract(prod.rhs());
199 typename ProductDerived::Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs().derived()) * RhsBlasTraits::extractScalarFactor(prod.rhs().derived());
201 internal::general_matrix_matrix_triangular_product<Index,
205 ::run(m_matrix.cols(), actualLhs.cols(),
206 &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &actualRhs.coeffRef(0,0), actualRhs.outerStride(),
207 const_cast<Scalar*
>(m_matrix.data()), m_matrix.outerStride(), actualAlpha);
214 #endif // EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H