jacobi.cpp
Go to the documentation of this file.
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();
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 
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 
79  }
80 }
s
RealScalar s
Definition: level1_cplx_impl.h:126
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
b
Scalar * b
Definition: benchVecAdd.cpp:17
Eigen::JacobiRotation
Rotation given by a cosine-sine pair.
Definition: ForwardDeclarations.h:283
rot
int EIGEN_BLAS_FUNC() rot(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps)
Definition: level1_real_impl.h:79
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(jacobi)
Definition: jacobi.cpp:60
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
Eigen::numext::q
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:1984
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
Eigen::g_repeat
static int g_repeat
Definition: main.h:169
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
conj
AnnoyingScalar conj(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:104
CALL_SUBTEST_6
#define CALL_SUBTEST_6(FUNC)
Definition: split_test_helper.h:34
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
jacobi
void jacobi(const MatrixType &m=MatrixType())
Definition: jacobi.cpp:15
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
main.h
EIGEN_TEST_MAX_SIZE
#define EIGEN_TEST_MAX_SIZE
Definition: boostmultiprec.cpp:16
p
float * p
Definition: Tutorial_Map_using.cpp:9
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
TEST_SET_BUT_UNUSED_VARIABLE
#define TEST_SET_BUT_UNUSED_VARIABLE(X)
Definition: main.h:121
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:55