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 }
Matrix3f m
int EIGEN_BLAS_FUNC() rot(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pc, RealScalar *ps)
#define CALL_SUBTEST_6(FUNC)
#define CALL_SUBTEST_4(FUNC)
Scalar * b
Definition: benchVecAdd.cpp:17
#define CALL_SUBTEST_3(FUNC)
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
Rotation given by a cosine-sine pair.
MatrixXf MatrixType
AnnoyingScalar conj(const AnnoyingScalar &x)
#define VERIFY_IS_APPROX(a, b)
#define CALL_SUBTEST_1(FUNC)
static int g_repeat
Definition: main.h:169
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Array< int, Dynamic, 1 > v
void jacobi(const MatrixType &m=MatrixType())
Definition: jacobi.cpp:15
RealScalar s
EIGEN_DEVICE_FUNC const Scalar & q
#define TEST_SET_BUT_UNUSED_VARIABLE(X)
Definition: main.h:121
#define CALL_SUBTEST_5(FUNC)
#define EIGEN_TEST_MAX_SIZE
float * p
EIGEN_DECLARE_TEST(jacobi)
Definition: jacobi.cpp:60
#define CALL_SUBTEST_2(FUNC)
The matrix class, also used for vectors and row-vectors.


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:27