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, 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typename MatrixType::Index i, typename MatrixType::Index j); 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute1x2offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typename MatrixType::Index i, typename MatrixType::Index j); 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute2x1offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typename MatrixType::Index i, typename MatrixType::Index j); 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute2x2offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 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, 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 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{ 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez result.resize(m_A.rows(), m_A.cols()); 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez computeDiagonalPartOfSqrt(result, m_A); 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez computeOffDiagonalPartOfSqrt(result, m_A); 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// pre: T is quasi-upper-triangular and sqrtT is a zero matrix of the same size 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// post: the diagonal blocks of sqrtT are the square roots of the diagonal blocks of T 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType>::computeDiagonalPartOfSqrt(MatrixType& sqrtT, 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& T) 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = m_A.rows(); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = 0; i < size; i++) { 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i == size - 1 || T.coeff(i+1, i) == 0) { 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(T(i,i) >= 0); 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez sqrtT.coeffRef(i,i) = sqrt(T.coeff(i,i)); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else { 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute2x2diagonalBlock(sqrtT, T, i); 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++i; 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// pre: T is quasi-upper-triangular and diagonal blocks of sqrtT are square root of diagonal blocks of T. 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// post: sqrtT is the square root of T. 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType>::computeOffDiagonalPartOfSqrt(MatrixType& sqrtT, 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& T) 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = m_A.rows(); 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = 1; j < size; j++) { 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (T.coeff(j, j-1) != 0) // if T(j-1:j, j-1:j) is a 2-by-2 block 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath continue; 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = j-1; i >= 0; i--) { 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i > 0 && T.coeff(i, i-1) != 0) // if T(i-1:i, i-1:i) is a 2-by-2 block 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath continue; 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool iBlockIs2x2 = (i < size - 1) && (T.coeff(i+1, i) != 0); 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool jBlockIs2x2 = (j < size - 1) && (T.coeff(j+1, j) != 0); 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (iBlockIs2x2 && jBlockIs2x2) 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute2x2offDiagonalBlock(sqrtT, T, i, j); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (iBlockIs2x2 && !jBlockIs2x2) 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute2x1offDiagonalBlock(sqrtT, T, i, j); 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (!iBlockIs2x2 && jBlockIs2x2) 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute1x2offDiagonalBlock(sqrtT, T, i, j); 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (!iBlockIs2x2 && !jBlockIs2x2) 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute1x1offDiagonalBlock(sqrtT, T, i, j); 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// pre: T.block(i,i,2,2) has complex conjugate eigenvalues 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// post: sqrtT.block(i,i,2,2) is square root of T.block(i,i,2,2) 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType> 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::compute2x2diagonalBlock(MatrixType& sqrtT, const MatrixType& T, typename MatrixType::Index i) 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO: This case (2-by-2 blocks with complex conjugate eigenvalues) is probably hidden somewhere 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // in EigenSolver. If we expose it, we could call it directly from here. 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> block = T.template block<2,2>(i,i); 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver<Matrix<Scalar,2,2> > es(block); 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sqrtT.template block<2,2>(i,i) 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath = (es.eigenvectors() * es.eigenvalues().cwiseSqrt().asDiagonal() * es.eigenvectors().inverse()).real(); 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// pre: block structure of T is such that (i,j) is a 1x1 block, 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// all blocks of sqrtT to left of and below (i,j) are correct 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// post: sqrtT(i,j) has the correct value 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType> 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::compute1x1offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Index i, typename MatrixType::Index j) 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar tmp = (sqrtT.row(i).segment(i+1,j-i-1) * sqrtT.col(j).segment(i+1,j-i-1)).value(); 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sqrtT.coeffRef(i,j) = (T.coeff(i,j) - tmp) / (sqrtT.coeff(i,i) + sqrtT.coeff(j,j)); 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// similar to compute1x1offDiagonalBlock() 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType> 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::compute1x2offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Index i, typename MatrixType::Index j) 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,1,2> rhs = T.template block<1,2>(i,j); 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (j-i > 1) 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs -= sqrtT.block(i, i+1, 1, j-i-1) * sqrtT.block(i+1, j, j-i-1, 2); 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> A = sqrtT.coeff(i,i) * Matrix<Scalar,2,2>::Identity(); 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath A += sqrtT.template block<2,2>(j,j).transpose(); 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sqrtT.template block<1,2>(i,j).transpose() = A.fullPivLu().solve(rhs.transpose()); 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// similar to compute1x1offDiagonalBlock() 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType> 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::compute2x1offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Index i, typename MatrixType::Index j) 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,1> rhs = T.template block<2,1>(i,j); 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (j-i > 2) 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 1); 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> A = sqrtT.coeff(j,j) * Matrix<Scalar,2,2>::Identity(); 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath A += sqrtT.template block<2,2>(i,i); 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sqrtT.template block<2,1>(i,j) = A.fullPivLu().solve(rhs); 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// similar to compute1x1offDiagonalBlock() 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType> 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::compute2x2offDiagonalBlock(MatrixType& sqrtT, const MatrixType& T, 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename MatrixType::Index i, typename MatrixType::Index j) 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> A = sqrtT.template block<2,2>(i,i); 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> B = sqrtT.template block<2,2>(j,j); 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> C = T.template block<2,2>(i,j); 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (j-i > 2) 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath C -= sqrtT.block(i, i+2, 2, j-i-2) * sqrtT.block(i+2, j, j-i-2, 2); 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,2,2> X; 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solveAuxiliaryEquation(X, A, B, C); 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sqrtT.template block<2,2>(i,j) = X; 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// solves the equation A X + X B = C where all matrices are 2-by-2 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename SmallMatrixType> 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootQuasiTriangular<MatrixType> 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ::solveAuxiliaryEquation(SmallMatrixType& X, const SmallMatrixType& A, 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const SmallMatrixType& B, const SmallMatrixType& C) 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_STATIC_ASSERT((internal::is_same<SmallMatrixType, Matrix<Scalar,2,2> >::value), 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_INTERNAL_ERROR_PLEASE_FILE_A_BUG_REPORT); 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,4,4> coeffMatrix = Matrix<Scalar,4,4>::Zero(); 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(0,0) = A.coeff(0,0) + B.coeff(0,0); 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(1,1) = A.coeff(0,0) + B.coeff(1,1); 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(2,2) = A.coeff(1,1) + B.coeff(0,0); 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(3,3) = A.coeff(1,1) + B.coeff(1,1); 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(0,1) = B.coeff(1,0); 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(0,2) = A.coeff(0,1); 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(1,0) = B.coeff(0,1); 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(1,3) = A.coeff(0,1); 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(2,0) = A.coeff(1,0); 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(2,3) = B.coeff(1,0); 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(3,1) = A.coeff(1,0); 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath coeffMatrix.coeffRef(3,2) = B.coeff(0,1); 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,4,1> rhs; 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs.coeffRef(0) = C.coeff(0,0); 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs.coeffRef(1) = C.coeff(0,1); 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs.coeffRef(2) = C.coeff(1,0); 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rhs.coeffRef(3) = C.coeff(1,1); 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,4,1> result; 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result = coeffMatrix.fullPivLu().solve(rhs); 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath X.coeffRef(0,0) = result.coeff(0); 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath X.coeffRef(0,1) = result.coeff(1); 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath X.coeffRef(1,0) = result.coeff(2); 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath X.coeffRef(1,1) = result.coeff(3); 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup MatrixFunctions_Module 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Class for computing matrix square roots of upper triangular matrices. 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType type of the argument of the matrix square root, 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * expected to be an instantiation of the Matrix class template. 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class computes the square root of the upper triangular matrix 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * stored in the upper triangular part (including the diagonal) of 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the matrix passed to the constructor. 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixSquareRoot, MatrixSquareRootQuasiTriangular 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass MatrixSquareRootTriangular 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRootTriangular(const MatrixType& A) 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_A(A) 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(A.rows() == A.cols()); 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Compute the matrix square root 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] result square root of \p A, as specified in the constructor. 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only the upper triangular part (including the diagonal) of 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p result is updated, the rest is not touched. See 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * MatrixBase::sqrt() for details on how this computation is 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * implemented. 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> void compute(ResultType &result); 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& m_A; 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename ResultType> 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixSquareRootTriangular<MatrixType>::compute(ResultType &result) 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 2847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Compute square root of m_A and store it in upper triangular part of result 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // This uses that the square root of triangular matrices can be computed directly. 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result.resize(m_A.rows(), m_A.cols()); 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = 0; i < m_A.rows(); i++) { 2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez result.coeffRef(i,i) = sqrt(m_A.coeff(i,i)); 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = 1; j < m_A.cols(); j++) { 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = j-1; i >= 0; i--) { 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // if i = j-1, then segment has length 0 so tmp = 0 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar tmp = (result.row(i).segment(i+1,j-i-1) * result.col(j).segment(i+1,j-i-1)).value(); 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // denominator may be zero if original matrix is singular 2997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez result.coeffRef(i,j) = (m_A.coeff(i,j) - tmp) / (result.coeff(i,i) + result.coeff(j,j)); 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup MatrixFunctions_Module 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Class for computing matrix square roots of general matrices. 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType type of the argument of the matrix square root, 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * expected to be an instantiation of the Matrix class template. 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixSquareRootTriangular, MatrixSquareRootQuasiTriangular, MatrixBase::sqrt() 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType, int IsComplex = NumTraits<typename internal::traits<MatrixType>::Scalar>::IsComplex> 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass MatrixSquareRoot 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor. 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] A matrix whose square root is to be computed. 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The class stores a reference to \p A, so it should not be 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * changed (or destroyed) before compute() is called. 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRoot(const MatrixType& A); 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Compute the matrix square root 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] result square root of \p A, as specified in the constructor. 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * See MatrixBase::sqrt() for details on how this computation is 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * implemented. 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> void compute(ResultType &result); 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// ********** Partial specialization for real matrices ********** 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass MatrixSquareRoot<MatrixType, 0> 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRoot(const MatrixType& A) 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_A(A) 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(A.rows() == A.cols()); 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> void compute(ResultType &result) 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute Schur decomposition of m_A 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealSchur<MatrixType> schurOfA(m_A); 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& T = schurOfA.matrixT(); 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& U = schurOfA.matrixU(); 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute square root of T 3587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType sqrtT = MatrixType::Zero(m_A.rows(), m_A.cols()); 3597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixSquareRootQuasiTriangular<MatrixType>(T).compute(sqrtT); 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute square root of m_A 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result = U * sqrtT * U.adjoint(); 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& m_A; 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// ********** Partial specialization for complex matrices ********** 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename MatrixType> 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass MatrixSquareRoot<MatrixType, 1> 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRoot(const MatrixType& A) 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_A(A) 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(A.rows() == A.cols()); 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> void compute(ResultType &result) 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute Schur decomposition of m_A 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const ComplexSchur<MatrixType> schurOfA(m_A); 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& T = schurOfA.matrixT(); 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& U = schurOfA.matrixU(); 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute square root of T 3917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType sqrtT; 3927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixSquareRootTriangular<MatrixType>(T).compute(sqrtT); 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute square root of m_A 3957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez result = U * (sqrtT.template triangularView<Upper>() * U.adjoint()); 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& m_A; 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup MatrixFunctions_Module 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Proxy for the matrix square root of some matrix (expression). 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam Derived Type of the argument to the matrix square root. 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class holds the argument to the matrix square root until it 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * is assigned or evaluated for some other reason (so the argument 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * should not be changed in the meantime). It is the return type of 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * MatrixBase::sqrt() and most of the time this is the only way it is 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * used. 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> class MatrixSquareRootReturnValue 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath: public ReturnByValue<MatrixSquareRootReturnValue<Derived> > 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Derived::Index Index; 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor. 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] src %Matrix (expression) forming the argument of the 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix square root. 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRootReturnValue(const Derived& src) : m_src(src) { } 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Compute the matrix square root. 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] result the matrix square root of \p src in the 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * constructor. 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline void evalTo(ResultType& result) const 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const typename Derived::PlainObject srcEvaluated = m_src.eval(); 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRoot<typename Derived::PlainObject> me(srcEvaluated); 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath me.compute(result); 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows() const { return m_src.rows(); } 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols() const { return m_src.cols(); } 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Derived& m_src; 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixSquareRootReturnValue& operator=(const MatrixSquareRootReturnValue&); 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<MatrixSquareRootReturnValue<Derived> > 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename Derived::PlainObject ReturnType; 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate <typename Derived> 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst MatrixSquareRootReturnValue<Derived> MatrixBase<Derived>::sqrt() const 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows() == cols()); 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return MatrixSquareRootReturnValue<Derived>(derived()); 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_MATRIX_FUNCTION 467