diff options
author | 2012-09-29 17:41:06 +0800 | |
---|---|---|
committer | 2012-09-29 17:41:06 +0800 | |
commit | 5814a5f1a00122f197ee017a62599ed8f1108e2a (patch) | |
tree | 881e393879cfce33146b9c4b004d4b965f2bb5fe /unsupported | |
parent | 067a5a98c8979bea06ce2a9caef92a37b78b48eb (diff) |
Abort the extension. MatrixSquareRootTriangular only takes upper triangular matrices.
Diffstat (limited to 'unsupported')
-rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h | 216 |
1 files changed, 51 insertions, 165 deletions
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h index e83d4abfd..9616659ca 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h @@ -82,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 - enum { maxPadeDegree = 7 }; + 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 - enum { maxPadeDegree = 8 }; + 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 - enum { maxPadeDegree = 10 }; + 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 - enum { maxPadeDegree = 10 }; + const int maxPadeDegree = 10; const double maxNormForPade[] = { 5.524506147036624377378713555116378e-5L /* degree = 3 */ , 6.640600568157479679823602193345995e-4L, 3.227716520106894279249709728084626e-3L, 9.619593944683432960546978734646284e-3L, 2.134595382433742403911124458161147e-2L, @@ -111,126 +111,7 @@ 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 { private: @@ -246,6 +127,7 @@ 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: @@ -253,14 +135,14 @@ class MatrixPowerTriangularAtomic void compute(MatrixType& res, RealScalar p) const; }; -template<typename MatrixType, unsigned int Mode> -MatrixPowerTriangularAtomic<MatrixType,Mode>::MatrixPowerTriangularAtomic(const MatrixType& T) : +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, unsigned int Mode> -void MatrixPowerTriangularAtomic<MatrixType,Mode>::compute(MatrixType& res, RealScalar p) const +template<typename MatrixType> +void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const { switch (m_T.rows()) { case 0: @@ -269,39 +151,65 @@ void MatrixPowerTriangularAtomic<MatrixType,Mode>::compute(MatrixType& res, Real res(0,0) = std::pow(m_T(0,0), p); break; case 2: - MatrixPowerTriangular2x2<MatrixType,Mode>(m_T).compute(res, p); + compute2x2(res, p); break; default: computeBig(res, p); } } -template<typename MatrixType, unsigned int Mode> -void MatrixPowerTriangularAtomic<MatrixType,Mode>::computePade(int degree, const MatrixType& IminusT, MatrixType& res, +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<Mode>().solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : + 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, unsigned int Mode> -void MatrixPowerTriangularAtomic<MatrixType,Mode>::computeBig(MatrixType& res, RealScalar p) const +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 { - enum { digits = std::numeric_limits<RealScalar>::digits }; + 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 - const MatrixPowerTriangular2x2<MatrixType,Mode> atomic2x2(m_T); MatrixType IminusT, sqrtT, T=m_T; RealScalar normIminusT; - int degree, degree2, numberOfSquareRoots=0, numberOfExtraSquareRoots=0; + int degree, degree2, numberOfSquareRoots=0; + bool hasExtraSquareRoot=false; while (true) { IminusT = MatrixType::Identity(m_T.rows(), m_T.cols()) - T; @@ -309,9 +217,9 @@ void MatrixPowerTriangularAtomic<MatrixType,Mode>::computeBig(MatrixType& res, R if (normIminusT < maxNormForPade) { degree = internal::matrix_power_get_pade_degree(normIminusT); degree2 = internal::matrix_power_get_pade_degree(normIminusT/2); - if (degree - degree2 <= 1 || numberOfExtraSquareRoots) + if (degree - degree2 <= 1 || hasExtraSquareRoot) break; - ++numberOfExtraSquareRoots; + hasExtraSquareRoot = true; } MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT); T = sqrtT; @@ -320,10 +228,10 @@ void MatrixPowerTriangularAtomic<MatrixType,Mode>::computeBig(MatrixType& res, R computePade(degree, IminusT, res, p); for (; numberOfSquareRoots; --numberOfSquareRoots) { - atomic2x2.compute(res, std::ldexp(p,-numberOfSquareRoots)); + compute2x2(res, std::ldexp(p,-numberOfSquareRoots)); res *= res; } - atomic2x2.compute(res, p); + compute2x2(res, p); } #define EIGEN_MATRIX_POWER_PUBLIC_INTERFACE(Derived) \ @@ -391,11 +299,6 @@ class MatrixPowerBase explicit MatrixPowerBase(const MatrixType& A, RealScalar cond); - template<typename OtherDerived> - explicit MatrixPowerBase(const MatrixBase<OtherDerived>& A, RealScalar cond); - - ~MatrixPowerBase(); - void compute(MatrixType& res, RealScalar p); template<typename OtherDerived, typename ResultType> @@ -411,31 +314,14 @@ class MatrixPowerBase 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; } + m_conditionNumber(cond) +{ eigen_assert(A.rows() == A.cols()); } template<typename Derived, typename MatrixType> void MatrixPowerBase<Derived,MatrixType>::compute(MatrixType& res, RealScalar p) |