1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. Eigen itself is part of the KDE project. 3 // 4 // Copyright (C) 2008 Gael Guennebaud <g.gael (at) free.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #include "main.h" 11 #include <Eigen/SVD> 12 13 template<typename MatrixType> void svd(const MatrixType& m) 14 { 15 /* this test covers the following files: 16 SVD.h 17 */ 18 int rows = m.rows(); 19 int cols = m.cols(); 20 21 typedef typename MatrixType::Scalar Scalar; 22 typedef typename NumTraits<Scalar>::Real RealScalar; 23 MatrixType a = MatrixType::Random(rows,cols); 24 Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> b = 25 Matrix<Scalar, MatrixType::RowsAtCompileTime, 1>::Random(rows,1); 26 Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> x(cols,1), x2(cols,1); 27 28 RealScalar largerEps = test_precision<RealScalar>(); 29 if (ei_is_same_type<RealScalar,float>::ret) 30 largerEps = 1e-3f; 31 32 { 33 SVD<MatrixType> svd(a); 34 MatrixType sigma = MatrixType::Zero(rows,cols); 35 MatrixType matU = MatrixType::Zero(rows,rows); 36 sigma.block(0,0,cols,cols) = svd.singularValues().asDiagonal(); 37 matU.block(0,0,rows,cols) = svd.matrixU(); 38 VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose()); 39 } 40 41 42 if (rows==cols) 43 { 44 if (ei_is_same_type<RealScalar,float>::ret) 45 { 46 MatrixType a1 = MatrixType::Random(rows,cols); 47 a += a * a.adjoint() + a1 * a1.adjoint(); 48 } 49 SVD<MatrixType> svd(a); 50 svd.solve(b, &x); 51 VERIFY_IS_APPROX(a * x,b); 52 } 53 54 55 if(rows==cols) 56 { 57 SVD<MatrixType> svd(a); 58 MatrixType unitary, positive; 59 svd.computeUnitaryPositive(&unitary, &positive); 60 VERIFY_IS_APPROX(unitary * unitary.adjoint(), MatrixType::Identity(unitary.rows(),unitary.rows())); 61 VERIFY_IS_APPROX(positive, positive.adjoint()); 62 for(int i = 0; i < rows; i++) VERIFY(positive.diagonal()[i] >= 0); // cheap necessary (not sufficient) condition for positivity 63 VERIFY_IS_APPROX(unitary*positive, a); 64 65 svd.computePositiveUnitary(&positive, &unitary); 66 VERIFY_IS_APPROX(unitary * unitary.adjoint(), MatrixType::Identity(unitary.rows(),unitary.rows())); 67 VERIFY_IS_APPROX(positive, positive.adjoint()); 68 for(int i = 0; i < rows; i++) VERIFY(positive.diagonal()[i] >= 0); // cheap necessary (not sufficient) condition for positivity 69 VERIFY_IS_APPROX(positive*unitary, a); 70 } 71 } 72 73 void test_eigen2_svd() 74 { 75 for(int i = 0; i < g_repeat; i++) { 76 CALL_SUBTEST_1( svd(Matrix3f()) ); 77 CALL_SUBTEST_2( svd(Matrix4d()) ); 78 CALL_SUBTEST_3( svd(MatrixXf(7,7)) ); 79 CALL_SUBTEST_4( svd(MatrixXd(14,7)) ); 80 // complex are not implemented yet 81 // CALL_SUBTEST( svd(MatrixXcd(6,6)) ); 82 // CALL_SUBTEST( svd(MatrixXcf(3,3)) ); 83 SVD<MatrixXf> s; 84 MatrixXf m = MatrixXf::Random(10,1); 85 s.compute(m); 86 } 87 } 88
