10 #ifndef EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H 11 #define EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H 15 template<
typename Scalar,
typename Index,
int StorageOrder,
int UpLo,
bool ConjLhs,
bool ConjRhs>
28 template<
typename LhsScalar,
typename RhsScalar,
typename Index,
int mr,
int nr,
bool ConjLhs,
bool ConjRhs,
int UpLo>
32 template <
typename Index,
33 typename LhsScalar,
int LhsStorageOrder,
bool ConjugateLhs,
34 typename RhsScalar,
int RhsStorageOrder,
bool ConjugateRhs,
35 int ResStorageOrder,
int UpLo,
int Version =
Specialized>
39 template <
typename Index,
typename LhsScalar,
int LhsStorageOrder,
bool ConjugateLhs,
40 typename RhsScalar,
int RhsStorageOrder,
bool ConjugateRhs,
int UpLo,
int Version>
45 const RhsScalar*
rhs, Index rhsStride, ResScalar* res, Index resStride,
const ResScalar& alpha)
51 ::run(size,depth,rhs,rhsStride,lhs,lhsStride,res,resStride,alpha);
55 template <
typename Index,
typename LhsScalar,
int LhsStorageOrder,
bool ConjugateLhs,
56 typename RhsScalar,
int RhsStorageOrder,
bool ConjugateRhs,
int UpLo,
int Version>
61 const RhsScalar* _rhs, Index rhsStride, ResScalar* res, Index resStride,
const ResScalar& alpha)
71 computeProductBlockingSizes<LhsScalar,RhsScalar>(kc, mc, nc);
74 mc = (mc/Traits::nr)*Traits::nr;
76 std::size_t sizeW = kc*Traits::WorkSpaceFactor;
77 std::size_t sizeB = sizeW + kc*size;
80 RhsScalar* blockB = allocatedBlockB + sizeW;
87 for(Index k2=0; k2<depth; k2+=kc)
89 const Index actual_kc = (std::min)(k2+kc,depth)-k2;
92 pack_rhs(blockB, &
rhs(k2,0), rhsStride, actual_kc, size);
94 for(Index i2=0; i2<size; i2+=mc)
96 const Index actual_mc = (std::min)(i2+mc,size)-i2;
98 pack_lhs(blockA, &lhs(i2, k2), lhsStride, actual_kc, actual_mc);
105 gebp(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, (std::min)(size,i2), alpha,
106 -1, -1, 0, 0, allocatedBlockB);
108 sybb(res+resStride*i2 + i2, resStride, blockA, blockB + actual_kc*i2, actual_mc, actual_kc, alpha, allocatedBlockB);
112 Index j2 = i2+actual_mc;
113 gebp(res+resStride*j2+i2, resStride, blockA, blockB+actual_kc*j2, actual_mc, actual_kc, (std::max)(Index(0), size-j2), alpha,
114 -1, -1, 0, 0, allocatedBlockB);
130 template<
typename LhsScalar,
typename RhsScalar,
typename Index,
int mr,
int nr,
bool ConjLhs,
bool ConjRhs,
int UpLo>
139 void operator()(ResScalar* res, Index resStride,
const LhsScalar* blockA,
const RhsScalar* blockB, Index size, Index depth,
const ResScalar& alpha, RhsScalar* workspace)
146 for (Index j=0; j<size; j+=BlockSize)
148 Index actualBlockSize = std::min<Index>(BlockSize,size - j);
149 const RhsScalar* actual_b = blockB+j*depth;
152 gebp_kernel(res+j*resStride, resStride, blockA, actual_b, j, depth, actualBlockSize, alpha,
153 -1, -1, 0, 0, workspace);
160 gebp_kernel(buffer.
data(), BlockSize, blockA+depth*i, actual_b, actualBlockSize, depth, actualBlockSize, alpha,
161 -1, -1, 0, 0, workspace);
163 for(Index j1=0; j1<actualBlockSize; ++j1)
165 ResScalar* r = res + (j+j1)*resStride + i;
166 for(Index i1=UpLo==
Lower ? j1 : 0;
167 UpLo==
Lower ? i1<actualBlockSize : i1<=j1; ++i1)
168 r[i1] += buffer(i1,j1);
174 Index i = j+actualBlockSize;
175 gebp_kernel(res+j*resStride+i, resStride, blockA+depth*i, actual_b, size-i, depth, actualBlockSize, alpha,
176 -1, -1, 0, 0, workspace);
186 template<
typename MatrixType,
typename ProductType,
int UpLo,
bool IsOuterProduct>
190 template<
typename MatrixType,
typename ProductType,
int UpLo>
193 static void run(MatrixType& mat,
const ProductType& prod,
const typename MatrixType::Scalar& alpha)
195 typedef typename MatrixType::Scalar Scalar;
196 typedef typename MatrixType::Index Index;
200 typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhs;
206 typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhs;
210 Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs().derived()) * RhsBlasTraits::extractScalarFactor(prod.rhs().derived());
214 UseLhsDirectly = _ActualLhs::InnerStrideAtCompileTime==1,
215 UseRhsDirectly = _ActualRhs::InnerStrideAtCompileTime==1
220 (UseLhsDirectly ?
const_cast<Scalar*
>(actualLhs.data()) : static_lhs.data()));
225 (UseRhsDirectly ?
const_cast<Scalar*
>(actualRhs.data()) : static_rhs.data()));
231 RhsBlasTraits::NeedToConjugate && NumTraits<Scalar>::IsComplex>
232 ::run(actualLhs.size(), mat.data(), mat.outerStride(), actualLhsPtr, actualRhsPtr, actualAlpha);
236 template<
typename MatrixType,
typename ProductType,
int UpLo>
239 static void run(MatrixType& mat,
const ProductType& prod,
const typename MatrixType::Scalar& alpha)
241 typedef typename MatrixType::Index Index;
245 typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhs;
251 typedef typename RhsBlasTraits::DirectLinearAccessType ActualRhs;
255 typename ProductType::Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs().derived()) * RhsBlasTraits::extractScalarFactor(prod.rhs().derived());
261 ::run(mat.cols(), actualLhs.cols(),
262 &actualLhs.coeffRef(0,0), actualLhs.outerStride(), &actualRhs.coeffRef(0,0), actualRhs.outerStride(),
263 mat.data(), mat.outerStride(), actualAlpha);
267 template<
typename MatrixType,
unsigned int UpLo>
268 template<
typename ProductDerived,
typename _Lhs,
typename _Rhs>
278 #endif // EIGEN_GENERAL_MATRIX_MATRIX_TRIANGULAR_H
#define EIGEN_STRONG_INLINE
Traits::ResScalar ResScalar
#define ei_declare_aligned_stack_constructed_variable(TYPE, NAME, SIZE, BUFFER)
A matrix or vector expression mapping an existing array of data.
static void run(MatrixType &mat, const ProductType &prod, const typename MatrixType::Scalar &alpha)
iterative scaling algorithm to equilibrate rows and column norms in matrices
void operator()(ResScalar *res, Index resStride, const LhsScalar *blockA, const RhsScalar *blockB, Index size, Index depth, const ResScalar &alpha, RhsScalar *workspace)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
const unsigned int RowMajorBit
static void run(MatrixType &mat, const ProductType &prod, const typename MatrixType::Scalar &alpha)
EIGEN_STRONG_INLINE TriangularView & assignProduct(const ProductBase< ProductDerived, Lhs, Rhs > &prod, const Scalar &alpha)
internal::traits< TriangularView >::Scalar Scalar
static EIGEN_STRONG_INLINE void run(Index size, Index depth, const LhsScalar *_lhs, Index lhsStride, const RhsScalar *_rhs, Index rhsStride, ResScalar *res, Index resStride, const ResScalar &alpha)
static EIGEN_STRONG_INLINE void run(Index size, Index depth, const LhsScalar *lhs, Index lhsStride, const RhsScalar *rhs, Index rhsStride, ResScalar *res, Index resStride, const ResScalar &alpha)
gebp_traits< LhsScalar, RhsScalar, ConjLhs, ConjRhs > Traits
Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > & setZero(Index size)
void rhs(const real_t *x, real_t *f)
scalar_product_traits< LhsScalar, RhsScalar >::ReturnType ResScalar
EIGEN_STRONG_INLINE const Scalar * data() const
scalar_product_traits< LhsScalar, RhsScalar >::ReturnType ResScalar
Base class for triangular part in a matrix.
The matrix class, also used for vectors and row-vectors.
#define EIGEN_PLAIN_ENUM_MAX(a, b)
scalar_product_traits< LhsScalar, RhsScalar >::ReturnType ResScalar