diff options
Diffstat (limited to 'unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h')
| -rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h | 359 |
1 files changed, 359 insertions, 0 deletions
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h new file mode 100644 index 000000000..9616659ca --- /dev/null +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h @@ -0,0 +1,359 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Chen-Pang He <jdh8@ms63.hinet.net> +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_MATRIX_POWER_BASE +#define EIGEN_MATRIX_POWER_BASE + +namespace Eigen { + +namespace internal { +template<int IsComplex> +struct recompose_complex_schur +{ + template<typename ResultType, typename MatrixType> + static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U) + { res.noalias() = U * (T.template triangularView<Upper>() * U.adjoint()); } +}; + +template<> +struct recompose_complex_schur<0> +{ + template<typename ResultType, typename MatrixType> + static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U) + { 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> +inline int binary_powering_cost(T p, int* squarings) +{ + int applyings=0, tmp; + + frexp(p, squarings); + --*squarings; + + while (std::frexp(p, &tmp), tmp > 0) { + p -= std::ldexp(static_cast<T>(0.5), tmp); + ++applyings; + } + return applyings; +} + +inline int matrix_power_get_pade_degree(float normIminusT) +{ + const float maxNormForPade[] = { 2.8064004e-1f /* degree = 3 */ , 4.3386528e-1f }; + int degree = 3; + for (; degree <= 4; ++degree) + if (normIminusT <= maxNormForPade[degree - 3]) + break; + return degree; +} + +inline int matrix_power_get_pade_degree(double normIminusT) +{ + const double maxNormForPade[] = { 1.884160592658218e-2 /* degree = 3 */ , 6.038881904059573e-2, 1.239917516308172e-1, + 1.999045567181744e-1, 2.789358995219730e-1 }; + int degree = 3; + for (; degree <= 7; ++degree) + if (normIminusT <= maxNormForPade[degree - 3]) + break; + return degree; +} + +inline int matrix_power_get_pade_degree(long double normIminusT) +{ +#if LDBL_MANT_DIG == 53 + const int 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; + 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; + 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; + const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ , + 6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L, + 9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L, + 3.908166513900489428442993794761185e-2L, 6.266780814639442865832535460550138e-2L, + 9.134603732914548552537150753385375e-2L }; +#endif + int degree = 3; + for (; degree <= maxPadeDegree; ++degree) + if (normIminusT <= maxNormForPade[degree - 3]) + break; + return degree; +} +} // namespace internal + +template<typename MatrixType> +class MatrixPowerTriangularAtomic +{ + 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; + 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: + explicit MatrixPowerTriangularAtomic(const MatrixType& T); + void compute(MatrixType& res, RealScalar p) const; +}; + +template<typename MatrixType> +MatrixPowerTriangularAtomic<MatrixType>::MatrixPowerTriangularAtomic(const MatrixType& T) : + m_T(T), + m_Id(MatrixType::Identity(T.rows(), T.cols())) +{ eigen_assert(T.rows() == T.cols()); } + +template<typename MatrixType> +void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const +{ + switch (m_T.rows()) { + case 0: + break; + case 1: + res(0,0) = std::pow(m_T(0,0), p); + break; + case 2: + compute2x2(res, p); + break; + default: + computeBig(res, p); + } +} + +template<typename MatrixType> +void MatrixPowerTriangularAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res, + RealScalar p) const +{ + int i = degree<<1; + res = (p-degree) / ((i-1)<<1) * IminusT; + for (--i; i; --i) { + res = (m_Id + res).template triangularView<Upper>().solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : + (p-(i>>1))/((i-1)<<1)) * IminusT).eval(); + } + res += m_Id; +} + +template<typename MatrixType> +void MatrixPowerTriangularAtomic<MatrixType>::compute2x2(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 MatrixPowerTriangularAtomic<MatrixType>::computeBig(MatrixType& res, RealScalar p) const +{ + const int 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-1L: // double-double + 9.134603732914548552537150753385375e-2L; // quadruple precision + MatrixType IminusT, sqrtT, T=m_T; + RealScalar normIminusT; + int degree, degree2, numberOfSquareRoots=0; + bool hasExtraSquareRoot=false; + + while (true) { + IminusT = MatrixType::Identity(m_T.rows(), m_T.cols()) - T; + normIminusT = IminusT.cwiseAbs().colwise().sum().maxCoeff(); + if (normIminusT < maxNormForPade) { + degree = internal::matrix_power_get_pade_degree(normIminusT); + degree2 = internal::matrix_power_get_pade_degree(normIminusT/2); + if (degree - degree2 <= 1 || hasExtraSquareRoot) + break; + hasExtraSquareRoot = true; + } + MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT); + T = sqrtT; + ++numberOfSquareRoots; + } + computePade(degree, IminusT, res, p); + + for (; numberOfSquareRoots; --numberOfSquareRoots) { + compute2x2(res, std::ldexp(p,-numberOfSquareRoots)); + res *= res; + } + compute2x2(res, p); +} + +#define EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(Derived) \ + typedef MatrixPowerBase<Derived, MatrixType> Base; \ + using Base::RowsAtCompileTime; \ + using Base::ColsAtCompileTime; \ + using Base::Options; \ + using Base::MaxRowsAtCompileTime; \ + using Base::MaxColsAtCompileTime; \ + typedef typename Base::Scalar Scalar; \ + typedef typename Base::RealScalar RealScalar; \ + 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; \ + EIGEN_DENSE_PUBLIC_INTERFACE(Derived) + +namespace internal { +template<typename Derived, typename _Lhs, typename _Rhs> +struct traits<MatrixPowerProductBase<Derived,_Lhs,_Rhs> > +{ + typedef MatrixXpr XprKind; + typedef typename remove_all<_Lhs>::type Lhs; + typedef typename remove_all<_Rhs>::type Rhs; + typedef typename remove_all<Derived>::type PlainObject; + typedef typename scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar; + typedef typename promote_storage_type<typename traits<Lhs>::StorageKind, + typename traits<Rhs>::StorageKind>::ret StorageKind; + typedef typename promote_index_type<typename traits<Lhs>::Index, + typename traits<Rhs>::Index>::type Index; + + enum { + RowsAtCompileTime = traits<Lhs>::RowsAtCompileTime, + ColsAtCompileTime = traits<Rhs>::ColsAtCompileTime, + MaxRowsAtCompileTime = traits<Lhs>::MaxRowsAtCompileTime, + MaxColsAtCompileTime = traits<Rhs>::MaxColsAtCompileTime, + Flags = (MaxRowsAtCompileTime==1 ? RowMajorBit : 0) + | EvalBeforeNestingBit | EvalBeforeAssigningBit | NestByRefBit, + CoeffReadCost = 0 + }; +}; +} // namespace internal + +template<typename Derived, typename MatrixType> +class MatrixPowerBase +{ + public: + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + Options = MatrixType::Options, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + + explicit MatrixPowerBase(const MatrixType& A, RealScalar cond); + + void compute(MatrixType& res, RealScalar p); + + template<typename OtherDerived, typename ResultType> + 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; +}; + +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) +{ eigen_assert(A.rows() == A.cols()); } + +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> +{ + public: + typedef MatrixBase<Derived> Base; + EIGEN_DENSE_PUBLIC_INTERFACE(MatrixPowerProductBase) + + inline Index rows() const { return derived().rows(); } + inline Index cols() const { return derived().cols(); } + + template<typename ResultType> + inline void evalTo(ResultType& res) const { derived().evalTo(res); } +}; + +template<typename Derived> +template<typename ProductDerived, typename Lhs, typename Rhs> +Derived& MatrixBase<Derived>::lazyAssign(const MatrixPowerProductBase<ProductDerived,Lhs,Rhs>& other) +{ + other.derived().evalTo(derived()); + return derived(); +} + +} // namespace Eigen + +#endif // EIGEN_MATRIX_POWER |