1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// discard stack allocation as that too bypasses malloc 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_STACK_ALLOCATION_LIMIT 0 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_RUNTIME_NO_MALLOC 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "main.h" 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <Eigen/SVD> 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int QRPreconditioner> 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid jacobisvd_check_full(const MatrixType& m, const JacobiSVD<MatrixType, QRPreconditioner>& svd) 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows = m.rows(); 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols = m.cols(); 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType; 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType; 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType sigma = MatrixType::Zero(rows,cols); 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sigma.diagonal() = svd.singularValues().template cast<Scalar>(); 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixUType u = svd.matrixU(); 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixVType v = svd.matrixV(); 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(m, u * sigma * v.adjoint()); 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_UNITARY(u); 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_UNITARY(v); 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int QRPreconditioner> 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid jacobisvd_compare_to_full(const MatrixType& m, 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath unsigned int computationOptions, 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const JacobiSVD<MatrixType, QRPreconditioner>& referenceSvd) 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows = m.rows(); 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols = m.cols(); 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index diagSize = (std::min)(rows, cols); 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions); 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(svd.singularValues(), referenceSvd.singularValues()); 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computationOptions & ComputeFullU) 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU()); 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computationOptions & ComputeThinU) 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize)); 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computationOptions & ComputeFullV) 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV()); 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computationOptions & ComputeThinV) 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize)); 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int QRPreconditioner> 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid jacobisvd_solve(const MatrixType& m, unsigned int computationOptions) 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::RealScalar RealScalar; 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows = m.rows(); 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols = m.cols(); 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, RowsAtCompileTime, Dynamic> RhsType; 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType; 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols)); 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions); 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(internal::is_same<RealScalar,double>::value) svd.setThreshold(1e-8); 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if(internal::is_same<RealScalar,float>::value) svd.setThreshold(1e-4); 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SolutionType x = svd.solve(rhs); 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar residual = (m*x-rhs).norm(); 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Check that there is no significantly better solution in the neighborhood of x 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(!test_isMuchSmallerThan(residual,rhs.norm())) 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // If the residual is very small, then we have an exact solution, so we are already good. 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int k=0;k<x.rows();++k) 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SolutionType y(x); 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y.row(k).array() += 2*NumTraits<RealScalar>::epsilon(); 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar residual_y = (m*y-rhs).norm(); 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); 1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y.row(k) = x.row(k).array() - 2*NumTraits<RealScalar>::epsilon(); 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez residual_y = (m*y-rhs).norm(); 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY( test_isApprox(residual_y,residual) || residual < residual_y ); 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // evaluate normal equation which works also for least-squares solutions 1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(internal::is_same<RealScalar,double>::value) 1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // This test is not stable with single precision. 1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // This is probably because squaring m signicantly affects the precision. 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs); 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check minimal norm solutions 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // generate a full-rank m x n problem with m<n 1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum { 1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RankAtCompileTime2 = ColsAtCompileTime==Dynamic ? Dynamic : (ColsAtCompileTime)/2+1, 1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RowsAtCompileTime3 = ColsAtCompileTime==Dynamic ? Dynamic : ColsAtCompileTime+1 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, RankAtCompileTime2, ColsAtCompileTime> MatrixType2; 1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, RankAtCompileTime2, 1> RhsType2; 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, ColsAtCompileTime, RankAtCompileTime2> MatrixType2T; 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index rank = RankAtCompileTime2==Dynamic ? internal::random<Index>(1,cols) : Index(RankAtCompileTime2); 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType2 m2(rank,cols); 1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int guard = 0; 1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez do { 1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m2.setRandom(); 1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } while(m2.jacobiSvd().setThreshold(test_precision<Scalar>()).rank()!=rank && (++guard)<10); 1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(guard<10); 1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RhsType2 rhs2 = RhsType2::Random(rank); 1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // use QR to find a reference minimal norm solution 1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez HouseholderQR<MatrixType2T> qr(m2.adjoint()); 1377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar,Dynamic,1> tmp = qr.matrixQR().topLeftCorner(rank,rank).template triangularView<Upper>().adjoint().solve(rhs2); 1387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp.conservativeResize(cols); 1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp.tail(cols-rank).setZero(); 1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SolutionType x21 = qr.householderQ() * tmp; 1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // now check with SVD 1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiSVD<MatrixType2, ColPivHouseholderQRPreconditioner> svd2(m2, computationOptions); 1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SolutionType x22 = svd2.solve(rhs2); 1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(m2*x21, rhs2); 1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(m2*x22, rhs2); 1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x21, x22); 1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Now check with a rank deficient matrix 1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, RowsAtCompileTime3, ColsAtCompileTime> MatrixType3; 1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, RowsAtCompileTime3, 1> RhsType3; 1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index rows3 = RowsAtCompileTime3==Dynamic ? internal::random<Index>(rank+1,2*cols) : Index(RowsAtCompileTime3); 1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar,RowsAtCompileTime3,Dynamic> C = Matrix<Scalar,RowsAtCompileTime3,Dynamic>::Random(rows3,rank); 1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType3 m3 = C * m2; 1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RhsType3 rhs3 = C * rhs2; 1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JacobiSVD<MatrixType3, ColPivHouseholderQRPreconditioner> svd3(m3, computationOptions); 1567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SolutionType x3 = svd3.solve(rhs3); 1577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(svd3.rank()!=rank) { 1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::cout << m3 << "\n\n"; 1597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::cout << svd3.singularValues().transpose() << "\n"; 1607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::cout << svd3.rank() << " == " << rank << "\n"; 1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::cout << x21.norm() << " == " << x3.norm() << "\n"; 1627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// VERIFY_IS_APPROX(m3*x3, rhs3); 1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(m3*x21, rhs3); 1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(m2*x3, rhs2); 1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x21, x3); 1687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int QRPreconditioner> 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid jacobisvd_test_all_computation_options(const MatrixType& m) 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols()) 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return; 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD<MatrixType, QRPreconditioner> fullSvd(m, ComputeFullU|ComputeFullV); 1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_check_full(m, fullSvd) )); 1787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV) )); 1797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez #if defined __INTEL_COMPILER 1817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // remark #111: statement is unreachable 1827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez #pragma warning disable 111 1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez #endif 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(QRPreconditioner == FullPivHouseholderQRPreconditioner) 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return; 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU, fullSvd) )); 1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullV, fullSvd) )); 1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_compare_to_full(m, 0, fullSvd) )); 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (MatrixType::ColsAtCompileTime == Dynamic) { 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // thin U/V are only available with dynamic number of columns 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd) )); 1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinV, fullSvd) )); 1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd) )); 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU , fullSvd) )); 1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd) )); 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV) )); 1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV) )); 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV) )); 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // test reconstruction 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index diagSize = (std::min)(m.rows(), m.cols()); 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD<MatrixType, QRPreconditioner> svd(m, ComputeThinU | ComputeThinV); 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(m, svd.matrixU().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint()); 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 2137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType m = a; 2147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(pickrandom) 2157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::RealScalar RealScalar; 2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index diagSize = (std::min)(a.rows(), a.cols()); 2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar s = std::numeric_limits<RealScalar>::max_exponent10/4; 2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez s = internal::random<RealScalar>(1,s); 2227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(diagSize); 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index k=0; k<diagSize; ++k) 2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s)); 2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m = Matrix<Scalar,Dynamic,Dynamic>::Random(a.rows(),diagSize) * d.asDiagonal() * Matrix<Scalar,Dynamic,Dynamic>::Random(diagSize,a.cols()); 2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // cancel some coeffs 2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index n = internal::random<Index>(0,m.size()-1); 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index i=0; i<n; ++i) 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = Scalar(0); 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m) )); 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m) )); 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m) )); 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m) )); 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> void jacobisvd_verify_assert(const MatrixType& m) 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows = m.rows(); 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols = m.cols(); 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, RowsAtCompileTime, 1> RhsType; 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RhsType rhs(rows); 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD<MatrixType> svd; 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.matrixU()) 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.singularValues()) 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.matrixV()) 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.solve(rhs)) 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType a = MatrixType::Zero(rows, cols); 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath a.setZero(); 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.compute(a, 0); 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.matrixU()) 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.matrixV()) 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.singularValues(); 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.solve(rhs)) 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (ColsAtCompileTime == Dynamic) 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.compute(a, ComputeThinU); 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.matrixU(); 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.matrixV()) 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.solve(rhs)) 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.compute(a, ComputeThinV); 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.matrixV(); 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.matrixU()) 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.solve(rhs)) 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner> svd_fullqr; 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeFullU|ComputeThinV)) 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeThinV)) 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeFullV)) 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU)) 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV)) 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid jacobisvd_method() 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { Size = MatrixType::RowsAtCompileTime }; 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<RealScalar, Size, 1> RealVecType; 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m = MatrixType::Identity(); 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(m.jacobiSvd().singularValues(), RealVecType::Ones()); 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(m.jacobiSvd().matrixU()); 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(m.jacobiSvd().matrixV()); 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(m.jacobiSvd(ComputeFullU|ComputeFullV).solve(m), m); 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// work around stupid msvc error when constructing at compile time an expression that involves 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// a division by zero, even if the numeric type has floating point 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathEIGEN_DONT_INLINE Scalar zero() { return Scalar(0); } 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// workaround aggressive optimization in ICC 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid jacobisvd_inf_nan() 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // all this function does is verify we don't iterate infinitely on nan/inf values 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD<MatrixType> svd; 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar some_inf = Scalar(1) / zero<Scalar>(); 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf)); 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV); 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar some_nan = zero<Scalar>() / zero<Scalar>(); 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY(some_nan != some_nan); 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV); 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m = MatrixType::Zero(10,10); 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf; 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.compute(m, ComputeFullU | ComputeFullV); 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m = MatrixType::Zero(10,10); 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_nan; 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.compute(m, ComputeFullU | ComputeFullV); 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Regression test for bug 286: JacobiSVD loops indefinitely with some 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// matrices containing denormal numbers. 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid jacobisvd_bug286() 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#if defined __INTEL_COMPILER 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// shut up warning #239: floating point underflow 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#pragma warning push 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#pragma warning disable 239 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix2d M; 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath M << -7.90884e-313, -4.94e-324, 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0, 5.60844e-313; 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#if defined __INTEL_COMPILER 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#pragma warning pop 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD<Matrix2d> svd; 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.compute(M); // just check we don't loop indefinitely 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid jacobisvd_preallocate() 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Vector3f v(3.f, 2.f, 1.f); 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXf m = v.asDiagonal(); 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::set_is_malloc_allowed(false); 3627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(VectorXf tmp(10);) 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD<MatrixXf> svd; 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::set_is_malloc_allowed(true); 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd.compute(m); 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(svd.singularValues(), v); 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD<MatrixXf> svd2(3,3); 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::set_is_malloc_allowed(false); 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd2.compute(m); 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::set_is_malloc_allowed(true); 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(svd2.singularValues(), v); 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd2.matrixU()); 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_RAISES_ASSERT(svd2.matrixV()); 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd2.compute(m, ComputeFullU | ComputeFullV); 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::set_is_malloc_allowed(false); 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd2.compute(m); 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::set_is_malloc_allowed(true); 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiSVD<MatrixXf> svd3(3,3,ComputeFullU|ComputeFullV); 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::set_is_malloc_allowed(false); 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd2.compute(m); 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::set_is_malloc_allowed(true); 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(svd2.singularValues(), v); 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::set_is_malloc_allowed(false); 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath svd2.compute(m, ComputeFullU|ComputeFullV); 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::set_is_malloc_allowed(true); 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid test_jacobisvd() 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) )); 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_4(( jacobisvd_verify_assert(Matrix4d()) )); 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_7(( jacobisvd_verify_assert(MatrixXf(10,12)) )); 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_8(( jacobisvd_verify_assert(MatrixXcd(7,5)) )); 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i = 0; i < g_repeat; i++) { 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix2cd m; 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m << 0, 1, 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0, 1; 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_1(( jacobisvd(m, false) )); 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m << 1, 0, 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1, 0; 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_1(( jacobisvd(m, false) )); 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix2d n; 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath n << 0, 0, 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0, 0; 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2(( jacobisvd(n, false) )); 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath n << 0, 0, 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0, 1; 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2(( jacobisvd(n, false) )); 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_3(( jacobisvd<Matrix3f>() )); 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_4(( jacobisvd<Matrix4d>() )); 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_5(( jacobisvd<Matrix<float,3,5> >() )); 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_6(( jacobisvd<Matrix<double,Dynamic,2> >(Matrix<double,Dynamic,2>(10,2)) )); 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int r = internal::random<int>(1, 30), 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath c = internal::random<int>(1, 30); 4257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez TEST_SET_BUT_UNUSED_VARIABLE(r) 4277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez TEST_SET_BUT_UNUSED_VARIABLE(c) 4287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST_10(( jacobisvd<MatrixXd>(MatrixXd(r,c)) )); 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(r,c)) )); 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(r,c)) )); 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (void) r; 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath (void) c; 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Test on inf/nan matrix 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_7( jacobisvd_inf_nan<MatrixXf>() ); 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3))) )); 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // test matrixbase method 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_1(( jacobisvd_method<Matrix2cd>() )); 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_3(( jacobisvd_method<Matrix3f>() )); 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Test problem size constructors 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_7( JacobiSVD<MatrixXf>(10,10) ); 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check that preallocation avoids subsequent mallocs 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_9( jacobisvd_preallocate() ); 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Regression check for bug 286 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST_2( jacobisvd_bug286() ); 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 455