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::Index Index; 18 Index rows = m.rows(); 19 Index cols = m.cols(); 20 21 enum { 22 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 23 ColsAtCompileTime = MatrixType::ColsAtCompileTime 24 }; 25 26 typedef Matrix<JacobiScalar, 2, 1> JacobiVector; 27 28 const MatrixType a(MatrixType::Random(rows, cols)); 29 30 JacobiVector v = JacobiVector::Random().normalized(); 31 JacobiScalar c = v.x(), s = v.y(); 32 JacobiRotation<JacobiScalar> rot(c, s); 33 34 { 35 Index p = internal::random<Index>(0, rows-1); 36 Index q; 37 do { 38 q = internal::random<Index>(0, rows-1); 39 } while (q == p); 40 41 MatrixType b = a; 42 b.applyOnTheLeft(p, q, rot); 43 VERIFY_IS_APPROX(b.row(p), c * a.row(p) + numext::conj(s) * a.row(q)); 44 VERIFY_IS_APPROX(b.row(q), -s * a.row(p) + numext::conj(c) * a.row(q)); 45 } 46 47 { 48 Index p = internal::random<Index>(0, cols-1); 49 Index q; 50 do { 51 q = internal::random<Index>(0, cols-1); 52 } while (q == p); 53 54 MatrixType b = a; 55 b.applyOnTheRight(p, q, rot); 56 VERIFY_IS_APPROX(b.col(p), c * a.col(p) - s * a.col(q)); 57 VERIFY_IS_APPROX(b.col(q), numext::conj(s) * a.col(p) + numext::conj(c) * a.col(q)); 58 } 59 } 60 61 void test_jacobi() 62 { 63 for(int i = 0; i < g_repeat; i++) { 64 CALL_SUBTEST_1(( jacobi<Matrix3f, float>() )); 65 CALL_SUBTEST_2(( jacobi<Matrix4d, double>() )); 66 CALL_SUBTEST_3(( jacobi<Matrix4cf, float>() )); 67 CALL_SUBTEST_3(( jacobi<Matrix4cf, std::complex<float> >() )); 68 69 int r = internal::random<int>(2, internal::random<int>(1,EIGEN_TEST_MAX_SIZE)/2), 70 c = internal::random<int>(2, internal::random<int>(1,EIGEN_TEST_MAX_SIZE)/2); 71 CALL_SUBTEST_4(( jacobi<MatrixXf, float>(MatrixXf(r,c)) )); 72 CALL_SUBTEST_5(( jacobi<MatrixXcd, double>(MatrixXcd(r,c)) )); 73 CALL_SUBTEST_5(( jacobi<MatrixXcd, std::complex<double> >(MatrixXcd(r,c)) )); 74 // complex<float> is really important to test as it is the only way to cover conjugation issues in certain unaligned paths 75 CALL_SUBTEST_6(( jacobi<MatrixXcf, float>(MatrixXcf(r,c)) )); 76 CALL_SUBTEST_6(( jacobi<MatrixXcf, std::complex<float> >(MatrixXcf(r,c)) )); 77 78 TEST_SET_BUT_UNUSED_VARIABLE(r); 79 TEST_SET_BUT_UNUSED_VARIABLE(c); 80 } 81 } 82
