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
14template<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 (internal::abs2(m1(i,i))<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 subsitution 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}
118
119
120template<typename MatrixType> void triangular_rect(const MatrixType& m)
121{
122  typedef const typename MatrixType::Index Index;
123  typedef typename MatrixType::Scalar Scalar;
124  typedef typename NumTraits<Scalar>::Real RealScalar;
125  enum { Rows =  MatrixType::RowsAtCompileTime, Cols =  MatrixType::ColsAtCompileTime };
126  typedef Matrix<Scalar, Rows, 1> VectorType;
127  typedef Matrix<Scalar, Rows, Rows> RMatrixType;
128
129
130  Index rows = m.rows();
131  Index cols = m.cols();
132
133  MatrixType m1 = MatrixType::Random(rows, cols),
134             m2 = MatrixType::Random(rows, cols),
135             m3(rows, cols),
136             m4(rows, cols),
137             r1(rows, cols),
138             r2(rows, cols);
139
140  MatrixType m1up = m1.template triangularView<Upper>();
141  MatrixType m2up = m2.template triangularView<Upper>();
142
143  if (rows>1 && cols>1)
144  {
145    VERIFY(m1up.isUpperTriangular());
146    VERIFY(m2up.transpose().isLowerTriangular());
147    VERIFY(!m2.isLowerTriangular());
148  }
149
150  // test overloaded operator+=
151  r1.setZero();
152  r2.setZero();
153  r1.template triangularView<Upper>() +=  m1;
154  r2 += m1up;
155  VERIFY_IS_APPROX(r1,r2);
156
157  // test overloaded operator=
158  m1.setZero();
159  m1.template triangularView<Upper>() = 3 * m2;
160  m3 = 3 * m2;
161  VERIFY_IS_APPROX(m3.template triangularView<Upper>().toDenseMatrix(), m1);
162
163
164  m1.setZero();
165  m1.template triangularView<Lower>() = 3 * m2;
166  VERIFY_IS_APPROX(m3.template triangularView<Lower>().toDenseMatrix(), m1);
167
168  m1.setZero();
169  m1.template triangularView<StrictlyUpper>() = 3 * m2;
170  VERIFY_IS_APPROX(m3.template triangularView<StrictlyUpper>().toDenseMatrix(), m1);
171
172
173  m1.setZero();
174  m1.template triangularView<StrictlyLower>() = 3 * m2;
175  VERIFY_IS_APPROX(m3.template triangularView<StrictlyLower>().toDenseMatrix(), m1);
176  m1.setRandom();
177  m2 = m1.template triangularView<Upper>();
178  VERIFY(m2.isUpperTriangular());
179  VERIFY(!m2.isLowerTriangular());
180  m2 = m1.template triangularView<StrictlyUpper>();
181  VERIFY(m2.isUpperTriangular());
182  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
183  m2 = m1.template triangularView<UnitUpper>();
184  VERIFY(m2.isUpperTriangular());
185  m2.diagonal().array() -= Scalar(1);
186  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
187  m2 = m1.template triangularView<Lower>();
188  VERIFY(m2.isLowerTriangular());
189  VERIFY(!m2.isUpperTriangular());
190  m2 = m1.template triangularView<StrictlyLower>();
191  VERIFY(m2.isLowerTriangular());
192  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
193  m2 = m1.template triangularView<UnitLower>();
194  VERIFY(m2.isLowerTriangular());
195  m2.diagonal().array() -= Scalar(1);
196  VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1)));
197  // test swap
198  m1.setOnes();
199  m2.setZero();
200  m2.template triangularView<Upper>().swap(m1);
201  m3.setZero();
202  m3.template triangularView<Upper>().setOnes();
203  VERIFY_IS_APPROX(m2,m3);
204}
205
206void bug_159()
207{
208  Matrix3d m = Matrix3d::Random().triangularView<Lower>();
209  EIGEN_UNUSED_VARIABLE(m)
210}
211
212void test_triangular()
213{
214  int maxsize = (std::min)(EIGEN_TEST_MAX_SIZE,20);
215  for(int i = 0; i < g_repeat ; i++)
216  {
217    int r = internal::random<int>(2,maxsize); EIGEN_UNUSED_VARIABLE(r);
218    int c = internal::random<int>(2,maxsize); EIGEN_UNUSED_VARIABLE(c);
219
220    CALL_SUBTEST_1( triangular_square(Matrix<float, 1, 1>()) );
221    CALL_SUBTEST_2( triangular_square(Matrix<float, 2, 2>()) );
222    CALL_SUBTEST_3( triangular_square(Matrix3d()) );
223    CALL_SUBTEST_4( triangular_square(Matrix<std::complex<float>,8, 8>()) );
224    CALL_SUBTEST_5( triangular_square(MatrixXcd(r,r)) );
225    CALL_SUBTEST_6( triangular_square(Matrix<float,Dynamic,Dynamic,RowMajor>(r, r)) );
226
227    CALL_SUBTEST_7( triangular_rect(Matrix<float, 4, 5>()) );
228    CALL_SUBTEST_8( triangular_rect(Matrix<double, 6, 2>()) );
229    CALL_SUBTEST_9( triangular_rect(MatrixXcf(r, c)) );
230    CALL_SUBTEST_5( triangular_rect(MatrixXcd(r, c)) );
231    CALL_SUBTEST_6( triangular_rect(Matrix<float,Dynamic,Dynamic,RowMajor>(r, c)) );
232  }
233
234  CALL_SUBTEST_1( bug_159() );
235}
236