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