diff options
author | Chen-Pang He <jdh8@ms63.hinet.net> | 2012-09-30 19:21:53 +0800 |
---|---|---|
committer | Chen-Pang He <jdh8@ms63.hinet.net> | 2012-09-30 19:21:53 +0800 |
commit | e92fe88159a762710651f1a8f9a428bc572a27df (patch) | |
tree | 36dc8b23012ab0f90fc29222ed78ebf7dc124b1e | |
parent | eb33d307af8cda6876b4eb334eaf258fbbfc8bff (diff) |
Add test for real MatrixPowerTriangular.
-rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixPower.h | 6 | ||||
-rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h | 69 | ||||
-rw-r--r-- | unsupported/test/matrix_power.cpp | 56 |
3 files changed, 88 insertions, 43 deletions
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index a618a536f..2cb1d95d9 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -20,7 +20,7 @@ namespace Eigen { * \tparam MatrixType type of the base, expected to be an instantiation * of the Matrix class template. * - * This class is capable of computing complex upper triangular matrices raised + * This class is capable of computing upper triangular matrices raised * to an arbitrary real power. */ template<typename MatrixType> @@ -37,7 +37,7 @@ class MatrixPowerTriangular : public MatrixPowerBase<MatrixPowerTriangular<Matri * The class stores a reference to A, so it should not be changed * (or destroyed) before evaluation. */ - explicit MatrixPowerTriangular(const MatrixType& A) : Base(A,0), m_T(Base::m_A) + explicit MatrixPowerTriangular(const MatrixType& A) : Base(A), m_T(Base::m_A) { } #ifdef EIGEN_PARSED_BY_DOXYGEN @@ -262,7 +262,7 @@ class MatrixPower : public MatrixPowerBase<MatrixPower<MatrixType>,MatrixType> * The class stores a reference to A, so it should not be changed * (or destroyed) before evaluation. */ - explicit MatrixPower(const MatrixType& A) : Base(A,0) + explicit MatrixPower(const MatrixType& A) : Base(A) { } #ifdef EIGEN_PARSED_BY_DOXYGEN diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h index b67939039..aa28b821b 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPowerBase.h @@ -75,10 +75,10 @@ class MatrixPowerBase typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; - explicit MatrixPowerBase(const MatrixType& A, RealScalar cond) : + explicit MatrixPowerBase(const MatrixType& A) : m_A(A), - m_Id(MatrixType::Identity(A.rows(),A.cols())), - m_conditionNumber(cond) + m_Id(MatrixType::Identity(A.rows(), A.cols())), + m_conditionNumber(0) { eigen_assert(A.rows() == A.cols()); } #ifndef EIGEN_PARSED_BY_DOXYGEN @@ -173,10 +173,8 @@ struct traits<MatrixPowerProduct<Derived,_Lhs,_Rhs> > }; }; -template<bool IsComplex> struct recompose_complex_schur; - -template<> -struct recompose_complex_schur<true> +template<int IsComplex> +struct recompose_complex_schur { template<typename ResultType, typename MatrixType> static inline void run(ResultType& res, const MatrixType& T, const MatrixType& U) @@ -184,14 +182,14 @@ struct recompose_complex_schur<true> }; template<> -struct recompose_complex_schur<false> +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> +template<typename Scalar, int IsComplex = NumTraits<Scalar>::IsComplex> struct matrix_power_unwinder { static inline Scalar run(const Scalar& eival, const Scalar& eival0, int unwindingNumber) @@ -274,11 +272,8 @@ inline int matrix_power_get_pade_degree(long double normIminusT) } // namespace internal -template<typename MatrixType, bool IsComplex=NumTraits<typename MatrixType::RealScalar>::IsComplex> -class MatrixPowerTriangularAtomic; - template<typename MatrixType> -class MatrixPowerTriangularAtomic<MatrixType,true> +class MatrixPowerTriangularAtomic { private: enum { @@ -289,7 +284,7 @@ class MatrixPowerTriangularAtomic<MatrixType,true> typedef typename MatrixType::RealScalar RealScalar; typedef Array<Scalar,RowsAtCompileTime,1,ColMajor,MaxRowsAtCompileTime> ArrayType; - const MatrixType& m_T; + const MatrixType& m_A; const MatrixType m_Id; void computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p) const; @@ -302,19 +297,19 @@ class MatrixPowerTriangularAtomic<MatrixType,true> }; template<typename MatrixType> -MatrixPowerTriangularAtomic<MatrixType,true>::MatrixPowerTriangularAtomic(const MatrixType& T) : - m_T(T), +MatrixPowerTriangularAtomic<MatrixType>::MatrixPowerTriangularAtomic(const MatrixType& T) : + m_A(T), m_Id(MatrixType::Identity(T.rows(), T.cols())) { eigen_assert(T.rows() == T.cols()); } template<typename MatrixType> -void MatrixPowerTriangularAtomic<MatrixType,true>::compute(MatrixType& res, RealScalar p) const +void MatrixPowerTriangularAtomic<MatrixType>::compute(MatrixType& res, RealScalar p) const { - switch (m_T.rows()) { + switch (m_A.rows()) { case 0: break; case 1: - res(0,0) = std::pow(m_T(0,0), p); + res(0,0) = std::pow(m_A(0,0), p); break; case 2: compute2x2(res, p); @@ -325,7 +320,7 @@ void MatrixPowerTriangularAtomic<MatrixType,true>::compute(MatrixType& res, Real } template<typename MatrixType> -void MatrixPowerTriangularAtomic<MatrixType,true>::computePade(int degree, const MatrixType& IminusT, MatrixType& res, +void MatrixPowerTriangularAtomic<MatrixType>::computePade(int degree, const MatrixType& IminusT, MatrixType& res, RealScalar p) const { int i = degree<<1; @@ -338,33 +333,33 @@ void MatrixPowerTriangularAtomic<MatrixType,true>::computePade(int degree, const } template<typename MatrixType> -void MatrixPowerTriangularAtomic<MatrixType,true>::compute2x2(MatrixType& res, RealScalar p) const +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); + ArrayType logTdiag = m_A.diagonal().array().log(); + res.coeffRef(0,0) = pow(m_A.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); + for (int i=1; i < m_A.cols(); ++i) { + res.coeffRef(i,i) = pow(m_A.coeff(i,i), p); + if (m_A.coeff(i-1,i-1) == m_A.coeff(i,i)) { + res.coeffRef(i-1,i) = p * pow(m_A.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 if (2*abs(m_A.coeff(i-1,i-1)) < abs(m_A.coeff(i,i)) || 2*abs(m_A.coeff(i,i)) < abs(m_A.coeff(i-1,i-1))) { + res.coeffRef(i-1,i) = m_A.coeff(i-1,i) * (res.coeff(i,i)-res.coeff(i-1,i-1)) / (m_A.coeff(i,i)-m_A.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)); + Scalar w = internal::matrix_power_unwinder<Scalar>::run(m_A.coeff(i,i), m_A.coeff(i-1,i-1), unwindingNumber); + res.coeffRef(i-1,i) = m_A.coeff(i-1,i) * RealScalar(2) * std::exp(RealScalar(0.5)*p*(logTdiag[i]+logTdiag[i-1])) * + std::sinh(p * w) / (m_A.coeff(i,i) - m_A.coeff(i-1,i-1)); } } } template<typename MatrixType> -void MatrixPowerTriangularAtomic<MatrixType,true>::computeBig(MatrixType& res, RealScalar p) const +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 @@ -372,10 +367,10 @@ void MatrixPowerTriangularAtomic<MatrixType,true>::computeBig(MatrixType& res, R 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; + MatrixType IminusT, sqrtT, T = m_A.template triangularView<Upper>(); RealScalar normIminusT; - int degree, degree2, numberOfSquareRoots=0; - bool hasExtraSquareRoot=false; + int degree, degree2, numberOfSquareRoots = 0; + bool hasExtraSquareRoot = false; while (true) { IminusT = m_Id - T; @@ -388,7 +383,7 @@ void MatrixPowerTriangularAtomic<MatrixType,true>::computeBig(MatrixType& res, R hasExtraSquareRoot = true; } MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT); - T = sqrtT; + T = sqrtT.template triangularView<Upper>(); ++numberOfSquareRoots; } computePade(degree, IminusT, res, p); diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index 95c63c574..b7b6423a8 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -9,6 +9,33 @@ #include "matrix_functions.h" +template <typename MatrixType, int IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> +struct generateTriangularMatrix; + +// for real matrices, make sure none of the eigenvalues are negative +template <typename MatrixType> +struct generateTriangularMatrix<MatrixType,0> +{ + static void run(MatrixType& result, typename MatrixType::Index size) + { + result.resize(size, size); + result.template triangularView<Upper>() = MatrixType::Random(size, size); + for (typename MatrixType::Index i = 0; i < size; ++i) + result.coeffRef(i,i) = std::abs(result.coeff(i,i)); + } +}; + +// for complex matrices, any matrix is fine +template <typename MatrixType> +struct generateTriangularMatrix<MatrixType,1> +{ + static void run(MatrixType& result, typename MatrixType::Index size) + { + result.resize(size, size); + result.template triangularView<Upper>() = MatrixType::Random(size, size); + } +}; + template<typename T> void test2dRotation(double tol) { @@ -59,7 +86,7 @@ void testExponentLaws(const MatrixType& m, double tol) MatrixType m1, m2, m3, m4, m5; RealScalar x, y; - for (int i=0; i<g_repeat; ++i) { + for (int i=0; i < g_repeat; ++i) { generateTestMatrix<MatrixType>::run(m1, m.rows()); MatrixPower<MatrixType> mpow(m1); @@ -90,7 +117,7 @@ void testProduct(const MatrixType& m, const VectorType& v, double tol) VectorType v1, v2, v3; RealScalar p; - for (int i=0; i<g_repeat; ++i) { + for (int i=0; i < g_repeat; ++i) { generateTestMatrix<MatrixType>::run(m1, m.rows()); MatrixPower<MatrixType> mpow(m1); @@ -99,7 +126,29 @@ void testProduct(const MatrixType& m, const VectorType& v, double tol) v2.noalias() = mpow(p) * v1; v3.noalias() = mpow(p).eval() * v1; - std::cout << "testMatrixVectorProduct: error powerm = " << relerr(v2, v3) << '\n'; + std::cout << "testProduct: error powerm = " << relerr(v2, v3) << '\n'; + VERIFY(v2.isApprox(v3, static_cast<RealScalar>(tol))); + } +} + +template<typename MatrixType, typename VectorType> +void testTriangularProduct(const MatrixType& m, const VectorType& v, double tol) +{ + typedef typename MatrixType::RealScalar RealScalar; + MatrixType m1; + VectorType v1, v2, v3; + RealScalar p; + + for (int i=0; i < g_repeat; ++i) { + generateTriangularMatrix<MatrixType>::run(m1, m.rows()); + MatrixPowerTriangular<MatrixType> mpow(m1); + + v1 = VectorType::Random(v.rows(), v.cols()); + p = internal::random<RealScalar>(); + + v2.noalias() = mpow(p) * v1; + v3.noalias() = mpow(p).eval() * v1; + std::cout << "testTriangularProduct: error powerm = " << relerr(v2, v3) << '\n'; VERIFY(v2.isApprox(v3, static_cast<RealScalar>(tol))); } } @@ -109,6 +158,7 @@ void testMatrixVector(const MatrixType& m, const VectorType& v, double tol) { testExponentLaws(m,tol); testProduct(m,v,tol); + testTriangularProduct(m,v,tol); } void test_matrix_power() |