diff options
author | 2012-09-29 02:02:12 +0800 | |
---|---|---|
committer | 2012-09-29 02:02:12 +0800 | |
commit | 067a5a98c8979bea06ce2a9caef92a37b78b48eb (patch) | |
tree | 2ffdcfc36f293e6119fa6280f543b01b9dd7ce38 /unsupported | |
parent | ed18d6f2adcadc521090cd392f22dcd715e1f95f (diff) |
Extend MatrixPowerTriangularAtomic for future implementation for triangular matrix power.
Diffstat (limited to 'unsupported')
-rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixPower.h | 24 | ||||
-rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h | 249 |
2 files changed, 205 insertions, 68 deletions
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<MatrixPower<MatrixType>,MatrixType> void compute(const Derived& b, ResultType& res, RealScalar p); private: + EIGEN_MATRIX_POWER_PROTECTED_MEMBERS(MatrixPower) + typedef Matrix<std::complex<RealScalar>, 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<MatrixPower<MatrixType>,MatrixType> template<typename MatrixType> template<typename MatrixExpression> MatrixPower<MatrixType>::MatrixPower(const MatrixExpression& A) : - Base(A), - m_init(false) -{ /* empty body */ } + Base(A, 0) +{ } template<typename MatrixType> void MatrixPower<MatrixType>::compute(MatrixType& res, RealScalar p) @@ -113,7 +112,7 @@ void MatrixPower<MatrixType>::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<MatrixType>::compute(const Derived& b, ResultType& res, RealSca template<typename MatrixType> typename MatrixPower<MatrixType>::Base::RealScalar MatrixPower<MatrixType>::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<MatrixType> 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<typename Derived, typename ResultType> void MatrixPower<MatrixType>::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<typename ResultType, typename MatrixType> static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U) - { res = U * (T.template triangularView<Upper>() * U.adjoint()); } + { res.noalias() = U * (T.template triangularView<Upper>() * U.adjoint()); } }; template<> @@ -26,7 +26,21 @@ struct recompose_complex_schur<0> { template<typename ResultType, typename MatrixType> static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U) - { res = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); } + { res.noalias() = (U * (T.template triangularView<Upper>() * U.adjoint())).real(); } +}; + +template<typename Scalar, int IsComplex=NumTraits<Scalar>::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<typename Scalar> +struct matrix_power_unwinder<Scalar,0> +{ + static inline Scalar run(Scalar eival, Scalar eival0, int) + { return internal::atanh2(eival-eival0, eival+eival0); } }; template<typename T> @@ -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<typename MatrixType> \ + class MatrixPowerTriangular2x2<MatrixType,Mode> \ + { \ + private: \ + enum { \ + RowsAtCompileTime = MatrixType::RowsAtCompileTime, \ + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime \ + }; \ + typedef typename MatrixType::Scalar Scalar; \ + typedef typename MatrixType::RealScalar RealScalar; \ + typedef Array<Scalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> ArrayType; \ + const MatrixType& m_T; \ + public: \ + explicit MatrixPowerTriangular2x2(const MatrixType& T) : m_T(T) { } \ + void compute(MatrixType& res, RealScalar p) const; \ + }; + +template<typename MatrixType, unsigned int Mode> +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<typename MatrixType> +void MatrixPowerTriangular2x2<MatrixType,Upper>::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<Scalar>::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<typename MatrixType> +void MatrixPowerTriangular2x2<MatrixType,Lower>::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<Scalar>::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<typename MatrixType> +void MatrixPowerTriangular2x2<MatrixType,UnitUpper>::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<typename MatrixType> +void MatrixPowerTriangular2x2<MatrixType,UnitLower>::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<typename MatrixType> +void MatrixPowerTriangular2x2<MatrixType,StrictlyUpper>::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<typename MatrixType> +void MatrixPowerTriangular2x2<MatrixType,StrictlyLower>::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<typename MatrixType, unsigned int Mode=Upper> 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<typename MatrixType, unsigned int Mode> MatrixPowerTriangularAtomic<MatrixType,Mode>::MatrixPowerTriangularAtomic(const MatrixType& T) : m_T(T), m_Id(MatrixType::Identity(T.rows(), T.cols())) -{ /* empty body */ } +{ } template<typename MatrixType, unsigned int Mode> void MatrixPowerTriangularAtomic<MatrixType,Mode>::compute(MatrixType& res, RealScalar p) const @@ -137,7 +269,7 @@ void MatrixPowerTriangularAtomic<MatrixType,Mode>::compute(MatrixType& res, Real res(0,0) = std::pow(m_T(0,0), p); break; case 2: - compute2x2(res, p); + MatrixPowerTriangular2x2<MatrixType,Mode>(m_T).compute(res, p); break; default: computeBig(res, p); @@ -158,41 +290,15 @@ void MatrixPowerTriangularAtomic<MatrixType,Mode>::computePade(int degree, const } template<typename MatrixType, unsigned int Mode> -void MatrixPowerTriangularAtomic<MatrixType,Mode>::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<typename MatrixType, unsigned int Mode> void MatrixPowerTriangularAtomic<MatrixType,Mode>::computeBig(MatrixType& res, RealScalar p) const { - const int digits = std::numeric_limits<RealScalar>::digits; + enum { digits = std::numeric_limits<RealScalar>::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<MatrixType,Mode> atomic2x2(m_T); MatrixType IminusT, sqrtT, T=m_T; RealScalar normIminusT; int degree, degree2, numberOfSquareRoots=0, numberOfExtraSquareRoots=0; @@ -214,14 +320,14 @@ void MatrixPowerTriangularAtomic<MatrixType,Mode>::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<Derived<MatrixType>, MatrixType> Base; \ + typedef MatrixPowerBase<Derived, MatrixType> Base; \ using Base::RowsAtCompileTime; \ using Base::ColsAtCompileTime; \ using Base::Options; \ @@ -229,8 +335,14 @@ void MatrixPowerTriangularAtomic<MatrixType,Mode>::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<Derived, Lhs, Rhs> 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<typename OtherDerived> - explicit MatrixPowerBase(const MatrixBase<OtherDerived>& A) - : m_A(*new MatrixType(A)), m_del(true) { } + explicit MatrixPowerBase(const MatrixBase<OtherDerived>& A, RealScalar cond); - ~MatrixPowerBase() - { if (m_del) delete &m_A; } + ~MatrixPowerBase(); - void compute(MatrixType& res, RealScalar p) - { static_cast<Derived*>(this)->compute(res,p); } + void compute(MatrixType& res, RealScalar p); template<typename OtherDerived, typename ResultType> - void compute(const OtherDerived& b, ResultType& res, RealScalar p) - { static_cast<Derived*>(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<RealScalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> 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<typename Derived, typename MatrixType> +MatrixPowerBase<Derived,MatrixType>::MatrixPowerBase(const MatrixType& A, RealScalar cond) : + m_A(A), + m_Id(MatrixType::Identity(A.rows(),A.cols())), + m_conditionNumber(cond), + m_del(false) +{ } + +template<typename Derived, typename MatrixType> +template<typename OtherDerived> +MatrixPowerBase<Derived,MatrixType>::MatrixPowerBase(const MatrixBase<OtherDerived>& A, RealScalar cond) : + m_A(*new MatrixType(A)), + m_Id(MatrixType::Identity(A.rows(),A.cols())), + m_conditionNumber(cond), + m_del(true) +{ } + +template<typename Derived, typename MatrixType> +MatrixPowerBase<Derived,MatrixType>::~MatrixPowerBase() +{ if (m_del) delete &m_A; } + +template<typename Derived, typename MatrixType> +void MatrixPowerBase<Derived,MatrixType>::compute(MatrixType& res, RealScalar p) +{ static_cast<Derived*>(this)->compute(res,p); } + +template<typename Derived, typename MatrixType> +template<typename OtherDerived, typename ResultType> +void MatrixPowerBase<Derived,MatrixType>::compute(const OtherDerived& b, ResultType& res, RealScalar p) +{ static_cast<Derived*>(this)->compute(b,res,p); } + template<typename Derived, typename Lhs, typename Rhs> class MatrixPowerProductBase : public MatrixBase<Derived> { |