1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1 (at) gmail.com> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #include "main.h" 12 #include <Eigen/SVD> 13 14 template<typename MatrixType, typename JacobiScalar> 15 void jacobi(const MatrixType& m = MatrixType()) 16 { 17 typedef typename MatrixType::Scalar Scalar; 18 typedef typename MatrixType::Index Index; 19 Index rows = m.rows(); 20 Index cols = m.cols(); 21 22 enum { 23 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 24 ColsAtCompileTime = MatrixType::ColsAtCompileTime 25 }; 26 27 typedef Matrix<JacobiScalar, 2, 1> JacobiVector; 28 29 const MatrixType a(MatrixType::Random(rows, cols)); 30 31 JacobiVector v = JacobiVector::Random().normalized(); 32 JacobiScalar c = v.x(), s = v.y(); 33 JacobiRotation<JacobiScalar> rot(c, s); 34 35 { 36 Index p = internal::random<Index>(0, rows-1); 37 Index q; 38 do { 39 q = internal::random<Index>(0, rows-1); 40 } while (q == p); 41 42 MatrixType b = a; 43 b.applyOnTheLeft(p, q, rot); 44 VERIFY_IS_APPROX(b.row(p), c * a.row(p) + internal::conj(s) * a.row(q)); 45 VERIFY_IS_APPROX(b.row(q), -s * a.row(p) + internal::conj(c) * a.row(q)); 46 } 47 48 { 49 Index p = internal::random<Index>(0, cols-1); 50 Index q; 51 do { 52 q = internal::random<Index>(0, cols-1); 53 } while (q == p); 54 55 MatrixType b = a; 56 b.applyOnTheRight(p, q, rot); 57 VERIFY_IS_APPROX(b.col(p), c * a.col(p) - s * a.col(q)); 58 VERIFY_IS_APPROX(b.col(q), internal::conj(s) * a.col(p) + internal::conj(c) * a.col(q)); 59 } 60 } 61 62 void test_jacobi() 63 { 64 for(int i = 0; i < g_repeat; i++) { 65 CALL_SUBTEST_1(( jacobi<Matrix3f, float>() )); 66 CALL_SUBTEST_2(( jacobi<Matrix4d, double>() )); 67 CALL_SUBTEST_3(( jacobi<Matrix4cf, float>() )); 68 CALL_SUBTEST_3(( jacobi<Matrix4cf, std::complex<float> >() )); 69 70 int r = internal::random<int>(2, internal::random<int>(1,EIGEN_TEST_MAX_SIZE)/2), 71 c = internal::random<int>(2, internal::random<int>(1,EIGEN_TEST_MAX_SIZE)/2); 72 CALL_SUBTEST_4(( jacobi<MatrixXf, float>(MatrixXf(r,c)) )); 73 CALL_SUBTEST_5(( jacobi<MatrixXcd, double>(MatrixXcd(r,c)) )); 74 CALL_SUBTEST_5(( jacobi<MatrixXcd, std::complex<double> >(MatrixXcd(r,c)) )); 75 // complex<float> is really important to test as it is the only way to cover conjugation issues in certain unaligned paths 76 CALL_SUBTEST_6(( jacobi<MatrixXcf, float>(MatrixXcf(r,c)) )); 77 CALL_SUBTEST_6(( jacobi<MatrixXcf, std::complex<float> >(MatrixXcf(r,c)) )); 78 (void) r; 79 (void) c; 80 } 81 } 82
