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@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
svd(const MatrixType & m)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
test_eigen2_svd()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