GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-05-16 17:56:44 +08:00
Code Issues Packages Projects Releases Wiki Activity

Decouple MatrixFunction and MatrixFunctionAtomic

Browse Source 
in preparation for implementation of matrix log.
This commit is contained in:
Jitse Niesen 2011-06-07 14:40:27 +01:00
parent 86ca35ccff
commit a6d42e28fe
1 changed files with 80 additions and 68 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

148
unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h
View File

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