triangular.cpp
Go to the documentation of this file.
1 // This file is triangularView of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifdef EIGEN_TEST_PART_100
11 # define EIGEN_NO_DEPRECATED_WARNING
12 #endif
13 
14 #include "main.h"
15 
16 
17 template<typename MatrixType> void triangular_deprecated(const MatrixType &m)
18 {
19  Index rows = m.rows();
20  Index cols = m.cols();
21  MatrixType m1, m2, m3, m4;
22  m1.setRandom(rows,cols);
23  m2.setRandom(rows,cols);
24  m3 = m1; m4 = m2;
25  // deprecated method:
26  m1.template triangularView<Eigen::Upper>().swap(m2);
27  // use this method instead:
28  m3.template triangularView<Eigen::Upper>().swap(m4.template triangularView<Eigen::Upper>());
29  VERIFY_IS_APPROX(m1,m3);
30  VERIFY_IS_APPROX(m2,m4);
31  // deprecated method:
32  m1.template triangularView<Eigen::Lower>().swap(m4);
33  // use this method instead:
34  m3.template triangularView<Eigen::Lower>().swap(m2.template triangularView<Eigen::Lower>());
35  VERIFY_IS_APPROX(m1,m3);
36  VERIFY_IS_APPROX(m2,m4);
37 }
38 
39 
40 template<typename MatrixType> void triangular_square(const MatrixType& m)
41 {
42  typedef typename MatrixType::Scalar Scalar;
43  typedef typename NumTraits<Scalar>::Real RealScalar;
44  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
45 
46  RealScalar largerEps = 10*test_precision<RealScalar>();
47 
48  Index rows = m.rows();
49  Index cols = m.cols();
50 
51  MatrixType m1 = MatrixType::Random(rows, cols),
52  m2 = MatrixType::Random(rows, cols),
53  m3(rows, cols),
54  m4(rows, cols),
55  r1(rows, cols),
56  r2(rows, cols);
57  VectorType v2 = VectorType::Random(rows);
58 
59  MatrixType m1up = m1.template triangularView<Upper>();
60  MatrixType m2up = m2.template triangularView<Upper>();
61 
62  if (rows*cols>1)
63  {
64  VERIFY(m1up.isUpperTriangular());
65  VERIFY(m2up.transpose().isLowerTriangular());
66  VERIFY(!m2.isLowerTriangular());
67  }
68 
69 // VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2);
70 
71  // test overloaded operator+=
72  r1.setZero();
73  r2.setZero();
74  r1.template triangularView<Upper>() += m1;
75  r2 += m1up;
76  VERIFY_IS_APPROX(r1,r2);
77 
78  // test overloaded operator=
79  m1.setZero();
80  m1.template triangularView<Upper>() = m2.transpose() + m2;
81  m3 = m2.transpose() + m2;
82  VERIFY_IS_APPROX(m3.template triangularView<Lower>().transpose().toDenseMatrix(), m1);
83 
84  // test overloaded operator=
85  m1.setZero();
86  m1.template triangularView<Lower>() = m2.transpose() + m2;
87  VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);
88 
89  VERIFY_IS_APPROX(m3.template triangularView<Lower>().conjugate().toDenseMatrix(),
90  m3.conjugate().template triangularView<Lower>().toDenseMatrix());
91 
92  m1 = MatrixType::Random(rows, cols);
93  for (int i=0; i<rows; ++i)
94  while (numext::abs2(m1(i,i))<RealScalar(1e-1)) m1(i,i) = internal::random<Scalar>();
95 
96  Transpose<MatrixType> trm4(m4);
97  // test back and forward substitution with a vector as the rhs
98  m3 = m1.template triangularView<Upper>();
99  VERIFY(v2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(v2)), largerEps));
100  m3 = m1.template triangularView<Lower>();
101  VERIFY(v2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(v2)), largerEps));
102  m3 = m1.template triangularView<Upper>();
103  VERIFY(v2.isApprox(m3 * (m1.template triangularView<Upper>().solve(v2)), largerEps));
104  m3 = m1.template triangularView<Lower>();
105  VERIFY(v2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(v2)), largerEps));
106 
107  // test back and forward substitution with a matrix as the rhs
108  m3 = m1.template triangularView<Upper>();
109  VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(m2)), largerEps));
110  m3 = m1.template triangularView<Lower>();
111  VERIFY(m2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(m2)), largerEps));
112  m3 = m1.template triangularView<Upper>();
113  VERIFY(m2.isApprox(m3 * (m1.template triangularView<Upper>().solve(m2)), largerEps));
114  m3 = m1.template triangularView<Lower>();
115  VERIFY(m2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(m2)), largerEps));
116 
117  // check M * inv(L) using in place API
118  m4 = m3;
119  m1.transpose().template triangularView<Eigen::Upper>().solveInPlace(trm4);
120  VERIFY_IS_APPROX(m4 * m1.template triangularView<Eigen::Lower>(), m3);
121 
122  // check M * inv(U) using in place API
123  m3 = m1.template triangularView<Upper>();
124  m4 = m3;
125  m3.transpose().template triangularView<Eigen::Lower>().solveInPlace(trm4);
126  VERIFY_IS_APPROX(m4 * m1.template triangularView<Eigen::Upper>(), m3);
127 
128  // check solve with unit diagonal
129  m3 = m1.template triangularView<UnitUpper>();
130  VERIFY(m2.isApprox(m3 * (m1.template triangularView<UnitUpper>().solve(m2)), largerEps));
131 
132 // VERIFY(( m1.template triangularView<Upper>()
133 // * m2.template triangularView<Upper>()).isUpperTriangular());
134 
135  // test swap
136  m1.setOnes();
137  m2.setZero();
138  m2.template triangularView<Upper>().swap(m1.template triangularView<Eigen::Upper>());
139  m3.setZero();
140  m3.template triangularView<Upper>().setOnes();
141  VERIFY_IS_APPROX(m2,m3);
142  VERIFY_RAISES_STATIC_ASSERT(m1.template triangularView<Eigen::Lower>().swap(m2.template triangularView<Eigen::Upper>()));
143 
144  m1.setRandom();
145  m3 = m1.template triangularView<Upper>();
146  Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic> m5(cols, internal::random<int>(1,20)); m5.setRandom();
147  Matrix<Scalar, Dynamic, MatrixType::RowsAtCompileTime> m6(internal::random<int>(1,20), rows); m6.setRandom();
148  VERIFY_IS_APPROX(m1.template triangularView<Upper>() * m5, m3*m5);
149  VERIFY_IS_APPROX(m6*m1.template triangularView<Upper>(), m6*m3);
150 
151  m1up = m1.template triangularView<Upper>();
152  VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up);
153  VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up);
154  VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(), m1up.adjoint());
155  VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(), m1up.adjoint());
156 
157  VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().diagonal(), m1.diagonal());
158 
159  m3.setRandom();
160  const MatrixType& m3c(m3);
161  VERIFY( is_same_type(m3c.template triangularView<Lower>(),m3.template triangularView<Lower>().template conjugateIf<false>()) );
162  VERIFY( is_same_type(m3c.template triangularView<Lower>().conjugate(),m3.template triangularView<Lower>().template conjugateIf<true>()) );
163  VERIFY_IS_APPROX(m3.template triangularView<Lower>().template conjugateIf<true>().toDenseMatrix(),
164  m3.conjugate().template triangularView<Lower>().toDenseMatrix());
165  VERIFY_IS_APPROX(m3.template triangularView<Lower>().template conjugateIf<false>().toDenseMatrix(),
166  m3.template triangularView<Lower>().toDenseMatrix());
167 
168  VERIFY( is_same_type(m3c.template selfadjointView<Lower>(),m3.template selfadjointView<Lower>().template conjugateIf<false>()) );
169  VERIFY( is_same_type(m3c.template selfadjointView<Lower>().conjugate(),m3.template selfadjointView<Lower>().template conjugateIf<true>()) );
170  VERIFY_IS_APPROX(m3.template selfadjointView<Lower>().template conjugateIf<true>().toDenseMatrix(),
171  m3.conjugate().template selfadjointView<Lower>().toDenseMatrix());
172  VERIFY_IS_APPROX(m3.template selfadjointView<Lower>().template conjugateIf<false>().toDenseMatrix(),
173  m3.template selfadjointView<Lower>().toDenseMatrix());
174 
175 }
176 
177 
178 template<typename MatrixType> void triangular_rect(const MatrixType& m)
179 {
180  typedef typename MatrixType::Scalar Scalar;
181  typedef typename NumTraits<Scalar>::Real RealScalar;
182  enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
183 
184  Index rows = m.rows();
185  Index cols = m.cols();
186 
187  MatrixType m1 = MatrixType::Random(rows, cols),
188  m2 = MatrixType::Random(rows, cols),
189  m3(rows, cols),
190  m4(rows, cols),
191  r1(rows, cols),
192  r2(rows, cols);
193 
194  MatrixType m1up = m1.template triangularView<Upper>();
195  MatrixType m2up = m2.template triangularView<Upper>();
196 
197  if (rows>1 && cols>1)
198  {
199  VERIFY(m1up.isUpperTriangular());
200  VERIFY(m2up.transpose().isLowerTriangular());
201  VERIFY(!m2.isLowerTriangular());
202  }
203 
204  // test overloaded operator+=
205  r1.setZero();
206  r2.setZero();
207  r1.template triangularView<Upper>() += m1;
208  r2 += m1up;
209  VERIFY_IS_APPROX(r1,r2);
210 
211  // test overloaded operator=
212  m1.setZero();
213  m1.template triangularView<Upper>() = 3 * m2;
214  m3 = 3 * m2;
215  VERIFY_IS_APPROX(m3.template triangularView<Upper>().toDenseMatrix(), m1);
216 
217 
218  m1.setZero();
219  m1.template triangularView<Lower>() = 3 * m2;
220  VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);
221 
222  m1.setZero();
223  m1.template triangularView<StrictlyUpper>() = 3 * m2;
224  VERIFY_IS_APPROX(m3.template triangularView<StrictlyUpper>().toDenseMatrix(), m1);
225 
226 
227  m1.setZero();
228  m1.template triangularView<StrictlyLower>() = 3 * m2;
229  VERIFY_IS_APPROX(m3.template triangularView<StrictlyLower>().toDenseMatrix(), m1);
230  m1.setRandom();
231  m2 = m1.template triangularView<Upper>();
232  VERIFY(m2.isUpperTriangular());
233  VERIFY(!m2.isLowerTriangular());
234  m2 = m1.template triangularView<StrictlyUpper>();
235  VERIFY(m2.isUpperTriangular());
236  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
237  m2 = m1.template triangularView<UnitUpper>();
238  VERIFY(m2.isUpperTriangular());
239  m2.diagonal().array() -= Scalar(1);
240  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
241  m2 = m1.template triangularView<Lower>();
242  VERIFY(m2.isLowerTriangular());
243  VERIFY(!m2.isUpperTriangular());
244  m2 = m1.template triangularView<StrictlyLower>();
245  VERIFY(m2.isLowerTriangular());
246  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
247  m2 = m1.template triangularView<UnitLower>();
248  VERIFY(m2.isLowerTriangular());
249  m2.diagonal().array() -= Scalar(1);
250  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
251  // test swap
252  m1.setOnes();
253  m2.setZero();
254  m2.template triangularView<Upper>().swap(m1.template triangularView<Eigen::Upper>());
255  m3.setZero();
256  m3.template triangularView<Upper>().setOnes();
257  VERIFY_IS_APPROX(m2,m3);
258 }
259 
260 void bug_159()
261 {
262  Matrix3d m = Matrix3d::Random().triangularView<Lower>();
263  EIGEN_UNUSED_VARIABLE(m)
264 }
265 
266 EIGEN_DECLARE_TEST(triangular)
267 {
268  int maxsize = (std::min)(EIGEN_TEST_MAX_SIZE,20);
269  for(int i = 0; i < g_repeat ; i++)
270  {
271  int r = internal::random<int>(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(r)
272  int c = internal::random<int>(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(c)
273 
274  CALL_SUBTEST_1( triangular_square(Matrix<float, 1, 1>()) );
275  CALL_SUBTEST_2( triangular_square(Matrix<float, 2, 2>()) );
276  CALL_SUBTEST_3( triangular_square(Matrix3d()) );
277  CALL_SUBTEST_4( triangular_square(Matrix<std::complex<float>,8, 8>()) );
278  CALL_SUBTEST_5( triangular_square(MatrixXcd(r,r)) );
279  CALL_SUBTEST_6( triangular_square(Matrix<float,Dynamic,Dynamic,RowMajor>(r, r)) );
280 
281  CALL_SUBTEST_7( triangular_rect(Matrix<float, 4, 5>()) );
282  CALL_SUBTEST_8( triangular_rect(Matrix<double, 6, 2>()) );
283  CALL_SUBTEST_9( triangular_rect(MatrixXcf(r, c)) );
284  CALL_SUBTEST_5( triangular_rect(MatrixXcd(r, c)) );
285  CALL_SUBTEST_6( triangular_rect(Matrix<float,Dynamic,Dynamic,RowMajor>(r, c)) );
286 
287  CALL_SUBTEST_100( triangular_deprecated(Matrix<float, 5, 7>()) );
288  CALL_SUBTEST_100( triangular_deprecated(MatrixXd(r,c)) );
289  }
290 
291  CALL_SUBTEST_1( bug_159() );
292 }
maxsize
const size_t maxsize
Definition: benchmark-blocking-sizes.cpp:49
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(triangular)
Definition: triangular.cpp:266
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::setRandom
Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > & setRandom(Index size)
Definition: Random.h:151
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
r2
static const double r2
Definition: testSmartRangeFactor.cpp:32
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
bug_159
void bug_159()
Definition: triangular.cpp:260
m1
Matrix3d m1
Definition: IOFormat.cpp:2
Eigen::is_same_type
internal::enable_if< internal::is_same< T1, T2 >::value, bool >::type is_same_type(const T1 &, const T2 &)
Definition: main.h:410
abs2
EIGEN_DEVICE_FUNC const EIGEN_STRONG_INLINE Abs2ReturnType abs2() const
Definition: ArrayCwiseUnaryOps.h:80
r1
static const double r1
Definition: testSmartRangeFactor.cpp:32
CALL_SUBTEST_9
#define CALL_SUBTEST_9(FUNC)
Definition: split_test_helper.h:52
triangular_square
void triangular_square(const MatrixType &m)
Definition: triangular.cpp:40
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
Eigen::Transpose
Expression of the transpose of a matrix.
Definition: Transpose.h:52
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
m2
MatrixType m2(n_dims)
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_UNUSED_VARIABLE
#define EIGEN_UNUSED_VARIABLE(var)
Definition: Macros.h:1076
CALL_SUBTEST_100
#define CALL_SUBTEST_100(FUNC)
Definition: split_test_helper.h:598
triangular_deprecated
void triangular_deprecated(const MatrixType &m)
Definition: triangular.cpp:17
triangular_rect
void triangular_rect(const MatrixType &m)
Definition: triangular.cpp:178
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
VERIFY_RAISES_STATIC_ASSERT
#define VERIFY_RAISES_STATIC_ASSERT(a)
Definition: main.h:344
CALL_SUBTEST_6
#define CALL_SUBTEST_6(FUNC)
Definition: split_test_helper.h:34
Eigen::Lower
@ Lower
Definition: Constants.h:209
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
main.h
EIGEN_TEST_MAX_SIZE
#define EIGEN_TEST_MAX_SIZE
Definition: boostmultiprec.cpp:16
v2
Vector v2
Definition: testSerializationBase.cpp:39
min
#define min(a, b)
Definition: datatypes.h:19
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
VectorType
Definition: FFTW.cpp:65
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
triangularView< Lower >
A triangularView< Lower >().adjoint().solveInPlace(B)
CALL_SUBTEST_7
#define CALL_SUBTEST_7(FUNC)
Definition: split_test_helper.h:40
CALL_SUBTEST_8
#define CALL_SUBTEST_8(FUNC)
Definition: split_test_helper.h:46
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
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
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
VERIFY
#define VERIFY(a)
Definition: main.h:380
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Cols
static const int Cols
Definition: testCallRecord.cpp:32

gtsam
Author(s):
autogenerated on Tue Jun 25 2024 03:07:37