1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> 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 EIGEN2_SVD_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN2_SVD_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup SVD_Module 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \nonstableyet 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class SVD 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Standard SVD decomposition of a matrix and associated features 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param MatrixType the type of the matrix of which we are computing the SVD decomposition 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with \c M \>= \c N. 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::SVD() 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> class SVD 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PacketSize = internal::packet_traits<Scalar>::size, 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath AlignmentMask = int(PacketSize)-1, 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MinSize = EIGEN_SIZE_MIN_PREFER_DYNAMIC(MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime) 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVector; 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> RowVector; 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MinSize> MatrixUType; 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixVType; 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, MinSize, 1> SingularValuesType; 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SVD() {} // a user who relied on compiler-generated default compiler reported problems with MSVC in 2.0.7 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SVD(const MatrixType& matrix) 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_matU(matrix.rows(), (std::min)(matrix.rows(), matrix.cols())), 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matV(matrix.cols(),matrix.cols()), 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma((std::min)(matrix.rows(),matrix.cols())) 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename OtherDerived, typename ResultType> 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool solve(const MatrixBase<OtherDerived> &b, ResultType* result) const; 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixUType& matrixU() const { return m_matU; } 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const SingularValuesType& singularValues() const { return m_sigma; } 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixVType& matrixV() const { return m_matV; } 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void compute(const MatrixType& matrix); 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SVD& sort(); 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename UnitaryType, typename PositiveType> 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void computeUnitaryPositive(UnitaryType *unitary, PositiveType *positive) const; 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename PositiveType, typename UnitaryType> 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void computePositiveUnitary(PositiveType *positive, UnitaryType *unitary) const; 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename RotationType, typename ScalingType> 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void computeRotationScaling(RotationType *unitary, ScalingType *positive) const; 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename ScalingType, typename RotationType> 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void computeScalingRotation(ScalingType *positive, RotationType *unitary) const; 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixUType m_matU; 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixVType m_matV; 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal */ 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SingularValuesType m_sigma; 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \note this code has been adapted from JAMA (public domain) 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid SVD<MatrixType>::compute(const MatrixType& matrix) 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int m = matrix.rows(); 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n = matrix.cols(); 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int nu = (std::min)(m,n); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ei_assert(m>=n && "In Eigen 2.0, SVD only works for MxN matrices with M>=N. Sorry!"); 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ei_assert(m>1 && "In Eigen 2.0, SVD doesn't work on 1x1 matrices"); 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU.resize(m, nu); 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU.setZero(); 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma.resize((std::min)(m,n)); 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matV.resize(n,n); 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowVector e(n); 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColVector work(m); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matA(matrix); 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const bool wantu = true; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const bool wantv = true; 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i=0, j=0, k=0; 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Reduce A to bidiagonal form, storing the diagonal elements 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // in s and the super-diagonal elements in e. 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int nct = (std::min)(m-1,n); 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int nrt = (std::max)(0,(std::min)(n-2,m)); 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (k = 0; k < (std::max)(nct,nrt); ++k) 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k < nct) 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the transformation for the k-th column and 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // place the k-th diagonal in m_sigma[k]. 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[k] = matA.col(k).end(m-k).norm(); 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_sigma[k] != 0.0) // FIXME 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (matA(k,k) < 0.0) 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[k] = -m_sigma[k]; 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA.col(k).end(m-k) /= m_sigma[k]; 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA(k,k) += 1.0; 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[k] = -m_sigma[k]; 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j = k+1; j < n; ++j) 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ((k < nct) && (m_sigma[k] != 0.0)) 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Apply the transformation. 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = matA.col(k).end(m-k).eigen2_dot(matA.col(j).end(m-k)); // FIXME dot product or cwise prod + .sum() ?? 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t = -t/matA(k,k); 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA.col(j).end(m-k) += t * matA.col(k).end(m-k); 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Place the k-th row of A into e for the 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // subsequent calculation of the row transformation. 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[j] = matA(k,j); 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Place the transformation in U for subsequent back multiplication. 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (wantu & (k < nct)) 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU.col(k).end(m-k) = matA.col(k).end(m-k); 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k < nrt) 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute the k-th row transformation and place the 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // k-th super-diagonal in e[k]. 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[k] = e.end(n-k-1).norm(); 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (e[k] != 0.0) 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (e[k+1] < 0.0) 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[k] = -e[k]; 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e.end(n-k-1) /= e[k]; 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[k+1] += 1.0; 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[k] = -e[k]; 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ((k+1 < m) & (e[k] != 0.0)) 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Apply the transformation. 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath work.end(m-k-1) = matA.corner(BottomRight,m-k-1,n-k-1) * e.end(n-k-1); 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j = k+1; j < n; ++j) 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA.col(j).end(m-k-1) += (-e[j]/e[k+1]) * work.end(m-k-1); 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Place the transformation in V for subsequent back multiplication. 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (wantv) 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matV.col(k).end(n-k-1) = e.end(n-k-1); 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Set up the final bidiagonal matrix or order p. 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int p = (std::min)(n,m+1); 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (nct < n) 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[nct] = matA(nct,nct); 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m < p) 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[p-1] = 0.0; 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (nrt+1 < p) 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[nrt] = matA(nrt,p-1); 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[p-1] = 0.0; 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // If required, generate U. 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (wantu) 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j = nct; j < nu; ++j) 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU.col(j).setZero(); 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU(j,j) = 1.0; 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (k = nct-1; k >= 0; k--) 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_sigma[k] != 0.0) 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j = k+1; j < nu; ++j) 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = m_matU.col(k).end(m-k).eigen2_dot(m_matU.col(j).end(m-k)); // FIXME is it really a dot product we want ? 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t = -t/m_matU(k,k); 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU.col(j).end(m-k) += t * m_matU.col(k).end(m-k); 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU.col(k).end(m-k) = - m_matU.col(k).end(m-k); 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU(k,k) = Scalar(1) + m_matU(k,k); 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k-1>0) 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU.col(k).start(k-1).setZero(); 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU.col(k).setZero(); 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU(k,k) = 1.0; 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // If required, generate V. 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (wantv) 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (k = n-1; k >= 0; k--) 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ((k < nrt) & (e[k] != 0.0)) 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j = k+1; j < nu; ++j) 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = m_matV.col(k).end(n-k-1).eigen2_dot(m_matV.col(j).end(n-k-1)); // FIXME is it really a dot product we want ? 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t = -t/m_matV(k+1,k); 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matV.col(j).end(n-k-1) += t * m_matV.col(k).end(n-k-1); 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matV.col(k).setZero(); 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matV(k,k) = 1.0; 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Main iteration loop for the singular values. 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int pp = p-1; 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int iter = 0; 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52)); 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (p > 0) 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int k=0; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int kase=0; 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Here is where a test for too many iterations would go. 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // This section of the program inspects for 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // negligible elements in the s and e arrays. On 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // completion the variables kase and k are set as follows. 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // kase = 1 if s(p) and e[k-1] are negligible and k<p 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // kase = 2 if s(k) is negligible and k<p 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // kase = 3 if e[k-1] is negligible, k<p, and 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // s(k), ..., s(p) are not negligible (qr step). 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // kase = 4 if e(p-1) is negligible (convergence). 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (k = p-2; k >= -1; --k) 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k == -1) 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ei_abs(e[k]) <= eps*(ei_abs(m_sigma[k]) + ei_abs(m_sigma[k+1]))) 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[k] = 0.0; 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k == p-2) 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath kase = 4; 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int ks; 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (ks = p-1; ks >= k; --ks) 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ks == k) 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = (ks != p ? ei_abs(e[ks]) : Scalar(0)) + (ks != k+1 ? ei_abs(e[ks-1]) : Scalar(0)); 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ei_abs(m_sigma[ks]) <= eps*t) 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[ks] = 0.0; 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ks == k) 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath kase = 3; 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (ks == p-1) 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath kase = 1; 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath kase = 2; 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath k = ks; 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++k; 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Perform the task indicated by kase. 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath switch (kase) 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Deflate negligible s(p). 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case 1: 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar f(e[p-2]); 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[p-2] = 0.0; 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j = p-2; j >= k; --j) 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar t(numext::hypot(m_sigma[j],f)); 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar cs(m_sigma[j]/t); 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar sn(f/t); 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[j] = t; 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (j != k) 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath f = -sn*e[j-1]; 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[j-1] = cs*e[j-1]; 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (wantv) 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (i = 0; i < n; ++i) 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t = cs*m_matV(i,j) + sn*m_matV(i,p-1); 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matV(i,p-1) = -sn*m_matV(i,j) + cs*m_matV(i,p-1); 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matV(i,j) = t; 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Split at negligible s(k). 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case 2: 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar f(e[k-1]); 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[k-1] = 0.0; 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j = k; j < p; ++j) 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar t(numext::hypot(m_sigma[j],f)); 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar cs( m_sigma[j]/t); 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar sn(f/t); 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[j] = t; 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath f = -sn*e[j]; 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[j] = cs*e[j]; 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (wantu) 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (i = 0; i < m; ++i) 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t = cs*m_matU(i,j) + sn*m_matU(i,k-1); 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU(i,k-1) = -sn*m_matU(i,j) + cs*m_matU(i,k-1); 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU(i,j) = t; 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Perform one qr step. 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case 3: 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Calculate the shift. 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar scale = (std::max)((std::max)((std::max)((std::max)( 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ei_abs(m_sigma[p-1]),ei_abs(m_sigma[p-2])),ei_abs(e[p-2])), 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ei_abs(m_sigma[k])),ei_abs(e[k])); 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar sp = m_sigma[p-1]/scale; 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar spm1 = m_sigma[p-2]/scale; 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar epm1 = e[p-2]/scale; 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar sk = m_sigma[k]/scale; 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar ek = e[k]/scale; 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/Scalar(2); 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar c = (sp*epm1)*(sp*epm1); 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar shift(0); 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ((b != 0.0) || (c != 0.0)) 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath shift = ei_sqrt(b*b + c); 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (b < 0.0) 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath shift = -shift; 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath shift = c/(b + shift); 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar f = (sk + sp)*(sk - sp) + shift; 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar g = sk*ek; 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Chase zeros. 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (j = k; j < p-1; ++j) 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar t = numext::hypot(f,g); 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar cs = f/t; 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar sn = g/t; 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (j != k) 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[j-1] = t; 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath f = cs*m_sigma[j] + sn*e[j]; 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[j] = cs*e[j] - sn*m_sigma[j]; 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath g = sn*m_sigma[j+1]; 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[j+1] = cs*m_sigma[j+1]; 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (wantv) 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (i = 0; i < n; ++i) 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t = cs*m_matV(i,j) + sn*m_matV(i,j+1); 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matV(i,j+1) = -sn*m_matV(i,j) + cs*m_matV(i,j+1); 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matV(i,j) = t; 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 4137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez t = numext::hypot(f,g); 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cs = f/t; 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sn = g/t; 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[j] = t; 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath f = cs*e[j] + sn*m_sigma[j+1]; 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[j+1] = -sn*e[j] + cs*m_sigma[j+1]; 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath g = sn*e[j+1]; 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[j+1] = cs*e[j+1]; 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (wantu && (j < m-1)) 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (i = 0; i < m; ++i) 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t = cs*m_matU(i,j) + sn*m_matU(i,j+1); 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU(i,j+1) = -sn*m_matU(i,j) + cs*m_matU(i,j+1); 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU(i,j) = t; 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath e[p-2] = f; 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iter = iter + 1; 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Convergence. 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath case 4: 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Make the singular values positive. 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_sigma[k] <= 0.0) 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[k] = m_sigma[k] < Scalar(0) ? -m_sigma[k] : Scalar(0); 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (wantv) 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matV.col(k).start(pp+1) = -m_matV.col(k).start(pp+1); 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Order the singular values. 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (k < pp) 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_sigma[k] >= m_sigma[k+1]) 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = m_sigma[k]; 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[k] = m_sigma[k+1]; 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma[k+1] = t; 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (wantv && (k < n-1)) 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matV.col(k).swap(m_matV.col(k+1)); 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (wantu && (k < m-1)) 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU.col(k).swap(m_matU.col(k+1)); 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++k; 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iter = 0; 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath p--; 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } // end big switch 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } // end iterations 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSVD<MatrixType>& SVD<MatrixType>::sort() 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int mu = m_matU.rows(); 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int mv = m_matV.rows(); 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int n = m_matU.cols(); 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i=0; i<n; ++i) 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int k = i; 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar p = m_sigma.coeff(i); 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int j=i+1; j<n; ++j) 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_sigma.coeff(j) > p) 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath k = j; 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath p = m_sigma.coeff(j); 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k != i) 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma.coeffRef(k) = m_sigma.coeff(i); // i.e. 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_sigma.coeffRef(i) = p; // swaps the i-th and the k-th elements 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int j = mu; 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int s=0; j!=0; ++s, --j) 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::swap(m_matU.coeffRef(s,i), m_matU.coeffRef(s,k)); 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath j = mv; 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int s=0; j!=0; ++s, --j) 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::swap(m_matV.coeffRef(s,i), m_matV.coeffRef(s,k)); 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \returns the solution of \f$ A x = b \f$ using the current SVD decomposition of A. 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The parts of the solution corresponding to zero singular values are ignored. 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::svd(), LU::solve(), LLT::solve() 510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename OtherDerived, typename ResultType> 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathbool SVD<MatrixType>::solve(const MatrixBase<OtherDerived> &b, ResultType* result) const 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 5157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ei_assert(b.rows() == m_matU.rows()); 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar maxVal = m_sigma.cwise().abs().maxCoeff(); 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int j=0; j<b.cols(); ++j) 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar,MatrixUType::RowsAtCompileTime,1> aux = m_matU.transpose() * b.col(j); 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i = 0; i <m_matU.cols(); ++i) 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar si = m_sigma.coeff(i); 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ei_isMuchSmallerThan(ei_abs(si),maxVal)) 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath aux.coeffRef(i) = 0; 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath aux.coeffRef(i) /= si; 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result->col(j) = m_matV * aux; 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Computes the polar decomposition of the matrix, as a product unitary x positive. 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * If either pointer is zero, the corresponding computation is skipped. 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only for square matrices. 541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa computePositiveUnitary(), computeRotationScaling() 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename UnitaryType, typename PositiveType> 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid SVD<MatrixType>::computeUnitaryPositive(UnitaryType *unitary, 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PositiveType *positive) const 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ei_assert(m_matU.cols() == m_matV.cols() && "Polar decomposition is only for square matrices"); 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(unitary) *unitary = m_matU * m_matV.adjoint(); 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(positive) *positive = m_matV * m_sigma.asDiagonal() * m_matV.adjoint(); 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Computes the polar decomposition of the matrix, as a product positive x unitary. 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * If either pointer is zero, the corresponding computation is skipped. 557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Only for square matrices. 559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa computeUnitaryPositive(), computeRotationScaling() 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename UnitaryType, typename PositiveType> 564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid SVD<MatrixType>::computePositiveUnitary(UnitaryType *positive, 565c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PositiveType *unitary) const 566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices"); 568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(unitary) *unitary = m_matU * m_matV.adjoint(); 569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(positive) *positive = m_matU * m_sigma.asDiagonal() * m_matU.adjoint(); 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** decomposes the matrix as a product rotation x scaling, the scaling being 573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * not necessarily positive. 574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * If either pointer is zero, the corresponding computation is skipped. 576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This method requires the Geometry module. 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa computeScalingRotation(), computeUnitaryPositive() 580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename RotationType, typename ScalingType> 583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid SVD<MatrixType>::computeRotationScaling(RotationType *rotation, ScalingType *scaling) const 584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices"); 586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar x = (m_matU * m_matV.adjoint()).determinant(); // so x has absolute value 1 587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> sv(m_sigma); 588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sv.coeffRef(0) *= x; 589c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(scaling) scaling->lazyAssign(m_matV * sv.asDiagonal() * m_matV.adjoint()); 590c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rotation) 591c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m(m_matU); 593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m.col(0) /= x; 594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rotation->lazyAssign(m * m_matV.adjoint()); 595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 597c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** decomposes the matrix as a product scaling x rotation, the scaling being 599c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * not necessarily positive. 600c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 601c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * If either pointer is zero, the corresponding computation is skipped. 602c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 603c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This method requires the Geometry module. 604c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 605c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa computeRotationScaling(), computeUnitaryPositive() 606c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 607c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 608c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename ScalingType, typename RotationType> 609c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid SVD<MatrixType>::computeScalingRotation(ScalingType *scaling, RotationType *rotation) const 610c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 611c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices"); 612c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar x = (m_matU * m_matV.adjoint()).determinant(); // so x has absolute value 1 613c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> sv(m_sigma); 614c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sv.coeffRef(0) *= x; 615c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(scaling) scaling->lazyAssign(m_matU * sv.asDiagonal() * m_matU.adjoint()); 616c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rotation) 617c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 618c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m(m_matU); 619c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m.col(0) /= x; 620c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rotation->lazyAssign(m * m_matV.adjoint()); 621c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 622c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 623c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 624c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 625c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \svd_module 626c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns the SVD decomposition of \c *this 627c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 628c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 629c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline SVD<typename MatrixBase<Derived>::PlainObject> 630c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixBase<Derived>::svd() const 631c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 632c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return SVD<PlainObject>(derived()); 633c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 634c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 635c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 636c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 637c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN2_SVD_H 638