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 #include "main.h"
11 
12 
13 
14 template<typename MatrixType> void triangular_square(const MatrixType& m)
15 {
16  typedef typename MatrixType::Scalar Scalar;
17  typedef typename NumTraits<Scalar>::Real RealScalar;
18  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
19 
20  RealScalar largerEps = 10*test_precision<RealScalar>();
21 
22  typename MatrixType::Index rows = m.rows();
23  typename MatrixType::Index cols = m.cols();
24 
25  MatrixType m1 = MatrixType::Random(rows, cols),
26  m2 = MatrixType::Random(rows, cols),
27  m3(rows, cols),
28  m4(rows, cols),
29  r1(rows, cols),
30  r2(rows, cols);
31  VectorType v2 = VectorType::Random(rows);
32 
33  MatrixType m1up = m1.template triangularView<Upper>();
34  MatrixType m2up = m2.template triangularView<Upper>();
35 
36  if (rows*cols>1)
37  {
38  VERIFY(m1up.isUpperTriangular());
39  VERIFY(m2up.transpose().isLowerTriangular());
40  VERIFY(!m2.isLowerTriangular());
41  }
42 
43 // VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2);
44 
45  // test overloaded operator+=
46  r1.setZero();
47  r2.setZero();
48  r1.template triangularView<Upper>() += m1;
49  r2 += m1up;
50  VERIFY_IS_APPROX(r1,r2);
51 
52  // test overloaded operator=
53  m1.setZero();
54  m1.template triangularView<Upper>() = m2.transpose() + m2;
55  m3 = m2.transpose() + m2;
56  VERIFY_IS_APPROX(m3.template triangularView<Lower>().transpose().toDenseMatrix(), m1);
57 
58  // test overloaded operator=
59  m1.setZero();
60  m1.template triangularView<Lower>() = m2.transpose() + m2;
61  VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);
62 
63  VERIFY_IS_APPROX(m3.template triangularView<Lower>().conjugate().toDenseMatrix(),
64  m3.conjugate().template triangularView<Lower>().toDenseMatrix());
65 
66  m1 = MatrixType::Random(rows, cols);
67  for (int i=0; i<rows; ++i)
68  while (numext::abs2(m1(i,i))<RealScalar(1e-1)) m1(i,i) = internal::random<Scalar>();
69 
70  Transpose<MatrixType> trm4(m4);
71  // test back and forward subsitution with a vector as the rhs
72  m3 = m1.template triangularView<Upper>();
73  VERIFY(v2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(v2)), largerEps));
74  m3 = m1.template triangularView<Lower>();
75  VERIFY(v2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(v2)), largerEps));
76  m3 = m1.template triangularView<Upper>();
77  VERIFY(v2.isApprox(m3 * (m1.template triangularView<Upper>().solve(v2)), largerEps));
78  m3 = m1.template triangularView<Lower>();
79  VERIFY(v2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(v2)), largerEps));
80 
81  // test back and forward substitution with a matrix as the rhs
82  m3 = m1.template triangularView<Upper>();
83  VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView<Lower>().solve(m2)), largerEps));
84  m3 = m1.template triangularView<Lower>();
85  VERIFY(m2.isApprox(m3.transpose() * (m1.transpose().template triangularView<Upper>().solve(m2)), largerEps));
86  m3 = m1.template triangularView<Upper>();
87  VERIFY(m2.isApprox(m3 * (m1.template triangularView<Upper>().solve(m2)), largerEps));
88  m3 = m1.template triangularView<Lower>();
89  VERIFY(m2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView<Lower>().solve(m2)), largerEps));
90 
91  // check M * inv(L) using in place API
92  m4 = m3;
93  m1.transpose().template triangularView<Eigen::Upper>().solveInPlace(trm4);
94  VERIFY_IS_APPROX(m4 * m1.template triangularView<Eigen::Lower>(), m3);
95 
96  // check M * inv(U) using in place API
97  m3 = m1.template triangularView<Upper>();
98  m4 = m3;
99  m3.transpose().template triangularView<Eigen::Lower>().solveInPlace(trm4);
100  VERIFY_IS_APPROX(m4 * m1.template triangularView<Eigen::Upper>(), m3);
101 
102  // check solve with unit diagonal
103  m3 = m1.template triangularView<UnitUpper>();
104  VERIFY(m2.isApprox(m3 * (m1.template triangularView<UnitUpper>().solve(m2)), largerEps));
105 
106 // VERIFY(( m1.template triangularView<Upper>()
107 // * m2.template triangularView<Upper>()).isUpperTriangular());
108 
109  // test swap
110  m1.setOnes();
111  m2.setZero();
112  m2.template triangularView<Upper>().swap(m1);
113  m3.setZero();
114  m3.template triangularView<Upper>().setOnes();
115  VERIFY_IS_APPROX(m2,m3);
116 
117  m1.setRandom();
118  m3 = m1.template triangularView<Upper>();
119  Matrix<Scalar, MatrixType::ColsAtCompileTime, Dynamic> m5(cols, internal::random<int>(1,20)); m5.setRandom();
120  Matrix<Scalar, Dynamic, MatrixType::RowsAtCompileTime> m6(internal::random<int>(1,20), rows); m6.setRandom();
121  VERIFY_IS_APPROX(m1.template triangularView<Upper>() * m5, m3*m5);
122  VERIFY_IS_APPROX(m6*m1.template triangularView<Upper>(), m6*m3);
123 
124  m1up = m1.template triangularView<Upper>();
125  VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up);
126  VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Upper>().toDenseMatrix(), m1up);
127  VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(), m1up.adjoint());
128  VERIFY_IS_APPROX(m1up.template selfadjointView<Upper>().template triangularView<Lower>().toDenseMatrix(), m1up.adjoint());
129 
130  VERIFY_IS_APPROX(m1.template selfadjointView<Upper>().diagonal(), m1.diagonal());
131 
132 }
133 
134 
135 template<typename MatrixType> void triangular_rect(const MatrixType& m)
136 {
137  typedef typename MatrixType::Scalar Scalar;
138  typedef typename NumTraits<Scalar>::Real RealScalar;
139  enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
140 
141  Index rows = m.rows();
142  Index cols = m.cols();
143 
144  MatrixType m1 = MatrixType::Random(rows, cols),
145  m2 = MatrixType::Random(rows, cols),
146  m3(rows, cols),
147  m4(rows, cols),
148  r1(rows, cols),
149  r2(rows, cols);
150 
151  MatrixType m1up = m1.template triangularView<Upper>();
152  MatrixType m2up = m2.template triangularView<Upper>();
153 
154  if (rows>1 && cols>1)
155  {
156  VERIFY(m1up.isUpperTriangular());
157  VERIFY(m2up.transpose().isLowerTriangular());
158  VERIFY(!m2.isLowerTriangular());
159  }
160 
161  // test overloaded operator+=
162  r1.setZero();
163  r2.setZero();
164  r1.template triangularView<Upper>() += m1;
165  r2 += m1up;
166  VERIFY_IS_APPROX(r1,r2);
167 
168  // test overloaded operator=
169  m1.setZero();
170  m1.template triangularView<Upper>() = 3 * m2;
171  m3 = 3 * m2;
172  VERIFY_IS_APPROX(m3.template triangularView<Upper>().toDenseMatrix(), m1);
173 
174 
175  m1.setZero();
176  m1.template triangularView<Lower>() = 3 * m2;
177  VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);
178 
179  m1.setZero();
180  m1.template triangularView<StrictlyUpper>() = 3 * m2;
181  VERIFY_IS_APPROX(m3.template triangularView<StrictlyUpper>().toDenseMatrix(), m1);
182 
183 
184  m1.setZero();
185  m1.template triangularView<StrictlyLower>() = 3 * m2;
186  VERIFY_IS_APPROX(m3.template triangularView<StrictlyLower>().toDenseMatrix(), m1);
187  m1.setRandom();
188  m2 = m1.template triangularView<Upper>();
189  VERIFY(m2.isUpperTriangular());
190  VERIFY(!m2.isLowerTriangular());
191  m2 = m1.template triangularView<StrictlyUpper>();
192  VERIFY(m2.isUpperTriangular());
193  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
194  m2 = m1.template triangularView<UnitUpper>();
195  VERIFY(m2.isUpperTriangular());
196  m2.diagonal().array() -= Scalar(1);
197  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
198  m2 = m1.template triangularView<Lower>();
199  VERIFY(m2.isLowerTriangular());
200  VERIFY(!m2.isUpperTriangular());
201  m2 = m1.template triangularView<StrictlyLower>();
202  VERIFY(m2.isLowerTriangular());
203  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
204  m2 = m1.template triangularView<UnitLower>();
205  VERIFY(m2.isLowerTriangular());
206  m2.diagonal().array() -= Scalar(1);
207  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
208  // test swap
209  m1.setOnes();
210  m2.setZero();
211  m2.template triangularView<Upper>().swap(m1);
212  m3.setZero();
213  m3.template triangularView<Upper>().setOnes();
214  VERIFY_IS_APPROX(m2,m3);
215 }
216 
217 void bug_159()
218 {
219  Matrix3d m = Matrix3d::Random().triangularView<Lower>();
220  EIGEN_UNUSED_VARIABLE(m)
221 }
222 
223 void test_triangular()
224 {
225  int maxsize = (std::min)(EIGEN_TEST_MAX_SIZE,20);
226  for(int i = 0; i < g_repeat ; i++)
227  {
228  int r = internal::random<int>(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(r)
229  int c = internal::random<int>(2,maxsize); TEST_SET_BUT_UNUSED_VARIABLE(c)
230 
231  CALL_SUBTEST_1( triangular_square(Matrix<float, 1, 1>()) );
232  CALL_SUBTEST_2( triangular_square(Matrix<float, 2, 2>()) );
233  CALL_SUBTEST_3( triangular_square(Matrix3d()) );
234  CALL_SUBTEST_4( triangular_square(Matrix<std::complex<float>,8, 8>()) );
235  CALL_SUBTEST_5( triangular_square(MatrixXcd(r,r)) );
236  CALL_SUBTEST_6( triangular_square(Matrix<float,Dynamic,Dynamic,RowMajor>(r, r)) );
237 
238  CALL_SUBTEST_7( triangular_rect(Matrix<float, 4, 5>()) );
239  CALL_SUBTEST_8( triangular_rect(Matrix<double, 6, 2>()) );
240  CALL_SUBTEST_9( triangular_rect(MatrixXcf(r, c)) );
241  CALL_SUBTEST_5( triangular_rect(MatrixXcd(r, c)) );
242  CALL_SUBTEST_6( triangular_rect(Matrix<float,Dynamic,Dynamic,RowMajor>(r, c)) );
243  }
244 
245  CALL_SUBTEST_1( bug_159() );
246 }
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
test_triangular
void test_triangular()
Definition: triangular.cpp:223
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:33
v2
Vector v2
Definition: testSerializationBase.cpp:36
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
min
#define min(a, b)
Definition: datatypes.h:19
Eigen::Transpose
Expression of the transpose of a matrix.
Definition: Transpose.h:52
m2
MatrixType m2(n_dims)
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
maxsize
const size_t maxsize
Definition: benchmark-blocking-sizes.cpp:49
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
swap
int EIGEN_BLAS_FUNC() swap(int *n, RealScalar *px, int *incx, RealScalar *py, int *incy)
Definition: level1_impl.h:152
Eigen::Lower
Definition: Constants.h:204
m1
Matrix3d m1
Definition: IOFormat.cpp:2
Eigen::g_repeat
static int g_repeat
Definition: main.h:144
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
setOnes
v setOnes(3)
e
Array< double, 1, 3 > e(1./3., 0.5, 2.)
r2
static const double r2
Definition: testSmartRangeFactor.cpp:32
triangular_square
void triangular_square(const MatrixType &m)
Definition: triangular.cpp:14
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
abs2
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Abs2ReturnType abs2() const
Definition: ArrayCwiseUnaryOps.h:71
Cols
static const int Cols
Definition: testCallRecord.cpp:32
TEST_SET_BUT_UNUSED_VARIABLE
#define TEST_SET_BUT_UNUSED_VARIABLE(X)
Definition: main.h:91
r1
static const double r1
Definition: testSmartRangeFactor.cpp:32
VERIFY
#define VERIFY(a)
Definition: main.h:325
EIGEN_TEST_MAX_SIZE
#define EIGEN_TEST_MAX_SIZE
Definition: boostmultiprec.cpp:16
triangularView< Lower >
A triangularView< Lower >().adjoint().solveInPlace(B)
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
bug_159
void bug_159()
Definition: triangular.cpp:217
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
triangular_rect
void triangular_rect(const MatrixType &m)
Definition: triangular.cpp:135
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::setRandom
Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > & setRandom(Index size)
EIGEN_UNUSED_VARIABLE
#define EIGEN_UNUSED_VARIABLE(var)
Definition: Macros.h:618

gtsam
Author(s):
autogenerated on Sat May 8 2021 02:51:10