devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Make use of TRMM (speed up), and remove useless condition (the triangular don't need LU)

Browse Source
This commit is contained in:
Chen-Pang He
2012-10-02 23:04:23 +08:00
parent 21c2b4e327
commit 4cfde4590f
2 changed files with 5 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h
View File

@@ -24,7 +24,6 @@ namespace Eigen {
typedef typename Base::RealArray RealArray;
#define EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(Derived) \
using Base::m_OKforLU; \
using Base::m_A; \
using Base::m_Id; \
using Base::m_tmp1; \
@@ -99,7 +98,6 @@ class MatrixPowerBase
protected:
typedef Array<RealScalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> RealArray;
static const bool m_OKforLU = RowsAtCompileTime == Dynamic || RowsAtCompileTime > 4;
const MatrixType& m_A;
const MatrixType m_Id;
@@ -284,6 +282,7 @@ class MatrixPowerTriangularAtomic
};
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
typedef Array<Scalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> ArrayType;
const MatrixType& m_A;
@@ -343,7 +342,7 @@ void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(MatrixType& res, RealSc
ArrayType logTdiag = m_A.diagonal().array().log();
res.coeffRef(0,0) = pow(m_A.coeff(0,0), p);
for (int i=1; i < m_A.cols(); ++i) {
for (Index i=1; i < m_A.cols(); ++i) {
res.coeffRef(i,i) = pow(m_A.coeff(i,i), p);
if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) {
res.coeffRef(i-1,i) = p * pow(m_A.coeff(i-1,i), p-1);
@@ -392,7 +391,7 @@ void MatrixPowerTriangularAtomic<MatrixType>::computeBig(MatrixType& res, RealSc
for (; numberOfSquareRoots; --numberOfSquareRoots) {
compute2x2(res, std::ldexp(p,-numberOfSquareRoots));
res *= res;
res = res.template triangularView<Upper>() * res;
}
compute2x2(res, p);
}