25 #ifndef EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H
26 #define EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H
40 template<
typename LhsScalar,
typename RhsScalar,
typename Index,
int mr,
int nr,
bool ConjLhs,
bool ConjRhs,
int UpLo>
44 template <
typename Index,
45 typename LhsScalar,
int LhsStorageOrder,
bool ConjugateLhs,
46 typename RhsScalar,
int RhsStorageOrder,
bool ConjugateRhs,
47 int ResStorageOrder,
int UpLo,
int Version =
Specialized>
51 template <
typename Index,
typename LhsScalar,
int LhsStorageOrder,
bool ConjugateLhs,
52 typename RhsScalar,
int RhsStorageOrder,
bool ConjugateRhs,
int UpLo,
int Version>
55 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
56 static EIGEN_STRONG_INLINE void run(Index size, Index depth,
const LhsScalar* lhs, Index lhsStride,
57 const RhsScalar* rhs, Index rhsStride, ResScalar* res, Index resStride, ResScalar alpha)
63 ::run(size,depth,rhs,rhsStride,lhs,lhsStride,res,resStride,alpha);
67 template <
typename Index,
typename LhsScalar,
int LhsStorageOrder,
bool ConjugateLhs,
68 typename RhsScalar,
int RhsStorageOrder,
bool ConjugateRhs,
int UpLo,
int Version>
71 typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
72 static EIGEN_STRONG_INLINE void run(Index size, Index depth,
const LhsScalar* _lhs, Index lhsStride,
73 const RhsScalar* _rhs, Index rhsStride, ResScalar* res, Index resStride, ResScalar alpha)
75 const_blas_data_mapper<LhsScalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
76 const_blas_data_mapper<RhsScalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
78 typedef gebp_traits<LhsScalar,RhsScalar> Traits;
83 computeProductBlockingSizes<LhsScalar,RhsScalar>(kc, mc, nc);
86 mc = (mc/Traits::nr)*Traits::nr;
88 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
89 std::size_t sizeB = sizeW + kc*size;
92 RhsScalar* blockB = allocatedBlockB + sizeW;
94 gemm_pack_lhs<LhsScalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
95 gemm_pack_rhs<RhsScalar, Index, Traits::nr, RhsStorageOrder> pack_rhs;
96 gebp_kernel <LhsScalar, RhsScalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp;
97 tribb_kernel<LhsScalar, RhsScalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs, UpLo> sybb;
99 for(Index k2=0; k2<depth; k2+=kc)
101 const Index actual_kc = (std::min)(k2+kc,depth)-k2;
104 pack_rhs(blockB, &rhs(k2,0), rhsStride, actual_kc, size);
106 for(Index i2=0; i2<size; i2+=mc)
108 const Index actual_mc = (std::min)(i2+mc,size)-i2;
110 pack_lhs(blockA, &lhs(i2, k2), lhsStride, actual_kc, actual_mc);
117 gebp(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, (std::min)(size,i2), alpha,
118 -1, -1, 0, 0, allocatedBlockB);
120 sybb(res+resStride*i2 + i2, resStride, blockA, blockB + actual_kc*i2, actual_mc, actual_kc, alpha, allocatedBlockB);
124 Index j2 = i2+actual_mc;
125 gebp(res+resStride*j2+i2, resStride, blockA, blockB+actual_kc*j2, actual_mc, actual_kc, (std::max)(Index(0), size-j2), alpha,
126 -1, -1, 0, 0, allocatedBlockB);
142 template<
typename LhsScalar,
typename RhsScalar,
typename Index,
int mr,
int nr,
bool ConjLhs,
bool ConjRhs,
int UpLo>
145 typedef gebp_traits<LhsScalar,RhsScalar,ConjLhs,ConjRhs> Traits;
146 typedef typename Traits::ResScalar ResScalar;
151 void operator()(ResScalar* res, Index resStride,
const LhsScalar* blockA,
const RhsScalar* blockB, Index size, Index depth, ResScalar alpha, RhsScalar* workspace)
153 gebp_kernel<LhsScalar, RhsScalar, Index, mr, nr, ConjLhs, ConjRhs> gebp_kernel;
154 Matrix<ResScalar,BlockSize,BlockSize,ColMajor> buffer;
158 for (Index j=0; j<size; j+=BlockSize)
160 Index actualBlockSize = std::min<Index>(BlockSize,size - j);
161 const RhsScalar* actual_b = blockB+j*depth;
164 gebp_kernel(res+j*resStride, resStride, blockA, actual_b, j, depth, actualBlockSize, alpha,
165 -1, -1, 0, 0, workspace);
172 gebp_kernel(buffer.data(), BlockSize, blockA+depth*i, actual_b, actualBlockSize, depth, actualBlockSize, alpha,
173 -1, -1, 0, 0, workspace);
175 for(Index j1=0; j1<actualBlockSize; ++j1)
177 ResScalar* r = res + (j+j1)*resStride + i;
178 for(Index i1=UpLo==
Lower ? j1 : 0;
179 UpLo==
Lower ? i1<actualBlockSize : i1<=j1; ++i1)
180 r[i1] += buffer(i1,j1);
186 Index i = j+actualBlockSize;
187 gebp_kernel(res+j*resStride+i, resStride, blockA+depth*i, actual_b, size-i, depth, actualBlockSize, alpha,
188 -1, -1, 0, 0, workspace);
198 template<
typename MatrixType,
unsigned int UpLo>
199 template<
typename ProductDerived,
typename _Lhs,
typename _Rhs>
202 typedef typename internal::remove_all<typename ProductDerived::LhsNested>::type Lhs;
203 typedef internal::blas_traits<Lhs> LhsBlasTraits;
204 typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhs;
205 typedef typename internal::remove_all<ActualLhs>::type _ActualLhs;
206 typename internal::add_const_on_value_type<ActualLhs>::type actualLhs = LhsBlasTraits::extract(prod.
lhs());
208 typedef typename internal::remove_all<typename ProductDerived::RhsNested>::type Rhs;
209 typedef internal::blas_traits<Rhs> RhsBlasTraits;
210 typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhs;
211 typedef typename internal::remove_all<ActualRhs>::type _ActualRhs;
212 typename internal::add_const_on_value_type<ActualRhs>::type actualRhs = RhsBlasTraits::extract(prod.
rhs());
214 typename ProductDerived::Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.
lhs().derived()) * RhsBlasTraits::extractScalarFactor(prod.
rhs().derived());
216 internal::general_matrix_matrix_triangular_product<
Index,
220 ::run(m_matrix.cols(), actualLhs.cols(),
221 &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &actualRhs.coeffRef(0,0), actualRhs.outerStride(),
222 const_cast<Scalar*
>(m_matrix.data()), m_matrix.outerStride(), actualAlpha);
229 #endif // EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H