ground/gcs/src/libs/eigen/test/jacobi.cpp

   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   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 }