From 75b3391e3fb047d028730cdff77baeb276ee3f99 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Thu, 4 Jul 2013 18:37:46 +0800 Subject: [PATCH 01/34] Enable singular matrix power using unitary similarities. --- .../Eigen/src/MatrixFunctions/MatrixPower.h | 129 ++++++++++-------- 1 file changed, 72 insertions(+), 57 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index c32437281..18f6703b6 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -49,14 +49,14 @@ class MatrixPowerAtomic typedef typename MatrixType::RealScalar RealScalar; typedef std::complex ComplexScalar; typedef typename MatrixType::Index Index; - typedef Array ArrayType; + typedef Block ResultType; const MatrixType& m_A; RealScalar m_p; - void computePade(int degree, const MatrixType& IminusT, MatrixType& res) const; - void compute2x2(MatrixType& res, RealScalar p) const; - void computeBig(MatrixType& res) const; + void computePade(int degree, const MatrixType& IminusT, ResultType& res) const; + void compute2x2(ResultType& res, RealScalar p) const; + void computeBig(ResultType& res) const; static int getPadeDegree(float normIminusT); static int getPadeDegree(double normIminusT); static int getPadeDegree(long double normIminusT); @@ -65,7 +65,7 @@ class MatrixPowerAtomic public: MatrixPowerAtomic(const MatrixType& T, RealScalar p); - void compute(MatrixType& res) const; + void compute(ResultType& res) const; }; template @@ -74,9 +74,8 @@ MatrixPowerAtomic::MatrixPowerAtomic(const MatrixType& T, RealScalar { eigen_assert(T.rows() == T.cols()); } template -void MatrixPowerAtomic::compute(MatrixType& res) const +void MatrixPowerAtomic::compute(ResultType& res) const { - res.resizeLike(m_A); switch (m_A.rows()) { case 0: break; @@ -92,25 +91,24 @@ void MatrixPowerAtomic::compute(MatrixType& res) const } template -void MatrixPowerAtomic::computePade(int degree, const MatrixType& IminusT, MatrixType& res) const +void MatrixPowerAtomic::computePade(int degree, const MatrixType& IminusT, ResultType& res) const { - int i = degree<<1; - res = (m_p-degree) / ((i-1)<<1) * IminusT; + int i = 2*degree; + res = (m_p-degree) / (2*i-2) * IminusT; for (--i; i; --i) { res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView() - .solve((i==1 ? -m_p : i&1 ? (-m_p-(i>>1))/(i<<1) : (m_p-(i>>1))/((i-1)<<1)) * IminusT).eval(); + .solve((i==1 ? -m_p : i&1 ? (-m_p-i/2)/(2*i) : (m_p-i/2)/(2*i-2)) * IminusT).eval(); } res += MatrixType::Identity(IminusT.rows(), IminusT.cols()); } // This function assumes that res has the correct size (see bug 614) template -void MatrixPowerAtomic::compute2x2(MatrixType& res, RealScalar p) const +void MatrixPowerAtomic::compute2x2(ResultType& res, RealScalar p) const { using std::abs; using std::pow; - ArrayType logTdiag = m_A.diagonal().array().log(); res.coeffRef(0,0) = pow(m_A.coeff(0,0), p); for (Index i=1; i < m_A.cols(); ++i) { @@ -126,7 +124,7 @@ void MatrixPowerAtomic::compute2x2(MatrixType& res, RealScalar p) co } template -void MatrixPowerAtomic::computeBig(MatrixType& res) const +void MatrixPowerAtomic::computeBig(ResultType& res) const { const int digits = std::numeric_limits::digits; const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision @@ -139,19 +137,6 @@ void MatrixPowerAtomic::computeBig(MatrixType& res) const int degree, degree2, numberOfSquareRoots = 0; bool hasExtraSquareRoot = false; - /* FIXME - * For singular T, norm(I - T) >= 1 but maxNormForPade < 1, leads to infinite - * loop. We should move 0 eigenvalues to bottom right corner. We need not - * worry about tiny values (e.g. 1e-300) because they will reach 1 if - * repetitively sqrt'ed. - * - * If the 0 eigenvalues are semisimple, they can form a 0 matrix at the - * bottom right corner. - * - * [ T A ]^p [ T^p (T^-1 T^p A) ] - * [ ] = [ ] - * [ 0 0 ] [ 0 0 ] - */ for (Index i=0; i < m_A.cols(); ++i) eigen_assert(m_A(i,i) != RealScalar(0)); @@ -275,12 +260,6 @@ template class MatrixPower { private: - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime, - MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, - MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime - }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; @@ -294,7 +273,11 @@ class MatrixPower * The class stores a reference to A, so it should not be changed * (or destroyed) before evaluation. */ - explicit MatrixPower(const MatrixType& A) : m_A(A), m_conditionNumber(0) + explicit MatrixPower(const MatrixType& A) : + m_A(A), + m_conditionNumber(0), + m_rank(A.cols()), + m_nulls(0) { eigen_assert(A.rows() == A.cols()); } /** @@ -322,15 +305,17 @@ class MatrixPower private: typedef std::complex ComplexScalar; - typedef Matrix ComplexMatrix; + typedef Matrix ComplexMatrix; typename MatrixType::Nested m_A; MatrixType m_tmp; ComplexMatrix m_T, m_U, m_fT; RealScalar m_conditionNumber; + Index m_rank, m_nulls; - RealScalar modfAndInit(RealScalar, RealScalar*); + RealScalar split(RealScalar, RealScalar*); + void initialize(); template void computeIntPower(ResultType&, RealScalar); @@ -362,37 +347,62 @@ void MatrixPower::compute(ResultType& res, RealScalar p) res(0,0) = std::pow(m_A.coeff(0,0), p); break; default: - RealScalar intpart, x = modfAndInit(p, &intpart); + RealScalar intpart, x = split(p, &intpart); computeIntPower(res, intpart); - computeFracPower(res, x); + if (x) computeFracPower(res, x); } } template -typename MatrixPower::RealScalar -MatrixPower::modfAndInit(RealScalar x, RealScalar* intpart) +typename MatrixPower::RealScalar MatrixPower::split(RealScalar x, RealScalar* intpart) { - typedef Array RealArray; - *intpart = std::floor(x); RealScalar res = x - *intpart; - if (!m_conditionNumber && res) { - const ComplexSchur schurOfA(m_A); - m_T = schurOfA.matrixT(); - m_U = schurOfA.matrixU(); - - const RealArray absTdiag = m_T.diagonal().array().abs(); - m_conditionNumber = absTdiag.maxCoeff() / absTdiag.minCoeff(); - } + if (!m_conditionNumber && res) + initialize(); - if (res>RealScalar(0.5) && res>(1-res)*std::pow(m_conditionNumber, res)) { + if (res > RealScalar(0.5) && res > (1-res) * std::pow(m_conditionNumber, res)) { --res; ++*intpart; } return res; } +template +void MatrixPower::initialize() +{ + const ComplexSchur schurOfA(m_A); + JacobiRotation rot; + ComplexScalar eigenvalue; + + m_fT.resizeLike(m_A); + m_T = schurOfA.matrixT(); + m_U = schurOfA.matrixU(); + m_conditionNumber = m_T.diagonal().array().abs().maxCoeff() / m_T.diagonal().array().abs().minCoeff(); + + for (Index i = cols()-1; i>=0; --i) { + if (m_rank <= 2) + return; + if (m_T.coeff(i,i) == RealScalar(0)) { + for (Index j=i+1; j < m_rank; ++j) { + eigenvalue = m_T.coeff(j,j); + rot.makeGivens(m_T.coeff(j-1,j), eigenvalue); + m_T.applyOnTheRight(j-1, j, rot); + m_T.applyOnTheLeft(j-1, j, rot.adjoint()); + m_T.coeffRef(j-1,j-1) = eigenvalue; + m_T.coeffRef(j,j) = RealScalar(0); + m_U.applyOnTheRight(j-1, j, rot); + } + --m_rank; + } + } + + m_nulls = rows() - m_rank; + if (m_nulls) + m_fT.bottomRows(m_nulls).fill(RealScalar(0)); +} + template template void MatrixPower::computeIntPower(ResultType& res, RealScalar p) @@ -415,12 +425,17 @@ template template void MatrixPower::computeFracPower(ResultType& res, RealScalar p) { - if (p) { - eigen_assert(m_conditionNumber); - MatrixPowerAtomic(m_T, p).compute(m_fT); - revertSchur(m_tmp, m_fT, m_U); - res = m_tmp * res; + Block blockTp(m_fT, 0, 0, m_rank, m_rank); + eigen_assert(m_conditionNumber); + eigen_assert(m_rank + m_nulls == rows()); + + MatrixPowerAtomic(m_T.topLeftCorner(m_rank, m_rank), p).compute(blockTp); + if (m_nulls) { + m_fT.topRightCorner(m_rank, m_nulls) = m_T.topLeftCorner(m_rank, m_rank).template triangularView() + .solve(blockTp * m_T.topRightCorner(m_rank, m_nulls)); } + revertSchur(m_tmp, m_fT, m_U); + res = m_tmp * res; } template From cce68d4e91a6b8f53bdea448144a300c60f0097b Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Thu, 4 Jul 2013 18:39:33 +0800 Subject: [PATCH 02/34] Remove unused inclusions. --- unsupported/Eigen/MatrixFunctions | 2 -- 1 file changed, 2 deletions(-) diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index 0991817d5..df49fdafd 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -13,8 +13,6 @@ #include #include -#include -#include #include #include From 4e26057f66d62d65739557bee72f4b9ad211e495 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Fri, 5 Jul 2013 00:08:11 +0800 Subject: [PATCH 03/34] Remove unused declarations for MatrixPowerProduct. --- Eigen/src/Core/MatrixBase.h | 3 --- Eigen/src/Core/util/ForwardDeclarations.h | 1 - 2 files changed, 4 deletions(-) diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 198e51084..84b83ffdf 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -162,9 +162,6 @@ template class MatrixBase #ifndef EIGEN_PARSED_BY_DOXYGEN template Derived& lazyAssign(const ProductBase& other); - - template - Derived& lazyAssign(const MatrixPowerProduct& other); #endif // not EIGEN_PARSED_BY_DOXYGEN template diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index d6a814586..7dc87968c 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -271,7 +271,6 @@ template class MatrixFunctionReturnValue; template class MatrixSquareRootReturnValue; template class MatrixLogarithmReturnValue; template class MatrixPowerReturnValue; -template class MatrixPowerProduct; namespace internal { template From 04a9ad6e10a5b216c868d1f182380ac079be9837 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Fri, 5 Jul 2013 00:28:28 +0800 Subject: [PATCH 04/34] Let complex power fall back to "log, scale, exp". --- Eigen/src/Core/MatrixBase.h | 1 + Eigen/src/Core/util/ForwardDeclarations.h | 1 + .../Eigen/src/MatrixFunctions/MatrixPower.h | 32 +++++++++++++++++++ 3 files changed, 34 insertions(+) diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 84b83ffdf..0628ebd1f 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -455,6 +455,7 @@ template class MatrixBase const MatrixSquareRootReturnValue sqrt() const; const MatrixLogarithmReturnValue log() const; const MatrixPowerReturnValue pow(const RealScalar& p) const; + const MatrixComplexPowerReturnValue pow(const std::complex& p) const; #ifdef EIGEN2_SUPPORT template diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 7dc87968c..cbedb7da4 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -271,6 +271,7 @@ template class MatrixFunctionReturnValue; template class MatrixSquareRootReturnValue; template class MatrixLogarithmReturnValue; template class MatrixPowerReturnValue; +template class MatrixComplexPowerReturnValue; namespace internal { template diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index 18f6703b6..a49db1916 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -503,6 +503,30 @@ class MatrixPowerReturnValue : public ReturnByValue< MatrixPowerReturnValue +class MatrixComplexPowerReturnValue : public ReturnByValue< MatrixComplexPowerReturnValue > +{ + public: + typedef typename Derived::PlainObject PlainObject; + typedef typename std::complex ComplexScalar; + typedef typename Derived::Index Index; + + MatrixComplexPowerReturnValue(const Derived& A, const ComplexScalar& p) : m_A(A), m_p(p) + { } + + template + inline void evalTo(ResultType& res) const + { res = (m_p * m_A.log()).exp(); } + + Index rows() const { return m_A.rows(); } + Index cols() const { return m_A.cols(); } + + private: + const Derived& m_A; + const ComplexScalar m_p; + MatrixComplexPowerReturnValue& operator=(const MatrixComplexPowerReturnValue&); +}; + namespace internal { template @@ -513,12 +537,20 @@ template struct traits< MatrixPowerReturnValue > { typedef typename Derived::PlainObject ReturnType; }; +template +struct traits< MatrixComplexPowerReturnValue > +{ typedef typename Derived::PlainObject ReturnType; }; + } template const MatrixPowerReturnValue MatrixBase::pow(const RealScalar& p) const { return MatrixPowerReturnValue(derived(), p); } +template +const MatrixComplexPowerReturnValue MatrixBase::pow(const std::complex& p) const +{ return MatrixComplexPowerReturnValue(derived(), p); } + } // namespace Eigen #endif // EIGEN_MATRIX_POWER From c273a6c37c74949ee62136b73290f0f9568bb5e1 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Fri, 5 Jul 2013 03:33:56 +0800 Subject: [PATCH 05/34] Avoid pow(Scalar, int) for C++11 conformance. --- .../src/MatrixFunctions/MatrixExponential.h | 44 ++++++++++++------- 1 file changed, 28 insertions(+), 16 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index 6825a7882..f8bbd59c1 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -157,6 +157,24 @@ class MatrixExponential { /** \brief L1 norm of m_M. */ RealScalar m_l1norm; + + /** \brief Scaling operator. */ + struct ScalingOp; +}; + +template +struct MatrixExponential::ScalingOp +{ + ScalingOp(int squarings) : m_squarings(squarings) { } + + inline const RealScalar operator() (const RealScalar& x) const + { return std::ldexp(x, -m_squarings); } + + inline const ComplexScalar operator() (const ComplexScalar& x) const + { return ComplexScalar(std::ldexp(x.real(), -m_squarings), std::ldexp(x.imag(), -m_squarings)); } + + private: + int m_squarings; }; template @@ -283,17 +301,15 @@ EIGEN_STRONG_INLINE void MatrixExponential::pade17(const MatrixType template void MatrixExponential::computeUV(float) { - using std::frexp; - using std::pow; if (m_l1norm < 4.258730016922831e-001) { pade3(m_M); } else if (m_l1norm < 1.880152677804762e+000) { pade5(m_M); } else { const float maxnorm = 3.925724783138660f; - frexp(m_l1norm / maxnorm, &m_squarings); + std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = m_M / pow(Scalar(2), m_squarings); + MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade7(A); } } @@ -301,8 +317,6 @@ void MatrixExponential::computeUV(float) template void MatrixExponential::computeUV(double) { - using std::frexp; - using std::pow; if (m_l1norm < 1.495585217958292e-002) { pade3(m_M); } else if (m_l1norm < 2.539398330063230e-001) { @@ -313,9 +327,9 @@ void MatrixExponential::computeUV(double) pade9(m_M); } else { const double maxnorm = 5.371920351148152; - frexp(m_l1norm / maxnorm, &m_squarings); + std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = m_M / pow(Scalar(2), m_squarings); + MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade13(A); } } @@ -323,8 +337,6 @@ void MatrixExponential::computeUV(double) template void MatrixExponential::computeUV(long double) { - using std::frexp; - using std::pow; #if LDBL_MANT_DIG == 53 // double precision computeUV(double()); #elif LDBL_MANT_DIG <= 64 // extended precision @@ -338,9 +350,9 @@ void MatrixExponential::computeUV(long double) pade9(m_M); } else { const long double maxnorm = 4.0246098906697353063L; - frexp(m_l1norm / maxnorm, &m_squarings); + std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = m_M / pow(Scalar(2), m_squarings); + MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade13(A); } #elif LDBL_MANT_DIG <= 106 // double-double @@ -356,9 +368,9 @@ void MatrixExponential::computeUV(long double) pade13(m_M); } else { const long double maxnorm = 3.2579440895405400856599663723517L; - frexp(m_l1norm / maxnorm, &m_squarings); + std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = m_M / pow(Scalar(2), m_squarings); + MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade17(A); } #elif LDBL_MANT_DIG <= 112 // quadruple precison @@ -374,9 +386,9 @@ void MatrixExponential::computeUV(long double) pade13(m_M); } else { const long double maxnorm = 2.884233277829519311757165057717815L; - frexp(m_l1norm / maxnorm, &m_squarings); + std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = m_M / pow(Scalar(2), m_squarings); + MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade17(A); } #else From 9e2b4eeac0dcf62f27e05bb0ae923d5c48ee98ed Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Fri, 5 Jul 2013 23:28:57 +0800 Subject: [PATCH 06/34] Const-correct the scaling functor. --- .../Eigen/src/MatrixFunctions/MatrixExponential.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index f8bbd59c1..ca4532357 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -309,7 +309,7 @@ void MatrixExponential::computeUV(float) const float maxnorm = 3.925724783138660f; std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = CwiseUnaryOp(m_M, m_squarings); + MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade7(A); } } @@ -329,7 +329,7 @@ void MatrixExponential::computeUV(double) const double maxnorm = 5.371920351148152; std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = CwiseUnaryOp(m_M, m_squarings); + MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade13(A); } } @@ -352,7 +352,7 @@ void MatrixExponential::computeUV(long double) const long double maxnorm = 4.0246098906697353063L; std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = CwiseUnaryOp(m_M, m_squarings); + MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade13(A); } #elif LDBL_MANT_DIG <= 106 // double-double @@ -370,7 +370,7 @@ void MatrixExponential::computeUV(long double) const long double maxnorm = 3.2579440895405400856599663723517L; std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = CwiseUnaryOp(m_M, m_squarings); + MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade17(A); } #elif LDBL_MANT_DIG <= 112 // quadruple precison @@ -388,7 +388,7 @@ void MatrixExponential::computeUV(long double) const long double maxnorm = 2.884233277829519311757165057717815L; std::frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; - MatrixType A = CwiseUnaryOp(m_M, m_squarings); + MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade17(A); } #else From 55ec3cc6d50807ae99dfe603e15ff5f7d075ecb3 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sun, 7 Jul 2013 19:34:13 +0800 Subject: [PATCH 07/34] Prevent copying with internal::noncopyable. --- .../Eigen/src/MatrixFunctions/MatrixExponential.h | 8 ++------ unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h | 10 +++------- .../Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h | 6 +----- .../Eigen/src/MatrixFunctions/MatrixLogarithm.h | 6 +----- unsupported/Eigen/src/MatrixFunctions/MatrixPower.h | 4 ++-- .../Eigen/src/MatrixFunctions/MatrixSquareRoot.h | 6 +++--- unsupported/Eigen/src/MatrixFunctions/StemFunction.h | 2 +- 7 files changed, 13 insertions(+), 29 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index ca4532357..e5d4e3ad2 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -21,8 +21,8 @@ namespace Eigen { * expected to be an instantiation of the Matrix class template. */ template -class MatrixExponential { - +class MatrixExponential : internal::noncopyable +{ public: /** \brief Constructor. @@ -43,10 +43,6 @@ class MatrixExponential { private: - // Prevent copying - MatrixExponential(const MatrixExponential&); - MatrixExponential& operator=(const MatrixExponential&); - /** \brief Compute the (3,3)-Padé approximant to the exponential. * * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h index 7d426640c..a5476ac26 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h @@ -34,7 +34,7 @@ namespace Eigen { template ::Scalar>::IsComplex> -class MatrixFunction +class MatrixFunction : internal::noncopyable { public: @@ -65,7 +65,7 @@ class MatrixFunction * \brief Partial specialization of MatrixFunction for real matrices */ template -class MatrixFunction +class MatrixFunction : internal::noncopyable { private: @@ -111,8 +111,6 @@ class MatrixFunction private: typename internal::nested::type m_A; /**< \brief Reference to argument of matrix function. */ AtomicType& m_atomic; /**< \brief Class for computing matrix function of atomic blocks. */ - - MatrixFunction& operator=(const MatrixFunction&); }; @@ -120,7 +118,7 @@ class MatrixFunction * \brief Partial specialization of MatrixFunction for complex matrices */ template -class MatrixFunction +class MatrixFunction : internal::noncopyable { private: @@ -176,8 +174,6 @@ class MatrixFunction * separation constant is set to 0.1, a value taken from the * paper by Davies and Higham. */ static const RealScalar separation() { return static_cast(0.1); } - - MatrixFunction& operator=(const MatrixFunction&); }; /** \brief Constructor. diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h index efe332c48..d6ff5f1ce 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h @@ -21,7 +21,7 @@ namespace Eigen { * entries are close to each other. */ template -class MatrixFunctionAtomic +class MatrixFunctionAtomic : internal::noncopyable { public: @@ -44,10 +44,6 @@ class MatrixFunctionAtomic private: - // Prevent copying - MatrixFunctionAtomic(const MatrixFunctionAtomic&); - MatrixFunctionAtomic& operator=(const MatrixFunctionAtomic&); - void computeMu(); bool taylorConverged(Index s, const MatrixType& F, const MatrixType& Fincr, const MatrixType& P); diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h index c744fc05f..586c034af 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h @@ -28,7 +28,7 @@ namespace Eigen { * \sa class MatrixFunctionAtomic, MatrixBase::log() */ template -class MatrixLogarithmAtomic +class MatrixLogarithmAtomic : internal::noncopyable { public: @@ -71,10 +71,6 @@ private: std::numeric_limits::digits<= 64? 8: // extended precision std::numeric_limits::digits<=106? 10: // double-double 11; // quadruple precision - - // Prevent copying - MatrixLogarithmAtomic(const MatrixLogarithmAtomic&); - MatrixLogarithmAtomic& operator=(const MatrixLogarithmAtomic&); }; /** \brief Compute logarithm of triangular matrix with clustered eigenvalues. */ diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index a49db1916..5d580ac9e 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -38,7 +38,7 @@ class MatrixPowerRetval : public ReturnByValue< MatrixPowerRetval > }; template -class MatrixPowerAtomic +class MatrixPowerAtomic : internal::noncopyable { private: enum { @@ -257,7 +257,7 @@ MatrixPowerAtomic::computeSuperDiag(RealScalar curr, RealScalar prev * Output: \verbinclude MatrixPower_optimal.out */ template -class MatrixPower +class MatrixPower : internal::noncopyable { private: typedef typename MatrixType::Scalar Scalar; diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h index b48ea9d46..f3bcef409 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h @@ -24,7 +24,7 @@ namespace Eigen { * \sa MatrixSquareRoot, MatrixSquareRootTriangular */ template -class MatrixSquareRootQuasiTriangular +class MatrixSquareRootQuasiTriangular : internal::noncopyable { public: @@ -253,7 +253,7 @@ void MatrixSquareRootQuasiTriangular * \sa MatrixSquareRoot, MatrixSquareRootQuasiTriangular */ template -class MatrixSquareRootTriangular +class MatrixSquareRootTriangular : internal::noncopyable { public: MatrixSquareRootTriangular(const MatrixType& A) @@ -370,7 +370,7 @@ class MatrixSquareRoot // ********** Partial specialization for complex matrices ********** template -class MatrixSquareRoot +class MatrixSquareRoot : internal::noncopyable { public: diff --git a/unsupported/Eigen/src/MatrixFunctions/StemFunction.h b/unsupported/Eigen/src/MatrixFunctions/StemFunction.h index 724e55c1d..0a815d834 100644 --- a/unsupported/Eigen/src/MatrixFunctions/StemFunction.h +++ b/unsupported/Eigen/src/MatrixFunctions/StemFunction.h @@ -16,7 +16,7 @@ namespace Eigen { * \brief Stem functions corresponding to standard mathematical functions. */ template -class StdStemFunctions +class StdStemFunctions : internal::noncopyable { public: From 00e30a5fc43bd3ff177c5e12fb53eaba8fb65ce8 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sun, 7 Jul 2013 19:57:23 +0800 Subject: [PATCH 08/34] We need not prohibit assignment here. Thanks to changeset 3edd4681f2f04c1164cb3805f1ac37fbf9a618c0 . --- unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h | 2 -- unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h | 2 -- unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h | 2 -- unsupported/Eigen/src/MatrixFunctions/MatrixPower.h | 3 --- unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h | 2 -- 5 files changed, 11 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index e5d4e3ad2..221ec3115 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -435,8 +435,6 @@ template struct MatrixExponentialReturnValue protected: const Derived& m_src; - private: - MatrixExponentialReturnValue& operator=(const MatrixExponentialReturnValue&); }; namespace internal { diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h index a5476ac26..ad54d8ed0 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h @@ -527,8 +527,6 @@ template class MatrixFunctionReturnValue private: typename internal::nested::type m_A; StemFunction *m_f; - - MatrixFunctionReturnValue& operator=(const MatrixFunctionReturnValue&); }; namespace internal { diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h index 586c034af..32ec385e0 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h @@ -454,8 +454,6 @@ public: private: typename internal::nested::type m_A; - - MatrixLogarithmReturnValue& operator=(const MatrixLogarithmReturnValue&); }; namespace internal { diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index 5d580ac9e..e0a687978 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -34,7 +34,6 @@ class MatrixPowerRetval : public ReturnByValue< MatrixPowerRetval > private: MatrixPower& m_pow; const RealScalar m_p; - MatrixPowerRetval& operator=(const MatrixPowerRetval&); }; template @@ -500,7 +499,6 @@ class MatrixPowerReturnValue : public ReturnByValue< MatrixPowerReturnValue @@ -524,7 +522,6 @@ class MatrixComplexPowerReturnValue : public ReturnByValue< MatrixComplexPowerRe private: const Derived& m_A; const ComplexScalar m_p; - MatrixComplexPowerReturnValue& operator=(const MatrixComplexPowerReturnValue&); }; namespace internal { diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h index f3bcef409..c405bfd61 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h @@ -442,8 +442,6 @@ template class MatrixSquareRootReturnValue protected: const Derived& m_src; - private: - MatrixSquareRootReturnValue& operator=(const MatrixSquareRootReturnValue&); }; namespace internal { From 04bd1e3fc003f4220e82a553c3496a4f8d509959 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Mon, 8 Jul 2013 16:49:27 +0800 Subject: [PATCH 09/34] Slightly optimize atanh2. --- Eigen/src/Core/MathFunctions.h | 52 ++++++++++++++++++++-------------- 1 file changed, 30 insertions(+), 22 deletions(-) diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 5df2d8bca..d415bc719 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -332,37 +332,45 @@ inline NewType cast(const OldType& x) * Implementation of atanh2 * ****************************************************************************/ -template -struct atanh2_default_impl +template +struct atanh2_impl { - typedef Scalar retval; - typedef typename NumTraits::Real RealScalar; - static inline Scalar run(const Scalar& x, const Scalar& y) + static inline Scalar run(const Scalar& x, const Scalar& r) + { + EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar) + #if __cplusplus >= 201103L + using std::log1p; + return log1p(2 * x / (r - x)) / 2; + #else + using std::abs; + using std::log; + using std::sqrt; + Scalar z = x / r; + if (r == 0 || abs(z) > sqrt(NumTraits::epsilon())) + return log((r + x) / (r - x)) / 2; + else + return z + z*z*z / 3; + #endif + } +}; + +template +struct atanh2_impl > +{ + typedef std::complex Scalar; + static inline Scalar run(const Scalar& x, const Scalar& r) { - using std::abs; using std::log; + using std::norm; using std::sqrt; - Scalar z = x / y; - if (y == Scalar(0) || abs(z) > sqrt(NumTraits::epsilon())) - return RealScalar(0.5) * log((y + x) / (y - x)); + Scalar z = x / r; + if (r == Scalar(0) || norm(z) > NumTraits::epsilon()) + return RealScalar(0.5) * log((r + x) / (r - x)); else return z + z*z*z / RealScalar(3); } }; -template -struct atanh2_default_impl -{ - static inline Scalar run(const Scalar&, const Scalar&) - { - EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar) - return Scalar(0); - } -}; - -template -struct atanh2_impl : atanh2_default_impl::IsInteger> {}; - template struct atanh2_retval { From 25544dbec3429848226c9a567ccd7e82973c04e7 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Wed, 10 Jul 2013 02:36:34 +0800 Subject: [PATCH 10/34] Add assertion against undefined matrix power. --- unsupported/Eigen/src/MatrixFunctions/MatrixPower.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index e0a687978..e4f13c993 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -398,8 +398,11 @@ void MatrixPower::initialize() } m_nulls = rows() - m_rank; - if (m_nulls) + if (m_nulls) { + eigen_assert(m_T.bottomRightCorner(m_nulls, m_nulls).isZero() + && "Base of matrix power should be invertible or with a semisimple zero eigenvalue."); m_fT.bottomRows(m_nulls).fill(RealScalar(0)); + } } template From 159a3bed9e26274ccc8da07a08ea394285d05bd3 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Wed, 10 Jul 2013 02:43:10 +0800 Subject: [PATCH 11/34] Write doc for complex power of a matrix. --- unsupported/Eigen/MatrixFunctions | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index df49fdafd..0bdd379d7 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -228,15 +228,16 @@ const MatrixPowerReturnValue MatrixBase::pow(RealScalar p) con \endcode \param[in] M base of the matrix power, should be a square matrix. -\param[in] p exponent of the matrix power, should be real. +\param[in] p exponent of the matrix power. The matrix power \f$ M^p \f$ is defined as \f$ \exp(p \log(M)) \f$, where exp denotes the matrix exponential, and log denotes the matrix logarithm. -The matrix \f$ M \f$ should meet the conditions to be an argument of -matrix logarithm. If \p p is not of the real scalar type of \p M, it -is casted into the real scalar type of \p M. +If \p p is complex, the scalar type of \p M should be the type of \p +p . \f$ M^p \f$ simply evaluates into \f$ \exp(p \log(M)) \f$. +Therefore, the matrix \f$ M \f$ should meet the conditions to be an +argument of matrix logarithm. This function computes the matrix power using the Schur-Padé algorithm as implemented by class MatrixPower. The exponent is split From c52cbd9de95553f2f8c7a9020903081db719573c Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Wed, 10 Jul 2013 02:44:38 +0800 Subject: [PATCH 12/34] Write doc for positive power of a matrix with a semisimple zero eigenvalue. --- unsupported/Eigen/MatrixFunctions | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index 0bdd379d7..41dfab390 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -239,12 +239,29 @@ p . \f$ M^p \f$ simply evaluates into \f$ \exp(p \log(M)) \f$. Therefore, the matrix \f$ M \f$ should meet the conditions to be an argument of matrix logarithm. -This function computes the matrix power using the Schur-Padé +If \p p is real, it is casted into the real scalar type of \p M. Then +this function computes the matrix power using the Schur-Padé algorithm as implemented by class MatrixPower. The exponent is split into integral part and fractional part, where the fractional part is in the interval \f$ (-1, 1) \f$. The main diagonal and the first super-diagonal is directly computed. +If \p M is singular with a semisimple zero eigenvalue and \p p is +positive, the Schur factor \f$ T \f$ is reordered with Givens +rotations, i.e. + +\f[ T = \left[ \begin{array}{cc} + T_1 & T_2 \\ + 0 & 0 + \end{array} \right] \f] + +where \f$ T_1 \f$ is invertible. Then \f$ T^p \f$ is given by + +\f[ T^p = \left[ \begin{array}{cc} + T_1^p & T_1^{-1} T_1^p T_2 \\ + 0 & 0 + \end{array}. \right] \f] + Details of the algorithm can be found in: Nicholas J. Higham and Lijing Lin, "A Schur-Padé algorithm for fractional powers of a matrix," SIAM J. %Matrix Anal. Applic., From d204bb57d04e1d106c33790407ee4727e0f9816b Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Wed, 10 Jul 2013 02:48:17 +0800 Subject: [PATCH 13/34] Remove unused struct definition in test. --- unsupported/test/matrix_power.cpp | 27 --------------------------- 1 file changed, 27 deletions(-) diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index b9d513b45..38b16fba9 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -9,33 +9,6 @@ #include "matrix_functions.h" -template ::IsComplex> -struct generateTriangularMatrix; - -// for real matrices, make sure none of the eigenvalues are negative -template -struct generateTriangularMatrix -{ - static void run(MatrixType& result, typename MatrixType::Index size) - { - result.resize(size, size); - result.template triangularView() = 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 -struct generateTriangularMatrix -{ - static void run(MatrixType& result, typename MatrixType::Index size) - { - result.resize(size, size); - result.template triangularView() = MatrixType::Random(size, size); - } -}; - template void test2dRotation(double tol) { From 639d03d900e844171af0cc3b98dfa98eca79dad5 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Wed, 10 Jul 2013 02:53:15 +0800 Subject: [PATCH 14/34] These casts are unnecessary because isApprox already casts them. --- unsupported/test/matrix_power.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index 38b16fba9..8d53ec94f 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -26,7 +26,7 @@ void test2dRotation(double tol) C = Apow(std::ldexp(angle,1) / M_PI); std::cout << "test2dRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n'; - VERIFY(C.isApprox(B, static_cast(tol))); + VERIFY(C.isApprox(B, tol)); } } @@ -48,7 +48,7 @@ void test2dHyperbolicRotation(double tol) C = Apow(angle); std::cout << "test2dHyperbolicRotation: i = " << i << " error powerm = " << relerr(C,B) << '\n'; - VERIFY(C.isApprox(B, static_cast(tol))); + VERIFY(C.isApprox(B, tol)); } } @@ -70,15 +70,15 @@ void testExponentLaws(const MatrixType& m, double tol) m4 = mpow(x+y); m5.noalias() = m2 * m3; - VERIFY(m4.isApprox(m5, static_cast(tol))); + VERIFY(m4.isApprox(m5, tol)); m4 = mpow(x*y); m5 = m2.pow(y); - VERIFY(m4.isApprox(m5, static_cast(tol))); + VERIFY(m4.isApprox(m5, tol)); m4 = (std::abs(x) * m1).pow(y); m5 = std::pow(std::abs(x), y) * m3; - VERIFY(m4.isApprox(m5, static_cast(tol))); + VERIFY(m4.isApprox(m5, tol)); } } From 5c95892b83f9c1b2fb7f5a236cfb965bf8e55c3f Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Wed, 10 Jul 2013 02:57:54 +0800 Subject: [PATCH 15/34] Test power of singular matrices. --- unsupported/test/matrix_power.cpp | 90 +++++++++++++++++++++++++++---- 1 file changed, 80 insertions(+), 10 deletions(-) diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index 8d53ec94f..a3055c81e 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -9,6 +9,36 @@ #include "matrix_functions.h" +// for complex matrices, any matrix is fine +template::Scalar>::IsComplex> +struct generateSingularMatrix +{ + static void run(MatrixType& result, typename MatrixType::Index size) + { + result = MatrixType::Random(size, size); + result.col(0).fill(0); + } +}; + +// for real matrices, make sure none of the eigenvalues are negative +template +struct generateSingularMatrix +{ + static void run(MatrixType& result, typename MatrixType::Index size) + { + MatrixType mat = MatrixType::Random(size, size); + mat.col(0).fill(0); + ComplexSchur schur(mat); + typename ComplexSchur::ComplexMatrixType T = schur.matrixT(); + + for (typename MatrixType::Index i = 0; i < size; ++i) { + if (T.coeff(i,i).imag() == 0 && T.coeff(i,i).real() < 0) + T.coeffRef(i,i) = -T.coeff(i,i); + } + result = (schur.matrixU() * (T.template triangularView() * schur.matrixU().adjoint())).real(); + } +}; + template void test2dRotation(double tol) { @@ -53,7 +83,7 @@ void test2dHyperbolicRotation(double tol) } template -void testExponentLaws(const MatrixType& m, double tol) +void testGeneral(const MatrixType& m, double tol) { typedef typename MatrixType::RealScalar RealScalar; MatrixType m1, m2, m3, m4, m5; @@ -82,6 +112,36 @@ void testExponentLaws(const MatrixType& m, double tol) } } +template +void testSingular(MatrixType m, double tol) +{ + typedef typename MatrixType::RealScalar RealScalar; + MatrixType m1, m2, m3, m4, m5; + RealScalar x, y; + + for (int i=0; i < g_repeat; ++i) { + generateTestMatrix::run(m1, m.rows()); + MatrixPower mpow(m1); + + x = internal::random(0, 1); + y = internal::random(0, 1); + m2 = mpow(x); + m3 = mpow(y); + + m4 = mpow(x+y); + m5.noalias() = m2 * m3; + VERIFY(m4.isApprox(m5, tol)); + + m4 = mpow(x*y); + m5 = m2.pow(y); + VERIFY(m4.isApprox(m5, tol)); + + m4 = (x * m1).pow(y); + m5 = std::pow(x, y) * m3; + VERIFY(m4.isApprox(m5, tol)); + } +} + typedef Matrix Matrix3dRowMajor; typedef Matrix MatrixXe; @@ -94,13 +154,23 @@ void test_matrix_power() CALL_SUBTEST_1(test2dHyperbolicRotation(1e-5)); CALL_SUBTEST_9(test2dHyperbolicRotation(1e-14)); - CALL_SUBTEST_2(testExponentLaws(Matrix2d(), 1e-13)); - CALL_SUBTEST_7(testExponentLaws(Matrix3dRowMajor(), 1e-13)); - CALL_SUBTEST_3(testExponentLaws(Matrix4cd(), 1e-13)); - CALL_SUBTEST_4(testExponentLaws(MatrixXd(8,8), 2e-12)); - CALL_SUBTEST_1(testExponentLaws(Matrix2f(), 1e-4)); - CALL_SUBTEST_5(testExponentLaws(Matrix3cf(), 1e-4)); - CALL_SUBTEST_8(testExponentLaws(Matrix4f(), 1e-4)); - CALL_SUBTEST_6(testExponentLaws(MatrixXf(2,2), 1e-3)); // see bug 614 - CALL_SUBTEST_9(testExponentLaws(MatrixXe(7,7), 1e-13)); + CALL_SUBTEST_2(testGeneral(Matrix2d(), 1e-13)); + CALL_SUBTEST_7(testGeneral(Matrix3dRowMajor(), 1e-13)); + CALL_SUBTEST_3(testGeneral(Matrix4cd(), 1e-13)); + CALL_SUBTEST_4(testGeneral(MatrixXd(8,8), 2e-12)); + CALL_SUBTEST_1(testGeneral(Matrix2f(), 1e-4)); + CALL_SUBTEST_5(testGeneral(Matrix3cf(), 1e-4)); + CALL_SUBTEST_8(testGeneral(Matrix4f(), 1e-4)); + CALL_SUBTEST_6(testGeneral(MatrixXf(2,2), 1e-3)); // see bug 614 + CALL_SUBTEST_9(testGeneral(MatrixXe(7,7), 1e-13)); + + CALL_SUBTEST_2(testSingular(Matrix2d(), 1e-13)); + CALL_SUBTEST_7(testSingular(Matrix3dRowMajor(), 1e-13)); + CALL_SUBTEST_3(testSingular(Matrix4cd(), 1e-13)); + CALL_SUBTEST_4(testSingular(MatrixXd(8,8), 2e-12)); + CALL_SUBTEST_1(testSingular(Matrix2f(), 1e-4)); + CALL_SUBTEST_5(testSingular(Matrix3cf(), 1e-4)); + CALL_SUBTEST_8(testSingular(Matrix4f(), 1e-4)); + CALL_SUBTEST_6(testSingular(MatrixXf(2,2), 1e-3)); + CALL_SUBTEST_9(testSingular(MatrixXe(7,7), 1e-13)); } From 4466875d544c7a279b84fabb3b67a73da6b15487 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Wed, 10 Jul 2013 02:59:16 +0800 Subject: [PATCH 16/34] The only(?) way to test complex matrix power. --- unsupported/test/matrix_power.cpp | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index a3055c81e..b9b72b528 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -142,6 +142,19 @@ void testSingular(MatrixType m, double tol) } } +template +void testLogThenExp(MatrixType m, double tol) +{ + typedef typename MatrixType::Scalar Scalar; + Scalar x; + + for (int i=0; i < g_repeat; ++i) { + generateTestMatrix::run(m, m.rows()); + x = internal::random(); + VERIFY(m.pow(x).isApprox((x * m.log()).exp(), tol)); + } +} + typedef Matrix Matrix3dRowMajor; typedef Matrix MatrixXe; @@ -173,4 +186,14 @@ void test_matrix_power() CALL_SUBTEST_8(testSingular(Matrix4f(), 1e-4)); CALL_SUBTEST_6(testSingular(MatrixXf(2,2), 1e-3)); CALL_SUBTEST_9(testSingular(MatrixXe(7,7), 1e-13)); + + CALL_SUBTEST_2(testLogThenExp(Matrix2d(), 1e-13)); + CALL_SUBTEST_7(testLogThenExp(Matrix3dRowMajor(), 1e-13)); + CALL_SUBTEST_3(testLogThenExp(Matrix4cd(), 1e-13)); + CALL_SUBTEST_4(testLogThenExp(MatrixXd(8,8), 2e-12)); + CALL_SUBTEST_1(testLogThenExp(Matrix2f(), 1e-4)); + CALL_SUBTEST_5(testLogThenExp(Matrix3cf(), 1e-4)); + CALL_SUBTEST_8(testLogThenExp(Matrix4f(), 1e-4)); + CALL_SUBTEST_6(testLogThenExp(MatrixXf(2,2), 1e-3)); + CALL_SUBTEST_9(testLogThenExp(MatrixXe(7,7), 1e-13)); } From a992fa74ebb560d6c25a2debf8a1460eb326b8f8 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Thu, 11 Jul 2013 02:31:13 +0800 Subject: [PATCH 17/34] Make non-conversion unary constructors explicit. --- .../Eigen/src/MatrixFunctions/MatrixExponential.h | 4 ++-- .../Eigen/src/MatrixFunctions/MatrixLogarithm.h | 2 +- .../Eigen/src/MatrixFunctions/MatrixSquareRoot.h | 12 ++++++------ 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index 221ec3115..750b0b989 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -32,7 +32,7 @@ class MatrixExponential : internal::noncopyable * * \param[in] M matrix whose exponential is to be computed. */ - MatrixExponential(const MatrixType &M); + explicit MatrixExponential(const MatrixType &M); /** \brief Computes the matrix exponential. * @@ -415,7 +415,7 @@ template struct MatrixExponentialReturnValue * \param[in] src %Matrix (expression) forming the argument of the * matrix exponential. */ - MatrixExponentialReturnValue(const Derived& src) : m_src(src) { } + explicit MatrixExponentialReturnValue(const Derived& src) : m_src(src) { } /** \brief Compute the matrix exponential. * diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h index 32ec385e0..33cfadfb4 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h @@ -425,7 +425,7 @@ public: * * \param[in] A %Matrix (expression) forming the argument of the matrix logarithm. */ - MatrixLogarithmReturnValue(const Derived& A) : m_A(A) { } + explicit MatrixLogarithmReturnValue(const Derived& A) : m_A(A) { } /** \brief Compute the matrix logarithm. * diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h index c405bfd61..0cd39ebe4 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h @@ -36,7 +36,7 @@ class MatrixSquareRootQuasiTriangular : internal::noncopyable * The class stores a reference to \p A, so it should not be * changed (or destroyed) before compute() is called. */ - MatrixSquareRootQuasiTriangular(const MatrixType& A) + explicit MatrixSquareRootQuasiTriangular(const MatrixType& A) : m_A(A) { eigen_assert(A.rows() == A.cols()); @@ -256,7 +256,7 @@ template class MatrixSquareRootTriangular : internal::noncopyable { public: - MatrixSquareRootTriangular(const MatrixType& A) + explicit MatrixSquareRootTriangular(const MatrixType& A) : m_A(A) { eigen_assert(A.rows() == A.cols()); @@ -321,7 +321,7 @@ class MatrixSquareRoot * The class stores a reference to \p A, so it should not be * changed (or destroyed) before compute() is called. */ - MatrixSquareRoot(const MatrixType& A); + explicit MatrixSquareRoot(const MatrixType& A); /** \brief Compute the matrix square root * @@ -341,7 +341,7 @@ class MatrixSquareRoot { public: - MatrixSquareRoot(const MatrixType& A) + explicit MatrixSquareRoot(const MatrixType& A) : m_A(A) { eigen_assert(A.rows() == A.cols()); @@ -374,7 +374,7 @@ class MatrixSquareRoot : internal::noncopyable { public: - MatrixSquareRoot(const MatrixType& A) + explicit MatrixSquareRoot(const MatrixType& A) : m_A(A) { eigen_assert(A.rows() == A.cols()); @@ -422,7 +422,7 @@ template class MatrixSquareRootReturnValue * \param[in] src %Matrix (expression) forming the argument of the * matrix square root. */ - MatrixSquareRootReturnValue(const Derived& src) : m_src(src) { } + explicit MatrixSquareRootReturnValue(const Derived& src) : m_src(src) { } /** \brief Compute the matrix square root. * From 738d75d3eb187bf1b3169d0de4328c6a2a5ef500 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sat, 13 Jul 2013 22:11:36 +0800 Subject: [PATCH 18/34] Document on the return type of MatrixPower::operator() --- .../Eigen/src/MatrixFunctions/MatrixPower.h | 42 ++++++++++++++++--- 1 file changed, 37 insertions(+), 5 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index e4f13c993..5ba6278ae 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -14,16 +14,48 @@ namespace Eigen { template class MatrixPower; +/** + * \ingroup MatrixFunctions_Module + * + * \brief Proxy for the matrix power of some matrix. + * + * \tparam MatrixType type of the base, a matrix. + * + * This class holds the arguments to the matrix power until it is + * assigned or evaluated for some other reason (so the argument + * should not be changed in the meantime). It is the return type of + * MatrixPower::operator() and related functions and most of the + * time this is the only way it is used. + */ +/* TODO This class is only used by MatrixPower, so it should be nested + * into MatrixPower, like MatrixPower::ReturnValue. However, my + * compiler complained about unused template parameter in the + * following declaration in namespace internal. + * + * template + * struct traits::ReturnValue>; + */ template -class MatrixPowerRetval : public ReturnByValue< MatrixPowerRetval > +class MatrixPowerParenthesesReturnValue : public ReturnByValue< MatrixPowerParenthesesReturnValue > { public: typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; - MatrixPowerRetval(MatrixPower& pow, RealScalar p) : m_pow(pow), m_p(p) + /** + * \brief Constructor. + * + * \param[in] pow %MatrixPower storing the base. + * \param[in] p scalar, the exponent of the matrix power. + */ + MatrixPowerParenthesesReturnValue(MatrixPower& pow, RealScalar p) : m_pow(pow), m_p(p) { } + /** + * \brief Compute the matrix power. + * + * \param[out] result + */ template inline void evalTo(ResultType& res) const { m_pow.compute(res, m_p); } @@ -286,8 +318,8 @@ class MatrixPower : internal::noncopyable * \return The expression \f$ A^p \f$, where A is specified in the * constructor. */ - const MatrixPowerRetval operator()(RealScalar p) - { return MatrixPowerRetval(*this, p); } + const MatrixPowerParenthesesReturnValue operator()(RealScalar p) + { return MatrixPowerParenthesesReturnValue(*this, p); } /** * \brief Compute the matrix power. @@ -530,7 +562,7 @@ class MatrixComplexPowerReturnValue : public ReturnByValue< MatrixComplexPowerRe namespace internal { template -struct traits< MatrixPowerRetval > +struct traits< MatrixPowerParenthesesReturnValue > { typedef typename MatrixPowerType::PlainObject ReturnType; }; template From 3c423ccfe215283c5a5159ac74a8a23603e3f9bd Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sat, 13 Jul 2013 22:12:09 +0800 Subject: [PATCH 19/34] Document on complex matrix power. --- .../Eigen/src/MatrixFunctions/MatrixPower.h | 30 ++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index 5ba6278ae..ffb344065 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -513,7 +513,7 @@ class MatrixPowerReturnValue : public ReturnByValue< MatrixPowerReturnValue class MatrixComplexPowerReturnValue : public ReturnByValue< MatrixComplexPowerReturnValue > { @@ -544,9 +557,24 @@ class MatrixComplexPowerReturnValue : public ReturnByValue< MatrixComplexPowerRe typedef typename std::complex ComplexScalar; typedef typename Derived::Index Index; + /** + * \brief Constructor. + * + * \param[in] A %Matrix (expression), the base of the matrix power. + * \param[in] p complex scalar, the exponent of the matrix power. + */ MatrixComplexPowerReturnValue(const Derived& A, const ComplexScalar& p) : m_A(A), m_p(p) { } + /** + * \brief Compute the matrix power. + * + * Because \p p is complex, \f$ A^p \f$ is simply evaluated as \f$ + * \exp(p \log(A)) \f$. + * + * \param[out] result \f$ A^p \f$ where \p A and \p p are as in the + * constructor. + */ template inline void evalTo(ResultType& res) const { res = (m_p * m_A.log()).exp(); } From d5501d3a90c954a901fe95e654f4195744179831 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sat, 13 Jul 2013 23:13:07 +0800 Subject: [PATCH 20/34] Document on MatrixPowerAtomic. --- .../Eigen/src/MatrixFunctions/MatrixPower.h | 38 ++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index ffb344065..e1de7f606 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -68,6 +68,21 @@ class MatrixPowerParenthesesReturnValue : public ReturnByValue< MatrixPowerParen const RealScalar m_p; }; +/** + * \ingroup MatrixFunctions_Module + * + * \brief Class for computing matrix powers. + * + * \tparam MatrixType type of the base, expected to be an instantiation + * of the Matrix class template. + * + * This class is capable of computing triangular real/complex matrices + * raised to a power in the interval \f$ (-1, 1) \f$. + * + * \note Currently this class is only used by MatrixPower. One may + * insist that this be nested into MatrixPower. This class is here to + * faciliate future development of triangular matrix functions. + */ template class MatrixPowerAtomic : internal::noncopyable { @@ -95,14 +110,35 @@ class MatrixPowerAtomic : internal::noncopyable static RealScalar computeSuperDiag(RealScalar, RealScalar, RealScalar p); public: + /** + * \brief Constructor. + * + * \param[in] T the base of the matrix power. + * \param[in] p the exponent of the matrix power, should be in + * \f$ (-1, 1) \f$. + * + * The class stores a reference to T, so it should not be changed + * (or destroyed) before evaluation. Only the upper triangular + * part of T is read. + */ MatrixPowerAtomic(const MatrixType& T, RealScalar p); + + /** + * \brief Compute the matrix power. + * + * \param[out] res \f$ A^p \f$ where A and p are specified in the + * constructor. + */ void compute(ResultType& res) const; }; template MatrixPowerAtomic::MatrixPowerAtomic(const MatrixType& T, RealScalar p) : m_A(T), m_p(p) -{ eigen_assert(T.rows() == T.cols()); } +{ + eigen_assert(T.rows() == T.cols()); + eigen_assert(p > -1 && p < 1); +} template void MatrixPowerAtomic::compute(ResultType& res) const From eeb744dc8d8608df22d3e8c512f42fb507f402ae Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sun, 14 Jul 2013 02:00:50 +0800 Subject: [PATCH 21/34] Add test3dRotation. --- unsupported/test/matrix_power.cpp | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index b9b72b528..8e33d6038 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -82,6 +82,20 @@ void test2dHyperbolicRotation(double tol) } } +template +void test3dRotation(double tol) +{ + Matrix v; + T angle; + + for (int i=0; i<=20; ++i) { + v = Matrix::Random(); + v.normalize(); + angle = pow(10, (i-10) / 5.); + VERIFY(AngleAxis(angle, v).matrix().isApprox(AngleAxis(1,v).matrix().pow(angle), tol)); + } +} + template void testGeneral(const MatrixType& m, double tol) { @@ -156,6 +170,7 @@ void testLogThenExp(MatrixType m, double tol) } typedef Matrix Matrix3dRowMajor; +typedef Matrix Matrix3e; typedef Matrix MatrixXe; void test_matrix_power() @@ -167,6 +182,10 @@ void test_matrix_power() CALL_SUBTEST_1(test2dHyperbolicRotation(1e-5)); CALL_SUBTEST_9(test2dHyperbolicRotation(1e-14)); + CALL_SUBTEST_10(test3dRotation(1e-13)); + CALL_SUBTEST_11(test3dRotation(1e-5)); + CALL_SUBTEST_12(test3dRotation(1e-13)); + CALL_SUBTEST_2(testGeneral(Matrix2d(), 1e-13)); CALL_SUBTEST_7(testGeneral(Matrix3dRowMajor(), 1e-13)); CALL_SUBTEST_3(testGeneral(Matrix4cd(), 1e-13)); @@ -176,6 +195,9 @@ void test_matrix_power() CALL_SUBTEST_8(testGeneral(Matrix4f(), 1e-4)); CALL_SUBTEST_6(testGeneral(MatrixXf(2,2), 1e-3)); // see bug 614 CALL_SUBTEST_9(testGeneral(MatrixXe(7,7), 1e-13)); + CALL_SUBTEST_10(testGeneral(Matrix3d(), 1e-13)); + CALL_SUBTEST_11(testGeneral(Matrix3f(), 1e-4)); + CALL_SUBTEST_12(testGeneral(Matrix3e(), 1e-13)); CALL_SUBTEST_2(testSingular(Matrix2d(), 1e-13)); CALL_SUBTEST_7(testSingular(Matrix3dRowMajor(), 1e-13)); @@ -186,6 +208,9 @@ void test_matrix_power() CALL_SUBTEST_8(testSingular(Matrix4f(), 1e-4)); CALL_SUBTEST_6(testSingular(MatrixXf(2,2), 1e-3)); CALL_SUBTEST_9(testSingular(MatrixXe(7,7), 1e-13)); + CALL_SUBTEST_10(testSingular(Matrix3d(), 1e-13)); + CALL_SUBTEST_11(testSingular(Matrix3f(), 1e-4)); + CALL_SUBTEST_12(testSingular(Matrix3e(), 1e-13)); CALL_SUBTEST_2(testLogThenExp(Matrix2d(), 1e-13)); CALL_SUBTEST_7(testLogThenExp(Matrix3dRowMajor(), 1e-13)); @@ -196,4 +221,7 @@ void test_matrix_power() CALL_SUBTEST_8(testLogThenExp(Matrix4f(), 1e-4)); CALL_SUBTEST_6(testLogThenExp(MatrixXf(2,2), 1e-3)); CALL_SUBTEST_9(testLogThenExp(MatrixXe(7,7), 1e-13)); + CALL_SUBTEST_10(testLogThenExp(Matrix3d(), 1e-13)); + CALL_SUBTEST_11(testLogThenExp(Matrix3f(), 1e-4)); + CALL_SUBTEST_12(testLogThenExp(Matrix3e(), 1e-13)); } From cbe92de2b57a7580e4e5a2e00ad8ea52dda707b8 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sun, 14 Jul 2013 17:27:44 +0800 Subject: [PATCH 22/34] Fix typo in testSingular. --- unsupported/test/matrix_power.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index 8e33d6038..223af60bc 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -134,7 +134,7 @@ void testSingular(MatrixType m, double tol) RealScalar x, y; for (int i=0; i < g_repeat; ++i) { - generateTestMatrix::run(m1, m.rows()); + generateSingularMatrix::run(m1, m.rows()); MatrixPower mpow(m1); x = internal::random(0, 1); From b8f0364a1c56784fa666c783e21c4bd1b218b9b0 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Mon, 15 Jul 2013 00:10:17 +0800 Subject: [PATCH 23/34] Test singular matrix power with square roots. Exponent laws are too unstable. --- unsupported/test/matrix_power.cpp | 96 +++++++++++++++---------------- 1 file changed, 47 insertions(+), 49 deletions(-) diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index 223af60bc..3ee19fc56 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2012 Chen-Pang He +// Copyright (C) 2012, 2013 Chen-Pang He // // 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 @@ -9,36 +9,6 @@ #include "matrix_functions.h" -// for complex matrices, any matrix is fine -template::Scalar>::IsComplex> -struct generateSingularMatrix -{ - static void run(MatrixType& result, typename MatrixType::Index size) - { - result = MatrixType::Random(size, size); - result.col(0).fill(0); - } -}; - -// for real matrices, make sure none of the eigenvalues are negative -template -struct generateSingularMatrix -{ - static void run(MatrixType& result, typename MatrixType::Index size) - { - MatrixType mat = MatrixType::Random(size, size); - mat.col(0).fill(0); - ComplexSchur schur(mat); - typename ComplexSchur::ComplexMatrixType T = schur.matrixT(); - - for (typename MatrixType::Index i = 0; i < size; ++i) { - if (T.coeff(i,i).imag() == 0 && T.coeff(i,i).real() < 0) - T.coeffRef(i,i) = -T.coeff(i,i); - } - result = (schur.matrixU() * (T.template triangularView() * schur.matrixU().adjoint())).real(); - } -}; - template void test2dRotation(double tol) { @@ -126,33 +96,61 @@ void testGeneral(const MatrixType& m, double tol) } } +// For complex matrices, any matrix is fine. +template::Scalar>::IsComplex> +struct processTriangularMatrix +{ + static void run(MatrixType&, MatrixType&, const MatrixType&) + { } +}; + +// For real matrices, make sure none of the eigenvalues are negative. +template +struct processTriangularMatrix +{ + static void run(MatrixType& m, MatrixType& T, const MatrixType& U) + { + typedef typename MatrixType::Index Index; + const Index size = m.cols(); + + for (Index i=0; i < size; ++i) { + if (i == size - 1 || T.coeff(i+1,i) == 0) + T.coeffRef(i,i) = std::abs(T.coeff(i,i)); + else + ++i; + } + m = U * T * U.adjoint(); + } +}; + template void testSingular(MatrixType m, double tol) { - typedef typename MatrixType::RealScalar RealScalar; - MatrixType m1, m2, m3, m4, m5; - RealScalar x, y; + const int IsComplex = NumTraits::Scalar>::IsComplex; + typedef typename internal::conditional< IsComplex, MatrixSquareRootTriangular, + MatrixSquareRootQuasiTriangular >::type SquareRootType; + typedef typename internal::conditional, const MatrixType&>::type TriangularType; + typename internal::conditional< IsComplex, ComplexSchur, RealSchur >::type schur; + MatrixType T; for (int i=0; i < g_repeat; ++i) { - generateSingularMatrix::run(m1, m.rows()); - MatrixPower mpow(m1); + m.setRandom(); + m.col(0).fill(0); - x = internal::random(0, 1); - y = internal::random(0, 1); - m2 = mpow(x); - m3 = mpow(y); + schur.compute(m); + T = schur.matrixT(); + const MatrixType& U = schur.matrixU(); + processTriangularMatrix::run(m, T, U); + MatrixPower mpow(m); - m4 = mpow(x+y); - m5.noalias() = m2 * m3; - VERIFY(m4.isApprox(m5, tol)); + SquareRootType(T).compute(T); + VERIFY(mpow(0.5).isApprox(U * (TriangularType(T) * U.adjoint()), tol)); - m4 = mpow(x*y); - m5 = m2.pow(y); - VERIFY(m4.isApprox(m5, tol)); + SquareRootType(T).compute(T); + VERIFY(mpow(0.25).isApprox(U * (TriangularType(T) * U.adjoint()), tol)); - m4 = (x * m1).pow(y); - m5 = std::pow(x, y) * m3; - VERIFY(m4.isApprox(m5, tol)); + SquareRootType(T).compute(T); + VERIFY(mpow(0.125).isApprox(U * (TriangularType(T) * U.adjoint()), tol)); } } From 9be658f7015161989b7ecccd70fd050ce563cad9 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Mon, 15 Jul 2013 00:43:14 +0800 Subject: [PATCH 24/34] generateTestMatrix can use processTriangularMatrix --- unsupported/test/matrix_functions.h | 41 ++++++++++++++++++++++------- unsupported/test/matrix_power.cpp | 27 ------------------- 2 files changed, 31 insertions(+), 37 deletions(-) diff --git a/unsupported/test/matrix_functions.h b/unsupported/test/matrix_functions.h index 5817caef6..295da16b6 100644 --- a/unsupported/test/matrix_functions.h +++ b/unsupported/test/matrix_functions.h @@ -10,27 +10,48 @@ #include "main.h" #include +// For complex matrices, any matrix is fine. +template::Scalar>::IsComplex> +struct processTriangularMatrix +{ + static void run(MatrixType&, MatrixType&, const MatrixType&) + { } +}; + +// For real matrices, make sure none of the eigenvalues are negative. +template +struct processTriangularMatrix +{ + static void run(MatrixType& m, MatrixType& T, const MatrixType& U) + { + typedef typename MatrixType::Index Index; + const Index size = m.cols(); + + for (Index i=0; i < size; ++i) { + if (i == size - 1 || T.coeff(i+1,i) == 0) + T.coeffRef(i,i) = std::abs(T.coeff(i,i)); + else + ++i; + } + m = U * T * U.transpose(); + } +}; + template ::Scalar>::IsComplex> struct generateTestMatrix; -// for real matrices, make sure none of the eigenvalues are negative template struct generateTestMatrix { static void run(MatrixType& result, typename MatrixType::Index size) { - MatrixType mat = MatrixType::Random(size, size); - EigenSolver es(mat); - typename EigenSolver::EigenvalueType eivals = es.eigenvalues(); - for (typename MatrixType::Index i = 0; i < size; ++i) { - if (eivals(i).imag() == 0 && eivals(i).real() < 0) - eivals(i) = -eivals(i); - } - result = (es.eigenvectors() * eivals.asDiagonal() * es.eigenvectors().inverse()).real(); + result = MatrixType::Random(size, size); + RealSchur schur(result); + MatrixType T = schur.matrixT(); + processTriangularMatrix::run(result, T, schur.matrixU()); } }; -// for complex matrices, any matrix is fine template struct generateTestMatrix { diff --git a/unsupported/test/matrix_power.cpp b/unsupported/test/matrix_power.cpp index 3ee19fc56..849e4287b 100644 --- a/unsupported/test/matrix_power.cpp +++ b/unsupported/test/matrix_power.cpp @@ -96,33 +96,6 @@ void testGeneral(const MatrixType& m, double tol) } } -// For complex matrices, any matrix is fine. -template::Scalar>::IsComplex> -struct processTriangularMatrix -{ - static void run(MatrixType&, MatrixType&, const MatrixType&) - { } -}; - -// For real matrices, make sure none of the eigenvalues are negative. -template -struct processTriangularMatrix -{ - static void run(MatrixType& m, MatrixType& T, const MatrixType& U) - { - typedef typename MatrixType::Index Index; - const Index size = m.cols(); - - for (Index i=0; i < size; ++i) { - if (i == size - 1 || T.coeff(i+1,i) == 0) - T.coeffRef(i,i) = std::abs(T.coeff(i,i)); - else - ++i; - } - m = U * T * U.adjoint(); - } -}; - template void testSingular(MatrixType m, double tol) { From 4b780553e05e35cb4813cd57eaf12befb6062891 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Mon, 15 Jul 2013 09:10:17 +0800 Subject: [PATCH 25/34] Eliminate unnecessary copying for sparse Kronecker product. --- .../KroneckerProduct/KroneckerTensorProduct.h | 106 ++++++++++-------- unsupported/test/kronecker_product.cpp | 21 ++-- 2 files changed, 75 insertions(+), 52 deletions(-) diff --git a/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h b/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h index 532896c3b..6ec8eb558 100644 --- a/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h +++ b/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h @@ -14,7 +14,42 @@ namespace Eigen { -template class SparseMatrix; +template +class KroneckerProductBase : public ReturnByValue +{ + private: + typedef typename internal::traits Traits; + typedef typename Traits::Lhs Lhs; + typedef typename Traits::Rhs Rhs; + typedef typename Traits::Scalar Scalar; + + protected: + typedef typename Traits::Index Index; + + public: + KroneckerProductBase(const Lhs& A, const Rhs& B) + : m_A(A), m_B(B) + {} + + inline Index rows() const { return m_A.rows() * m_B.rows(); } + inline Index cols() const { return m_A.cols() * m_B.cols(); } + + Scalar coeff(Index row, Index col) const + { + return m_A.coeff(row / m_B.rows(), col / m_B.cols()) * + m_B.coeff(row % m_B.rows(), col % m_B.cols()); + } + + Scalar coeff(Index i) const + { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived); + return m_A.coeff(i / m_A.size()) * m_B.coeff(i % m_A.size()); + } + + protected: + typename Lhs::Nested m_A; + typename Rhs::Nested m_B; +}; /*! * \brief Kronecker tensor product helper class for dense matrices @@ -27,40 +62,21 @@ template class SparseMatrix; * \tparam Rhs Type of the rignt-hand side, a matrix expression. */ template -class KroneckerProduct : public ReturnByValue > +class KroneckerProduct : public KroneckerProductBase > { private: - typedef ReturnByValue Base; - typedef typename Base::Scalar Scalar; - typedef typename Base::Index Index; + typedef KroneckerProductBase Base; + using Base::m_A; + using Base::m_B; public: /*! \brief Constructor. */ KroneckerProduct(const Lhs& A, const Rhs& B) - : m_A(A), m_B(B) + : Base(A, B) {} /*! \brief Evaluate the Kronecker tensor product. */ template void evalTo(Dest& dst) const; - - inline Index rows() const { return m_A.rows() * m_B.rows(); } - inline Index cols() const { return m_A.cols() * m_B.cols(); } - - Scalar coeff(Index row, Index col) const - { - return m_A.coeff(row / m_B.rows(), col / m_B.cols()) * - m_B.coeff(row % m_B.rows(), col % m_B.cols()); - } - - Scalar coeff(Index i) const - { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(KroneckerProduct); - return m_A.coeff(i / m_A.size()) * m_B.coeff(i % m_A.size()); - } - - private: - typename Lhs::Nested m_A; - typename Rhs::Nested m_B; }; /*! @@ -77,40 +93,28 @@ class KroneckerProduct : public ReturnByValue > * \tparam Rhs Type of the rignt-hand side, a matrix expression. */ template -class KroneckerProductSparse : public EigenBase > +class KroneckerProductSparse : public KroneckerProductBase > { private: - typedef typename internal::traits::Index Index; + typedef KroneckerProductBase Base; + using Base::m_A; + using Base::m_B; public: /*! \brief Constructor. */ KroneckerProductSparse(const Lhs& A, const Rhs& B) - : m_A(A), m_B(B) + : Base(A, B) {} /*! \brief Evaluate the Kronecker tensor product. */ template void evalTo(Dest& dst) const; - - inline Index rows() const { return m_A.rows() * m_B.rows(); } - inline Index cols() const { return m_A.cols() * m_B.cols(); } - - template - operator SparseMatrix() - { - SparseMatrix result; - evalTo(result.derived()); - return result; - } - - private: - typename Lhs::Nested m_A; - typename Rhs::Nested m_B; }; template template void KroneckerProduct::evalTo(Dest& dst) const { + typedef typename Base::Index Index; const int BlockRows = Rhs::RowsAtCompileTime, BlockCols = Rhs::ColsAtCompileTime; const Index Br = m_B.rows(), @@ -124,9 +128,10 @@ template template void KroneckerProductSparse::evalTo(Dest& dst) const { + typedef typename Base::Index Index; const Index Br = m_B.rows(), Bc = m_B.cols(); - dst.resize(rows(),cols()); + dst.resize(this->rows(), this->cols()); dst.resizeNonZeros(0); dst.reserve(m_A.nonZeros() * m_B.nonZeros()); @@ -155,6 +160,7 @@ struct traits > typedef typename remove_all<_Lhs>::type Lhs; typedef typename remove_all<_Rhs>::type Rhs; typedef typename scalar_product_traits::ReturnType Scalar; + typedef typename promote_index_type::type Index; enum { Rows = size_at_compile_time::RowsAtCompileTime, traits::RowsAtCompileTime>::ret, @@ -193,6 +199,8 @@ struct traits > | EvalBeforeNestingBit | EvalBeforeAssigningBit, CoeffReadCost = Dynamic }; + + typedef SparseMatrix ReturnType; }; } // end namespace internal @@ -228,6 +236,16 @@ KroneckerProduct kroneckerProduct(const MatrixBase& a, const MatrixBase< * Computes Kronecker tensor product of two matrices, at least one of * which is sparse * + * \warning If you want to replace a matrix by its Kronecker product + * with some matrix, do \b NOT do this: + * \code + * A = kroneckerProduct(A,B); // bug!!! caused by aliasing effect + * \endcode + * instead, use eval() to work around this: + * \code + * A = kroneckerProduct(A,B).eval(); + * \endcode + * * \param a Dense/sparse matrix a * \param b Dense/sparse matrix b * \return Kronecker tensor product of a and b, stored in a sparse diff --git a/unsupported/test/kronecker_product.cpp b/unsupported/test/kronecker_product.cpp index 8ddc6ec28..c68a07de8 100644 --- a/unsupported/test/kronecker_product.cpp +++ b/unsupported/test/kronecker_product.cpp @@ -107,31 +107,34 @@ void test_kronecker_product() SparseMatrix SM_row_a(SM_a), SM_row_b(SM_b); - // test kroneckerProduct(DM_block,DM,DM_fixedSize) + // test DM_fixedSize = kroneckerProduct(DM_block,DM) Matrix DM_fix_ab = kroneckerProduct(DM_a.topLeftCorner<2,3>(),DM_b); CALL_SUBTEST(check_kronecker_product(DM_fix_ab)); + CALL_SUBTEST(check_kronecker_product(kroneckerProduct(DM_a.topLeftCorner<2,3>(),DM_b))); for(int i=0;i(2,5) = kroneckerProduct(DM_a,DM_b); CALL_SUBTEST(check_kronecker_product(DM_block_ab.block<6,6>(2,5))); - // test kroneckerProduct(DM,DM,DM) + // test DM = kroneckerProduct(DM,DM) MatrixXd DM_ab = kroneckerProduct(DM_a,DM_b); CALL_SUBTEST(check_kronecker_product(DM_ab)); + CALL_SUBTEST(check_kronecker_product(kroneckerProduct(DM_a,DM_b))); - // test kroneckerProduct(SM,DM,SM) + // test SM = kroneckerProduct(SM,DM) SparseMatrix SM_ab = kroneckerProduct(SM_a,DM_b); CALL_SUBTEST(check_kronecker_product(SM_ab)); SparseMatrix SM_ab2 = kroneckerProduct(SM_a,DM_b); CALL_SUBTEST(check_kronecker_product(SM_ab2)); + CALL_SUBTEST(check_kronecker_product(kroneckerProduct(SM_a,DM_b))); - // test kroneckerProduct(DM,SM,SM) + // test SM = kroneckerProduct(DM,SM) SM_ab.setZero(); SM_ab.insert(0,0)=37.0; SM_ab = kroneckerProduct(DM_a,SM_b); @@ -140,8 +143,9 @@ void test_kronecker_product() SM_ab2.insert(0,0)=37.0; SM_ab2 = kroneckerProduct(DM_a,SM_b); CALL_SUBTEST(check_kronecker_product(SM_ab2)); + CALL_SUBTEST(check_kronecker_product(kroneckerProduct(DM_a,SM_b))); - // test kroneckerProduct(SM,SM,SM) + // test SM = kroneckerProduct(SM,SM) SM_ab.resize(2,33); SM_ab.insert(0,0)=37.0; SM_ab = kroneckerProduct(SM_a,SM_b); @@ -150,8 +154,9 @@ void test_kronecker_product() SM_ab2.insert(0,0)=37.0; SM_ab2 = kroneckerProduct(SM_a,SM_b); CALL_SUBTEST(check_kronecker_product(SM_ab2)); + CALL_SUBTEST(check_kronecker_product(kroneckerProduct(SM_a,SM_b))); - // test kroneckerProduct(SM,SM,SM) with sparse pattern + // test SM = kroneckerProduct(SM,SM) with sparse pattern SM_a.resize(4,5); SM_b.resize(3,2); SM_a.resizeNonZeros(0); @@ -169,7 +174,7 @@ void test_kronecker_product() SM_ab = kroneckerProduct(SM_a,SM_b); CALL_SUBTEST(check_sparse_kronecker_product(SM_ab)); - // test dimension of result of kroneckerProduct(DM,DM,DM) + // test dimension of result of DM = kroneckerProduct(DM,DM) MatrixXd DM_a2(2,1); MatrixXd DM_b2(5,4); MatrixXd DM_ab2 = kroneckerProduct(DM_a2,DM_b2); From c587e63631b6c37e842033c2b0438f269865dd0a Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sat, 20 Jul 2013 17:49:38 +0800 Subject: [PATCH 26/34] Simplify MatrixPower::split --- .../Eigen/src/MatrixFunctions/MatrixPower.h | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index e1de7f606..b43789bd0 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -381,7 +381,7 @@ class MatrixPower : internal::noncopyable RealScalar m_conditionNumber; Index m_rank, m_nulls; - RealScalar split(RealScalar, RealScalar*); + void split(RealScalar&, RealScalar&); void initialize(); template @@ -414,26 +414,26 @@ void MatrixPower::compute(ResultType& res, RealScalar p) res(0,0) = std::pow(m_A.coeff(0,0), p); break; default: - RealScalar intpart, x = split(p, &intpart); + RealScalar intpart; + split(p, intpart); computeIntPower(res, intpart); - if (x) computeFracPower(res, x); + if (p) computeFracPower(res, p); } } template -typename MatrixPower::RealScalar MatrixPower::split(RealScalar x, RealScalar* intpart) +void MatrixPower::split(RealScalar& p, RealScalar& intpart) { - *intpart = std::floor(x); - RealScalar res = x - *intpart; + intpart = std::floor(p); + p -= intpart; - if (!m_conditionNumber && res) + if (!m_conditionNumber && p) initialize(); - if (res > RealScalar(0.5) && res > (1-res) * std::pow(m_conditionNumber, res)) { - --res; - ++*intpart; + if (p > RealScalar(0.5) && p > (1-p) * std::pow(m_conditionNumber, p)) { + --p; + ++intpart; } - return res; } template From 2320073e4153b60405e5460eab397c93c51de7e9 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sat, 20 Jul 2013 17:58:12 +0800 Subject: [PATCH 27/34] Comment on private members of MatrixPower. --- .../Eigen/src/MatrixFunctions/MatrixPower.h | 49 ++++++++++++++++--- 1 file changed, 42 insertions(+), 7 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index b43789bd0..be13095e5 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -375,20 +375,51 @@ class MatrixPower : internal::noncopyable typedef Matrix ComplexMatrix; + /** \brief Reference to the base of matrix power. */ typename MatrixType::Nested m_A; - MatrixType m_tmp; - ComplexMatrix m_T, m_U, m_fT; - RealScalar m_conditionNumber; - Index m_rank, m_nulls; - void split(RealScalar&, RealScalar&); + /** \brief Temporary storage for computing integral power. */ + MatrixType m_tmp; + + /** \brief Store the result of Schur decomposition. */ + ComplexMatrix m_T, m_U; + + /** \brief Store fractional power of m_T. */ + ComplexMatrix m_fT; + + /** + * \brief Condition number of m_A. + * + * It is initialized as 0 to avoid performing unnecessary Schur + * decomposition, which is the bottleneck. + */ + RealScalar m_conditionNumber; + + /** \brief Rank of m_A. */ + Index m_rank; + + /** \brief Rank deficiency of m_A. */ + Index m_nulls; + + /** + * \brief Split p into integral part and fractional part. + * + * \param[in] p The exponent. + * \param[out] p The fractional part ranging in \f$ (-1, 1) \f$. + * \param[out] intpart The integral part. + * + * Only if the fractional part is nonzero, it calls initialize(). + */ + void split(RealScalar& p, RealScalar& intpart); + + /** \brief Perform Schur decomposition for fractional power. */ void initialize(); template - void computeIntPower(ResultType&, RealScalar); + void computeIntPower(ResultType& res, RealScalar p); template - void computeFracPower(ResultType&, RealScalar); + void computeFracPower(ResultType& res, RealScalar p); template static void revertSchur( @@ -427,9 +458,12 @@ void MatrixPower::split(RealScalar& p, RealScalar& intpart) intpart = std::floor(p); p -= intpart; + // Perform Schur decomposition if it is not yet performed and the power is + // not an integer. if (!m_conditionNumber && p) initialize(); + // Choose the more stable of intpart = floor(p) and intpart = ceil(p). if (p > RealScalar(0.5) && p > (1-p) * std::pow(m_conditionNumber, p)) { --p; ++intpart; @@ -448,6 +482,7 @@ void MatrixPower::initialize() m_U = schurOfA.matrixU(); m_conditionNumber = m_T.diagonal().array().abs().maxCoeff() / m_T.diagonal().array().abs().minCoeff(); + // Move zero eigenvalues to the bottom right corner. for (Index i = cols()-1; i>=0; --i) { if (m_rank <= 2) return; From dda869051db084d22dcfc46b7732b04d77ff9b4b Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sat, 20 Jul 2013 18:47:54 +0800 Subject: [PATCH 28/34] Optimize MatrixPower::computeIntPower --- .../Eigen/src/MatrixFunctions/MatrixPower.h | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index be13095e5..7124874f7 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -447,6 +447,8 @@ void MatrixPower::compute(ResultType& res, RealScalar p) default: RealScalar intpart; split(p, intpart); + + res = MatrixType::Identity(rows(), cols()); computeIntPower(res, intpart); if (p) computeFracPower(res, p); } @@ -514,15 +516,18 @@ void MatrixPower::computeIntPower(ResultType& res, RealScalar p) { RealScalar pp = std::abs(p); - if (p<0) m_tmp = m_A.inverse(); - else m_tmp = m_A; + if (p<0) + m_tmp = m_A.inverse(); + else + m_tmp = m_A; - res = MatrixType::Identity(rows(), cols()); - while (pp >= 1) { + while (true) { if (std::fmod(pp, 2) >= 1) res = m_tmp * res; - m_tmp *= m_tmp; pp /= 2; + if (pp < 1) + break; + m_tmp *= m_tmp; } } From ede27f5780880e5fb89662d11db5a0e3b5586766 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sat, 20 Jul 2013 23:30:37 +0800 Subject: [PATCH 29/34] Apply argument-dependent lookup on user-defined types. (using std::) --- .../src/MatrixFunctions/MatrixExponential.h | 23 +++++++--- .../Eigen/src/MatrixFunctions/MatrixPower.h | 43 +++++++++++++------ 2 files changed, 46 insertions(+), 20 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index 750b0b989..810f426f9 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -164,10 +164,16 @@ struct MatrixExponential::ScalingOp ScalingOp(int squarings) : m_squarings(squarings) { } inline const RealScalar operator() (const RealScalar& x) const - { return std::ldexp(x, -m_squarings); } + { + using std::ldexp; + return ldexp(x, -m_squarings); + } inline const ComplexScalar operator() (const ComplexScalar& x) const - { return ComplexScalar(std::ldexp(x.real(), -m_squarings), std::ldexp(x.imag(), -m_squarings)); } + { + using std::ldexp; + return ComplexScalar(ldexp(x.real(), -m_squarings), ldexp(x.imag(), -m_squarings)); + } private: int m_squarings; @@ -297,13 +303,14 @@ EIGEN_STRONG_INLINE void MatrixExponential::pade17(const MatrixType template void MatrixExponential::computeUV(float) { + using std::frexp; if (m_l1norm < 4.258730016922831e-001) { pade3(m_M); } else if (m_l1norm < 1.880152677804762e+000) { pade5(m_M); } else { const float maxnorm = 3.925724783138660f; - std::frexp(m_l1norm / maxnorm, &m_squarings); + frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade7(A); @@ -313,6 +320,7 @@ void MatrixExponential::computeUV(float) template void MatrixExponential::computeUV(double) { + using std::frexp; if (m_l1norm < 1.495585217958292e-002) { pade3(m_M); } else if (m_l1norm < 2.539398330063230e-001) { @@ -323,7 +331,7 @@ void MatrixExponential::computeUV(double) pade9(m_M); } else { const double maxnorm = 5.371920351148152; - std::frexp(m_l1norm / maxnorm, &m_squarings); + frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade13(A); @@ -333,6 +341,7 @@ void MatrixExponential::computeUV(double) template void MatrixExponential::computeUV(long double) { + using std::frexp; #if LDBL_MANT_DIG == 53 // double precision computeUV(double()); #elif LDBL_MANT_DIG <= 64 // extended precision @@ -346,7 +355,7 @@ void MatrixExponential::computeUV(long double) pade9(m_M); } else { const long double maxnorm = 4.0246098906697353063L; - std::frexp(m_l1norm / maxnorm, &m_squarings); + frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade13(A); @@ -364,7 +373,7 @@ void MatrixExponential::computeUV(long double) pade13(m_M); } else { const long double maxnorm = 3.2579440895405400856599663723517L; - std::frexp(m_l1norm / maxnorm, &m_squarings); + frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade17(A); @@ -382,7 +391,7 @@ void MatrixExponential::computeUV(long double) pade13(m_M); } else { const long double maxnorm = 2.884233277829519311757165057717815L; - std::frexp(m_l1norm / maxnorm, &m_squarings); + frexp(m_l1norm / maxnorm, &m_squarings); if (m_squarings < 0) m_squarings = 0; MatrixType A = CwiseUnaryOp(m_M, m_squarings); pade17(A); diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index 7124874f7..b419e245b 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -143,11 +143,12 @@ MatrixPowerAtomic::MatrixPowerAtomic(const MatrixType& T, RealScalar template void MatrixPowerAtomic::compute(ResultType& res) const { + using std::pow; switch (m_A.rows()) { case 0: break; case 1: - res(0,0) = std::pow(m_A(0,0), m_p); + res(0,0) = pow(m_A(0,0), m_p); break; case 2: compute2x2(res, m_p); @@ -162,6 +163,7 @@ void MatrixPowerAtomic::computePade(int degree, const MatrixType& Im { int i = 2*degree; res = (m_p-degree) / (2*i-2) * IminusT; + for (--i; i; --i) { res = (MatrixType::Identity(IminusT.rows(), IminusT.cols()) + res).template triangularView() .solve((i==1 ? -m_p : i&1 ? (-m_p-i/2)/(2*i) : (m_p-i/2)/(2*i-2)) * IminusT).eval(); @@ -175,7 +177,6 @@ void MatrixPowerAtomic::compute2x2(ResultType& res, RealScalar p) co { using std::abs; using std::pow; - res.coeffRef(0,0) = pow(m_A.coeff(0,0), p); for (Index i=1; i < m_A.cols(); ++i) { @@ -193,6 +194,7 @@ void MatrixPowerAtomic::compute2x2(ResultType& res, RealScalar p) co template void MatrixPowerAtomic::computeBig(ResultType& res) const { + using std::ldexp; const int digits = std::numeric_limits::digits; const RealScalar maxNormForPade = digits <= 24? 4.3386528e-1f: // sigle precision digits <= 53? 2.789358995219730e-1: // double precision @@ -224,7 +226,7 @@ void MatrixPowerAtomic::computeBig(ResultType& res) const computePade(degree, IminusT, res); for (; numberOfSquareRoots; --numberOfSquareRoots) { - compute2x2(res, std::ldexp(m_p, -numberOfSquareRoots)); + compute2x2(res, ldexp(m_p, -numberOfSquareRoots)); res = res.template triangularView() * res; } compute2x2(res, m_p); @@ -289,19 +291,28 @@ template inline typename MatrixPowerAtomic::ComplexScalar MatrixPowerAtomic::computeSuperDiag(const ComplexScalar& curr, const ComplexScalar& prev, RealScalar p) { - ComplexScalar logCurr = std::log(curr); - ComplexScalar logPrev = std::log(prev); - int unwindingNumber = std::ceil((numext::imag(logCurr - logPrev) - M_PI) / (2*M_PI)); + using std::ceil; + using std::exp; + using std::log; + using std::sinh; + + ComplexScalar logCurr = log(curr); + ComplexScalar logPrev = log(prev); + int unwindingNumber = ceil((numext::imag(logCurr - logPrev) - M_PI) / (2*M_PI)); ComplexScalar w = numext::atanh2(curr - prev, curr + prev) + ComplexScalar(0, M_PI*unwindingNumber); - return RealScalar(2) * std::exp(RealScalar(0.5) * p * (logCurr + logPrev)) * std::sinh(p * w) / (curr - prev); + return RealScalar(2) * exp(RealScalar(0.5) * p * (logCurr + logPrev)) * sinh(p * w) / (curr - prev); } template inline typename MatrixPowerAtomic::RealScalar MatrixPowerAtomic::computeSuperDiag(RealScalar curr, RealScalar prev, RealScalar p) { + using std::exp; + using std::log; + using std::sinh; + RealScalar w = numext::atanh2(curr - prev, curr + prev); - return 2 * std::exp(p * (std::log(curr) + std::log(prev)) / 2) * std::sinh(p * w) / (curr - prev); + return 2 * exp(p * (log(curr) + log(prev)) / 2) * sinh(p * w) / (curr - prev); } /** @@ -438,11 +449,12 @@ template template void MatrixPower::compute(ResultType& res, RealScalar p) { + using std::pow; switch (cols()) { case 0: break; case 1: - res(0,0) = std::pow(m_A.coeff(0,0), p); + res(0,0) = pow(m_A.coeff(0,0), p); break; default: RealScalar intpart; @@ -457,7 +469,10 @@ void MatrixPower::compute(ResultType& res, RealScalar p) template void MatrixPower::split(RealScalar& p, RealScalar& intpart) { - intpart = std::floor(p); + using std::floor; + using std::pow; + + intpart = floor(p); p -= intpart; // Perform Schur decomposition if it is not yet performed and the power is @@ -466,7 +481,7 @@ void MatrixPower::split(RealScalar& p, RealScalar& intpart) initialize(); // Choose the more stable of intpart = floor(p) and intpart = ceil(p). - if (p > RealScalar(0.5) && p > (1-p) * std::pow(m_conditionNumber, p)) { + if (p > RealScalar(0.5) && p > (1-p) * pow(m_conditionNumber, p)) { --p; ++intpart; } @@ -514,7 +529,9 @@ template template void MatrixPower::computeIntPower(ResultType& res, RealScalar p) { - RealScalar pp = std::abs(p); + using std::abs; + using std::fmod; + RealScalar pp = abs(p); if (p<0) m_tmp = m_A.inverse(); @@ -522,7 +539,7 @@ void MatrixPower::computeIntPower(ResultType& res, RealScalar p) m_tmp = m_A; while (true) { - if (std::fmod(pp, 2) >= 1) + if (fmod(pp, 2) >= 1) res = m_tmp * res; pp /= 2; if (pp < 1) From 3d94ed9fa0bc06771c39970bcae6705ca2524cad Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sun, 21 Jul 2013 00:18:19 +0800 Subject: [PATCH 30/34] Document on MatrixExponential::ScalingOp --- .../src/MatrixFunctions/MatrixExponential.h | 20 ++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index 810f426f9..30b36a179 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -126,6 +126,7 @@ class MatrixExponential : internal::noncopyable */ void computeUV(long double); + struct ScalingOp; typedef typename internal::traits::Scalar Scalar; typedef typename NumTraits::Real RealScalar; typedef typename std::complex ComplexScalar; @@ -153,22 +154,35 @@ class MatrixExponential : internal::noncopyable /** \brief L1 norm of m_M. */ RealScalar m_l1norm; - - /** \brief Scaling operator. */ - struct ScalingOp; }; +/** \brief Scaling operator. + * + * This struct is used by CwiseUnaryOp to scale a matrix by \f$ 2^{-s} \f$. + */ template struct MatrixExponential::ScalingOp { + /** \brief Constructor. + * + * \param[in] squarings The integer \f$ s \f$ in this document. + */ ScalingOp(int squarings) : m_squarings(squarings) { } + /** \brief Scale a matrix coefficient. + * + * \param[in,out] x The scalar to be scaled, becoming \f$ 2^{-s} x \f$. + */ inline const RealScalar operator() (const RealScalar& x) const { using std::ldexp; return ldexp(x, -m_squarings); } + /** \brief Scale a matrix coefficient. + * + * \param[in,out] x The scalar to be scaled, becoming \f$ 2^{-s} x \f$. + */ inline const ComplexScalar operator() (const ComplexScalar& x) const { using std::ldexp; From 1191949e8725106035d0f7276857217585a115ef Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sun, 21 Jul 2013 00:19:46 +0800 Subject: [PATCH 31/34] Improve documentation on Kronecker product module. --- .../KroneckerProduct/KroneckerTensorProduct.h | 22 ++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h b/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h index 6ec8eb558..80be9fd7a 100644 --- a/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h +++ b/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h @@ -12,8 +12,15 @@ #ifndef KRONECKER_TENSOR_PRODUCT_H #define KRONECKER_TENSOR_PRODUCT_H -namespace Eigen { +namespace Eigen { +/*! + * \ingroup KroneckerProduct_Module + * + * \brief The base class of dense and sparse Kronecker product. + * + * \tparam Derived is the derived type. + */ template class KroneckerProductBase : public ReturnByValue { @@ -27,6 +34,7 @@ class KroneckerProductBase : public ReturnByValue typedef typename Traits::Index Index; public: + /*! \brief Constructor. */ KroneckerProductBase(const Lhs& A, const Rhs& B) : m_A(A), m_B(B) {} @@ -34,12 +42,20 @@ class KroneckerProductBase : public ReturnByValue inline Index rows() const { return m_A.rows() * m_B.rows(); } inline Index cols() const { return m_A.cols() * m_B.cols(); } + /*! + * This overrides ReturnByValue::coeff because this function is + * efficient enough. + */ Scalar coeff(Index row, Index col) const { return m_A.coeff(row / m_B.rows(), col / m_B.cols()) * m_B.coeff(row % m_B.rows(), col % m_B.cols()); } + /*! + * This overrides ReturnByValue::coeff because this function is + * efficient enough. + */ Scalar coeff(Index i) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived); @@ -52,6 +68,8 @@ class KroneckerProductBase : public ReturnByValue }; /*! + * \ingroup KroneckerProduct_Module + * * \brief Kronecker tensor product helper class for dense matrices * * This class is the return value of kroneckerProduct(MatrixBase, @@ -80,6 +98,8 @@ class KroneckerProduct : public KroneckerProductBase > }; /*! + * \ingroup KroneckerProduct_Module + * * \brief Kronecker tensor product helper class for sparse matrices * * If at least one of the operands is a sparse matrix expression, From 51573da3a4407288b900e8764e20de4fc09018d3 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sun, 21 Jul 2013 01:00:36 +0800 Subject: [PATCH 32/34] Warn about power of a matrix with non-semisimple 0 eigenvalue. --- unsupported/Eigen/MatrixFunctions | 4 ++++ .../doc/examples/MatrixPower_failure.cpp | 20 +++++++++++++++++++ 2 files changed, 24 insertions(+) create mode 100644 unsupported/doc/examples/MatrixPower_failure.cpp diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index 41dfab390..ff466a519 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -262,6 +262,10 @@ where \f$ T_1 \f$ is invertible. Then \f$ T^p \f$ is given by 0 & 0 \end{array}. \right] \f] +\warning Fractional power of a matrix with a non-semisimple zero +eigenvalue is not well-defined. We introduce an assertion failure +against inaccurate result, e.g. \include MatrixPower_failure.cpp + Details of the algorithm can be found in: Nicholas J. Higham and Lijing Lin, "A Schur-Padé algorithm for fractional powers of a matrix," SIAM J. %Matrix Anal. Applic., diff --git a/unsupported/doc/examples/MatrixPower_failure.cpp b/unsupported/doc/examples/MatrixPower_failure.cpp new file mode 100644 index 000000000..d20de78a0 --- /dev/null +++ b/unsupported/doc/examples/MatrixPower_failure.cpp @@ -0,0 +1,20 @@ +#include +#include + +int main() +{ + Eigen::Matrix4d A; + A << 0, 0, 2, 3, + 0, 0, 4, 5, + 0, 0, 6, 7, + 0, 0, 8, 9; + std::cout << A.pow(0.37) << std::endl; + + // The 1 makes eigenvalue 0 non-semisimple. + A.coeffRef(0, 1) = 1; + + // This fails if EIGEN_NO_DEBUG is undefined. + std::cout << A.pow(0.37) << std::endl; + + return 0; +} From c00f688c643038650f941f786227a76cd6fd2310 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sun, 21 Jul 2013 05:40:56 +0800 Subject: [PATCH 33/34] Fix doc. (It is also used by computeFracPower) --- unsupported/Eigen/src/MatrixFunctions/MatrixPower.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index b419e245b..5548bd95c 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -389,7 +389,7 @@ class MatrixPower : internal::noncopyable /** \brief Reference to the base of matrix power. */ typename MatrixType::Nested m_A; - /** \brief Temporary storage for computing integral power. */ + /** \brief Temporary storage. */ MatrixType m_tmp; /** \brief Store the result of Schur decomposition. */ From 01190b3544cd0a674be6475185d5dd8e4b7890c5 Mon Sep 17 00:00:00 2001 From: Chen-Pang He Date: Sun, 21 Jul 2013 18:09:11 +0800 Subject: [PATCH 34/34] Directly code failing example, or it breaks `make doc`. --- unsupported/Eigen/MatrixFunctions | 23 ++++++++++++++++++- .../doc/examples/MatrixPower_failure.cpp | 20 ---------------- 2 files changed, 22 insertions(+), 21 deletions(-) delete mode 100644 unsupported/doc/examples/MatrixPower_failure.cpp diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index ff466a519..0b12aaffb 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -264,7 +264,28 @@ where \f$ T_1 \f$ is invertible. Then \f$ T^p \f$ is given by \warning Fractional power of a matrix with a non-semisimple zero eigenvalue is not well-defined. We introduce an assertion failure -against inaccurate result, e.g. \include MatrixPower_failure.cpp +against inaccurate result, e.g. \code +#include +#include + +int main() +{ + Eigen::Matrix4d A; + A << 0, 0, 2, 3, + 0, 0, 4, 5, + 0, 0, 6, 7, + 0, 0, 8, 9; + std::cout << A.pow(0.37) << std::endl; + + // The 1 makes eigenvalue 0 non-semisimple. + A.coeffRef(0, 1) = 1; + + // This fails if EIGEN_NO_DEBUG is undefined. + std::cout << A.pow(0.37) << std::endl; + + return 0; +} +\endcode Details of the algorithm can be found in: Nicholas J. Higham and Lijing Lin, "A Schur-Padé algorithm for fractional powers of a diff --git a/unsupported/doc/examples/MatrixPower_failure.cpp b/unsupported/doc/examples/MatrixPower_failure.cpp deleted file mode 100644 index d20de78a0..000000000 --- a/unsupported/doc/examples/MatrixPower_failure.cpp +++ /dev/null @@ -1,20 +0,0 @@ -#include -#include - -int main() -{ - Eigen::Matrix4d A; - A << 0, 0, 2, 3, - 0, 0, 4, 5, - 0, 0, 6, 7, - 0, 0, 8, 9; - std::cout << A.pow(0.37) << std::endl; - - // The 1 makes eigenvalue 0 non-semisimple. - A.coeffRef(0, 1) = 1; - - // This fails if EIGEN_NO_DEBUG is undefined. - std::cout << A.pow(0.37) << std::endl; - - return 0; -}