JacobiSVD.h revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba
1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com> 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_JACOBISVD_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_JACOBISVD_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// forward declaration (needed by ICC) 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// the empty body is required by MSVC 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int QRPreconditioner, 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct svd_precondition_2x2_block_to_be_real {}; 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*** QR preconditioners (R-SVD) 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *** 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *** Their role is to reduce the problem of computing the SVD to the case of a square matrix. 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *** JacobiSVD which by itself is only able to work on square matrices. 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ***/ 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathenum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols }; 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int QRPreconditioner, int Case> 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct qr_preconditioner_should_do_anything 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { a = MatrixType::RowsAtCompileTime != Dynamic && 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType::ColsAtCompileTime != Dynamic && 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime, 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath b = MatrixType::RowsAtCompileTime != Dynamic && 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType::ColsAtCompileTime != Dynamic && 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime, 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ret = !( (QRPreconditioner == NoQRPreconditioner) || 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (Case == PreconditionIfMoreColsThanRows && bool(a)) || 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (Case == PreconditionIfMoreRowsThanCols && bool(b)) ) 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int QRPreconditioner, int Case, 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath> struct qr_preconditioner_impl {}; 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int QRPreconditioner, int Case> 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false> 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {} 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&) 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return false; 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*** preconditioner using FullPivHouseholderQR ***/ 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType; 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd) 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.~QRType(); 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::new (&m_qr) QRType(svd.rows(), svd.cols()); 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(matrix.rows() > matrix.cols()) 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_qr.compute(matrix); 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace); 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return false; 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprivate: 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef FullPivHouseholderQR<MatrixType> QRType; 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QRType m_qr; 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath WorkspaceType m_workspace; 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime, 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Options = MatrixType::Options 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime> 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TransposeTypeWithSameStorageOrder; 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd) 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.~QRType(); 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::new (&m_qr) QRType(svd.cols(), svd.rows()); 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_adjoint.resize(svd.cols(), svd.rows()); 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(matrix.cols() > matrix.rows()) 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_adjoint = matrix.adjoint(); 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_qr.compute(m_adjoint); 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace); 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else return false; 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprivate: 1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QRType m_qr; 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TransposeTypeWithSameStorageOrder m_adjoint; 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::plain_row_type<MatrixType>::type m_workspace; 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*** preconditioner using ColPivHouseholderQR ***/ 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd) 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.~QRType(); 1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::new (&m_qr) QRType(svd.rows(), svd.cols()); 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(matrix.rows() > matrix.cols()) 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_qr.compute(matrix); 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(svd.m_computeThinU) 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return false; 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprivate: 1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef ColPivHouseholderQR<MatrixType> QRType; 1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QRType m_qr; 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::plain_col_type<MatrixType>::type m_workspace; 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime, 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Options = MatrixType::Options 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime> 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TransposeTypeWithSameStorageOrder; 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd) 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.~QRType(); 2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::new (&m_qr) QRType(svd.cols(), svd.rows()); 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_adjoint.resize(svd.cols(), svd.rows()); 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(matrix.cols() > matrix.rows()) 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_adjoint = matrix.adjoint(); 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_qr.compute(m_adjoint); 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(svd.m_computeThinV) 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else return false; 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprivate: 2477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 2487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QRType m_qr; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TransposeTypeWithSameStorageOrder m_adjoint; 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::plain_row_type<MatrixType>::type m_workspace; 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*** preconditioner using HouseholderQR ***/ 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd) 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 2657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.~QRType(); 2667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::new (&m_qr) QRType(svd.rows(), svd.cols()); 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(matrix.rows() > matrix.cols()) 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_qr.compute(matrix); 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(svd.m_computeThinU) 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols()); 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return false; 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprivate: 2907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef HouseholderQR<MatrixType> QRType; 2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QRType m_qr; 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::plain_col_type<MatrixType>::type m_workspace; 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathclass qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathpublic: 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime, 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Options = MatrixType::Options 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime> 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TransposeTypeWithSameStorageOrder; 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd) 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.~QRType(); 3187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::new (&m_qr) QRType(svd.cols(), svd.rows()); 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_adjoint.resize(svd.cols(), svd.rows()); 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(matrix.cols() > matrix.rows()) 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_adjoint = matrix.adjoint(); 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_qr.compute(m_adjoint); 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if(svd.m_computeThinV) 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows()); 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else return false; 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathprivate: 3467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 3477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QRType m_qr; 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TransposeTypeWithSameStorageOrder m_adjoint; 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typename internal::plain_row_type<MatrixType>::type m_workspace; 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/*** 2x2 SVD implementation 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *** 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *** JacobiSVD consists in performing a series of 2x2 SVD subproblems 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ***/ 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int QRPreconditioner> 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false> 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename SVD::Index Index; 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {} 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int QRPreconditioner> 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true> 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename SVD::Index Index; 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q) 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar z; 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation<Scalar> rot; 3777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p))); 3787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(n==0) 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath work_matrix.row(p) *= z; 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); 3847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(work_matrix.coeff(q,q)!=Scalar(0)) 3857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); 3877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez work_matrix.row(q) *= z; 3887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); 3897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // otherwise the second row is already zero, so we have nothing to do. 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rot.c() = conj(work_matrix.coeff(p,p)) / n; 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rot.s() = work_matrix.coeff(q,p) / n; 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath work_matrix.applyOnTheLeft(p,q,rot); 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(work_matrix.coeff(p,q) != Scalar(0)) 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath work_matrix.col(q) *= z; 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.computeV()) svd.m_matrixV.col(q) *= z; 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(work_matrix.coeff(q,q) != Scalar(0)) 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath work_matrix.row(q) *= z; 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, typename RealScalar, typename Index> 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation<RealScalar> *j_left, 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation<RealScalar> *j_right) 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 4197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 4207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<RealScalar,2,2> m; 4227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), 4237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation<RealScalar> rot1; 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar t = m.coeff(0,0) + m.coeff(1,1); 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar d = m.coeff(1,0) - m.coeff(0,1); 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(t == RealScalar(0)) 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rot1.c() = RealScalar(0); 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1); 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 4347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar t2d2 = numext::hypot(t,d); 4357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rot1.c() = abs(t)/t2d2; 4367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rot1.s() = d/t2d2; 4377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(t<RealScalar(0)) 4387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rot1.s() = -rot1.s(); 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m.applyOnTheLeft(0,1,rot1); 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath j_right->makeJacobi(m,0,1); 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *j_left = rot1 * j_right->transpose(); 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup SVD_Module 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class JacobiSVD 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param MatrixType the type of the matrix of which we are computing the SVD decomposition 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param QRPreconditioner this optional parameter allows to specify the type of QR decomposition that will be used internally 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * for the R-SVD step for non-square matrices. See discussion of possible values below. 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f[ A = U S V^* \f] 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * where \a U is a n-by-n unitary, \a V is a p-by-p unitary, and \a S is a n-by-p real positive matrix which is zero outside of its main diagonal; 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the diagonal entries of S are known as the \em singular \em values of \a A and the columns of \a U and \a V are known as the left 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and right \em singular \em vectors of \a A respectively. 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Singular values are always sorted in decreasing order. 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This JacobiSVD decomposition computes only the singular values by default. If you want \a U or \a V, you need to ask for them explicitly. 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * You can ask for only \em thin \a U or \a V to be computed, meaning the following. In case of a rectangular n-by-p matrix, letting \a m be the 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * smaller value among \a n and \a p, there are only \a m singular vectors; the remaining columns of \a U and \a V do not correspond to actual 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * singular vectors. Asking for \em thin \a U or \a V means asking for only their \a m first columns to be formed. So \a U is then a n-by-m matrix, 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and \a V is then a p-by-m matrix. Notice that thin \a U and \a V are all you need for (least squares) solving. 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Here's an example demonstrating basic usage: 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \include JacobiSVD_basic.cpp 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude JacobiSVD_basic.out 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This JacobiSVD class is a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The downside is that it's slower than 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * bidiagonalizing SVD algorithms for large square matrices; however its complexity is still \f$ O(n^2p) \f$ where \a n is the smaller dimension and 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \a p is the greater dimension, meaning that it is still of the same order of complexity as the faster bidiagonalizing R-SVD algorithms. 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * In particular, like any R-SVD, it takes advantage of non-squareness in that its complexity is only linear in the greater dimension. 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * terminate in finite (and reasonable) time. 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The possible values for QRPreconditioner are: 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR. 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR. 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Contrary to other QRs, it doesn't allow computing thin unitaries. 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \li HouseholderQRPreconditioner is the fastest, and less safe and accurate than the pivoting variants. It uses non-pivoting QR. 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is very similar in safety and accuracy to the bidiagonalization process used by bidiagonalizing SVD algorithms (since bidiagonalization 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * is inherently non-pivoting). However the resulting SVD is still more reliable than bidiagonalizing SVDs because the Jacobi-based iterarive 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * process is more reliable than the optimized bidiagonal SVD iterations. 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \li NoQRPreconditioner allows not to use a QR preconditioner at all. This is useful if you know that you will only be computing 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * JacobiSVD decompositions of square matrices. Non-square matrices require a QR preconditioner. Using this option will result in 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * faster compilation and smaller executable code. It won't significantly speed up computation, since JacobiSVD is always checking 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * if QR preconditioning is needed before applying it anyway. 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::jacobiSvd() 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int QRPreconditioner> class JacobiSVD 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime, 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime), 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime), 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixOptions = MatrixType::Options 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime> 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixUType; 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixVType; 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType; 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::plain_row_type<MatrixType>::type RowType; 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::plain_col_type<MatrixType>::type ColType; 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime> 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath WorkMatrixType; 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default Constructor. 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The default constructor is useful in cases in which the user intends to 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * perform decompositions via JacobiSVD::compute(const MatrixType&). 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD() 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_isInitialized(false), 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isAllocated(false), 5397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_usePrescribedThreshold(false), 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_computationOptions(0), 5417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_rows(-1), m_cols(-1), m_diagSize(0) 542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default Constructor with memory preallocation 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Like the default constructor but with preallocation of the internal data 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * according to the specified problem size. 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa JacobiSVD() 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0) 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_isInitialized(false), 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isAllocated(false), 5547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_usePrescribedThreshold(false), 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_computationOptions(0), 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_rows(-1), m_cols(-1) 557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath allocate(rows, cols, computationOptions); 559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor performing the decomposition of given matrix. 562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param matrix the matrix to decompose 564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. 565c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU, 566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * #ComputeFullV, #ComputeThinV. 567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not 569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * available with the (non-default) FullPivHouseholderQR preconditioner. 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0) 572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_isInitialized(false), 573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isAllocated(false), 5747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_usePrescribedThreshold(false), 575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_computationOptions(0), 576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_rows(-1), m_cols(-1) 577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix, computationOptions); 579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Method performing the decomposition of given matrix using custom options. 582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param matrix the matrix to decompose 584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. 585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU, 586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * #ComputeFullV, #ComputeThinV. 587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not 589c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * available with the (non-default) FullPivHouseholderQR preconditioner. 590c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 591c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions); 592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Method performing the decomposition of given matrix using current options. 594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param matrix the matrix to decompose 596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 597c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int). 598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 599c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD& compute(const MatrixType& matrix) 600c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 601c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return compute(matrix, m_computationOptions); 602c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 603c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 604c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the \a U matrix. 605c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 606c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, 607c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the U matrix is n-by-n if you asked for #ComputeFullU, and is n-by-m if you asked for #ComputeThinU. 608c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 609c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The \a m first columns of \a U are the left singular vectors of the matrix being decomposed. 610c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 611c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This method asserts that you asked for \a U to be computed. 612c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 613c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixUType& matrixU() const 614c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 615c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 616c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(computeU() && "This JacobiSVD decomposition didn't compute U. Did you ask for it?"); 617c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_matrixU; 618c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 619c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 620c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the \a V matrix. 621c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 622c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, 623c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the V matrix is p-by-p if you asked for #ComputeFullV, and is p-by-m if you asked for ComputeThinV. 624c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 625c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The \a m first columns of \a V are the right singular vectors of the matrix being decomposed. 626c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 627c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This method asserts that you asked for \a V to be computed. 628c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 629c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixVType& matrixV() const 630c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 631c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 632c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(computeV() && "This JacobiSVD decomposition didn't compute V. Did you ask for it?"); 633c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_matrixV; 634c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 635c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 636c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the vector of singular values. 637c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 638c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, the 639c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * returned vector has size \a m. Singular values are always sorted in decreasing order. 640c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 641c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const SingularValuesType& singularValues() const 642c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 643c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 644c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_singularValues; 645c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 646c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 647c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns true if \a U (full or thin) is asked for in this SVD decomposition */ 648c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline bool computeU() const { return m_computeFullU || m_computeThinU; } 649c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns true if \a V (full or thin) is asked for in this SVD decomposition */ 650c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline bool computeV() const { return m_computeFullV || m_computeThinV; } 651c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 652c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A. 653c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 654c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param b the right-hand-side of the equation to solve. 655c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 656c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \note Solving requires both U and V to be computed. Thin U and V are enough, there is no need for full U or V. 657c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 658c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving. 659c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$. 660c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 661c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs> 662c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const internal::solve_retval<JacobiSVD, Rhs> 663c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solve(const MatrixBase<Rhs>& b) const 664c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 665c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 666c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(computeU() && computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice)."); 667c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived()); 668c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 669c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 670c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the number of singular values that are not exactly 0 */ 671c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index nonzeroSingularValues() const 672c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 673c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 674c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_nonzeroSingularValues; 675c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 6767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns the rank of the matrix of which \c *this is the SVD. 6787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 6797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \note This method has to determine which singular values should be considered nonzero. 6807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * For that, it uses the threshold value that you can control by calling 6817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * setThreshold(const RealScalar&). 6827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 6837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline Index rank() const 6847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 6867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); 6877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(m_singularValues.size()==0) return 0; 6887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar premultiplied_threshold = m_singularValues.coeff(0) * threshold(); 6897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index i = m_nonzeroSingularValues-1; 6907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i; 6917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return i+1; 6927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** Allows to prescribe a threshold to be used by certain methods, such as rank() and solve(), 6957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * which need to determine when singular values are to be considered nonzero. 6967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This is not used for the SVD decomposition itself. 6977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 6987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * When it needs to get the threshold value, Eigen calls threshold(). 6997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The default is \c NumTraits<Scalar>::epsilon() 7007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 7017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param threshold The new value to use as the threshold. 7027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 7037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * A singular value will be considered nonzero if its value is strictly greater than 7047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \f$ \vert singular value \vert \leqslant threshold \times \vert max singular value \vert \f$. 7057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 7067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * If you want to come back to the default behavior, call setThreshold(Default_t) 7077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 7087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiSVD& setThreshold(const RealScalar& threshold) 7097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_usePrescribedThreshold = true; 7117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_prescribedThreshold = threshold; 7127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 7137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** Allows to come back to the default behavior, letting Eigen use its default formula for 7167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * determining the threshold. 7177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 7187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * You should pass the special object Eigen::Default as parameter here. 7197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \code svd.setThreshold(Eigen::Default); \endcode 7207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 7217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * See the documentation of setThreshold(const RealScalar&). 7227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 7237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiSVD& setThreshold(Default_t) 7247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_usePrescribedThreshold = false; 7267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 7277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** Returns the threshold that will be used by certain methods such as rank(). 7307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 7317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * See the documentation of setThreshold(const RealScalar&). 7327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 7337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar threshold() const 7347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized || m_usePrescribedThreshold); 7367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_usePrescribedThreshold ? m_prescribedThreshold 7377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : (std::max<Index>)(1,m_diagSize)*NumTraits<Scalar>::epsilon(); 7387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 739c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 740c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index rows() const { return m_rows; } 741c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index cols() const { return m_cols; } 742c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 743c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 744c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void allocate(Index rows, Index cols, unsigned int computationOptions); 745c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 746c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 747c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixUType m_matrixU; 748c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixVType m_matrixV; 749c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SingularValuesType m_singularValues; 750c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath WorkMatrixType m_workMatrix; 7517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold; 752c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_computeFullU, m_computeThinU; 753c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_computeFullV, m_computeThinV; 754c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath unsigned int m_computationOptions; 755c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize; 7567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar m_prescribedThreshold; 757c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 758c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex> 759c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath friend struct internal::svd_precondition_2x2_block_to_be_real; 760c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything> 761c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath friend struct internal::qr_preconditioner_impl; 762c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 763c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols; 764c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows; 765c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 766c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 767c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int QRPreconditioner> 768c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions) 769c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 770c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rows >= 0 && cols >= 0); 771c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 772c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_isAllocated && 773c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rows == m_rows && 774c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cols == m_cols && 775c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath computationOptions == m_computationOptions) 776c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 777c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return; 778c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 779c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 780c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_rows = rows; 781c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_cols = cols; 782c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = false; 783c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isAllocated = true; 784c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_computationOptions = computationOptions; 785c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_computeFullU = (computationOptions & ComputeFullU) != 0; 786c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_computeThinU = (computationOptions & ComputeThinU) != 0; 787c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_computeFullV = (computationOptions & ComputeFullV) != 0; 788c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_computeThinV = (computationOptions & ComputeThinV) != 0; 789c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U"); 790c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V"); 791c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) && 792c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns."); 793c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (QRPreconditioner == FullPivHouseholderQRPreconditioner) 794c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 795c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(!(m_computeThinU || m_computeThinV) && 796c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " 797c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath "Use the ColPivHouseholderQR preconditioner instead."); 798c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 799c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_diagSize = (std::min)(m_rows, m_cols); 800c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_singularValues.resize(m_diagSize); 8017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(RowsAtCompileTime==Dynamic) 8027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matrixU.resize(m_rows, m_computeFullU ? m_rows 8037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_computeThinU ? m_diagSize 8047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : 0); 8057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(ColsAtCompileTime==Dynamic) 8067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matrixV.resize(m_cols, m_computeFullV ? m_cols 8077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_computeThinV ? m_diagSize 8087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : 0); 809c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_workMatrix.resize(m_diagSize, m_diagSize); 810c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 811c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this); 812c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this); 813c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 814c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 815c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int QRPreconditioner> 816c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathJacobiSVD<MatrixType, QRPreconditioner>& 817c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathJacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions) 818c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 8197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 820c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath allocate(matrix.rows(), matrix.cols(), computationOptions); 821c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 822c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations, 823c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // only worsening the precision of U and V as we accumulate more rotations 824c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon(); 825c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 826c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286) 827c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min(); 828c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 829c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ 830c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 831c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix)) 832c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 833c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize); 834c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows); 835c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize); 836c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols); 837c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize); 838c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 8397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 8407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Scaling factor to reduce over/under-flows 8417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar scale = m_workMatrix.cwiseAbs().maxCoeff(); 8427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(scale==RealScalar(0)) scale = RealScalar(1); 8437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_workMatrix /= scale; 844c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 845c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /*** step 2. The main Jacobi SVD iteration. ***/ 846c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 847c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool finished = false; 848c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while(!finished) 849c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 850c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath finished = true; 851c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 852c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix 853c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 854c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index p = 1; p < m_diagSize; ++p) 855c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 856c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index q = 0; q < p; ++q) 857c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 858c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // if this 2x2 sub-matrix is not diagonal already... 859c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't 860c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // keep us iterating forever. Similarly, small denormal numbers are considered zero. 861c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using std::max; 8627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)), 8637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez abs(m_workMatrix.coeff(q,q)))); 8647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold) 865c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 866c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath finished = false; 867c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 868c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal 869c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q); 870c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation<RealScalar> j_left, j_right; 871c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); 872c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 873c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // accumulate resulting Jacobi rotations 874c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_workMatrix.applyOnTheLeft(p,q,j_left); 875c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); 876c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 877c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_workMatrix.applyOnTheRight(p,q,j_right); 878c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right); 879c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 880c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 881c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 882c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 883c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 884c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/ 885c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 886c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i = 0; i < m_diagSize; ++i) 887c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 8887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar a = abs(m_workMatrix.coeff(i,i)); 889c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_singularValues.coeffRef(i) = a; 890c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a; 891c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 892c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 893c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/ 894c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 895c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_nonzeroSingularValues = m_diagSize; 896c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i = 0; i < m_diagSize; i++) 897c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 898c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index pos; 899c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos); 900c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(maxRemainingSingularValue == RealScalar(0)) 901c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 902c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_nonzeroSingularValues = i; 903c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 904c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 905c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(pos) 906c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 907c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath pos += i; 908c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); 909c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU()) m_matrixU.col(pos).swap(m_matrixU.col(i)); 910c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i)); 911c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 912c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 9137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 9147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_singularValues *= scale; 915c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 916c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 917c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 918c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 919c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 920c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 921c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int QRPreconditioner, typename Rhs> 922c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> 923c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> 924c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 925c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType; 926c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs) 927c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 928c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void evalTo(Dest& dst) const 929c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 930c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rhs().rows() == dec().rows()); 931c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 932c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // A = U S V^* 933c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // So A^{-1} = V S^{-1} U^* 934c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 9357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar, Dynamic, Rhs::ColsAtCompileTime, 0, _MatrixType::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime> tmp; 9367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index rank = dec().rank(); 9377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 9387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp.noalias() = dec().matrixU().leftCols(rank).adjoint() * rhs(); 9397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp = dec().singularValues().head(rank).asDiagonal().inverse() * tmp; 9407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dst = dec().matrixV().leftCols(rank) * tmp; 941c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 942c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 943c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 944c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 945c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \svd_module 946c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 947c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \return the singular value decomposition of \c *this computed by two-sided 948c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Jacobi transformations. 949c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 950c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class JacobiSVD 951c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 952c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 953c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathJacobiSVD<typename MatrixBase<Derived>::PlainObject> 954c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const 955c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 956c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return JacobiSVD<PlainObject>(*this, computationOptions); 957c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 958c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 959c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 960c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 961c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_JACOBISVD_H 962