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