diff options
author | Jitse Niesen <jitse@maths.leeds.ac.uk> | 2011-05-07 22:57:46 +0100 |
---|---|---|
committer | Jitse Niesen <jitse@maths.leeds.ac.uk> | 2011-05-07 22:57:46 +0100 |
commit | 6b4e215710dd5c12ad1fe8e820875674bdd849c8 (patch) | |
tree | f0aabd83f0302f1fda83da172cad4ce02f7020ea /unsupported/Eigen | |
parent | 0896c6d97dc587f76b92687d87ae3e21d0c5aab4 (diff) |
Implement matrix square root for complex matrices.
I hope to implement the real case soon, but it's a bit more
complicated due to the 2-by-2 blocks in the real Schur decomposition.
Diffstat (limited to 'unsupported/Eigen')
-rw-r--r-- | unsupported/Eigen/MatrixFunctions | 1 | ||||
-rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h | 8 | ||||
-rw-r--r-- | unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h | 135 |
3 files changed, 140 insertions, 4 deletions
diff --git a/unsupported/Eigen/MatrixFunctions b/unsupported/Eigen/MatrixFunctions index d39c49e53..79340e206 100644 --- a/unsupported/Eigen/MatrixFunctions +++ b/unsupported/Eigen/MatrixFunctions @@ -69,6 +69,7 @@ namespace Eigen { #include "src/MatrixFunctions/MatrixExponential.h" #include "src/MatrixFunctions/MatrixFunction.h" +#include "src/MatrixFunctions/MatrixSquareRoot.h" diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h index 4b9d8a102..b343f38bc 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixFunction.h @@ -68,8 +68,8 @@ class MatrixFunction }; -/** \ingroup MatrixFunctions_Module - * \brief Partial specialization of MatrixFunction for real matrices \internal +/** \internal \ingroup MatrixFunctions_Module + * \brief Partial specialization of MatrixFunction for real matrices */ template <typename MatrixType> class MatrixFunction<MatrixType, 0> @@ -124,8 +124,8 @@ class MatrixFunction<MatrixType, 0> }; -/** \ingroup MatrixFunctions_Module - * \brief Partial specialization of MatrixFunction for complex matrices \internal +/** \internal \ingroup MatrixFunctions_Module + * \brief Partial specialization of MatrixFunction for complex matrices */ template <typename MatrixType> class MatrixFunction<MatrixType, 1> diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h new file mode 100644 index 000000000..5eeda11ec --- /dev/null +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixSquareRoot.h @@ -0,0 +1,135 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2011 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_SQUARE_ROOT +#define EIGEN_MATRIX_SQUARE_ROOT + +/** \ingroup MatrixFunctions_Module + * \brief Class for computing matrix square roots. + * \tparam MatrixType type of the argument of the matrix square root, + * expected to be an instantiation of the Matrix class template. + */ +template <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> +class MatrixSquareRoot +{ + public: + + /** \brief Constructor. + * + * \param[in] A matrix whose square root is to be computed. + * + * The class stores a reference to \p A, so it should not be + * changed (or destroyed) before compute() is called. + */ + MatrixSquareRoot(const MatrixType& A); + + /** \brief Compute the matrix square root + * + * \param[out] result square root of \p A, as specified in the constructor. + * + * See MatrixBase::sqrt() for details on how this computation + * is implemented. + */ + template <typename ResultType> + void compute(ResultType &result); +}; + + +// ********** Partial specialization for real matrices ********** + +template <typename MatrixType> +class MatrixSquareRoot<MatrixType, 0> +{ + public: + MatrixSquareRoot(const MatrixType& A) + : m_A(A) + { + eigen_assert(A.rows() == A.cols()); + } + + template <typename ResultType> void compute(ResultType &result); + + private: + const MatrixType& m_A; +}; + +template <typename MatrixType> +template <typename ResultType> +void MatrixSquareRoot<MatrixType, 0>::compute(ResultType &result) +{ + eigen_assert("Square root of real matrices is not implemented!"); +} + + +// ********** Partial specialization for complex matrices ********** + +template <typename MatrixType> +class MatrixSquareRoot<MatrixType, 1> +{ + public: + MatrixSquareRoot(const MatrixType& A) + : m_A(A) + { + eigen_assert(A.rows() == A.cols()); + } + + template <typename ResultType> void compute(ResultType &result); + + private: + const MatrixType& m_A; +}; + +template <typename MatrixType> +template <typename ResultType> +void MatrixSquareRoot<MatrixType, 1>::compute(ResultType &result) +{ + // Compute Schur decomposition of m_A + const ComplexSchur<MatrixType> schurOfA(m_A); + const MatrixType& T = schurOfA.matrixT(); + const MatrixType& U = schurOfA.matrixU(); + + // Compute square root of T and store it in upper triangular part of result + // This uses that the square root of triangular matrices can be computed directly. + result.resize(m_A.rows(), m_A.cols()); + typedef typename MatrixType::Index Index; + for (Index i = 0; i < m_A.rows(); i++) { + result.coeffRef(i,i) = internal::sqrt(T.coeff(i,i)); + } + for (Index j = 1; j < m_A.cols(); j++) { + for (Index i = j-1; i >= 0; i--) { + typedef typename MatrixType::Scalar Scalar; + // if i = j-1, then segment has length 0 so tmp = 0 + Scalar tmp = result.row(i).segment(i+1,j-i-1) * result.col(j).segment(i+1,j-i-1); + // denominator may be zero if original matrix is singular + result.coeffRef(i,j) = (T.coeff(i,j) - tmp) / (result.coeff(i,i) + result.coeff(j,j)); + } + } + + // Compute square root of m_A as U * result * U.adjoint() + MatrixType tmp; + tmp.noalias() = U * result.template triangularView<Upper>(); + result.noalias() = tmp * U.adjoint(); +} + +#endif // EIGEN_MATRIX_FUNCTION |