1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2011 Jitse Niesen <jitse@maths.leeds.ac.uk> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_MATRIX_SQUARE_ROOT 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_MATRIX_SQUARE_ROOT 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup MatrixFunctions_Module 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Class for computing matrix square roots of upper quasi-triangular matrices. 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType type of the argument of the matrix square root, 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * expected to be an instantiation of the Matrix class template. 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class computes the square root of the upper quasi-triangular 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix stored in the upper Hessenberg part of the matrix passed to 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the constructor. 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixSquareRoot, MatrixSquareRootTriangular 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass MatrixSquareRootQuasiTriangular 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor. 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] A upper quasi-triangular matrix whose square root 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * is to be computed. 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The class stores a reference to \p A, so it should not be 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * changed (or destroyed) before compute() is called. 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRootQuasiTriangular(const MatrixType& A) 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_A(A) 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(A.rows() == A.cols()); 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Compute the matrix square root 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] result square root of \p A, as specified in the constructor. 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only the upper Hessenberg part of \p result is updated, the 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * rest is not touched. See MatrixBase::sqrt() for details on 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * how this computation is implemented. 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> void compute(ResultType &result); 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void computeDiagonalPartOfSqrt(MatrixType& sqrtT, const MatrixType& T); 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void computeOffDiagonalPartOfSqrt(MatrixType& sqrtT, const MatrixType& T); 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute2x2diagonalBlock(MatrixType& sqrtT, const MatrixType& T, typename MatrixType::Index i); 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute1x1offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Index i, typename MatrixType::Index j); 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute1x2offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Index i, typename MatrixType::Index j); 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute2x1offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Index i, typename MatrixType::Index j); 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute2x2offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Index i, typename MatrixType::Index j); 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename SmallMatrixType> 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void solveAuxiliaryEquation(SmallMatrixType& X, const SmallMatrixType& A, 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const SmallMatrixType& B, const SmallMatrixType& C); 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& m_A; 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename ResultType> 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType>::compute(ResultType &result) 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute Schur decomposition of m_A 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealSchur<MatrixType> schurOfA(m_A); 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& T = schurOfA.matrixT(); 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& U = schurOfA.matrixU(); 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute square root of T 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType sqrtT = MatrixType::Zero(m_A.rows(), m_A.rows()); 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath computeDiagonalPartOfSqrt(sqrtT, T); 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath computeOffDiagonalPartOfSqrt(sqrtT, T); 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute square root of m_A 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result = U * sqrtT * U.adjoint(); 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// pre: T is quasi-upper-triangular and sqrtT is a zero matrix of the same size 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// post: the diagonal blocks of sqrtT are the square roots of the diagonal blocks of T 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType>::computeDiagonalPartOfSqrt(MatrixType& sqrtT, 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& T) 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = m_A.rows(); 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = 0; i < size; i++) { 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i == size - 1 || T.coeff(i+1, i) == 0) { 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(T(i,i) > 0); 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sqrtT.coeffRef(i,i) = internal::sqrt(T.coeff(i,i)); 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else { 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute2x2diagonalBlock(sqrtT, T, i); 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++i; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// pre: T is quasi-upper-triangular and diagonal blocks of sqrtT are square root of diagonal blocks of T. 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// post: sqrtT is the square root of T. 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType>::computeOffDiagonalPartOfSqrt(MatrixType& sqrtT, 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& T) 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = m_A.rows(); 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = 1; j < size; j++) { 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (T.coeff(j, j-1) != 0) // if T(j-1:j, j-1:j) is a 2-by-2 block 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath continue; 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = j-1; i >= 0; i--) { 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i > 0 && T.coeff(i, i-1) != 0) // if T(i-1:i, i-1:i) is a 2-by-2 block 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath continue; 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool iBlockIs2x2 = (i < size - 1) && (T.coeff(i+1, i) != 0); 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool jBlockIs2x2 = (j < size - 1) && (T.coeff(j+1, j) != 0); 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (iBlockIs2x2 && jBlockIs2x2) 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute2x2offDiagonalBlock(sqrtT, T, i, j); 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (iBlockIs2x2 && !jBlockIs2x2) 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute2x1offDiagonalBlock(sqrtT, T, i, j); 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (!iBlockIs2x2 && jBlockIs2x2) 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute1x2offDiagonalBlock(sqrtT, T, i, j); 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (!iBlockIs2x2 && !jBlockIs2x2) 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute1x1offDiagonalBlock(sqrtT, T, i, j); 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// pre: T.block(i,i,2,2) has complex conjugate eigenvalues 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// post: sqrtT.block(i,i,2,2) is square root of T.block(i,i,2,2) 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType> 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::compute2x2diagonalBlock(MatrixType& sqrtT, const MatrixType& T, typename MatrixType::Index i) 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO: This case (2-by-2 blocks with complex conjugate eigenvalues) is probably hidden somewhere 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // in EigenSolver. If we expose it, we could call it directly from here. 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> block = T.template block<2,2>(i,i); 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver<Matrix<Scalar,2,2> > es(block); 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sqrtT.template block<2,2>(i,i) 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath = (es.eigenvectors() * es.eigenvalues().cwiseSqrt().asDiagonal() * es.eigenvectors().inverse()).real(); 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// pre: block structure of T is such that (i,j) is a 1x1 block, 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// all blocks of sqrtT to left of and below (i,j) are correct 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// post: sqrtT(i,j) has the correct value 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType> 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::compute1x1offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Index i, typename MatrixType::Index j) 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar tmp = (sqrtT.row(i).segment(i+1,j-i-1) * sqrtT.col(j).segment(i+1,j-i-1)).value(); 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sqrtT.coeffRef(i,j) = (T.coeff(i,j) - tmp) / (sqrtT.coeff(i,i) + sqrtT.coeff(j,j)); 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// similar to compute1x1offDiagonalBlock() 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType> 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::compute1x2offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Index i, typename MatrixType::Index j) 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,1,2> rhs = T.template block<1,2>(i,j); 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (j-i > 1) 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs -= sqrtT.block(i, i+1, 1, j-i-1) * sqrtT.block(i+1, j, j-i-1, 2); 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> A = sqrtT.coeff(i,i) * Matrix<Scalar,2,2>::Identity(); 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath A += sqrtT.template block<2,2>(j,j).transpose(); 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sqrtT.template block<1,2>(i,j).transpose() = A.fullPivLu().solve(rhs.transpose()); 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// similar to compute1x1offDiagonalBlock() 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType> 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::compute2x1offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Index i, typename MatrixType::Index j) 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,1> rhs = T.template block<2,1>(i,j); 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (j-i > 2) 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 1); 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> A = sqrtT.coeff(j,j) * Matrix<Scalar,2,2>::Identity(); 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath A += sqrtT.template block<2,2>(i,i); 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sqrtT.template block<2,1>(i,j) = A.fullPivLu().solve(rhs); 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// similar to compute1x1offDiagonalBlock() 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType> 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::compute2x2offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Index i, typename MatrixType::Index j) 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> A = sqrtT.template block<2,2>(i,i); 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> B = sqrtT.template block<2,2>(j,j); 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> C = T.template block<2,2>(i,j); 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (j-i > 2) 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath C -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 2); 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> X; 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solveAuxiliaryEquation(X, A, B, C); 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sqrtT.template block<2,2>(i,j) = X; 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// solves the equation A X + X B = C where all matrices are 2-by-2 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename SmallMatrixType> 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType> 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::solveAuxiliaryEquation(SmallMatrixType& X, const SmallMatrixType& A, 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const SmallMatrixType& B, const SmallMatrixType& C) 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_STATIC_ASSERT((internal::is_same<SmallMatrixType, Matrix<Scalar,2,2> >::value), 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_INTERNAL_ERROR_PLEASE_FILE_A_BUG_REPORT); 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,4,4> coeffMatrix = Matrix<Scalar,4,4>::Zero(); 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(0,0) = A.coeff(0,0) + B.coeff(0,0); 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(1,1) = A.coeff(0,0) + B.coeff(1,1); 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(2,2) = A.coeff(1,1) + B.coeff(0,0); 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(3,3) = A.coeff(1,1) + B.coeff(1,1); 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(0,1) = B.coeff(1,0); 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(0,2) = A.coeff(0,1); 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(1,0) = B.coeff(0,1); 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(1,3) = A.coeff(0,1); 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(2,0) = A.coeff(1,0); 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(2,3) = B.coeff(1,0); 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(3,1) = A.coeff(1,0); 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(3,2) = B.coeff(0,1); 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,4,1> rhs; 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs.coeffRef(0) = C.coeff(0,0); 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs.coeffRef(1) = C.coeff(0,1); 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs.coeffRef(2) = C.coeff(1,0); 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs.coeffRef(3) = C.coeff(1,1); 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,4,1> result; 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result = coeffMatrix.fullPivLu().solve(rhs); 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath X.coeffRef(0,0) = result.coeff(0); 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath X.coeffRef(0,1) = result.coeff(1); 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath X.coeffRef(1,0) = result.coeff(2); 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath X.coeffRef(1,1) = result.coeff(3); 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup MatrixFunctions_Module 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Class for computing matrix square roots of upper triangular matrices. 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType type of the argument of the matrix square root, 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * expected to be an instantiation of the Matrix class template. 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class computes the square root of the upper triangular matrix 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * stored in the upper triangular part (including the diagonal) of 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the matrix passed to the constructor. 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixSquareRoot, MatrixSquareRootQuasiTriangular 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass MatrixSquareRootTriangular 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRootTriangular(const MatrixType& A) 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_A(A) 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(A.rows() == A.cols()); 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Compute the matrix square root 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] result square root of \p A, as specified in the constructor. 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only the upper triangular part (including the diagonal) of 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p result is updated, the rest is not touched. See 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * MatrixBase::sqrt() for details on how this computation is 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * implemented. 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> void compute(ResultType &result); 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& m_A; 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename ResultType> 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootTriangular<MatrixType>::compute(ResultType &result) 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute Schur decomposition of m_A 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const ComplexSchur<MatrixType> schurOfA(m_A); 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& T = schurOfA.matrixT(); 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& U = schurOfA.matrixU(); 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute square root of T and store it in upper triangular part of result 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // This uses that the square root of triangular matrices can be computed directly. 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result.resize(m_A.rows(), m_A.cols()); 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = 0; i < m_A.rows(); i++) { 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result.coeffRef(i,i) = internal::sqrt(T.coeff(i,i)); 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = 1; j < m_A.cols(); j++) { 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = j-1; i >= 0; i--) { 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // if i = j-1, then segment has length 0 so tmp = 0 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar tmp = (result.row(i).segment(i+1,j-i-1) * result.col(j).segment(i+1,j-i-1)).value(); 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // denominator may be zero if original matrix is singular 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result.coeffRef(i,j) = (T.coeff(i,j) - tmp) / (result.coeff(i,i) + result.coeff(j,j)); 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute square root of m_A as U * result * U.adjoint() 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType tmp; 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp.noalias() = U * result.template triangularView<Upper>(); 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result.noalias() = tmp * U.adjoint(); 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup MatrixFunctions_Module 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Class for computing matrix square roots of general matrices. 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType type of the argument of the matrix square root, 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * expected to be an instantiation of the Matrix class template. 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixSquareRootTriangular, MatrixSquareRootQuasiTriangular, MatrixBase::sqrt() 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass MatrixSquareRoot 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor. 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] A matrix whose square root is to be computed. 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The class stores a reference to \p A, so it should not be 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * changed (or destroyed) before compute() is called. 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRoot(const MatrixType& A); 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Compute the matrix square root 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] result square root of \p A, as specified in the constructor. 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * See MatrixBase::sqrt() for details on how this computation is 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * implemented. 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> void compute(ResultType &result); 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// ********** Partial specialization for real matrices ********** 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass MatrixSquareRoot<MatrixType, 0> 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRoot(const MatrixType& A) 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_A(A) 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(A.rows() == A.cols()); 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> void compute(ResultType &result) 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute Schur decomposition of m_A 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealSchur<MatrixType> schurOfA(m_A); 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& T = schurOfA.matrixT(); 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& U = schurOfA.matrixU(); 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute square root of T 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRootQuasiTriangular<MatrixType> tmp(T); 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType sqrtT = MatrixType::Zero(m_A.rows(), m_A.rows()); 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp.compute(sqrtT); 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute square root of m_A 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result = U * sqrtT * U.adjoint(); 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& m_A; 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// ********** Partial specialization for complex matrices ********** 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass MatrixSquareRoot<MatrixType, 1> 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRoot(const MatrixType& A) 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_A(A) 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(A.rows() == A.cols()); 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> void compute(ResultType &result) 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute Schur decomposition of m_A 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const ComplexSchur<MatrixType> schurOfA(m_A); 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& T = schurOfA.matrixT(); 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& U = schurOfA.matrixU(); 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute square root of T 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRootTriangular<MatrixType> tmp(T); 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType sqrtT = MatrixType::Zero(m_A.rows(), m_A.rows()); 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp.compute(sqrtT); 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute square root of m_A 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result = U * sqrtT * U.adjoint(); 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& m_A; 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup MatrixFunctions_Module 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Proxy for the matrix square root of some matrix (expression). 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam Derived Type of the argument to the matrix square root. 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class holds the argument to the matrix square root until it 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * is assigned or evaluated for some other reason (so the argument 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * should not be changed in the meantime). It is the return type of 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * MatrixBase::sqrt() and most of the time this is the only way it is 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * used. 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> class MatrixSquareRootReturnValue 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath: public ReturnByValue<MatrixSquareRootReturnValue<Derived> > 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Derived::Index Index; 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor. 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] src %Matrix (expression) forming the argument of the 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix square root. 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRootReturnValue(const Derived& src) : m_src(src) { } 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Compute the matrix square root. 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] result the matrix square root of \p src in the 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * constructor. 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline void evalTo(ResultType& result) const 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const typename Derived::PlainObject srcEvaluated = m_src.eval(); 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRoot<typename Derived::PlainObject> me(srcEvaluated); 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath me.compute(result); 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows() const { return m_src.rows(); } 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols() const { return m_src.cols(); } 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Derived& m_src; 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRootReturnValue& operator=(const MatrixSquareRootReturnValue&); 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<MatrixSquareRootReturnValue<Derived> > 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Derived::PlainObject ReturnType; 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename Derived> 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt() const 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows() == cols()); 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return MatrixSquareRootReturnValue<Derived>(derived()); 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_MATRIX_FUNCTION 485