17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com> 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr> 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr> 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr> 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla 137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed 147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// discard stack allocation as that too bypasses malloc 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define EIGEN_STACK_ALLOCATION_LIMIT 0 187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define EIGEN_RUNTIME_NO_MALLOC 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include "main.h" 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include <unsupported/Eigen/SVD> 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include <Eigen/LU> 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// check if "svd" is the good image of "m" 267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, typename SVD> 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid svd_check_full(const MatrixType& m, const SVD& svd) 287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index rows = m.rows(); 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index cols = m.cols(); 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum { 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RowsAtCompileTime = MatrixType::RowsAtCompileTime, 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ColsAtCompileTime = MatrixType::ColsAtCompileTime 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType; 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType; 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType sigma = MatrixType::Zero(rows, cols); 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez sigma.diagonal() = svd.singularValues().template cast<Scalar>(); 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixUType u = svd.matrixU(); 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixVType v = svd.matrixV(); 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(m, u * sigma * v.adjoint()); 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_UNITARY(u); 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_UNITARY(v); 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end svd_check_full 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Compare to a reference value 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, typename SVD> 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid svd_compare_to_full(const MatrixType& m, 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez unsigned int computationOptions, 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const SVD& referenceSvd) 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index rows = m.rows(); 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index cols = m.cols(); 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index diagSize = (std::min)(rows, cols); 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SVD svd(m, computationOptions); 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(svd.singularValues(), referenceSvd.singularValues()); 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(computationOptions & ComputeFullU) 687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU()); 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(computationOptions & ComputeThinU) 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize)); 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(computationOptions & ComputeFullV) 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV()); 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(computationOptions & ComputeThinV) 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize)); 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end svd_compare_to_full 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, typename SVD> 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid svd_solve(const MatrixType& m, unsigned int computationOptions) 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index rows = m.rows(); 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index cols = m.cols(); 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum { 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RowsAtCompileTime = MatrixType::RowsAtCompileTime, 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ColsAtCompileTime = MatrixType::ColsAtCompileTime 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, RowsAtCompileTime, Dynamic> RhsType; 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType; 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols)); 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SVD svd(m, computationOptions); 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SolutionType x = svd.solve(rhs); 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // evaluate normal equation which works also for least-squares solutions 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs); 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end svd_solve 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// test computations options 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 2 functions because Jacobisvd can return before the second function 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, typename SVD> 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid svd_test_computation_options_1(const MatrixType& m, const SVD& fullSvd) 1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_check_full< MatrixType, SVD >(m, fullSvd); 1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_solve< MatrixType, SVD >(m, ComputeFullU | ComputeFullV); 1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, typename SVD> 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid svd_test_computation_options_2(const MatrixType& m, const SVD& fullSvd) 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_compare_to_full< MatrixType, SVD >(m, ComputeFullU, fullSvd); 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_compare_to_full< MatrixType, SVD >(m, ComputeFullV, fullSvd); 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_compare_to_full< MatrixType, SVD >(m, 0, fullSvd); 1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (MatrixType::ColsAtCompileTime == Dynamic) { 1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // thin U/V are only available with dynamic number of columns 1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_compare_to_full< MatrixType, SVD >(m, ComputeFullU|ComputeThinV, fullSvd); 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_compare_to_full< MatrixType, SVD >(m, ComputeThinV, fullSvd); 1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_compare_to_full< MatrixType, SVD >(m, ComputeThinU|ComputeFullV, fullSvd); 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_compare_to_full< MatrixType, SVD >(m, ComputeThinU , fullSvd); 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_compare_to_full< MatrixType, SVD >(m, ComputeThinU|ComputeThinV, fullSvd); 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_solve<MatrixType, SVD>(m, ComputeFullU | ComputeThinV); 1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_solve<MatrixType, SVD>(m, ComputeThinU | ComputeFullV); 1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd_solve<MatrixType, SVD>(m, ComputeThinU | ComputeThinV); 1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index diagSize = (std::min)(m.rows(), m.cols()); 1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SVD svd(m, ComputeThinU | ComputeThinV); 1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(m, svd.matrixU().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint()); 1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 1387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, typename SVD> 1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid svd_verify_assert(const MatrixType& m) 1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index rows = m.rows(); 1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index cols = m.cols(); 1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum { 1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RowsAtCompileTime = MatrixType::RowsAtCompileTime, 1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ColsAtCompileTime = MatrixType::ColsAtCompileTime 1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, RowsAtCompileTime, 1> RhsType; 1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RhsType rhs(rows); 1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SVD svd; 1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.matrixU()) 1567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.singularValues()) 1577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.matrixV()) 1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.solve(rhs)) 1597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType a = MatrixType::Zero(rows, cols); 1607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez a.setZero(); 1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.compute(a, 0); 1627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.matrixU()) 1637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.matrixV()) 1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.singularValues(); 1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.solve(rhs)) 1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (ColsAtCompileTime == Dynamic) 1687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.compute(a, ComputeThinU); 1707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.matrixU(); 1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.matrixV()) 1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.solve(rhs)) 1737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.compute(a, ComputeThinV); 1747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.matrixV(); 1757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.matrixU()) 1767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.solve(rhs)) 1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 1797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU)) 1817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV)) 1827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 1847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// work around stupid msvc error when constructing at compile time an expression that involves 1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// a division by zero, even if the numeric type has floating point 1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Scalar> 1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezEIGEN_DONT_INLINE Scalar zero() { return Scalar(0); } 1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// workaround aggressive optimization in ICC 1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } 1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType, typename SVD> 1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid svd_inf_nan() 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // all this function does is verify we don't iterate infinitely on nan/inf values 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SVD svd; 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar some_inf = Scalar(1) / zero<Scalar>(); 2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf)); 2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV); 2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar some_nan = zero<Scalar> () / zero<Scalar> (); 2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(some_nan != some_nan); 2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV); 2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType m = MatrixType::Zero(10,10); 2107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf; 2117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.compute(m, ComputeFullU | ComputeFullV); 2127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m = MatrixType::Zero(10,10); 2147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_nan; 2157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.compute(m, ComputeFullU | ComputeFullV); 2167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename SVD> 2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid svd_preallocate() 2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 2227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Vector3f v(3.f, 2.f, 1.f); 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixXf m = v.asDiagonal(); 2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::set_is_malloc_allowed(false); 2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(VectorXf v(10);) 2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SVD svd; 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::set_is_malloc_allowed(true); 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd.compute(m); 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(svd.singularValues(), v); 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SVD svd2(3,3); 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::set_is_malloc_allowed(false); 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd2.compute(m); 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::set_is_malloc_allowed(true); 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(svd2.singularValues(), v); 2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd2.matrixU()); 2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_RAISES_ASSERT(svd2.matrixV()); 2397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd2.compute(m, ComputeFullU | ComputeFullV); 2407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); 2417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); 2427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::set_is_malloc_allowed(false); 2437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd2.compute(m); 2447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::set_is_malloc_allowed(true); 2457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SVD svd3(3,3,ComputeFullU|ComputeFullV); 2477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::set_is_malloc_allowed(false); 2487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd2.compute(m); 2497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::set_is_malloc_allowed(true); 2507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(svd2.singularValues(), v); 2517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); 2527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); 2537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::set_is_malloc_allowed(false); 2547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez svd2.compute(m, ComputeFullU|ComputeFullV); 2557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::set_is_malloc_allowed(true); 2567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 2577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 262