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
13template<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
73void 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