From a3a55357db7394281c872e911f13d69aba510aec Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Thu, 25 Jul 2013 15:08:53 +0100 Subject: Clean up MatrixFunction and MatrixLogarithm. --- .../Eigen/src/MatrixFunctions/MatrixFunction.h | 678 ++++++++++----------- .../src/MatrixFunctions/MatrixFunctionAtomic.h | 127 ---- .../Eigen/src/MatrixFunctions/MatrixLogarithm.h | 456 ++++++-------- .../Eigen/src/MatrixFunctions/MatrixSquareRoot.h | 2 +- 4 files changed, 509 insertions(+), 754 deletions(-) delete mode 100644 unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h index ad54d8ed0..12e28793d 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-2011 Jitse Niesen +// Copyright (C) 2009-2011, 2013 Jitse Niesen // // 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 @@ -11,394 +11,244 @@ #define EIGEN_MATRIX_FUNCTION #include "StemFunction.h" -#include "MatrixFunctionAtomic.h" namespace Eigen { +namespace internal { + +/** \brief Maximum distance allowed between eigenvalues to be considered "close". */ +static const float matrix_function_separation = 0.1; + /** \ingroup MatrixFunctions_Module - * \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. + * \class MatrixFunctionAtomic + * \brief Helper class for computing matrix functions of atomic matrices. * - * 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 + * Here, an atomic matrix is a triangular matrix whose diagonal entries are close to each other. */ -template ::Scalar>::IsComplex> -class MatrixFunction : internal::noncopyable -{ +template +class MatrixFunctionAtomic +{ public: - /** \brief Constructor. - * - * \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 references to \p A and \p atomic, so they should not be - * changed (or destroyed) before compute() is called. - */ - MatrixFunction(const MatrixType& A, AtomicType& atomic); - - /** \brief Compute the matrix function. - * - * \param[out] result the function \p f applied to \p A, as - * specified in the constructor. - * - * See MatrixBase::matrixFunction() for details on how this computation - * is implemented. - */ - template - void compute(ResultType &result); -}; - - -/** \internal \ingroup MatrixFunctions_Module - * \brief Partial specialization of MatrixFunction for real matrices - */ -template -class MatrixFunction : internal::noncopyable -{ - private: - - typedef internal::traits Traits; - typedef typename Traits::Scalar Scalar; - static const int Rows = Traits::RowsAtCompileTime; - static const int Cols = Traits::ColsAtCompileTime; - static const int Options = MatrixType::Options; - static const int MaxRows = Traits::MaxRowsAtCompileTime; - static const int MaxCols = Traits::MaxColsAtCompileTime; - - typedef std::complex ComplexScalar; - typedef Matrix ComplexMatrix; - - public: + typedef typename MatrixType::Scalar Scalar; + typedef typename stem_function::type StemFunction; - /** \brief Constructor. - * - * \param[in] A argument of matrix function, should be a square matrix. - * \param[in] atomic class for computing matrix function of atomic blocks. + /** \brief Constructor + * \param[in] f matrix function to compute. */ - MatrixFunction(const MatrixType& A, AtomicType& atomic) : m_A(A), m_atomic(atomic) { } + MatrixFunctionAtomic(StemFunction f) : m_f(f) { } - /** \brief Compute the matrix function. - * - * \param[out] result the function \p f applied to \p A, as - * specified in the constructor. - * - * This function converts the real matrix \c A to a complex matrix, - * uses MatrixFunction and then converts the result back to - * a real matrix. + /** \brief Compute matrix function of atomic matrix + * \param[in] A argument of matrix function, should be upper triangular and atomic + * \returns f(A), the matrix function evaluated at the given matrix */ - template - void compute(ResultType& result) - { - ComplexMatrix CA = m_A.template cast(); - ComplexMatrix Cresult; - 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. */ - AtomicType& m_atomic; /**< \brief Class for computing matrix function of atomic blocks. */ -}; - - -/** \internal \ingroup MatrixFunctions_Module - * \brief Partial specialization of MatrixFunction for complex matrices - */ -template -class MatrixFunction : internal::noncopyable -{ - private: - - typedef internal::traits Traits; - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::Index Index; - static const int RowsAtCompileTime = Traits::RowsAtCompileTime; - static const int ColsAtCompileTime = Traits::ColsAtCompileTime; - static const int Options = MatrixType::Options; - typedef typename NumTraits::Real RealScalar; - typedef Matrix VectorType; - typedef Matrix IntVectorType; - typedef Matrix DynamicIntVectorType; - typedef std::list Cluster; - typedef std::list ListOfClusters; - typedef Matrix DynMatrixType; - - public: - - MatrixFunction(const MatrixType& A, AtomicType& atomic); - template void compute(ResultType& result); + MatrixType compute(const MatrixType& A); private: - - void computeSchurDecomposition(); - void partitionEigenvalues(); - typename ListOfClusters::iterator findCluster(Scalar key); - void computeClusterSize(); - void computeBlockStart(); - void constructPermutation(); - void permuteSchur(); - void swapEntriesInSchur(Index index); - void computeBlockAtomic(); - Block block(MatrixType& A, Index i, Index j); - void computeOffDiagonal(); - DynMatrixType solveTriangularSylvester(const DynMatrixType& A, const DynMatrixType& B, const DynMatrixType& C); - - 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. */ - 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 */ - ListOfClusters m_clusters; /**< \brief Partition of eigenvalues into clusters of ei'vals "close" to each other */ - DynamicIntVectorType m_eivalToCluster; /**< \brief m_eivalToCluster[i] = j means i-th ei'val is in j-th cluster */ - DynamicIntVectorType m_clusterSize; /**< \brief Number of eigenvalues in each clusters */ - DynamicIntVectorType m_blockStart; /**< \brief Row index at which block corresponding to i-th cluster starts */ - IntVectorType m_permutation; /**< \brief Permutation which groups ei'vals in the same cluster together */ - - /** \brief Maximum distance allowed between eigenvalues to be considered "close". - * - * This is morally a \c static \c const \c Scalar, but only - * integers can be static constant class members in C++. The - * 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); } + StemFunction* m_f; }; -/** \brief Constructor. - * - * \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, AtomicType& atomic) - : m_A(A), m_atomic(atomic) +template +typename NumTraits::Real matrix_function_compute_mu(const MatrixType& A) { - /* empty body */ + typedef typename plain_col_type::type VectorType; + typename MatrixType::Index rows = A.rows(); + const MatrixType N = MatrixType::Identity(rows, rows) - A; + VectorType e = VectorType::Ones(rows); + N.template triangularView().solveInPlace(e); + return e.cwiseAbs().maxCoeff(); } -/** \brief Compute the matrix function. - * - * \param[out] result the function \p f applied to \p A, as - * specified in the constructor. - */ -template -template -void MatrixFunction::compute(ResultType& result) +template +MatrixType MatrixFunctionAtomic::compute(const MatrixType& A) { - computeSchurDecomposition(); - partitionEigenvalues(); - computeClusterSize(); - computeBlockStart(); - constructPermutation(); - permuteSchur(); - computeBlockAtomic(); - computeOffDiagonal(); - result = m_U * (m_fT.template triangularView() * m_U.adjoint()); + // TODO: Use that A is upper triangular + typedef typename NumTraits::Real RealScalar; + typedef typename MatrixType::Index Index; + Index rows = A.rows(); + Scalar avgEival = A.trace() / Scalar(RealScalar(rows)); + MatrixType Ashifted = A - avgEival * MatrixType::Identity(rows, rows); + RealScalar mu = matrix_function_compute_mu(Ashifted); + MatrixType F = m_f(avgEival, 0) * MatrixType::Identity(rows, rows); + MatrixType P = Ashifted; + MatrixType Fincr; + for (Index s = 1; s < 1.1 * rows + 10; s++) { // upper limit is fairly arbitrary + Fincr = m_f(avgEival, static_cast(s)) * P; + F += Fincr; + P = Scalar(RealScalar(1.0/(s + 1))) * P * Ashifted; + + // test whether Taylor series converged + const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff(); + const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff(); + if (Fincr_norm < NumTraits::epsilon() * F_norm) { + RealScalar delta = 0; + RealScalar rfactorial = 1; + for (Index r = 0; r < rows; r++) { + RealScalar mx = 0; + for (Index i = 0; i < rows; i++) + mx = (std::max)(mx, std::abs(m_f(Ashifted(i, i) + avgEival, static_cast(s+r)))); + if (r != 0) + rfactorial *= RealScalar(r); + delta = (std::max)(delta, mx / rfactorial); + } + const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff(); + if (mu * delta * P_norm < NumTraits::epsilon() * F_norm) // series converged + break; + } + } + return F; } -/** \brief Store the Schur decomposition of #m_A in #m_T and #m_U */ -template -void MatrixFunction::computeSchurDecomposition() +/** \brief Find cluster in \p clusters containing some value + * \param[in] key Value to find + * \returns Iterator to cluster containing \p key, or \c clusters.end() if no cluster in \p m_clusters + * contains \p key. + */ +template +typename ListOfClusters::iterator matrix_function_find_cluster(Scalar key, ListOfClusters& clusters) { - const ComplexSchur schurOfA(m_A); - m_T = schurOfA.matrixT(); - m_U = schurOfA.matrixU(); + typename std::list::iterator j; + for (typename ListOfClusters::iterator i = clusters.begin(); i != clusters.end(); ++i) { + j = std::find(i->begin(), i->end(), key); + if (j != i->end()) + return i; + } + return clusters.end(); } /** \brief Partition eigenvalues in clusters of ei'vals close to each other * - * This function computes #m_clusters. This is a partition of the - * eigenvalues of #m_T in clusters, such that - * # Any eigenvalue in a certain cluster is at most separation() away - * from another eigenvalue in the same cluster. - * # The distance between two eigenvalues in different clusters is - * more than separation(). - * The implementation follows Algorithm 4.1 in the paper of Davies - * and Higham. + * \param[in] eivals Eigenvalues + * \param[out] clusters Resulting partition of eigenvalues + * + * The partition satisfies the following two properties: + * # Any eigenvalue in a certain cluster is at most matrix_function_separation() away from another eigenvalue + * in the same cluster. + * # The distance between two eigenvalues in different clusters is more than matrix_function_separation(). + * The implementation follows Algorithm 4.1 in the paper of Davies and Higham. */ -template -void MatrixFunction::partitionEigenvalues() +template +void matrix_function_partition_eigenvalues(const EivalsType& eivals, std::list& clusters) { - using std::abs; - const Index rows = m_T.rows(); - VectorType diag = m_T.diagonal(); // contains eigenvalues of A - - for (Index i=0; i::iterator qi = matrix_function_find_cluster(eivals(i), clusters); + if (qi == clusters.end()) { Cluster l; - l.push_back(diag(i)); - m_clusters.push_back(l); - qi = m_clusters.end(); + l.push_back(eivals(i)); + clusters.push_back(l); + qi = clusters.end(); --qi; } // Look for other element to add to the set - for (Index j=i+1; jbegin(), qi->end(), diag(j)) == qi->end()) { - typename ListOfClusters::iterator qj = findCluster(diag(j)); - if (qj == m_clusters.end()) { - qi->push_back(diag(j)); + for (Index j=i+1; jbegin(), qi->end(), eivals(j)) == qi->end()) { + typename std::list::iterator qj = matrix_function_find_cluster(eivals(j), clusters); + if (qj == clusters.end()) { + qi->push_back(eivals(j)); } else { qi->insert(qi->end(), qj->begin(), qj->end()); - m_clusters.erase(qj); + clusters.erase(qj); } } } } } -/** \brief Find cluster in #m_clusters containing some value - * \param[in] key Value to find - * \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) +/** \brief Compute size of each cluster given a partitioning */ +template +void matrix_function_compute_cluster_size(const ListOfClusters& clusters, Matrix& clusterSize) { - typename Cluster::iterator j; - for (typename ListOfClusters::iterator i = m_clusters.begin(); i != m_clusters.end(); ++i) { - j = std::find(i->begin(), i->end(), key); - if (j != i->end()) - return i; + const Index numClusters = static_cast(clusters.size()); + clusterSize.setZero(numClusters); + Index clusterIndex = 0; + for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) { + clusterSize[clusterIndex] = cluster->size(); + ++clusterIndex; } - return m_clusters.end(); } -/** \brief Compute #m_clusterSize and #m_eivalToCluster using #m_clusters */ -template -void MatrixFunction::computeClusterSize() +/** \brief Compute start of each block using clusterSize */ +template +void matrix_function_compute_block_start(const VectorType& clusterSize, VectorType& blockStart) { - const Index rows = m_T.rows(); - VectorType diag = m_T.diagonal(); - const Index numClusters = static_cast(m_clusters.size()); + blockStart.resize(clusterSize.rows()); + blockStart(0) = 0; + for (typename VectorType::Index i = 1; i < clusterSize.rows(); i++) { + blockStart(i) = blockStart(i-1) + clusterSize(i-1); + } +} - m_clusterSize.setZero(numClusters); - m_eivalToCluster.resize(rows); +/** \brief Compute mapping of eigenvalue indices to cluster indices */ +template +void matrix_function_compute_map(const EivalsType& eivals, const ListOfClusters& clusters, VectorType& eivalToCluster) +{ + typedef typename EivalsType::Index Index; + eivalToCluster.resize(eivals.rows()); Index clusterIndex = 0; - for (typename ListOfClusters::const_iterator cluster = m_clusters.begin(); cluster != m_clusters.end(); ++cluster) { - for (Index i = 0; i < diag.rows(); ++i) { - if (std::find(cluster->begin(), cluster->end(), diag(i)) != cluster->end()) { - ++m_clusterSize[clusterIndex]; - m_eivalToCluster[i] = clusterIndex; + for (typename ListOfClusters::const_iterator cluster = clusters.begin(); cluster != clusters.end(); ++cluster) { + for (Index i = 0; i < eivals.rows(); ++i) { + if (std::find(cluster->begin(), cluster->end(), eivals(i)) != cluster->end()) { + eivalToCluster[i] = clusterIndex; } } ++clusterIndex; } } -/** \brief Compute #m_blockStart using #m_clusterSize */ -template -void MatrixFunction::computeBlockStart() +/** \brief Compute permutation which groups ei'vals in same cluster together */ +template +void matrix_function_compute_permutation(const DynVectorType& blockStart, const DynVectorType& eivalToCluster, VectorType& permutation) { - m_blockStart.resize(m_clusterSize.rows()); - m_blockStart(0) = 0; - for (Index i = 1; i < m_clusterSize.rows(); i++) { - m_blockStart(i) = m_blockStart(i-1) + m_clusterSize(i-1); - } -} - -/** \brief Compute #m_permutation using #m_eivalToCluster and #m_blockStart */ -template -void MatrixFunction::constructPermutation() -{ - DynamicIntVectorType indexNextEntry = m_blockStart; - m_permutation.resize(m_T.rows()); - for (Index i = 0; i < m_T.rows(); i++) { - Index cluster = m_eivalToCluster[i]; - m_permutation[i] = indexNextEntry[cluster]; + typedef typename VectorType::Index Index; + DynVectorType indexNextEntry = blockStart; + permutation.resize(eivalToCluster.rows()); + for (Index i = 0; i < eivalToCluster.rows(); i++) { + Index cluster = eivalToCluster[i]; + permutation[i] = indexNextEntry[cluster]; ++indexNextEntry[cluster]; } } -/** \brief Permute Schur decomposition in #m_U and #m_T according to #m_permutation */ -template -void MatrixFunction::permuteSchur() +/** \brief Permute Schur decomposition in U and T according to permutation */ +template +void matrix_function_permute_schur(VectorType& permutation, MatrixType& U, MatrixType& T) { - IntVectorType p = m_permutation; - for (Index i = 0; i < p.rows() - 1; i++) { + typedef typename VectorType::Index Index; + for (Index i = 0; i < permutation.rows() - 1; i++) { Index j; - for (j = i; j < p.rows(); j++) { - if (p(j) == i) break; + for (j = i; j < permutation.rows(); j++) { + if (permutation(j) == i) break; } - eigen_assert(p(j) == i); + eigen_assert(permutation(j) == i); for (Index k = j-1; k >= i; k--) { - swapEntriesInSchur(k); - std::swap(p.coeffRef(k), p.coeffRef(k+1)); + JacobiRotation rotation; + rotation.makeGivens(T(k, k+1), T(k+1, k+1) - T(k, k)); + T.applyOnTheLeft(k, k+1, rotation.adjoint()); + T.applyOnTheRight(k, k+1, rotation); + U.applyOnTheRight(k, k+1, rotation); + std::swap(permutation.coeffRef(k), permutation.coeffRef(k+1)); } } } -/** \brief Swap rows \a index and \a index+1 in Schur decomposition in #m_U and #m_T */ -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)); - m_T.applyOnTheLeft(index, index+1, rotation.adjoint()); - m_T.applyOnTheRight(index, index+1, rotation); - m_U.applyOnTheRight(index, index+1, rotation); -} - -/** \brief Compute block diagonal part of #m_fT. - * - * 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() -{ - m_fT.resize(m_T.rows(), m_T.cols()); - m_fT.setZero(); - for (Index i = 0; i < m_clusterSize.rows(); ++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) -{ - return A.block(m_blockStart(i), m_blockStart(j), m_clusterSize(i), m_clusterSize(j)); -} - -/** \brief Compute part of #m_fT above block diagonal. +/** \brief Compute block diagonal part of matrix function. * - * This routine assumes that the block diagonal part of #m_fT (which - * 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. + * This routine computes the matrix function applied to the block diagonal part of \p T (which should be + * upper triangular), with the blocking given by \p blockStart and \p clusterSize. The matrix function of + * each diagonal block is computed by \p atomic. The off-diagonal parts of \p fT are set to zero. */ -template -void MatrixFunction::computeOffDiagonal() +template +void matrix_function_compute_block_atomic(const MatrixType& T, AtomicType& atomic, const VectorType& blockStart, const VectorType& clusterSize, MatrixType& fT) { - for (Index diagIndex = 1; diagIndex < m_clusterSize.rows(); diagIndex++) { - for (Index blockIndex = 0; blockIndex < m_clusterSize.rows() - diagIndex; blockIndex++) { - // compute (blockIndex, blockIndex+diagIndex) block - DynMatrixType A = block(m_T, blockIndex, blockIndex); - DynMatrixType B = -block(m_T, blockIndex+diagIndex, blockIndex+diagIndex); - DynMatrixType C = block(m_fT, blockIndex, blockIndex) * block(m_T, blockIndex, blockIndex+diagIndex); - C -= block(m_T, blockIndex, blockIndex+diagIndex) * block(m_fT, blockIndex+diagIndex, blockIndex+diagIndex); - for (Index k = blockIndex + 1; k < blockIndex + diagIndex; k++) { - C += block(m_fT, blockIndex, k) * block(m_T, k, blockIndex+diagIndex); - C -= block(m_T, blockIndex, k) * block(m_fT, k, blockIndex+diagIndex); - } - block(m_fT, blockIndex, blockIndex+diagIndex) = solveTriangularSylvester(A, B, C); - } + fT.setZero(T.rows(), T.cols()); + for (typename VectorType::Index i = 0; i < clusterSize.rows(); ++i) { + fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)) + = atomic.compute(T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i))); } } @@ -410,8 +260,8 @@ void MatrixFunction::computeOffDiagonal() * * \returns the solution X. * - * If A is m-by-m and B is n-by-n, then both C and X are m-by-n. - * The (i,j)-th component of the Sylvester equation is + * If A is m-by-m and B is n-by-n, then both C and X are m-by-n. The (i,j)-th component of the Sylvester + * equation is * \f[ * \sum_{k=i}^m A_{ik} X_{kj} + \sum_{k=1}^j X_{ik} B_{kj} = C_{ij}. * \f] @@ -420,16 +270,12 @@ void MatrixFunction::computeOffDiagonal() * X_{ij} = \frac{1}{A_{ii} + B_{jj}} \Bigl( C_{ij} * - \sum_{k=i+1}^m A_{ik} X_{kj} - \sum_{k=1}^{j-1} X_{ik} B_{kj} \Bigr). * \f] - * It is assumed that A and B are such that the numerator is never - * zero (otherwise the Sylvester equation does not have a unique - * 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$. + * It is assumed that A and B are such that the numerator is never zero (otherwise the Sylvester equation + * does not have a unique 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( - const DynMatrixType& A, - const DynMatrixType& B, - const DynMatrixType& C) +template +MatrixType matrix_function_solve_triangular_sylvester(const MatrixType& A, const MatrixType& B, const MatrixType& C) { eigen_assert(A.rows() == A.cols()); eigen_assert(A.isUpperTriangular()); @@ -438,9 +284,12 @@ typename MatrixFunction::DynMatrixType MatrixFunction= 0; --i) { for (Index j = 0; j < n; ++j) { @@ -469,17 +318,164 @@ typename MatrixFunction::DynMatrixType MatrixFunction +void matrix_function_compute_above_diagonal(const MatrixType& T, const VectorType& blockStart, const VectorType& clusterSize, MatrixType& fT) +{ + typedef internal::traits Traits; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::Index Index; + static const int RowsAtCompileTime = Traits::RowsAtCompileTime; + static const int ColsAtCompileTime = Traits::ColsAtCompileTime; + static const int Options = MatrixType::Options; + typedef Matrix DynMatrixType; + + for (Index k = 1; k < clusterSize.rows(); k++) { + for (Index i = 0; i < clusterSize.rows() - k; i++) { + // compute (i, i+k) block + DynMatrixType A = T.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)); + DynMatrixType B = -T.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k)); + DynMatrixType C = fT.block(blockStart(i), blockStart(i), clusterSize(i), clusterSize(i)) + * T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k)); + C -= T.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k)) + * fT.block(blockStart(i+k), blockStart(i+k), clusterSize(i+k), clusterSize(i+k)); + for (Index m = i + 1; m < i + k; m++) { + C += fT.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m)) + * T.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k)); + C -= T.block(blockStart(i), blockStart(m), clusterSize(i), clusterSize(m)) + * fT.block(blockStart(m), blockStart(i+k), clusterSize(m), clusterSize(i+k)); + } + fT.block(blockStart(i), blockStart(i+k), clusterSize(i), clusterSize(i+k)) + = matrix_function_solve_triangular_sylvester(A, B, C); + } + } +} + +/** \ingroup MatrixFunctions_Module + * \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> +struct matrix_function_compute +{ + /** \brief Compute the matrix function. + * + * \param[in] A argument of matrix function, should be a square matrix. + * \param[in] atomic class for computing matrix function of atomic blocks. + * \param[out] result the function \p f applied to \p A, as + * specified in the constructor. + * + * See MatrixBase::matrixFunction() for details on how this computation + * is implemented. + */ + template + static void run(const MatrixType& A, AtomicType& atomic, ResultType &result); +}; + +/** \internal \ingroup MatrixFunctions_Module + * \brief Partial specialization of MatrixFunction for real matrices + * + * This converts the real matrix to a complex matrix, compute the matrix function of that matrix, and then + * converts the result back to a real matrix. + */ +template +struct matrix_function_compute +{ + template + static void run(const MatrixType& A, AtomicType& atomic, ResultType &result) + { + typedef internal::traits Traits; + typedef typename Traits::Scalar Scalar; + static const int Rows = Traits::RowsAtCompileTime, Cols = Traits::ColsAtCompileTime; + static const int Options = MatrixType::Options; + static const int MaxRows = Traits::MaxRowsAtCompileTime, MaxCols = Traits::MaxColsAtCompileTime; + + typedef std::complex ComplexScalar; + typedef Matrix ComplexMatrix; + + ComplexMatrix CA = A.template cast(); + ComplexMatrix Cresult; + matrix_function_compute::run(CA, atomic, Cresult); + result = Cresult.real(); + } +}; + +/** \internal \ingroup MatrixFunctions_Module + * \brief Partial specialization of MatrixFunction for complex matrices + */ +template +struct matrix_function_compute +{ + template + static void run(const MatrixType& A, AtomicType& atomic, ResultType &result) + { + typedef internal::traits Traits; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::Index Index; + + // compute Schur decomposition of A + const ComplexSchur schurOfA(A); + MatrixType T = schurOfA.matrixT(); + MatrixType U = schurOfA.matrixU(); + + // partition eigenvalues into clusters of ei'vals "close" to each other + std::list > clusters; + matrix_function_partition_eigenvalues(T.diagonal(), clusters); + + // compute size of each cluster + Matrix clusterSize; + matrix_function_compute_cluster_size(clusters, clusterSize); + + // blockStart[i] is row index at which block corresponding to i-th cluster starts + Matrix blockStart; + matrix_function_compute_block_start(clusterSize, blockStart); + + // compute map so that eivalToCluster[i] = j means that ei'val T(i,i) is in j-th cluster + Matrix eivalToCluster; + matrix_function_compute_map(T.diagonal(), clusters, eivalToCluster); + + // compute permutation which groups ei'vals in same cluster together + Matrix permutation; + matrix_function_compute_permutation(blockStart, eivalToCluster, permutation); + + // permute Schur decomposition + matrix_function_permute_schur(permutation, U, T); + + // compute result + MatrixType fT; // matrix function applied to T + matrix_function_compute_block_atomic(T, atomic, blockStart, clusterSize, fT); + matrix_function_compute_above_diagonal(T, blockStart, clusterSize, fT); + result = U * (fT.template triangularView() * U.adjoint()); + } +}; + +} // end of namespace internal + /** \ingroup MatrixFunctions_Module * * \brief Proxy for the matrix function of some matrix (expression). * * \tparam Derived Type of the argument to the matrix function. * - * This class holds the argument to the matrix function 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 - * matrixBase::matrixFunction() and related functions and most of the - * time this is the only way it is used. + * This class holds the argument to the matrix function 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 + * matrixBase::matrixFunction() and related functions and most of the time this is the only way it is used. */ template class MatrixFunctionReturnValue : public ReturnByValue > @@ -490,18 +486,16 @@ template class MatrixFunctionReturnValue typedef typename Derived::Index Index; typedef typename internal::stem_function::type StemFunction; - /** \brief Constructor. + /** \brief Constructor. * - * \param[in] A %Matrix (expression) forming the argument of the - * matrix function. + * \param[in] A %Matrix (expression) forming the argument of the matrix function. * \param[in] f Stem function for matrix function under consideration. */ MatrixFunctionReturnValue(const Derived& A, StemFunction f) : m_A(A), m_f(f) { } /** \brief Compute the matrix function. * - * \param[out] result \p f applied to \p A, where \p f and \p A - * are as in the constructor. + * \param[out] result \p f applied to \p A, where \p f and \p A are as in the constructor. */ template inline void evalTo(ResultType& result) const @@ -513,12 +507,12 @@ template class MatrixFunctionReturnValue static const int Options = PlainObject::Options; typedef std::complex::Real> ComplexScalar; typedef Matrix DynMatrixType; - typedef MatrixFunctionAtomic AtomicType; + + typedef internal::MatrixFunctionAtomic AtomicType; AtomicType atomic(m_f); const PlainObject Aevaluated = m_A.eval(); - MatrixFunction mf(Aevaluated, atomic); - mf.compute(result); + internal::matrix_function_compute::run(Aevaluated, atomic, result); } Index rows() const { return m_A.rows(); } diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h deleted file mode 100644 index d6ff5f1ce..000000000 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h +++ /dev/null @@ -1,127 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009 Jitse Niesen -// -// 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 -// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#ifndef EIGEN_MATRIX_FUNCTION_ATOMIC -#define EIGEN_MATRIX_FUNCTION_ATOMIC - -namespace Eigen { - -/** \ingroup MatrixFunctions_Module - * \class MatrixFunctionAtomic - * \brief Helper class for computing matrix functions of atomic matrices. - * - * \internal - * Here, an atomic matrix is a triangular matrix whose diagonal - * entries are close to each other. - */ -template -class MatrixFunctionAtomic : internal::noncopyable -{ - public: - - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::Index Index; - typedef typename NumTraits::Real RealScalar; - typedef typename internal::stem_function::type StemFunction; - typedef Matrix VectorType; - - /** \brief Constructor - * \param[in] f matrix function to compute. - */ - MatrixFunctionAtomic(StemFunction f) : m_f(f) { } - - /** \brief Compute matrix function of atomic matrix - * \param[in] A argument of matrix function, should be upper triangular and atomic - * \returns f(A), the matrix function evaluated at the given matrix - */ - MatrixType compute(const MatrixType& A); - - private: - - void computeMu(); - bool taylorConverged(Index s, const MatrixType& F, const MatrixType& Fincr, const MatrixType& P); - - /** \brief Pointer to scalar function */ - StemFunction* m_f; - - /** \brief Size of matrix function */ - Index m_Arows; - - /** \brief Mean of eigenvalues */ - Scalar m_avgEival; - - /** \brief Argument shifted by mean of eigenvalues */ - MatrixType m_Ashifted; - - /** \brief Constant used to determine whether Taylor series has converged */ - RealScalar m_mu; -}; - -template -MatrixType MatrixFunctionAtomic::compute(const MatrixType& A) -{ - // TODO: Use that A is upper triangular - m_Arows = A.rows(); - m_avgEival = A.trace() / Scalar(RealScalar(m_Arows)); - m_Ashifted = A - m_avgEival * MatrixType::Identity(m_Arows, m_Arows); - computeMu(); - MatrixType F = m_f(m_avgEival, 0) * MatrixType::Identity(m_Arows, m_Arows); - MatrixType P = m_Ashifted; - MatrixType Fincr; - for (Index s = 1; s < 1.1 * m_Arows + 10; s++) { // upper limit is fairly arbitrary - Fincr = m_f(m_avgEival, static_cast(s)) * P; - F += Fincr; - P = Scalar(RealScalar(1.0/(s + 1))) * P * m_Ashifted; - if (taylorConverged(s, F, Fincr, P)) { - return F; - } - } - eigen_assert("Taylor series does not converge" && 0); - return F; -} - -/** \brief Compute \c m_mu. */ -template -void MatrixFunctionAtomic::computeMu() -{ - const MatrixType N = MatrixType::Identity(m_Arows, m_Arows) - m_Ashifted; - VectorType e = VectorType::Ones(m_Arows); - N.template triangularView().solveInPlace(e); - m_mu = e.cwiseAbs().maxCoeff(); -} - -/** \brief Determine whether Taylor series has converged */ -template -bool MatrixFunctionAtomic::taylorConverged(Index s, const MatrixType& F, - const MatrixType& Fincr, const MatrixType& P) -{ - const Index n = F.rows(); - const RealScalar F_norm = F.cwiseAbs().rowwise().sum().maxCoeff(); - const RealScalar Fincr_norm = Fincr.cwiseAbs().rowwise().sum().maxCoeff(); - if (Fincr_norm < NumTraits::epsilon() * F_norm) { - RealScalar delta = 0; - RealScalar rfactorial = 1; - for (Index r = 0; r < n; r++) { - RealScalar mx = 0; - for (Index i = 0; i < n; i++) - mx = (std::max)(mx, std::abs(m_f(m_Ashifted(i, i) + m_avgEival, static_cast(s+r)))); - if (r != 0) - rfactorial *= RealScalar(r); - delta = (std::max)(delta, mx / rfactorial); - } - const RealScalar P_norm = P.cwiseAbs().rowwise().sum().maxCoeff(); - if (m_mu * delta * P_norm < NumTraits::epsilon() * F_norm) - return true; - } - return false; -} - -} // end namespace Eigen - -#endif // EIGEN_MATRIX_FUNCTION_ATOMIC diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h index 4b1eb5a34..6f84a31bd 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h @@ -17,81 +17,30 @@ namespace Eigen { -/** \ingroup MatrixFunctions_Module - * \class MatrixLogarithmAtomic - * \brief Helper class for computing matrix logarithm of atomic matrices. - * - * \internal - * Here, an atomic matrix is a triangular matrix whose diagonal - * entries are close to each other. - * - * \sa class MatrixFunctionAtomic, MatrixBase::log() - */ -template -class MatrixLogarithmAtomic : internal::noncopyable -{ -public: - - typedef typename MatrixType::Scalar Scalar; - // typedef typename MatrixType::Index Index; - typedef typename NumTraits::Real RealScalar; - // typedef typename internal::stem_function::type StemFunction; - // typedef Matrix VectorType; - - /** \brief Constructor. */ - MatrixLogarithmAtomic() { } - - /** \brief Compute matrix logarithm of atomic matrix - * \param[in] A argument of matrix logarithm, should be upper triangular and atomic - * \returns The logarithm of \p A. - */ - MatrixType compute(const MatrixType& A); +namespace internal { -private: - - void compute2x2(const MatrixType& A, MatrixType& result); - void computeBig(const MatrixType& A, MatrixType& result); - int getPadeDegree(float normTminusI); - int getPadeDegree(double normTminusI); - int getPadeDegree(long double normTminusI); - void computePade(MatrixType& result, const MatrixType& T, int degree); - void computePade3(MatrixType& result, const MatrixType& T); - void computePade4(MatrixType& result, const MatrixType& T); - void computePade5(MatrixType& result, const MatrixType& T); - void computePade6(MatrixType& result, const MatrixType& T); - void computePade7(MatrixType& result, const MatrixType& T); - void computePade8(MatrixType& result, const MatrixType& T); - void computePade9(MatrixType& result, const MatrixType& T); - void computePade10(MatrixType& result, const MatrixType& T); - void computePade11(MatrixType& result, const MatrixType& T); - - static const int minPadeDegree = 3; - static const int maxPadeDegree = std::numeric_limits::digits<= 24? 5: // single precision - std::numeric_limits::digits<= 53? 7: // double precision - std::numeric_limits::digits<= 64? 8: // extended precision - std::numeric_limits::digits<=106? 10: // double-double - 11; // quadruple precision +template +struct matrix_log_min_pade_degree +{ + static const int value = 3; }; -/** \brief Compute logarithm of triangular matrix with clustered eigenvalues. */ -template -MatrixType MatrixLogarithmAtomic::compute(const MatrixType& A) +template +struct matrix_log_max_pade_degree { - using std::log; - MatrixType result(A.rows(), A.rows()); - if (A.rows() == 1) - result(0,0) = log(A(0,0)); - else if (A.rows() == 2) - compute2x2(A, result); - else - computeBig(A, result); - return result; -} + typedef typename NumTraits::Real RealScalar; + static const int value = std::numeric_limits::digits<= 24? 5: // single precision + std::numeric_limits::digits<= 53? 7: // double precision + std::numeric_limits::digits<= 64? 8: // extended precision + std::numeric_limits::digits<=106? 10: // double-double + 11; // quadruple precision +}; /** \brief Compute logarithm of 2x2 triangular matrix. */ template -void MatrixLogarithmAtomic::compute2x2(const MatrixType& A, MatrixType& result) +void matrix_log_compute_2x2(const MatrixType& A, MatrixType& result) { + typedef typename MatrixType::Scalar Scalar; using std::abs; using std::ceil; using std::imag; @@ -116,47 +65,14 @@ void MatrixLogarithmAtomic::compute2x2(const MatrixType& A, MatrixTy } } -/** \brief Compute logarithm of triangular matrices with size > 2. - * \details This uses a inverse scale-and-square algorithm. */ -template -void MatrixLogarithmAtomic::computeBig(const MatrixType& A, MatrixType& result) -{ - using std::pow; - int numberOfSquareRoots = 0; - int numberOfExtraSquareRoots = 0; - int degree; - MatrixType T = A, sqrtT; - const RealScalar maxNormForPade = maxPadeDegree<= 5? 5.3149729967117310e-1: // single precision - maxPadeDegree<= 7? 2.6429608311114350e-1: // double precision - maxPadeDegree<= 8? 2.32777776523703892094e-1L: // extended precision - maxPadeDegree<=10? 1.05026503471351080481093652651105e-1L: // double-double - 1.1880960220216759245467951592883642e-1L; // quadruple precision - - while (true) { - RealScalar normTminusI = (T - MatrixType::Identity(T.rows(), T.rows())).cwiseAbs().colwise().sum().maxCoeff(); - if (normTminusI < maxNormForPade) { - degree = getPadeDegree(normTminusI); - int degree2 = getPadeDegree(normTminusI / RealScalar(2)); - if ((degree - degree2 <= 1) || (numberOfExtraSquareRoots == 1)) - break; - ++numberOfExtraSquareRoots; - } - matrix_sqrt_triangular(T, sqrtT); - T = sqrtT.template triangularView(); - ++numberOfSquareRoots; - } - - computePade(result, T, degree); - result *= pow(RealScalar(2), numberOfSquareRoots); -} - /* \brief Get suitable degree for Pade approximation. (specialized for RealScalar = float) */ -template -int MatrixLogarithmAtomic::getPadeDegree(float normTminusI) +int matrix_log_get_pade_degree(float normTminusI) { const float maxNormForPade[] = { 2.5111573934555054e-1 /* degree = 3 */ , 4.0535837411880493e-1, 5.3149729967117310e-1 }; - int degree = 3; + const int minPadeDegree = matrix_log_min_pade_degree::value; + const int maxPadeDegree = matrix_log_max_pade_degree::value; + int degree = minPadeDegree; for (; degree <= maxPadeDegree; ++degree) if (normTminusI <= maxNormForPade[degree - minPadeDegree]) break; @@ -164,12 +80,13 @@ int MatrixLogarithmAtomic::getPadeDegree(float normTminusI) } /* \brief Get suitable degree for Pade approximation. (specialized for RealScalar = double) */ -template -int MatrixLogarithmAtomic::getPadeDegree(double normTminusI) +int matrix_log_get_pade_degree(double normTminusI) { const double maxNormForPade[] = { 1.6206284795015624e-2 /* degree = 3 */ , 5.3873532631381171e-2, 1.1352802267628681e-1, 1.8662860613541288e-1, 2.642960831111435e-1 }; - int degree = 3; + const int minPadeDegree = matrix_log_min_pade_degree::value; + const int maxPadeDegree = matrix_log_max_pade_degree::value; + int degree = minPadeDegree; for (; degree <= maxPadeDegree; ++degree) if (normTminusI <= maxNormForPade[degree - minPadeDegree]) break; @@ -177,8 +94,7 @@ int MatrixLogarithmAtomic::getPadeDegree(double normTminusI) } /* \brief Get suitable degree for Pade approximation. (specialized for RealScalar = long double) */ -template -int MatrixLogarithmAtomic::getPadeDegree(long double normTminusI) +int matrix_log_get_pade_degree(long double normTminusI) { #if LDBL_MANT_DIG == 53 // double precision const long double maxNormForPade[] = { 1.6206284795015624e-2L /* degree = 3 */ , 5.3873532631381171e-2L, @@ -200,7 +116,9 @@ int MatrixLogarithmAtomic::getPadeDegree(long double normTminusI) 3.6688019729653446926585242192447447e-2L, 5.9290962294020186998954055264528393e-2L, 8.6998436081634343903250580992127677e-2L, 1.1880960220216759245467951592883642e-1L }; #endif - int degree = 3; + const int minPadeDegree = matrix_log_min_pade_degree::value; + const int maxPadeDegree = matrix_log_max_pade_degree::value; + int degree = minPadeDegree; for (; degree <= maxPadeDegree; ++degree) if (normTminusI <= maxNormForPade[degree - minPadeDegree]) break; @@ -209,197 +127,168 @@ int MatrixLogarithmAtomic::getPadeDegree(long double normTminusI) /* \brief Compute Pade approximation to matrix logarithm */ template -void MatrixLogarithmAtomic::computePade(MatrixType& result, const MatrixType& T, int degree) +void matrix_log_compute_pade(MatrixType& result, const MatrixType& T, int degree) { - switch (degree) { - case 3: computePade3(result, T); break; - case 4: computePade4(result, T); break; - case 5: computePade5(result, T); break; - case 6: computePade6(result, T); break; - case 7: computePade7(result, T); break; - case 8: computePade8(result, T); break; - case 9: computePade9(result, T); break; - case 10: computePade10(result, T); break; - case 11: computePade11(result, T); break; - default: assert(false); // should never happen - } -} + typedef typename NumTraits::Real RealScalar; + const int minPadeDegree = 3; + const int maxPadeDegree = 11; + assert(degree >= minPadeDegree && degree <= maxPadeDegree); + + const RealScalar nodes[][maxPadeDegree] = { + { 0.1127016653792583114820734600217600L, 0.5000000000000000000000000000000000L, // degree 3 + 0.8872983346207416885179265399782400L }, + { 0.0694318442029737123880267555535953L, 0.3300094782075718675986671204483777L, // degree 4 + 0.6699905217924281324013328795516223L, 0.9305681557970262876119732444464048L }, + { 0.0469100770306680036011865608503035L, 0.2307653449471584544818427896498956L, // degree 5 + 0.5000000000000000000000000000000000L, 0.7692346550528415455181572103501044L, + 0.9530899229693319963988134391496965L }, + { 0.0337652428984239860938492227530027L, 0.1693953067668677431693002024900473L, // degree 6 + 0.3806904069584015456847491391596440L, 0.6193095930415984543152508608403560L, + 0.8306046932331322568306997975099527L, 0.9662347571015760139061507772469973L }, + { 0.0254460438286207377369051579760744L, 0.1292344072003027800680676133596058L, // degree 7 + 0.2970774243113014165466967939615193L, 0.5000000000000000000000000000000000L, + 0.7029225756886985834533032060384807L, 0.8707655927996972199319323866403942L, + 0.9745539561713792622630948420239256L }, + { 0.0198550717512318841582195657152635L, 0.1016667612931866302042230317620848L, // degree 8 + 0.2372337950418355070911304754053768L, 0.4082826787521750975302619288199080L, + 0.5917173212478249024697380711800920L, 0.7627662049581644929088695245946232L, + 0.8983332387068133697957769682379152L, 0.9801449282487681158417804342847365L }, + { 0.0159198802461869550822118985481636L, 0.0819844463366821028502851059651326L, // degree 9 + 0.1933142836497048013456489803292629L, 0.3378732882980955354807309926783317L, + 0.5000000000000000000000000000000000L, 0.6621267117019044645192690073216683L, + 0.8066857163502951986543510196707371L, 0.9180155536633178971497148940348674L, + 0.9840801197538130449177881014518364L }, + { 0.0130467357414141399610179939577740L, 0.0674683166555077446339516557882535L, // degree 10 + 0.1602952158504877968828363174425632L, 0.2833023029353764046003670284171079L, + 0.4255628305091843945575869994351400L, 0.5744371694908156054424130005648600L, + 0.7166976970646235953996329715828921L, 0.8397047841495122031171636825574368L, + 0.9325316833444922553660483442117465L, 0.9869532642585858600389820060422260L }, + { 0.0108856709269715035980309994385713L, 0.0564687001159523504624211153480364L, // degree 11 + 0.1349239972129753379532918739844233L, 0.2404519353965940920371371652706952L, + 0.3652284220238275138342340072995692L, 0.5000000000000000000000000000000000L, + 0.6347715779761724861657659927004308L, 0.7595480646034059079628628347293048L, + 0.8650760027870246620467081260155767L, 0.9435312998840476495375788846519636L, + 0.9891143290730284964019690005614287L } }; + + const RealScalar weights[][maxPadeDegree] = { + { 0.2777777777777777777777777777777778L, 0.4444444444444444444444444444444444L, // degree 3 + 0.2777777777777777777777777777777778L }, + { 0.1739274225687269286865319746109997L, 0.3260725774312730713134680253890003L, // degree 4 + 0.3260725774312730713134680253890003L, 0.1739274225687269286865319746109997L }, + { 0.1184634425280945437571320203599587L, 0.2393143352496832340206457574178191L, // degree 5 + 0.2844444444444444444444444444444444L, 0.2393143352496832340206457574178191L, + 0.1184634425280945437571320203599587L }, + { 0.0856622461895851725201480710863665L, 0.1803807865240693037849167569188581L, // degree 6 + 0.2339569672863455236949351719947755L, 0.2339569672863455236949351719947755L, + 0.1803807865240693037849167569188581L, 0.0856622461895851725201480710863665L }, + { 0.0647424830844348466353057163395410L, 0.1398526957446383339507338857118898L, // degree 7 + 0.1909150252525594724751848877444876L, 0.2089795918367346938775510204081633L, + 0.1909150252525594724751848877444876L, 0.1398526957446383339507338857118898L, + 0.0647424830844348466353057163395410L }, + { 0.0506142681451881295762656771549811L, 0.1111905172266872352721779972131204L, // degree 8 + 0.1568533229389436436689811009933007L, 0.1813418916891809914825752246385978L, + 0.1813418916891809914825752246385978L, 0.1568533229389436436689811009933007L, + 0.1111905172266872352721779972131204L, 0.0506142681451881295762656771549811L }, + { 0.0406371941807872059859460790552618L, 0.0903240803474287020292360156214564L, // degree 9 + 0.1303053482014677311593714347093164L, 0.1561735385200014200343152032922218L, + 0.1651196775006298815822625346434870L, 0.1561735385200014200343152032922218L, + 0.1303053482014677311593714347093164L, 0.0903240803474287020292360156214564L, + 0.0406371941807872059859460790552618L }, + { 0.0333356721543440687967844049466659L, 0.0747256745752902965728881698288487L, // degree 10 + 0.1095431812579910219977674671140816L, 0.1346333596549981775456134607847347L, + 0.1477621123573764350869464973256692L, 0.1477621123573764350869464973256692L, + 0.1346333596549981775456134607847347L, 0.1095431812579910219977674671140816L, + 0.0747256745752902965728881698288487L, 0.0333356721543440687967844049466659L }, + { 0.0278342835580868332413768602212743L, 0.0627901847324523123173471496119701L, // degree 11 + 0.0931451054638671257130488207158280L, 0.1165968822959952399592618524215876L, + 0.1314022722551233310903444349452546L, 0.1364625433889503153572417641681711L, + 0.1314022722551233310903444349452546L, 0.1165968822959952399592618524215876L, + 0.0931451054638671257130488207158280L, 0.0627901847324523123173471496119701L, + 0.0278342835580868332413768602212743L } }; -template -void MatrixLogarithmAtomic::computePade3(MatrixType& result, const MatrixType& T) -{ - const int degree = 3; - const RealScalar nodes[] = { 0.1127016653792583114820734600217600L, 0.5000000000000000000000000000000000L, - 0.8872983346207416885179265399782400L }; - const RealScalar weights[] = { 0.2777777777777777777777777777777778L, 0.4444444444444444444444444444444444L, - 0.2777777777777777777777777777777778L }; - eigen_assert(degree <= maxPadeDegree); MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows()); result.setZero(T.rows(), T.rows()); - for (int k = 0; k < degree; ++k) - result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI) - .template triangularView().solve(TminusI); -} + for (int k = 0; k < degree; ++k) { + RealScalar weight = weights[degree-minPadeDegree][k]; + RealScalar node = nodes[degree-minPadeDegree][k]; + result += weight * (MatrixType::Identity(T.rows(), T.rows()) + node * TminusI) + .template triangularView().solve(TminusI); + } +} +/** \brief Compute logarithm of triangular matrices with size > 2. + * \details This uses a inverse scale-and-square algorithm. */ template -void MatrixLogarithmAtomic::computePade4(MatrixType& result, const MatrixType& T) +void matrix_log_compute_big(const MatrixType& A, MatrixType& result) { - const int degree = 4; - const RealScalar nodes[] = { 0.0694318442029737123880267555535953L, 0.3300094782075718675986671204483777L, - 0.6699905217924281324013328795516223L, 0.9305681557970262876119732444464048L }; - const RealScalar weights[] = { 0.1739274225687269286865319746109997L, 0.3260725774312730713134680253890003L, - 0.3260725774312730713134680253890003L, 0.1739274225687269286865319746109997L }; - eigen_assert(degree <= maxPadeDegree); - MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows()); - result.setZero(T.rows(), T.rows()); - for (int k = 0; k < degree; ++k) - result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI) - .template triangularView().solve(TminusI); -} + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + using std::pow; -template -void MatrixLogarithmAtomic::computePade5(MatrixType& result, const MatrixType& T) -{ - const int degree = 5; - const RealScalar nodes[] = { 0.0469100770306680036011865608503035L, 0.2307653449471584544818427896498956L, - 0.5000000000000000000000000000000000L, 0.7692346550528415455181572103501044L, - 0.9530899229693319963988134391496965L }; - const RealScalar weights[] = { 0.1184634425280945437571320203599587L, 0.2393143352496832340206457574178191L, - 0.2844444444444444444444444444444444L, 0.2393143352496832340206457574178191L, - 0.1184634425280945437571320203599587L }; - eigen_assert(degree <= maxPadeDegree); - MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows()); - result.setZero(T.rows(), T.rows()); - for (int k = 0; k < degree; ++k) - result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI) - .template triangularView().solve(TminusI); -} + int numberOfSquareRoots = 0; + int numberOfExtraSquareRoots = 0; + int degree; + MatrixType T = A, sqrtT; -template -void MatrixLogarithmAtomic::computePade6(MatrixType& result, const MatrixType& T) -{ - const int degree = 6; - const RealScalar nodes[] = { 0.0337652428984239860938492227530027L, 0.1693953067668677431693002024900473L, - 0.3806904069584015456847491391596440L, 0.6193095930415984543152508608403560L, - 0.8306046932331322568306997975099527L, 0.9662347571015760139061507772469973L }; - const RealScalar weights[] = { 0.0856622461895851725201480710863665L, 0.1803807865240693037849167569188581L, - 0.2339569672863455236949351719947755L, 0.2339569672863455236949351719947755L, - 0.1803807865240693037849167569188581L, 0.0856622461895851725201480710863665L }; - eigen_assert(degree <= maxPadeDegree); - MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows()); - result.setZero(T.rows(), T.rows()); - for (int k = 0; k < degree; ++k) - result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI) - .template triangularView().solve(TminusI); -} + int maxPadeDegree = matrix_log_max_pade_degree::value; + const RealScalar maxNormForPade = maxPadeDegree<= 5? 5.3149729967117310e-1: // single precision + maxPadeDegree<= 7? 2.6429608311114350e-1: // double precision + maxPadeDegree<= 8? 2.32777776523703892094e-1L: // extended precision + maxPadeDegree<=10? 1.05026503471351080481093652651105e-1L: // double-double + 1.1880960220216759245467951592883642e-1L; // quadruple precision -template -void MatrixLogarithmAtomic::computePade7(MatrixType& result, const MatrixType& T) -{ - const int degree = 7; - const RealScalar nodes[] = { 0.0254460438286207377369051579760744L, 0.1292344072003027800680676133596058L, - 0.2970774243113014165466967939615193L, 0.5000000000000000000000000000000000L, - 0.7029225756886985834533032060384807L, 0.8707655927996972199319323866403942L, - 0.9745539561713792622630948420239256L }; - const RealScalar weights[] = { 0.0647424830844348466353057163395410L, 0.1398526957446383339507338857118898L, - 0.1909150252525594724751848877444876L, 0.2089795918367346938775510204081633L, - 0.1909150252525594724751848877444876L, 0.1398526957446383339507338857118898L, - 0.0647424830844348466353057163395410L }; - eigen_assert(degree <= maxPadeDegree); - MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows()); - result.setZero(T.rows(), T.rows()); - for (int k = 0; k < degree; ++k) - result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI) - .template triangularView().solve(TminusI); -} + while (true) { + RealScalar normTminusI = (T - MatrixType::Identity(T.rows(), T.rows())).cwiseAbs().colwise().sum().maxCoeff(); + if (normTminusI < maxNormForPade) { + degree = matrix_log_get_pade_degree(normTminusI); + int degree2 = matrix_log_get_pade_degree(normTminusI / RealScalar(2)); + if ((degree - degree2 <= 1) || (numberOfExtraSquareRoots == 1)) + break; + ++numberOfExtraSquareRoots; + } + matrix_sqrt_triangular(T, sqrtT); + T = sqrtT.template triangularView(); + ++numberOfSquareRoots; + } -template -void MatrixLogarithmAtomic::computePade8(MatrixType& result, const MatrixType& T) -{ - const int degree = 8; - const RealScalar nodes[] = { 0.0198550717512318841582195657152635L, 0.1016667612931866302042230317620848L, - 0.2372337950418355070911304754053768L, 0.4082826787521750975302619288199080L, - 0.5917173212478249024697380711800920L, 0.7627662049581644929088695245946232L, - 0.8983332387068133697957769682379152L, 0.9801449282487681158417804342847365L }; - const RealScalar weights[] = { 0.0506142681451881295762656771549811L, 0.1111905172266872352721779972131204L, - 0.1568533229389436436689811009933007L, 0.1813418916891809914825752246385978L, - 0.1813418916891809914825752246385978L, 0.1568533229389436436689811009933007L, - 0.1111905172266872352721779972131204L, 0.0506142681451881295762656771549811L }; - eigen_assert(degree <= maxPadeDegree); - MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows()); - result.setZero(T.rows(), T.rows()); - for (int k = 0; k < degree; ++k) - result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI) - .template triangularView().solve(TminusI); + matrix_log_compute_pade(result, T, degree); + result *= pow(RealScalar(2), numberOfSquareRoots); } +/** \ingroup MatrixFunctions_Module + * \class MatrixLogarithmAtomic + * \brief Helper class for computing matrix logarithm of atomic matrices. + * + * Here, an atomic matrix is a triangular matrix whose diagonal entries are close to each other. + * + * \sa class MatrixFunctionAtomic, MatrixBase::log() + */ template -void MatrixLogarithmAtomic::computePade9(MatrixType& result, const MatrixType& T) +class MatrixLogarithmAtomic { - const int degree = 9; - const RealScalar nodes[] = { 0.0159198802461869550822118985481636L, 0.0819844463366821028502851059651326L, - 0.1933142836497048013456489803292629L, 0.3378732882980955354807309926783317L, - 0.5000000000000000000000000000000000L, 0.6621267117019044645192690073216683L, - 0.8066857163502951986543510196707371L, 0.9180155536633178971497148940348674L, - 0.9840801197538130449177881014518364L }; - const RealScalar weights[] = { 0.0406371941807872059859460790552618L, 0.0903240803474287020292360156214564L, - 0.1303053482014677311593714347093164L, 0.1561735385200014200343152032922218L, - 0.1651196775006298815822625346434870L, 0.1561735385200014200343152032922218L, - 0.1303053482014677311593714347093164L, 0.0903240803474287020292360156214564L, - 0.0406371941807872059859460790552618L }; - eigen_assert(degree <= maxPadeDegree); - MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows()); - result.setZero(T.rows(), T.rows()); - for (int k = 0; k < degree; ++k) - result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI) - .template triangularView().solve(TminusI); -} +public: + /** \brief Compute matrix logarithm of atomic matrix + * \param[in] A argument of matrix logarithm, should be upper triangular and atomic + * \returns The logarithm of \p A. + */ + MatrixType compute(const MatrixType& A); +}; template -void MatrixLogarithmAtomic::computePade10(MatrixType& result, const MatrixType& T) +MatrixType MatrixLogarithmAtomic::compute(const MatrixType& A) { - const int degree = 10; - const RealScalar nodes[] = { 0.0130467357414141399610179939577740L, 0.0674683166555077446339516557882535L, - 0.1602952158504877968828363174425632L, 0.2833023029353764046003670284171079L, - 0.4255628305091843945575869994351400L, 0.5744371694908156054424130005648600L, - 0.7166976970646235953996329715828921L, 0.8397047841495122031171636825574368L, - 0.9325316833444922553660483442117465L, 0.9869532642585858600389820060422260L }; - const RealScalar weights[] = { 0.0333356721543440687967844049466659L, 0.0747256745752902965728881698288487L, - 0.1095431812579910219977674671140816L, 0.1346333596549981775456134607847347L, - 0.1477621123573764350869464973256692L, 0.1477621123573764350869464973256692L, - 0.1346333596549981775456134607847347L, 0.1095431812579910219977674671140816L, - 0.0747256745752902965728881698288487L, 0.0333356721543440687967844049466659L }; - eigen_assert(degree <= maxPadeDegree); - MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows()); - result.setZero(T.rows(), T.rows()); - for (int k = 0; k < degree; ++k) - result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI) - .template triangularView().solve(TminusI); + using std::log; + MatrixType result(A.rows(), A.rows()); + if (A.rows() == 1) + result(0,0) = log(A(0,0)); + else if (A.rows() == 2) + matrix_log_compute_2x2(A, result); + else + matrix_log_compute_big(A, result); + return result; } -template -void MatrixLogarithmAtomic::computePade11(MatrixType& result, const MatrixType& T) -{ - const int degree = 11; - const RealScalar nodes[] = { 0.0108856709269715035980309994385713L, 0.0564687001159523504624211153480364L, - 0.1349239972129753379532918739844233L, 0.2404519353965940920371371652706952L, - 0.3652284220238275138342340072995692L, 0.5000000000000000000000000000000000L, - 0.6347715779761724861657659927004308L, 0.7595480646034059079628628347293048L, - 0.8650760027870246620467081260155767L, 0.9435312998840476495375788846519636L, - 0.9891143290730284964019690005614287L }; - const RealScalar weights[] = { 0.0278342835580868332413768602212743L, 0.0627901847324523123173471496119701L, - 0.0931451054638671257130488207158280L, 0.1165968822959952399592618524215876L, - 0.1314022722551233310903444349452546L, 0.1364625433889503153572417641681711L, - 0.1314022722551233310903444349452546L, 0.1165968822959952399592618524215876L, - 0.0931451054638671257130488207158280L, 0.0627901847324523123173471496119701L, - 0.0278342835580868332413768602212743L }; - eigen_assert(degree <= maxPadeDegree); - MatrixType TminusI = T - MatrixType::Identity(T.rows(), T.rows()); - result.setZero(T.rows(), T.rows()); - for (int k = 0; k < degree; ++k) - result += weights[k] * (MatrixType::Identity(T.rows(), T.rows()) + nodes[k] * TminusI) - .template triangularView().solve(TminusI); -} +} // end of namespace internal /** \ingroup MatrixFunctions_Module * @@ -441,12 +330,11 @@ public: static const int Options = PlainObject::Options; typedef std::complex::Real> ComplexScalar; typedef Matrix DynMatrixType; - typedef MatrixLogarithmAtomic AtomicType; + typedef internal::MatrixLogarithmAtomic AtomicType; AtomicType atomic; const PlainObject Aevaluated = m_A.eval(); - MatrixFunction mf(Aevaluated, atomic); - mf.compute(result); + internal::matrix_function_compute::run(Aevaluated, atomic, result); } Index rows() const { return m_A.rows(); } diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h index 0261d4aa9..314b3f38e 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2011 Jitse Niesen +// Copyright (C) 2011, 2013 Jitse Niesen // // 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 -- cgit v1.2.3