diff options
| -rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h | 75 | ||||
| -rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h | 142 |
2 files changed, 149 insertions, 68 deletions
diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h index 963244771..334b94336 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h @@ -78,6 +78,7 @@ EIGEN_STRONG_INLINE void ei_matrix_function(const MatrixBase<Derived>& M, typename ei_stem_function<typename ei_traits<Derived>::Scalar>::type f, typename MatrixBase<Derived>::PlainMatrixType* result); +#include "MatrixFunctionAtomic.h" /** \ingroup MatrixFunctions_Module * \class MatrixFunction @@ -169,15 +170,10 @@ class MatrixFunction<MatrixType, 1, 1> void computeTriangular(const MatrixType& T, MatrixType& result, const IntVectorType& blockSize); void computeBlockAtomic(const MatrixType& T, MatrixType& result, const IntVectorType& blockSize); MatrixType solveTriangularSylvester(const MatrixType& A, const MatrixType& B, const MatrixType& C); - MatrixType computeAtomic(const MatrixType& T); void divideInBlocks(const VectorType& v, listOfLists* result); void constructPermutation(const VectorType& diag, const listOfLists& blocks, IntVectorType& blockSize, IntVectorType& permutation); - RealScalar computeMu(const MatrixType& M); - bool taylorConverged(const MatrixType& T, int s, const MatrixType& F, - const MatrixType& Fincr, const MatrixType& P, RealScalar mu); - static const RealScalar separation() { return static_cast<RealScalar>(0.01); } StemFunction *m_f; }; @@ -271,11 +267,11 @@ void MatrixFunction<MatrixType,1,1>::computeTriangular(const MatrixType& T, Matr /** \brief Solve a triangular Sylvester equation AX + XB = C * - * \param[in] A The matrix A; should be square and upper triangular - * \param[in] B The matrix B; should be square and upper triangular - * \param[in] C The matrix C; should have correct size. + * \param[in] A the matrix A; should be square and upper triangular + * \param[in] B the matrix B; should be square and upper triangular + * \param[in] C the matrix C; should have correct size. * - * \returns The solution X. + * \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 @@ -346,8 +342,9 @@ void MatrixFunction<MatrixType,1,1>::computeBlockAtomic(const MatrixType& T, Mat result.resize(T.rows(), T.cols()); result.setZero(); for (int i = 0; i < blockSize.rows(); i++) { + MatrixFunctionAtomic<MatrixType> mfa(m_f); result.block(blockStart, blockStart, blockSize(i), blockSize(i)) - = computeAtomic(T.block(blockStart, blockStart, blockSize(i), blockSize(i))); + = mfa.compute(T.block(blockStart, blockStart, blockSize(i), blockSize(i))); blockStart += blockSize(i); } } @@ -434,64 +431,6 @@ void MatrixFunction<MatrixType,1,1>::constructPermutation(const VectorType& diag } } -template <typename MatrixType> -MatrixType MatrixFunction<MatrixType,1,1>::computeAtomic(const MatrixType& T) -{ - // TODO: Use that T is upper triangular - const int n = T.rows(); - const Scalar sigma = T.trace() / Scalar(n); - const MatrixType M = T - sigma * MatrixType::Identity(n, n); - const RealScalar mu = computeMu(M); - MatrixType F = m_f(sigma, 0) * MatrixType::Identity(n, n); - MatrixType P = M; - MatrixType Fincr; - for (int s = 1; s < 1.1*n + 10; s++) { // upper limit is fairly arbitrary - Fincr = m_f(sigma, s) * P; - F += Fincr; - P = (1/(s + 1.0)) * P * M; - if (taylorConverged(T, s, F, Fincr, P, mu)) { - return F; - } - } - ei_assert("Taylor series does not converge" && 0); - return F; -} - -template <typename MatrixType> -typename MatrixFunction<MatrixType,1,1>::RealScalar MatrixFunction<MatrixType,1,1>::computeMu(const MatrixType& M) -{ - const int n = M.rows(); - const MatrixType N = MatrixType::Identity(n, n) - M; - VectorType e = VectorType::Ones(n); - N.template triangularView<UpperTriangular>().solveInPlace(e); - return e.cwise().abs().maxCoeff(); -} - -template <typename MatrixType> -bool MatrixFunction<MatrixType,1,1>::taylorConverged(const MatrixType& T, int s, const MatrixType& F, - const MatrixType& Fincr, const MatrixType& P, RealScalar mu) -{ - const int n = F.rows(); - const RealScalar F_norm = F.cwise().abs().rowwise().sum().maxCoeff(); - const RealScalar Fincr_norm = Fincr.cwise().abs().rowwise().sum().maxCoeff(); - if (Fincr_norm < epsilon<Scalar>() * F_norm) { - RealScalar delta = 0; - RealScalar rfactorial = 1; - for (int r = 0; r < n; r++) { - RealScalar mx = 0; - for (int i = 0; i < n; i++) - mx = std::max(mx, std::abs(m_f(T(i, i), s+r))); - if (r != 0) - rfactorial *= r; - delta = std::max(delta, mx / rfactorial); - } - const RealScalar P_norm = P.cwise().abs().rowwise().sum().maxCoeff(); - if (mu * delta * P_norm < epsilon<Scalar>() * F_norm) - return true; - } - return false; -} - template <typename Derived> EIGEN_STRONG_INLINE void ei_matrix_function(const MatrixBase<Derived>& M, typename ei_stem_function<typename ei_traits<Derived>::Scalar>::type f, diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h new file mode 100644 index 000000000..117ee82d7 --- /dev/null +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixFunctionAtomic.h @@ -0,0 +1,142 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Jitse Niesen <jitse@maths.leeds.ac.uk> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_MATRIX_FUNCTION_ATOMIC +#define EIGEN_MATRIX_FUNCTION_ATOMIC + +/** \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 <typename MatrixType> +class MatrixFunctionAtomic +{ + public: + + typedef ei_traits<MatrixType> Traits; + typedef typename Traits::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef typename ei_stem_function<Scalar>::type StemFunction; + typedef Matrix<Scalar, Traits::RowsAtCompileTime, 1> 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: + + // Prevent copying + MatrixFunctionAtomic(const MatrixFunctionAtomic&); + MatrixFunctionAtomic& operator=(const MatrixFunctionAtomic&); + + void computeMu(); + bool taylorConverged(int s, const MatrixType& F, const MatrixType& Fincr, const MatrixType& P); + + /** \brief Pointer to scalar function */ + StemFunction* m_f; + + /** \brief Size of matrix function */ + int 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 <typename MatrixType> +MatrixType MatrixFunctionAtomic<MatrixType>::compute(const MatrixType& A) +{ + // TODO: Use that A is upper triangular + m_Arows = A.rows(); + m_avgEival = A.trace() / Scalar(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 (int s = 1; s < 1.1 * m_Arows + 10; s++) { // upper limit is fairly arbitrary + Fincr = m_f(m_avgEival, s) * P; + F += Fincr; + P = (1/(s + 1.0)) * P * m_Ashifted; + if (taylorConverged(s, F, Fincr, P)) { + return F; + } + } + ei_assert("Taylor series does not converge" && 0); + return F; +} + +/** \brief Compute \c m_mu. */ +template <typename MatrixType> +void MatrixFunctionAtomic<MatrixType>::computeMu() +{ + const MatrixType N = MatrixType::Identity(m_Arows, m_Arows) - m_Ashifted; + VectorType e = VectorType::Ones(m_Arows); + N.template triangularView<UpperTriangular>().solveInPlace(e); + m_mu = e.cwise().abs().maxCoeff(); +} + +/** \brief Determine whether Taylor series has converged */ +template <typename MatrixType> +bool MatrixFunctionAtomic<MatrixType>::taylorConverged(int s, const MatrixType& F, + const MatrixType& Fincr, const MatrixType& P) +{ + const int n = F.rows(); + const RealScalar F_norm = F.cwise().abs().rowwise().sum().maxCoeff(); + const RealScalar Fincr_norm = Fincr.cwise().abs().rowwise().sum().maxCoeff(); + if (Fincr_norm < epsilon<Scalar>() * F_norm) { + RealScalar delta = 0; + RealScalar rfactorial = 1; + for (int r = 0; r < n; r++) { + RealScalar mx = 0; + for (int i = 0; i < n; i++) + mx = std::max(mx, std::abs(m_f(m_Ashifted(i, i) + m_avgEival, s+r))); + if (r != 0) + rfactorial *= r; + delta = std::max(delta, mx / rfactorial); + } + const RealScalar P_norm = P.cwise().abs().rowwise().sum().maxCoeff(); + if (m_mu * delta * P_norm < epsilon<Scalar>() * F_norm) + return true; + } + return false; +} + +#endif // EIGEN_MATRIX_FUNCTION_ATOMIC |