diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h index b343f38bc2..6970b3a15c 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2009, 2010 Jitse Niesen +// Copyright (C) 2009-2011 Jitse Niesen // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -30,30 +30,36 @@ /** \ingroup MatrixFunctions_Module - * \brief Class for computing matrix exponentials. - * \tparam MatrixType type of the argument of the matrix function, - * expected to be an instantiation of the Matrix class template. + * \brief Class for computing matrix functions. + * \tparam MatrixType type of the argument of the matrix function, + * expected to be an instantiation of the Matrix class template. + * \tparam AtomicType type for computing matrix function of atomic blocks. + * \tparam IsComplex used internally to select correct specialization. + * + * This class implements the Schur-Parlett algorithm for computing matrix functions. The spectrum of the + * matrix is divided in clustered of eigenvalues that lies close together. This class delegates the + * computation of the matrix function on every block corresponding to these clusters to an object of type + * \p AtomicType and uses these results to compute the matrix function of the whole matrix. The class + * \p AtomicType should have a \p compute() member function for computing the matrix function of a block. + * + * \sa class MatrixFunctionAtomic, class MatrixLogarithmAtomic */ -template ::Scalar>::IsComplex> +template ::Scalar>::IsComplex> class MatrixFunction { - private: - - typedef typename internal::traits::Index Index; - typedef typename internal::traits::Scalar Scalar; - typedef typename internal::stem_function::type StemFunction; - public: /** \brief Constructor. * - * \param[in] A argument of matrix function, should be a square matrix. - * \param[in] f an entire function; \c f(x,n) should compute the n-th derivative of f at x. + * \param[in] A argument of matrix function, should be a square matrix. + * \param[in] atomic class for computing matrix function of atomic blocks. * - * The class stores a reference to \p A, so it should not be + * The class stores references to \p A and \p atomic, so they should not be * changed (or destroyed) before compute() is called. */ - MatrixFunction(const MatrixType& A, StemFunction f); + MatrixFunction(const MatrixType& A, AtomicType& atomic); /** \brief Compute the matrix function. * @@ -71,8 +77,8 @@ class MatrixFunction /** \internal \ingroup MatrixFunctions_Module * \brief Partial specialization of MatrixFunction for real matrices */ -template -class MatrixFunction +template +class MatrixFunction { private: @@ -86,16 +92,15 @@ class MatrixFunction typedef std::complex ComplexScalar; typedef Matrix ComplexMatrix; - typedef typename internal::stem_function::type StemFunction; public: /** \brief Constructor. * - * \param[in] A argument of matrix function, should be a square matrix. - * \param[in] f an entire function; \c f(x,n) should compute the n-th derivative of f at x. + * \param[in] A argument of matrix function, should be a square matrix. + * \param[in] atomic class for computing matrix function of atomic blocks. */ - MatrixFunction(const MatrixType& A, StemFunction f) : m_A(A), m_f(f) { } + MatrixFunction(const MatrixType& A, AtomicType& atomic) : m_A(A), m_atomic(atomic) { } /** \brief Compute the matrix function. * @@ -111,14 +116,14 @@ class MatrixFunction { ComplexMatrix CA = m_A.template cast(); ComplexMatrix Cresult; - MatrixFunction mf(CA, m_f); + MatrixFunction mf(CA, m_atomic); mf.compute(Cresult); result = Cresult.real(); } private: typename internal::nested::type m_A; /**< \brief Reference to argument of matrix function. */ - StemFunction *m_f; /**< \brief Stem function for matrix function under consideration */ + AtomicType& m_atomic; /**< \brief Class for computing matrix function of atomic blocks. */ MatrixFunction& operator=(const MatrixFunction&); }; @@ -127,8 +132,8 @@ class MatrixFunction /** \internal \ingroup MatrixFunctions_Module * \brief Partial specialization of MatrixFunction for complex matrices */ -template -class MatrixFunction +template +class MatrixFunction { private: @@ -139,7 +144,6 @@ class MatrixFunction static const int ColsAtCompileTime = Traits::ColsAtCompileTime; static const int Options = MatrixType::Options; typedef typename NumTraits::Real RealScalar; - typedef typename internal::stem_function::type StemFunction; typedef Matrix VectorType; typedef Matrix IntVectorType; typedef Matrix DynamicIntVectorType; @@ -149,7 +153,7 @@ class MatrixFunction public: - MatrixFunction(const MatrixType& A, StemFunction f); + MatrixFunction(const MatrixType& A, AtomicType& atomic); template void compute(ResultType& result); private: @@ -168,7 +172,7 @@ class MatrixFunction DynMatrixType solveTriangularSylvester(const DynMatrixType& A, const DynMatrixType& B, const DynMatrixType& C); typename internal::nested::type m_A; /**< \brief Reference to argument of matrix function. */ - StemFunction *m_f; /**< \brief Stem function for matrix function under consideration */ + AtomicType& m_atomic; /**< \brief Class for computing matrix function of atomic blocks. */ MatrixType m_T; /**< \brief Triangular part of Schur decomposition */ MatrixType m_U; /**< \brief Unitary part of Schur decomposition */ MatrixType m_fT; /**< \brief %Matrix function applied to #m_T */ @@ -191,12 +195,12 @@ class MatrixFunction /** \brief Constructor. * - * \param[in] A argument of matrix function, should be a square matrix. - * \param[in] f an entire function; \c f(x,n) should compute the n-th derivative of f at x. + * \param[in] A argument of matrix function, should be a square matrix. + * \param[in] atomic class for computing matrix function of atomic blocks. */ -template -MatrixFunction::MatrixFunction(const MatrixType& A, StemFunction f) : - m_A(A), m_f(f) +template +MatrixFunction::MatrixFunction(const MatrixType& A, AtomicType& atomic) + : m_A(A), m_atomic(atomic) { /* empty body */ } @@ -206,9 +210,9 @@ MatrixFunction::MatrixFunction(const MatrixType& A, StemFunction f * \param[out] result the function \p f applied to \p A, as * specified in the constructor. */ -template +template template -void MatrixFunction::compute(ResultType& result) +void MatrixFunction::compute(ResultType& result) { computeSchurDecomposition(); partitionEigenvalues(); @@ -222,8 +226,8 @@ void MatrixFunction::compute(ResultType& result) } /** \brief Store the Schur decomposition of #m_A in #m_T and #m_U */ -template -void MatrixFunction::computeSchurDecomposition() +template +void MatrixFunction::computeSchurDecomposition() { const ComplexSchur schurOfA(m_A); m_T = schurOfA.matrixT(); @@ -241,8 +245,8 @@ void MatrixFunction::computeSchurDecomposition() * The implementation follows Algorithm 4.1 in the paper of Davies * and Higham. */ -template -void MatrixFunction::partitionEigenvalues() +template +void MatrixFunction::partitionEigenvalues() { const Index rows = m_T.rows(); VectorType diag = m_T.diagonal(); // contains eigenvalues of A @@ -278,8 +282,8 @@ void MatrixFunction::partitionEigenvalues() * \returns Iterator to cluster containing \c key, or * \c m_clusters.end() if no cluster in m_clusters contains \c key. */ -template -typename MatrixFunction::ListOfClusters::iterator MatrixFunction::findCluster(Scalar key) +template +typename MatrixFunction::ListOfClusters::iterator MatrixFunction::findCluster(Scalar key) { typename Cluster::iterator j; for (typename ListOfClusters::iterator i = m_clusters.begin(); i != m_clusters.end(); ++i) { @@ -291,8 +295,8 @@ typename MatrixFunction::ListOfClusters::iterator MatrixFunction -void MatrixFunction::computeClusterSize() +template +void MatrixFunction::computeClusterSize() { const Index rows = m_T.rows(); VectorType diag = m_T.diagonal(); @@ -313,8 +317,8 @@ void MatrixFunction::computeClusterSize() } /** \brief Compute #m_blockStart using #m_clusterSize */ -template -void MatrixFunction::computeBlockStart() +template +void MatrixFunction::computeBlockStart() { m_blockStart.resize(m_clusterSize.rows()); m_blockStart(0) = 0; @@ -324,8 +328,8 @@ void MatrixFunction::computeBlockStart() } /** \brief Compute #m_permutation using #m_eivalToCluster and #m_blockStart */ -template -void MatrixFunction::constructPermutation() +template +void MatrixFunction::constructPermutation() { DynamicIntVectorType indexNextEntry = m_blockStart; m_permutation.resize(m_T.rows()); @@ -337,8 +341,8 @@ void MatrixFunction::constructPermutation() } /** \brief Permute Schur decomposition in #m_U and #m_T according to #m_permutation */ -template -void MatrixFunction::permuteSchur() +template +void MatrixFunction::permuteSchur() { IntVectorType p = m_permutation; for (Index i = 0; i < p.rows() - 1; i++) { @@ -355,8 +359,8 @@ void MatrixFunction::permuteSchur() } /** \brief Swap rows \a index and \a index+1 in Schur decomposition in #m_U and #m_T */ -template -void MatrixFunction::swapEntriesInSchur(Index index) +template +void MatrixFunction::swapEntriesInSchur(Index index) { JacobiRotation rotation; rotation.makeGivens(m_T(index, index+1), m_T(index+1, index+1) - m_T(index, index)); @@ -367,25 +371,23 @@ void MatrixFunction::swapEntriesInSchur(Index index) /** \brief Compute block diagonal part of #m_fT. * - * This routine computes the matrix function #m_f applied to the block - * diagonal part of #m_T, with the blocking given by #m_blockStart. The - * result is stored in #m_fT. The off-diagonal parts of #m_fT are set - * to zero. + * This routine computes the matrix function applied to the block diagonal part of #m_T, with the blocking + * given by #m_blockStart. The matrix function of each diagonal block is computed by #m_atomic. The + * off-diagonal parts of #m_fT are set to zero. */ -template -void MatrixFunction::computeBlockAtomic() +template +void MatrixFunction::computeBlockAtomic() { m_fT.resize(m_T.rows(), m_T.cols()); m_fT.setZero(); - MatrixFunctionAtomic mfa(m_f); for (Index i = 0; i < m_clusterSize.rows(); ++i) { - block(m_fT, i, i) = mfa.compute(block(m_T, i, i)); + block(m_fT, i, i) = m_atomic.compute(block(m_T, i, i)); } } /** \brief Return block of matrix according to blocking given by #m_blockStart */ -template -Block MatrixFunction::block(MatrixType& A, Index i, Index j) +template +Block MatrixFunction::block(MatrixType& A, Index i, Index j) { return A.block(m_blockStart(i), m_blockStart(j), m_clusterSize(i), m_clusterSize(j)); } @@ -393,12 +395,12 @@ Block MatrixFunction::block(MatrixType& A, Index i, In /** \brief Compute part of #m_fT above block diagonal. * * This routine assumes that the block diagonal part of #m_fT (which - * equals #m_f applied to #m_T) has already been computed and computes + * equals the matrix function applied to #m_T) has already been computed and computes * the part above the block diagonal. The part below the diagonal is * zero, because #m_T is upper triangular. */ -template -void MatrixFunction::computeOffDiagonal() +template +void MatrixFunction::computeOffDiagonal() { for (Index diagIndex = 1; diagIndex < m_clusterSize.rows(); diagIndex++) { for (Index blockIndex = 0; blockIndex < m_clusterSize.rows() - diagIndex; blockIndex++) { @@ -439,8 +441,8 @@ void MatrixFunction::computeOffDiagonal() * solution). In that case, these equations can be evaluated in the * order \f$ i=m,\ldots,1 \f$ and \f$ j=1,\ldots,n \f$. */ -template -typename MatrixFunction::DynMatrixType MatrixFunction::solveTriangularSylvester( +template +typename MatrixFunction::DynMatrixType MatrixFunction::solveTriangularSylvester( const DynMatrixType& A, const DynMatrixType& B, const DynMatrixType& C) @@ -520,8 +522,18 @@ template class MatrixFunctionReturnValue template inline void evalTo(ResultType& result) const { - const typename Derived::PlainObject Aevaluated = m_A.eval(); - MatrixFunction mf(Aevaluated, m_f); + typedef typename Derived::PlainObject PlainObject; + typedef internal::traits Traits; + static const int RowsAtCompileTime = Traits::RowsAtCompileTime; + static const int ColsAtCompileTime = Traits::ColsAtCompileTime; + static const int Options = PlainObject::Options; + typedef std::complex::Real> ComplexScalar; + typedef Matrix DynMatrixType; + typedef MatrixFunctionAtomic AtomicType; + AtomicType atomic(m_f); + + const PlainObject Aevaluated = m_A.eval(); + MatrixFunction mf(Aevaluated, atomic); mf.compute(result); }