xref: /aosp_15_r20/external/eigen/test/triangular.cpp (revision bf2c37156dfe67e5dfebd6d394bad8b2ab5804d4) 
1 // This file is triangularView of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Gael Guennebaud <[email protected]>
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 
triangular_deprecated(const MatrixType & m)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 
triangular_square(const MatrixType & m)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 
triangular_rect(const MatrixType & m)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 
bug_159()260 void bug_159()
261 {
262   Matrix3d m = Matrix3d::Random().triangularView<Lower>();
263   EIGEN_UNUSED_VARIABLE(m)
264 }
265 
EIGEN_DECLARE_TEST(triangular)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 }
293