JacobiSVD.h revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba
17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2009-2010 Benoit Jacob <jacob.benoit.1@gmail.com> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#ifndef EIGEN_JACOBISVD_H 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define EIGEN_JACOBISVD_H 127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace Eigen { 147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace internal { 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// forward declaration (needed by ICC) 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// the empty body is required by MSVC 187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, int QRPreconditioner, 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct svd_precondition_2x2_block_to_be_real {}; 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/*** QR preconditioners (R-SVD) 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *** 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *** Their role is to reduce the problem of computing the SVD to the case of a square matrix. 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for 267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *** JacobiSVD which by itself is only able to work on square matrices. 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ***/ 287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezenum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols }; 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, int QRPreconditioner, int Case> 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct qr_preconditioner_should_do_anything 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum { a = MatrixType::RowsAtCompileTime != Dynamic && 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType::ColsAtCompileTime != Dynamic && 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime, 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez b = MatrixType::RowsAtCompileTime != Dynamic && 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType::ColsAtCompileTime != Dynamic && 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime, 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ret = !( (QRPreconditioner == NoQRPreconditioner) || 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez (Case == PreconditionIfMoreColsThanRows && bool(a)) || 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez (Case == PreconditionIfMoreRowsThanCols && bool(b)) ) 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, int QRPreconditioner, int Case, 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez> struct qr_preconditioner_impl {}; 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, int QRPreconditioner, int Case> 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false> 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezpublic: 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {} 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&) 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return false; 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/*** preconditioner using FullPivHouseholderQR ***/ 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType> 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezpublic: 687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RowsAtCompileTime = MatrixType::RowsAtCompileTime, 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType; 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd) 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.~QRType(); 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::new (&m_qr) QRType(svd.rows(), svd.cols()); 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(matrix.rows() > matrix.cols()) 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.compute(matrix); 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.m_computeFullU) m_qr.matrixQ().evalTo(svd.m_matrixU, m_workspace); 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return true; 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return false; 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezprivate: 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef FullPivHouseholderQR<MatrixType> QRType; 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QRType m_qr; 1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez WorkspaceType m_workspace; 1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType> 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezpublic: 1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum 1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RowsAtCompileTime = MatrixType::RowsAtCompileTime, 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ColsAtCompileTime = MatrixType::ColsAtCompileTime, 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Options = MatrixType::Options 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime> 1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez TransposeTypeWithSameStorageOrder; 1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd) 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.~QRType(); 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::new (&m_qr) QRType(svd.cols(), svd.rows()); 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_adjoint.resize(svd.cols(), svd.rows()); 1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(matrix.cols() > matrix.rows()) 1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_adjoint = matrix.adjoint(); 1387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.compute(m_adjoint); 1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.m_computeFullV) m_qr.matrixQ().evalTo(svd.m_matrixV, m_workspace); 1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); 1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return true; 1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else return false; 1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezprivate: 1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QRType m_qr; 1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez TransposeTypeWithSameStorageOrder m_adjoint; 1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typename internal::plain_row_type<MatrixType>::type m_workspace; 1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/*** preconditioner using ColPivHouseholderQR ***/ 1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType> 1567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 1577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezpublic: 1597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 1607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd) 1627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.~QRType(); 1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::new (&m_qr) QRType(svd.rows(), svd.cols()); 1677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 1697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); 1707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 1737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(matrix.rows() > matrix.cols()) 1757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.compute(matrix); 1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 1787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); 1797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if(svd.m_computeThinU) 1807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); 1827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); 1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); 1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return true; 1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return false; 1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezprivate: 1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef ColPivHouseholderQR<MatrixType> QRType; 1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QRType m_qr; 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typename internal::plain_col_type<MatrixType>::type m_workspace; 1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType> 1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezpublic: 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum 2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RowsAtCompileTime = MatrixType::RowsAtCompileTime, 2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ColsAtCompileTime = MatrixType::ColsAtCompileTime, 2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Options = MatrixType::Options 2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 2107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime> 2127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez TransposeTypeWithSameStorageOrder; 2137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd) 2157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.~QRType(); 2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::new (&m_qr) QRType(svd.cols(), svd.rows()); 2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 2227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_adjoint.resize(svd.cols(), svd.rows()); 2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) 2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(matrix.cols() > matrix.rows()) 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_adjoint = matrix.adjoint(); 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.compute(m_adjoint); 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if(svd.m_computeThinV) 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); 2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); 2397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); 2417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return true; 2427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else return false; 2447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezprivate: 2477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 2487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QRType m_qr; 2497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez TransposeTypeWithSameStorageOrder m_adjoint; 2507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typename internal::plain_row_type<MatrixType>::type m_workspace; 2517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 2527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/*** preconditioner using HouseholderQR ***/ 2547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType> 2567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> 2577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 2587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezpublic: 2597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 2607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd) 2627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) 2647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.~QRType(); 2667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::new (&m_qr) QRType(svd.rows(), svd.cols()); 2677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (svd.m_computeFullU) m_workspace.resize(svd.rows()); 2697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); 2707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) 2737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(matrix.rows() > matrix.cols()) 2757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.compute(matrix); 2777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); 2787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.m_computeFullU) m_qr.householderQ().evalTo(svd.m_matrixU, m_workspace); 2797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if(svd.m_computeThinU) 2807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); 2827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU, m_workspace); 2837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols()); 2857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return true; 2867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return false; 2887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezprivate: 2907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef HouseholderQR<MatrixType> QRType; 2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QRType m_qr; 2927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typename internal::plain_col_type<MatrixType>::type m_workspace; 2937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 2947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType> 2967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> 2977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 2987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezpublic: 2997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 3007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 3017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum 3027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RowsAtCompileTime = MatrixType::RowsAtCompileTime, 3047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ColsAtCompileTime = MatrixType::ColsAtCompileTime, 3057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 3067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 3077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Options = MatrixType::Options 3087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 3097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime> 3117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez TransposeTypeWithSameStorageOrder; 3127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd) 3147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) 3167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.~QRType(); 3187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ::new (&m_qr) QRType(svd.cols(), svd.rows()); 3197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (svd.m_computeFullV) m_workspace.resize(svd.cols()); 3217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (svd.m_computeThinV) m_workspace.resize(svd.rows()); 3227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_adjoint.resize(svd.cols(), svd.rows()); 3237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) 3267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(matrix.cols() > matrix.rows()) 3287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_adjoint = matrix.adjoint(); 3307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.compute(m_adjoint); 3317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); 3337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.m_computeFullV) m_qr.householderQ().evalTo(svd.m_matrixV, m_workspace); 3347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if(svd.m_computeThinV) 3357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); 3377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV, m_workspace); 3387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows()); 3407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return true; 3417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else return false; 3437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezprivate: 3467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef HouseholderQR<TransposeTypeWithSameStorageOrder> QRType; 3477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QRType m_qr; 3487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez TransposeTypeWithSameStorageOrder m_adjoint; 3497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typename internal::plain_row_type<MatrixType>::type m_workspace; 3507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 3517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/*** 2x2 SVD implementation 3537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *** 3547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *** JacobiSVD consists in performing a series of 2x2 SVD subproblems 3557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ***/ 3567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, int QRPreconditioner> 3587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false> 3597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 3607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; 3617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename SVD::Index Index; 3627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {} 3637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 3647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, int QRPreconditioner> 3667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true> 3677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 3687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef JacobiSVD<MatrixType, QRPreconditioner> SVD; 3697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 3707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::RealScalar RealScalar; 3717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename SVD::Index Index; 3727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q) 3737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 3757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar z; 3767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiRotation<Scalar> rot; 3777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p))); 3787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(n==0) 3797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); 3817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez work_matrix.row(p) *= z; 3827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); 3837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); 3847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez work_matrix.row(q) *= z; 3857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); 3867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 3887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rot.c() = conj(work_matrix.coeff(p,p)) / n; 3907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rot.s() = work_matrix.coeff(q,p) / n; 3917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez work_matrix.applyOnTheLeft(p,q,rot); 3927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.computeU()) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); 3937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(work_matrix.coeff(p,q) != Scalar(0)) 3947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); 3967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez work_matrix.col(q) *= z; 3977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.computeV()) svd.m_matrixV.col(q) *= z; 3987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(work_matrix.coeff(q,q) != Scalar(0)) 4007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); 4027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez work_matrix.row(q) *= z; 4037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); 4047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 4087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, typename RealScalar, typename Index> 4107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, 4117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiRotation<RealScalar> *j_left, 4127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiRotation<RealScalar> *j_right) 4137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 4147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 4157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<RealScalar,2,2> m; 4167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), 4177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); 4187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiRotation<RealScalar> rot1; 4197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar t = m.coeff(0,0) + m.coeff(1,1); 4207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar d = m.coeff(1,0) - m.coeff(0,1); 4217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(t == RealScalar(0)) 4227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rot1.c() = RealScalar(0); 4247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1); 4257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 4277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar u = d / t; 4297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u)); 4307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez rot1.s() = rot1.c() * u; 4317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m.applyOnTheLeft(0,1,rot1); 4337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez j_right->makeJacobi(m,0,1); 4347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez *j_left = rot1 * j_right->transpose(); 4357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 4367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace internal 4387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** \ingroup SVD_Module 4407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \class JacobiSVD 4437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix 4457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param MatrixType the type of the matrix of which we are computing the SVD decomposition 4477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param QRPreconditioner this optional parameter allows to specify the type of QR decomposition that will be used internally 4487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * for the R-SVD step for non-square matrices. See discussion of possible values below. 4497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product 4517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \f[ A = U S V^* \f] 4527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 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; 4537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 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 4547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * and right \em singular \em vectors of \a A respectively. 4557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Singular values are always sorted in decreasing order. 4577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 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. 4597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 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 4617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 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 4627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 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, 4637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 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. 4647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Here's an example demonstrating basic usage: 4667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \include JacobiSVD_basic.cpp 4677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Output: \verbinclude JacobiSVD_basic.out 4687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This JacobiSVD class is a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The downside is that it's slower than 4707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 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 4717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \a p is the greater dimension, meaning that it is still of the same order of complexity as the faster bidiagonalizing R-SVD algorithms. 4727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * In particular, like any R-SVD, it takes advantage of non-squareness in that its complexity is only linear in the greater dimension. 4737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to 4757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * terminate in finite (and reasonable) time. 4767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The possible values for QRPreconditioner are: 4787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR. 4797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR. 4807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Contrary to other QRs, it doesn't allow computing thin unitaries. 4817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \li HouseholderQRPreconditioner is the fastest, and less safe and accurate than the pivoting variants. It uses non-pivoting QR. 4827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This is very similar in safety and accuracy to the bidiagonalization process used by bidiagonalizing SVD algorithms (since bidiagonalization 4837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * is inherently non-pivoting). However the resulting SVD is still more reliable than bidiagonalizing SVDs because the Jacobi-based iterarive 4847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * process is more reliable than the optimized bidiagonal SVD iterations. 4857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \li NoQRPreconditioner allows not to use a QR preconditioner at all. This is useful if you know that you will only be computing 4867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * JacobiSVD decompositions of square matrices. Non-square matrices require a QR preconditioner. Using this option will result in 4877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * faster compilation and smaller executable code. It won't significantly speed up computation, since JacobiSVD is always checking 4887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * if QR preconditioning is needed before applying it anyway. 4897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 4907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa MatrixBase::jacobiSvd() 4917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 4927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType, int QRPreconditioner> 4937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezclass JacobiSVD : public SVDBase<_MatrixType> 4947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 4957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez public: 4967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef _MatrixType MatrixType; 4987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 4997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 5007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 5017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum { 5027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RowsAtCompileTime = MatrixType::RowsAtCompileTime, 5037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ColsAtCompileTime = MatrixType::ColsAtCompileTime, 5047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez DiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_DYNAMIC(RowsAtCompileTime,ColsAtCompileTime), 5057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 5067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 5077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(MaxRowsAtCompileTime,MaxColsAtCompileTime), 5087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixOptions = MatrixType::Options 5097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 5107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, 5127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime> 5137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixUType; 5147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, 5157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> 5167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixVType; 5177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType; 5187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename internal::plain_row_type<MatrixType>::type RowType; 5197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename internal::plain_col_type<MatrixType>::type ColType; 5207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, 5217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime> 5227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez WorkMatrixType; 5237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Default Constructor. 5257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 5267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The default constructor is useful in cases in which the user intends to 5277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * perform decompositions via JacobiSVD::compute(const MatrixType&). 5287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 5297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiSVD() 5307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : SVDBase<_MatrixType>::SVDBase() 5317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez {} 5327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Default Constructor with memory preallocation 5357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 5367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Like the default constructor but with preallocation of the internal data 5377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * according to the specified problem size. 5387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa JacobiSVD() 5397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 5407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0) 5417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : SVDBase<_MatrixType>::SVDBase() 5427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez allocate(rows, cols, computationOptions); 5447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Constructor performing the decomposition of given matrix. 5477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 5487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param matrix the matrix to decompose 5497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. 5507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU, 5517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * #ComputeFullV, #ComputeThinV. 5527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 5537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not 5547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * available with the (non-default) FullPivHouseholderQR preconditioner. 5557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 5567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0) 5577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : SVDBase<_MatrixType>::SVDBase() 5587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez compute(matrix, computationOptions); 5607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Method performing the decomposition of given matrix using custom options. 5637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 5647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param matrix the matrix to decompose 5657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. 5667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU, 5677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * #ComputeFullV, #ComputeThinV. 5687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 5697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not 5707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * available with the (non-default) FullPivHouseholderQR preconditioner. 5717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 5727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SVDBase<MatrixType>& compute(const MatrixType& matrix, unsigned int computationOptions); 5737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Method performing the decomposition of given matrix using current options. 5757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 5767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param matrix the matrix to decompose 5777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 5787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int). 5797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 5807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SVDBase<MatrixType>& compute(const MatrixType& matrix) 5817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return compute(matrix, this->m_computationOptions); 5837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A. 5867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 5877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param b the right-hand-side of the equation to solve. 5887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 5897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \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. 5907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 5917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving. 5927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$. 5937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 5947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Rhs> 5957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline const internal::solve_retval<JacobiSVD, Rhs> 5967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez solve(const MatrixBase<Rhs>& b) const 5977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(this->m_isInitialized && "JacobiSVD is not initialized."); 5997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(SVDBase<MatrixType>::computeU() && SVDBase<MatrixType>::computeV() && "JacobiSVD::solve() requires both unitaries U and V to be computed (thin unitaries suffice)."); 6007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return internal::solve_retval<JacobiSVD, Rhs>(*this, b.derived()); 6017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez private: 6067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void allocate(Index rows, Index cols, unsigned int computationOptions); 6077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez protected: 6097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez WorkMatrixType m_workMatrix; 6107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex> 6127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez friend struct internal::svd_precondition_2x2_block_to_be_real; 6137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything> 6147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez friend struct internal::qr_preconditioner_impl; 6157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols; 6177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows; 6187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 6197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, int QRPreconditioner> 6217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions) 6227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 6237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (SVDBase<MatrixType>::allocate(rows, cols, computationOptions)) return; 6247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (QRPreconditioner == FullPivHouseholderQRPreconditioner) 6267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(!(this->m_computeThinU || this->m_computeThinV) && 6287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " 6297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez "Use the ColPivHouseholderQR preconditioner instead."); 6307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_workMatrix.resize(this->m_diagSize, this->m_diagSize); 6337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(this->m_cols>this->m_rows) m_qr_precond_morecols.allocate(*this); 6357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(this->m_rows>this->m_cols) m_qr_precond_morerows.allocate(*this); 6367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 6377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, int QRPreconditioner> 6397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezSVDBase<MatrixType>& 6407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezJacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions) 6417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 6427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 6437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez allocate(matrix.rows(), matrix.cols(), computationOptions); 6447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations, 6467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // only worsening the precision of U and V as we accumulate more rotations 6477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon(); 6487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286) 6507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min(); 6517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ 6537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix)) 6557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_workMatrix = matrix.block(0,0,this->m_diagSize,this->m_diagSize); 6577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(this->m_computeFullU) this->m_matrixU.setIdentity(this->m_rows,this->m_rows); 6587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(this->m_computeThinU) this->m_matrixU.setIdentity(this->m_rows,this->m_diagSize); 6597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(this->m_computeFullV) this->m_matrixV.setIdentity(this->m_cols,this->m_cols); 6607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(this->m_computeThinV) this->m_matrixV.setIdentity(this->m_cols, this->m_diagSize); 6617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /*** step 2. The main Jacobi SVD iteration. ***/ 6647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool finished = false; 6667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez while(!finished) 6677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez finished = true; 6697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix 6717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index p = 1; p < this->m_diagSize; ++p) 6737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index q = 0; q < p; ++q) 6757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // if this 2x2 sub-matrix is not diagonal already... 6777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't 6787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // keep us iterating forever. Similarly, small denormal numbers are considered zero. 6797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::max; 6807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)), 6817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez abs(m_workMatrix.coeff(q,q)))); 6827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold) 6837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez finished = false; 6857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal 6877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q); 6887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiRotation<RealScalar> j_left, j_right; 6897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); 6907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // accumulate resulting Jacobi rotations 6927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_workMatrix.applyOnTheLeft(p,q,j_left); 6937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(SVDBase<MatrixType>::computeU()) this->m_matrixU.applyOnTheRight(p,q,j_left.transpose()); 6947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_workMatrix.applyOnTheRight(p,q,j_right); 6967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(SVDBase<MatrixType>::computeV()) this->m_matrixV.applyOnTheRight(p,q,j_right); 6977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/ 7037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index i = 0; i < this->m_diagSize; ++i) 7057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar a = abs(m_workMatrix.coeff(i,i)); 7077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez this->m_singularValues.coeffRef(i) = a; 7087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(SVDBase<MatrixType>::computeU() && (a!=RealScalar(0))) this->m_matrixU.col(i) *= this->m_workMatrix.coeff(i,i)/a; 7097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/ 7127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez this->m_nonzeroSingularValues = this->m_diagSize; 7147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index i = 0; i < this->m_diagSize; i++) 7157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index pos; 7177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar maxRemainingSingularValue = this->m_singularValues.tail(this->m_diagSize-i).maxCoeff(&pos); 7187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(maxRemainingSingularValue == RealScalar(0)) 7197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez this->m_nonzeroSingularValues = i; 7217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez break; 7227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(pos) 7247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez pos += i; 7267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::swap(this->m_singularValues.coeffRef(i), this->m_singularValues.coeffRef(pos)); 7277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(SVDBase<MatrixType>::computeU()) this->m_matrixU.col(pos).swap(this->m_matrixU.col(i)); 7287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(SVDBase<MatrixType>::computeV()) this->m_matrixV.col(pos).swap(this->m_matrixV.col(i)); 7297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez this->m_isInitialized = true; 7337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 7347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 7357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace internal { 7377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType, int QRPreconditioner, typename Rhs> 7387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstruct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> 7397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : solve_retval_base<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> 7407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef JacobiSVD<_MatrixType, QRPreconditioner> JacobiSVDType; 7427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EIGEN_MAKE_SOLVE_HELPERS(JacobiSVDType,Rhs) 7437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename Dest> void evalTo(Dest& dst) const 7457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 7467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(rhs().rows() == dec().rows()); 7477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // A = U S V^* 7497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // So A^{-1} = V S^{-1} U^* 7507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index diagSize = (std::min)(dec().rows(), dec().cols()); 7527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize); 7537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index nonzeroSingVals = dec().nonzeroSingularValues(); 7557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse(); 7567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez invertedSingVals.tail(diagSize - nonzeroSingVals).setZero(); 7577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dst = dec().matrixV().leftCols(diagSize) 7597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * invertedSingVals.asDiagonal() 7607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * dec().matrixU().leftCols(diagSize).adjoint() 7617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * rhs(); 7627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 7637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 7647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace internal 7657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** \svd_module 7677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 7687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \return the singular value decomposition of \c *this computed by two-sided 7697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Jacobi transformations. 7707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 7717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa class JacobiSVD 7727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 7737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Derived> 7747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezJacobiSVD<typename MatrixBase<Derived>::PlainObject> 7757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezMatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const 7767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 7777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return JacobiSVD<PlainObject>(*this, computationOptions); 7787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 7797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace Eigen 7817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 7827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif // EIGEN_JACOBISVD_H 783