GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-04 12:15:11 +08:00
Code Issues Packages Projects Releases Wiki Activity

Make MatrixPowerTriangularAtomic::computePade static because it should be.

Browse Source
This commit is contained in:
Chen-Pang He 2012-10-07 02:25:00 +08:00
parent a5d348e30a
commit 0017d8c58f
2 changed files with 13 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
unsupported/Eigen/src/MatrixFunctions/MatrixPower.h
View File

@ -102,7 +102,7 @@ void MatrixPowerTriangular<MatrixType>::compute(MatrixType& res, RealScalar p)
break; break;
default: default:
RealScalar intpart, x = modfAndInit(p, &intpart); RealScalar intpart, x = modfAndInit(p, &intpart);
res = m_Id; res = MatrixType::Identity(m_A.rows(), m_A.cols());
computeIntPower(res, intpart); computeIntPower(res, intpart);
computeFracPower(res, x); computeFracPower(res, x);
} }
@ -162,7 +162,7 @@ void MatrixPowerTriangular<MatrixType>::computeIntPower(ResultType& res, RealSca
{ {
RealScalar pp = std::abs(p); RealScalar pp = std::abs(p);
if (p<0) m_tmp1 = m_T.solve(m_Id); if (p<0) m_tmp1 = m_T.solve(MatrixType::Identity(m_A.rows(), m_A.cols()));
else m_tmp1 = m_T; else m_tmp1 = m_T;
while (pp >= 1) { while (pp >= 1) {
@ -178,7 +178,7 @@ template<typename Derived, typename ResultType>
void MatrixPowerTriangular<MatrixType>::computeIntPower(const Derived& b, ResultType& res, RealScalar p) void MatrixPowerTriangular<MatrixType>::computeIntPower(const Derived& b, ResultType& res, RealScalar p)
{ {
if (b.cols() >= m_A.cols()) { if (b.cols() >= m_A.cols()) {
m_tmp2 = m_Id; m_tmp2 = MatrixType::Identity(m_A.rows(), m_A.cols());
computeIntPower(m_tmp2, p); computeIntPower(m_tmp2, p);
res.noalias() = m_tmp2.template triangularView<Upper>() * b; res.noalias() = m_tmp2.template triangularView<Upper>() * b;
} }
@ -201,7 +201,7 @@ void MatrixPowerTriangular<MatrixType>::computeIntPower(const Derived& b, Result
return; return;
} }
else { else {
m_tmp1 = m_T.solve(m_Id); m_tmp1 = m_T.solve(MatrixType::Identity(m_A.rows(), m_A.cols()));
} }
while (b.cols()*(pp-applyings) > m_A.cols()*squarings) { while (b.cols()*(pp-applyings) > m_A.cols()*squarings) {
@ -330,7 +330,7 @@ void MatrixPower<MatrixType>::compute(MatrixType& res, RealScalar p)
break; break;
default: default:
RealScalar intpart, x = modfAndInit(p, &intpart); RealScalar intpart, x = modfAndInit(p, &intpart);
res = m_Id; res = MatrixType::Identity(m_A.rows(), m_A.cols());
computeIntPower(res, intpart); computeIntPower(res, intpart);
computeFracPower(res, x); computeFracPower(res, x);
} }
@ -409,7 +409,7 @@ template<typename Derived, typename ResultType>
void MatrixPower<MatrixType>::computeIntPower(const Derived& b, ResultType& res, RealScalar p) void MatrixPower<MatrixType>::computeIntPower(const Derived& b, ResultType& res, RealScalar p)
{ {
if (b.cols() >= m_A.cols()) { if (b.cols() >= m_A.cols()) {
m_tmp2 = m_Id; m_tmp2 = MatrixType::Identity(m_A.rows(), m_A.cols());
computeIntPower(m_tmp2, p); computeIntPower(m_tmp2, p);
res.noalias() = m_tmp2 * b; res.noalias() = m_tmp2 * b;
} }

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

@ -25,7 +25,6 @@ namespace Eigen {
#define EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(Derived) \ #define EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(Derived) \
using Base::m_A; \ using Base::m_A; \
using Base::m_Id; \
using Base::m_tmp1; \ using Base::m_tmp1; \
using Base::m_tmp2; \ using Base::m_tmp2; \
using Base::m_conditionNumber; using Base::m_conditionNumber;
@ -77,7 +76,6 @@ class MatrixPowerBase
explicit MatrixPowerBase(const MatrixType& A) : explicit MatrixPowerBase(const MatrixType& A) :
m_A(A), m_A(A),
m_Id(MatrixType::Identity(A.rows(), A.cols())),
m_conditionNumber(0) m_conditionNumber(0)
{ eigen_assert(A.rows() == A.cols()); } { eigen_assert(A.rows() == A.cols()); }
@ -100,7 +98,6 @@ class MatrixPowerBase
typedef Array<RealScalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> RealArray; typedef Array<RealScalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> RealArray;
const MatrixType& m_A; const MatrixType& m_A;
const MatrixType m_Id;
MatrixType m_tmp1, m_tmp2; MatrixType m_tmp1, m_tmp2;
RealScalar m_conditionNumber; RealScalar m_conditionNumber;
}; };
@ -286,9 +283,8 @@ class MatrixPowerTriangularAtomic
typedef Array<Scalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> ArrayType; typedef Array<Scalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> ArrayType;
const MatrixType& m_A; const MatrixType& m_A;
const MatrixType m_Id;
void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p) const; static void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p);
void compute2x2(MatrixType& res, RealScalar p) const; void compute2x2(MatrixType& res, RealScalar p) const;
void computeBig(MatrixType& res, RealScalar p) const; void computeBig(MatrixType& res, RealScalar p) const;
@ -299,8 +295,7 @@ class MatrixPowerTriangularAtomic
template<typename MatrixType> template<typename MatrixType>
MatrixPowerTriangularAtomic<MatrixType>::MatrixPowerTriangularAtomic(const MatrixType& T) : MatrixPowerTriangularAtomic<MatrixType>::MatrixPowerTriangularAtomic(const MatrixType& T) :
m_A(T), m_A(T)
m_Id(MatrixType::Identity(T.rows(), T.cols()))
{ eigen_assert(T.rows() == T.cols()); } { eigen_assert(T.rows() == T.cols()); }
template<typename MatrixType> template<typename MatrixType>
@ -321,16 +316,15 @@ void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScala
} }
template<typename MatrixType> template<typename MatrixType>
void MatrixPowerTriangularAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res, void MatrixPowerTriangularAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p)
RealScalar p) const
{ {
int i = degree<<1; int i = degree<<1;
res = (p-degree) / ((i-1)<<1) * IminusT; res = (p-degree) / ((i-1)<<1) * IminusT;
for (--i; i; --i) { for (--i; i; --i) {
res = (m_Id + res).template triangularView<Upper>().solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView<Upper>()
(p-(i>>1))/((i-1)<<1)) * IminusT).eval(); .solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : (p-(i>>1))/((i-1)<<1)) * IminusT).eval();
} }
res += m_Id; res += MatrixType::Identity(IminusT.rows(), IminusT.cols());
} }
template<typename MatrixType> template<typename MatrixType>
@ -374,7 +368,7 @@ void MatrixPowerTriangularAtomic<MatrixType>::computeBig(MatrixType& res, RealSc
bool hasExtraSquareRoot = false; bool hasExtraSquareRoot = false;
while (true) { while (true) {
IminusT = m_Id - T; IminusT = MatrixType::Identity(m_A.rows(), m_A.cols()) - T;
normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff(); normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff();
if (normIminusT < maxNormForPade) { if (normIminusT < maxNormForPade) {
degree = internal::matrix_power_get_pade_degree(normIminusT); degree = internal::matrix_power_get_pade_degree(normIminusT);