From 067a5a98c8979bea06ce2a9caef92a37b78b48eb Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sat, 29 Sep 2012 02:02:12 +0800 Subject: Extend MatrixPowerTriangularAtomic for future implementation for triangular matrix power. --- .../Eigen/src/MatrixFunctions/MatrixPower.h | 24 +- .../Eigen/src/MatrixFunctions/MatrixPowerBase.h | 249 ++++++++++++++++----- 2 files changed, 205 insertions(+), 68 deletions(-) (limited to 'unsupported') diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index 64afdaab8..7beb209eb 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -74,11 +74,11 @@ class MatrixPower : public MatrixPowerBase,MatrixType> void compute(const Derived& b, ResultType& res, RealScalar p); private: + EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(MatrixPower) + typedef Matrix, RowsAtCompileTime, ColsAtCompileTime, Options,MaxRowsAtCompileTime,MaxColsAtCompileTime> ComplexMatrix; - MatrixType m_tmp1, m_tmp2; ComplexMatrix m_T, m_U, m_fT; - bool m_init; RealScalar modfAndInit(RealScalar, RealScalar*); @@ -98,9 +98,8 @@ class MatrixPower : public MatrixPowerBase,MatrixType> template template MatrixPower::MatrixPower(const MatrixExpression& A) : - Base(A), - m_init(false) -{ /* empty body */ } + Base(A, 0) +{ } template void MatrixPower::compute(MatrixType& res, RealScalar p) @@ -113,7 +112,7 @@ void MatrixPower::compute(MatrixType& res, RealScalar p) break; default: RealScalar intpart, x = modfAndInit(p, &intpart); - res = MatrixType::Identity(m_A.rows(), m_A.cols()); + res = m_Id; computeIntPower(res, intpart); computeFracPower(res, x); } @@ -139,22 +138,19 @@ void MatrixPower::compute(const Derived& b, ResultType& res, RealSca template typename MatrixPower::Base::RealScalar MatrixPower::modfAndInit(RealScalar x, RealScalar* intpart) { - static RealScalar maxAbsEival, minAbsEival; *intpart = std::floor(x); RealScalar res = x - *intpart; - if (!m_init && res) { + if (!m_conditionNumber && res) { const ComplexSchur schurOfA(m_A); m_T = schurOfA.matrixT(); m_U = schurOfA.matrixU(); - m_init = true; - + const RealArray absTdiag = m_T.diagonal().array().abs(); - maxAbsEival = absTdiag.maxCoeff(); - minAbsEival = absTdiag.minCoeff(); + m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff(); } - if (res > RealScalar(0.5) && res > (1-res) * std::pow(maxAbsEival/minAbsEival, res)) { + if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber, res)) { --res; ++*intpart; } @@ -195,7 +191,7 @@ template void MatrixPower::computeIntPower(const Derived& b, ResultType& res, RealScalar p) { if (b.cols() >= m_A.cols()) { - m_tmp2 = MatrixType::Identity(m_A.rows(), m_A.cols()); + m_tmp2 = m_Id; computeIntPower(m_tmp2, p); res.noalias() = m_tmp2 * b; } diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h index b5e8ed7ed..e83d4abfd 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h @@ -18,7 +18,7 @@ struct recompose_complex_schur { template static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U) - { res = U * (T.template triangularView() * U.adjoint()); } + { res.noalias() = U * (T.template triangularView() * U.adjoint()); } }; template<> @@ -26,7 +26,21 @@ struct recompose_complex_schur<0> { template static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U) - { res = (U * (T.template triangularView() * U.adjoint())).real(); } + { res.noalias() = (U * (T.template triangularView() * U.adjoint())).real(); } +}; + +template::IsComplex> +struct matrix_power_unwinder +{ + static inline Scalar run(const Scalar& eival, const Scalar& eival0, int unwindingNumber) + { return internal::atanh2(eival-eival0, eival+eival0) + Scalar(0, M_PI*unwindingNumber); } +}; + +template +struct matrix_power_unwinder +{ + static inline Scalar run(Scalar eival, Scalar eival0, int) + { return internal::atanh2(eival-eival0, eival+eival0); } }; template @@ -68,21 +82,21 @@ inline int matrix_power_get_pade_degree(double normIminusT) inline int matrix_power_get_pade_degree(long double normIminusT) { #if LDBL_MANT_DIG == 53 - const int maxPadeDegree = 7; + enum { maxPadeDegree = 7 }; const double maxNormForPade[] = { 1.884160592658218e-2L /* degree = 3 */ , 6.038881904059573e-2L, 1.239917516308172e-1L, 1.999045567181744e-1L, 2.789358995219730e-1L }; #elif LDBL_MANT_DIG <= 64 - const int maxPadeDegree = 8; + enum { maxPadeDegree = 8 }; const double maxNormForPade[] = { 6.3854693117491799460e-3L /* degree = 3 */ , 2.6394893435456973676e-2L, 6.4216043030404063729e-2L, 1.1701165502926694307e-1L, 1.7904284231268670284e-1L, 2.4471944416607995472e-1L }; #elif LDBL_MANT_DIG <= 106 - const int maxPadeDegree = 10; + enum { maxPadeDegree = 10 }; const double maxNormForPade[] = { 1.0007161601787493236741409687186e-4L /* degree = 3 */ , 1.0007161601787493236741409687186e-3L, 4.7069769360887572939882574746264e-3L, 1.3220386624169159689406653101695e-2L, 2.8063482381631737920612944054906e-2L, 4.9625993951953473052385361085058e-2L, 7.7367040706027886224557538328171e-2L, 1.1016843812851143391275867258512e-1L }; #else - const int maxPadeDegree = 10; + enum { maxPadeDegree = 10 }; const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ , 6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L, 9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L, @@ -97,6 +111,125 @@ inline int matrix_power_get_pade_degree(long double normIminusT) } } // namespace internal +#define MATRIX_POWER_TRIANGULAR_2x2_SPECIALIZATION(Mode) \ + template \ + class MatrixPowerTriangular2x2 \ + { \ + private: \ + enum { \ + RowsAtCompileTime = MatrixType::RowsAtCompileTime, \ + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime \ + }; \ + typedef typename MatrixType::Scalar Scalar; \ + typedef typename MatrixType::RealScalar RealScalar; \ + typedef Array ArrayType; \ + const MatrixType& m_T; \ + public: \ + explicit MatrixPowerTriangular2x2(const MatrixType& T) : m_T(T) { } \ + void compute(MatrixType& res, RealScalar p) const; \ + }; + +template +class MatrixPowerTriangular2x2; + +MATRIX_POWER_TRIANGULAR_2x2_SPECIALIZATION(Upper) +MATRIX_POWER_TRIANGULAR_2x2_SPECIALIZATION(Lower) +MATRIX_POWER_TRIANGULAR_2x2_SPECIALIZATION(UnitUpper) +MATRIX_POWER_TRIANGULAR_2x2_SPECIALIZATION(UnitLower) +MATRIX_POWER_TRIANGULAR_2x2_SPECIALIZATION(StrictlyUpper) +MATRIX_POWER_TRIANGULAR_2x2_SPECIALIZATION(StrictlyLower) + +template +void MatrixPowerTriangular2x2::compute(MatrixType& res, RealScalar p) const +{ + using std::abs; + using std::pow; + + ArrayType logTdiag = m_T.diagonal().array().log(); + res.coeffRef(0,0) = pow(m_T.coeff(0,0), p); + + for (int i=1; i < m_T.cols(); ++i) { + res.coeffRef(i,i) = pow(m_T.coeff(i,i), p); + if (m_T.coeff(i-1,i-1) == m_T.coeff(i,i)) { + res.coeffRef(i-1,i) = p * pow(m_T.coeff(i-1,i), p-1); + } + else if (2*abs(m_T.coeff(i-1,i-1)) < abs(m_T.coeff(i,i)) || 2*abs(m_T.coeff(i,i)) < abs(m_T.coeff(i-1,i-1))) { + res.coeffRef(i-1,i) = m_T.coeff(i-1,i) * (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_T.coeff(i,i)-m_T.coeff(i-1,i-1)); + } + else { + int unwindingNumber = std::ceil((internal::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI)); + Scalar w = internal::matrix_power_unwinder::run(m_T.coeff(i,i), m_T.coeff(i-1,i-1), unwindingNumber); + res.coeffRef(i-1,i) = m_T.coeff(i-1,i) * RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) * + std::sinh(p * w) / (m_T.coeff(i,i) - m_T.coeff(i-1,i-1)); + } + } +} + +template +void MatrixPowerTriangular2x2::compute(MatrixType& res, RealScalar p) const +{ + using std::abs; + using std::pow; + + ArrayType logTdiag = m_T.diagonal().array().log(); + res.coeffRef(0,0) = pow(m_T.coeff(0,0), p); + + for (int i=1; i < m_T.cols(); ++i) { + res.coeffRef(i,i) = pow(m_T.coeff(i,i), p); + if (m_T.coeff(i-1,i-1) == m_T.coeff(i,i)) { + res.coeffRef(i,i-1) = p * pow(m_T.coeff(i,i-1), p-1); + } + else if (2*abs(m_T.coeff(i-1,i-1)) < abs(m_T.coeff(i,i)) || 2*abs(m_T.coeff(i,i)) < abs(m_T.coeff(i-1,i-1))) { + res.coeffRef(i,i-1) = m_T.coeff(i,i-1) * (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_T.coeff(i,i)-m_T.coeff(i-1,i-1)); + } + else { + int unwindingNumber = std::ceil((internal::imag(logTdiag[i]-logTdiag[i-1]) - M_PI) / (2*M_PI)); + Scalar w = internal::matrix_power_unwinder::run(m_T.coeff(i,i), m_T.coeff(i-1,i-1), unwindingNumber); + res.coeffRef(i,i-1) = m_T.coeff(i,i-1) * RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) * + std::sinh(p * w) / (m_T.coeff(i,i) - m_T.coeff(i-1,i-1)); + } + } +} + +template +void MatrixPowerTriangular2x2::compute(MatrixType& res, RealScalar p) const +{ + for (int i=1; i < m_T.cols(); ++i) + res.coeffRef(i-1,i) = p * std::pow(m_T.coeff(i-1,i), p-1); +} + +template +void MatrixPowerTriangular2x2::compute(MatrixType& res, RealScalar p) const +{ + for (int i=1; i < m_T.cols(); ++i) { + res.coeffRef(i,i-1) = p * std::pow(m_T.coeff(i,i-1), p-1); + } +} + +template +void MatrixPowerTriangular2x2::compute(MatrixType& res, RealScalar p) const +{ + RealScalar diag = !p ? 1 : 0; + res.coeffRef(0,0) = diag; + + for (int i=1; i < m_T.cols(); ++i) { + res.coeffRef(i,i) = diag; + res.coeffRef(i-1,i) = p * std::pow(m_T.coeff(i-1,i), p-1); + } +} + +template +void MatrixPowerTriangular2x2::compute(MatrixType& res, RealScalar p) const +{ + RealScalar diag = !p ? 1 : 0; + res.coeffRef(0,0) = diag; + + for (int i=1; i < m_T.cols(); ++i) { + res.coeffRef(i,i) = diag; + res.coeffRef(i,i-1) = p * std::pow(m_T.coeff(i,i-1), p-1); + } +} + template class MatrixPowerTriangularAtomic { @@ -113,7 +246,6 @@ class MatrixPowerTriangularAtomic const MatrixType m_Id; void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p) const; - void compute2x2(MatrixType& res, RealScalar p) const; void computeBig(MatrixType& res, RealScalar p) const; public: @@ -125,7 +257,7 @@ template MatrixPowerTriangularAtomic::MatrixPowerTriangularAtomic(const MatrixType& T) : m_T(T), m_Id(MatrixType::Identity(T.rows(), T.cols())) -{ /* empty body */ } +{ } template void MatrixPowerTriangularAtomic::compute(MatrixType& res, RealScalar p) const @@ -137,7 +269,7 @@ void MatrixPowerTriangularAtomic::compute(MatrixType& res, Real res(0,0) = std::pow(m_T(0,0), p); break; case 2: - compute2x2(res, p); + MatrixPowerTriangular2x2(m_T).compute(res, p); break; default: computeBig(res, p); @@ -157,42 +289,16 @@ void MatrixPowerTriangularAtomic::computePade(int degree, const res += m_Id; } -template -void MatrixPowerTriangularAtomic::compute2x2(MatrixType& res, RealScalar p) const -{ - using std::abs; - using std::pow; - - ArrayType logTdiag = m_T.diagonal().array().log(); - res(0,0) = pow(m_T(0,0), p); - - for (int i=1; i < m_T.cols(); ++i) { - res(i,i) = pow(m_T(i,i), p); - if (m_T(i-1,i-1) == m_T(i,i)) { - res(i-1,i) = p * pow(m_T(i-1,i), p-1); - } - else if (2*abs(m_T(i-1,i-1)) < abs(m_T(i,i)) || 2*abs(m_T(i,i)) < abs(m_T(i-1,i-1))) { - res(i-1,i) = m_T(i-1,i) * (res(i,i)-res(i-1,i-1)) / (m_T(i,i)-m_T(i-1,i-1)); - } - else { - // computation in previous branch is inaccurate if abs(m_T(i,i)) \approx abs(m_T(i-1,i-1)) - int unwindingNumber = std::ceil(((logTdiag[i]-logTdiag[i-1]).imag() - M_PI) / (2*M_PI)); - Scalar w = internal::atanh2(m_T(i,i)-m_T(i-1,i-1), m_T(i,i)+m_T(i-1,i-1)) + Scalar(0, M_PI*unwindingNumber); - res(i-1,i) = m_T(i-1,i) * RealScalar(2) * std::exp(RealScalar(0.5) * p * (logTdiag[i]+logTdiag[i-1])) * - std::sinh(p * w) / (m_T(i,i) - m_T(i-1,i-1)); - } - } -} - template void MatrixPowerTriangularAtomic::computeBig(MatrixType& res, RealScalar p) const { - const int digits = std::numeric_limits::digits; + enum { digits = std::numeric_limits::digits }; const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision digits <= 53? 2.789358995219730e-1: // double precision digits <= 64? 2.4471944416607995472e-1L: // extended precision - digits <= 106? 1.1016843812851143391275867258512e-01: // double-double - 9.134603732914548552537150753385375e-02; // quadruple precision + digits <= 106? 1.1016843812851143391275867258512e-1L: // double-double + 9.134603732914548552537150753385375e-2L; // quadruple precision + const MatrixPowerTriangular2x2 atomic2x2(m_T); MatrixType IminusT, sqrtT, T=m_T; RealScalar normIminusT; int degree, degree2, numberOfSquareRoots=0, numberOfExtraSquareRoots=0; @@ -214,14 +320,14 @@ void MatrixPowerTriangularAtomic::computeBig(MatrixType& res, R computePade(degree, IminusT, res, p); for (; numberOfSquareRoots; --numberOfSquareRoots) { - compute2x2(res, std::ldexp(p,-numberOfSquareRoots)); + atomic2x2.compute(res, std::ldexp(p,-numberOfSquareRoots)); res *= res; } - compute2x2(res, p); + atomic2x2.compute(res, p); } #define EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(Derived) \ - typedef MatrixPowerBase, MatrixType> Base; \ + typedef MatrixPowerBase Base; \ using Base::RowsAtCompileTime; \ using Base::ColsAtCompileTime; \ using Base::Options; \ @@ -229,8 +335,14 @@ void MatrixPowerTriangularAtomic::computeBig(MatrixType& res, R using Base::MaxColsAtCompileTime; \ typedef typename Base::Scalar Scalar; \ typedef typename Base::RealScalar RealScalar; \ - typedef typename Base::RealArray RealArray; \ - using Base::m_A; + typedef typename Base::RealArray RealArray; + +#define EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(Derived) \ + using Base::m_A; \ + using Base::m_Id; \ + using Base::m_tmp1; \ + using Base::m_tmp2; \ + using Base::m_conditionNumber; #define EIGEN_MATRIX_POWER_PRODUCT_PUBLIC_INTERFACE(Derived) \ typedef MatrixPowerProductBase Base; \ @@ -277,34 +389,63 @@ class MatrixPowerBase typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; - explicit MatrixPowerBase(const MatrixType& A) - : m_A(A), m_del(false) { } + explicit MatrixPowerBase(const MatrixType& A, RealScalar cond); template - explicit MatrixPowerBase(const MatrixBase& A) - : m_A(*new MatrixType(A)), m_del(true) { } + explicit MatrixPowerBase(const MatrixBase& A, RealScalar cond); - ~MatrixPowerBase() - { if (m_del) delete &m_A; } + ~MatrixPowerBase(); - void compute(MatrixType& res, RealScalar p) - { static_cast(this)->compute(res,p); } + void compute(MatrixType& res, RealScalar p); template - void compute(const OtherDerived& b, ResultType& res, RealScalar p) - { static_cast(this)->compute(b,res,p); } + void compute(const OtherDerived& b, ResultType& res, RealScalar p); Index rows() const { return m_A.rows(); } Index cols() const { return m_A.cols(); } protected: typedef Array RealArray; + const MatrixType& m_A; + const MatrixType m_Id; + MatrixType m_tmp1, m_tmp2; + RealScalar m_conditionNumber; private: const bool m_del; // whether to delete the pointer at destruction }; +template +MatrixPowerBase::MatrixPowerBase(const MatrixType& A, RealScalar cond) : + m_A(A), + m_Id(MatrixType::Identity(A.rows(),A.cols())), + m_conditionNumber(cond), + m_del(false) +{ } + +template +template +MatrixPowerBase::MatrixPowerBase(const MatrixBase& A, RealScalar cond) : + m_A(*new MatrixType(A)), + m_Id(MatrixType::Identity(A.rows(),A.cols())), + m_conditionNumber(cond), + m_del(true) +{ } + +template +MatrixPowerBase::~MatrixPowerBase() +{ if (m_del) delete &m_A; } + +template +void MatrixPowerBase::compute(MatrixType& res, RealScalar p) +{ static_cast(this)->compute(res,p); } + +template +template +void MatrixPowerBase::compute(const OtherDerived& b, ResultType& res, RealScalar p) +{ static_cast(this)->compute(b,res,p); } + template class MatrixPowerProductBase : public MatrixBase { -- cgit v1.2.3