From 5814a5f1a00122f197ee017a62599ed8f1108e2a Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sat, 29 Sep 2012 17:41:06 +0800 Subject: [PATCH] Abort the extension. MatrixSquareRootTriangular only takes upper triangular matrices. --- .../src/MatrixFunctions/MatrixPowerBase.h | 258 +++++------------- 1 file changed, 72 insertions(+), 186 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,36 +111,68 @@ 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 MatrixPowerTriangularAtomic +{ + private: + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime + }; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Array ArrayType; -template -class MatrixPowerTriangular2x2; + const MatrixType& m_T; + const MatrixType m_Id; -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) + 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 -void MatrixPowerTriangular2x2::compute(MatrixType& res, RealScalar p) const +MatrixPowerTriangularAtomic::MatrixPowerTriangularAtomic(const MatrixType& T) : + m_T(T), + m_Id(MatrixType::Identity(T.rows(), T.cols())) +{ eigen_assert(T.rows() == T.cols()); } + +template +void MatrixPowerTriangularAtomic::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 +void MatrixPowerTriangularAtomic::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().solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : + (p-(i>>1))/((i-1)<<1)) * IminusT).eval(); + } + res += m_Id; +} + +template +void MatrixPowerTriangularAtomic::compute2x2(MatrixType& res, RealScalar p) const { using std::abs; using std::pow; @@ -166,142 +198,18 @@ void MatrixPowerTriangular2x2::compute(MatrixType& res, RealSc } template -void MatrixPowerTriangular2x2::compute(MatrixType& res, RealScalar p) const +void MatrixPowerTriangularAtomic::computeBig(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 -{ - 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; - const MatrixType m_Id; - - void computePade(int degree, const MatrixType& IminusT, 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 -MatrixPowerTriangularAtomic::MatrixPowerTriangularAtomic(const MatrixType& T) : - m_T(T), - m_Id(MatrixType::Identity(T.rows(), T.cols())) -{ } - -template -void MatrixPowerTriangularAtomic::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: - MatrixPowerTriangular2x2(m_T).compute(res, p); - break; - default: - computeBig(res, p); - } -} - -template -void MatrixPowerTriangularAtomic::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().solve((i==1 ? -p : i&1 ? (-p-(i>>1))/(i<<1) : - (p-(i>>1))/((i-1)<<1)) * IminusT).eval(); - } - res += m_Id; -} - -template -void MatrixPowerTriangularAtomic::computeBig(MatrixType& res, RealScalar p) const -{ - enum { digits = std::numeric_limits::digits }; + const int 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-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; + 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::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(T).compute(sqrtT); T = sqrtT; @@ -320,10 +228,10 @@ void MatrixPowerTriangularAtomic::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 - explicit MatrixPowerBase(const MatrixBase& A, RealScalar cond); - - ~MatrixPowerBase(); - void compute(MatrixType& res, RealScalar p); template @@ -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 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; } + m_conditionNumber(cond) +{ eigen_assert(A.rows() == A.cols()); } template void MatrixPowerBase::compute(MatrixType& res, RealScalar p)