From a20a744adc5414f9a79ae0637de9c3227c3dd420 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Thu, 19 Nov 2009 17:07:55 -0500 Subject: TriangularMatrix: extend to rectangular matrices --- test/triangular.cpp | 183 ++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 143 insertions(+), 40 deletions(-) (limited to 'test/triangular.cpp') diff --git a/test/triangular.cpp b/test/triangular.cpp index ee02c0022..3c89a8abb 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -24,7 +24,9 @@ #include "main.h" -template void triangular(const MatrixType& m) + + +template void triangular_square(const MatrixType& m) { typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; @@ -51,8 +53,8 @@ template void triangular(const MatrixType& m) v2 = VectorType::Random(rows), vzero = VectorType::Zero(rows); - MatrixType m1up = m1.template triangularView(); - MatrixType m2up = m2.template triangularView(); + MatrixType m1up = m1.template triangularView(); + MatrixType m2up = m2.template triangularView(); if (rows*cols>1) { @@ -66,20 +68,20 @@ template void triangular(const MatrixType& m) // test overloaded operator+= r1.setZero(); r2.setZero(); - r1.template triangularView() += m1; + r1.template triangularView() += m1; r2 += m1up; VERIFY_IS_APPROX(r1,r2); // test overloaded operator= m1.setZero(); - m1.template triangularView() = m2.transpose() + m2; + m1.template triangularView() = m2.transpose() + m2; m3 = m2.transpose() + m2; - VERIFY_IS_APPROX(m3.template triangularView().transpose().toDense(), m1); + VERIFY_IS_APPROX(m3.template triangularView().transpose().toDenseMatrix(), m1); // test overloaded operator= m1.setZero(); - m1.template triangularView() = m2.transpose() + m2; - VERIFY_IS_APPROX(m3.template triangularView().toDense(), m1); + m1.template triangularView() = m2.transpose() + m2; + VERIFY_IS_APPROX(m3.template triangularView().toDenseMatrix(), m1); m1 = MatrixType::Random(rows, cols); for (int i=0; i void triangular(const MatrixType& m) Transpose trm4(m4); // test back and forward subsitution with a vector as the rhs - m3 = m1.template triangularView(); - VERIFY(v2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView().solve(v2)), largerEps)); - m3 = m1.template triangularView(); - VERIFY(v2.isApprox(m3.transpose() * (m1.transpose().template triangularView().solve(v2)), largerEps)); - m3 = m1.template triangularView(); - VERIFY(v2.isApprox(m3 * (m1.template triangularView().solve(v2)), largerEps)); - m3 = m1.template triangularView(); - VERIFY(v2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView().solve(v2)), largerEps)); + m3 = m1.template triangularView(); + VERIFY(v2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView().solve(v2)), largerEps)); + m3 = m1.template triangularView(); + VERIFY(v2.isApprox(m3.transpose() * (m1.transpose().template triangularView().solve(v2)), largerEps)); + m3 = m1.template triangularView(); + VERIFY(v2.isApprox(m3 * (m1.template triangularView().solve(v2)), largerEps)); + m3 = m1.template triangularView(); + VERIFY(v2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView().solve(v2)), largerEps)); // test back and forward subsitution with a matrix as the rhs - m3 = m1.template triangularView(); - VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView().solve(m2)), largerEps)); - m3 = m1.template triangularView(); - VERIFY(m2.isApprox(m3.transpose() * (m1.transpose().template triangularView().solve(m2)), largerEps)); - m3 = m1.template triangularView(); - VERIFY(m2.isApprox(m3 * (m1.template triangularView().solve(m2)), largerEps)); - m3 = m1.template triangularView(); - VERIFY(m2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView().solve(m2)), largerEps)); + m3 = m1.template triangularView(); + VERIFY(m2.isApprox(m3.adjoint() * (m1.adjoint().template triangularView().solve(m2)), largerEps)); + m3 = m1.template triangularView(); + VERIFY(m2.isApprox(m3.transpose() * (m1.transpose().template triangularView().solve(m2)), largerEps)); + m3 = m1.template triangularView(); + VERIFY(m2.isApprox(m3 * (m1.template triangularView().solve(m2)), largerEps)); + m3 = m1.template triangularView(); + VERIFY(m2.isApprox(m3.conjugate() * (m1.conjugate().template triangularView().solve(m2)), largerEps)); // check M * inv(L) using in place API m4 = m3; - m3.transpose().template triangularView().solveInPlace(trm4); + m3.transpose().template triangularView().solveInPlace(trm4); VERIFY(m4.cwise().abs().isIdentity(test_precision())); // check M * inv(U) using in place API - m3 = m1.template triangularView(); + m3 = m1.template triangularView(); m4 = m3; - m3.transpose().template triangularView().solveInPlace(trm4); + m3.transpose().template triangularView().solveInPlace(trm4); VERIFY(m4.cwise().abs().isIdentity(test_precision())); // check solve with unit diagonal - m3 = m1.template triangularView(); - VERIFY(m2.isApprox(m3 * (m1.template triangularView().solve(m2)), largerEps)); + m3 = m1.template triangularView(); + VERIFY(m2.isApprox(m3 * (m1.template triangularView().solve(m2)), largerEps)); + +// VERIFY(( m1.template triangularView() +// * m2.template triangularView()).isUpperTriangular()); + + // test swap + m1.setOnes(); + m2.setZero(); + m2.template triangularView().swap(m1); + m3.setZero(); + m3.template triangularView().setOnes(); + VERIFY_IS_APPROX(m2,m3); + +} + + +template void triangular_rect(const MatrixType& m) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix VectorType; + + int rows = m.rows(); + int cols = m.cols(); + + MatrixType m1 = MatrixType::Random(rows, cols), + m2 = MatrixType::Random(rows, cols), + m3(rows, cols), + m4(rows, cols), + r1(rows, cols), + r2(rows, cols), + mzero = MatrixType::Zero(rows, cols), + mones = MatrixType::Ones(rows, cols), + identity = Matrix + ::Identity(rows, rows), + square = Matrix + ::Random(rows, rows); + VectorType v1 = VectorType::Random(rows), + v2 = VectorType::Random(rows), + vzero = VectorType::Zero(rows); + + MatrixType m1up = m1.template triangularView(); + MatrixType m2up = m2.template triangularView(); + + if (rows*cols>1) + { + VERIFY(m1up.isUpperTriangular()); + VERIFY(m2up.transpose().isLowerTriangular()); + VERIFY(!m2.isLowerTriangular()); + } + +// VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2); -// VERIFY(( m1.template triangularView() -// * m2.template triangularView()).isUpperTriangular()); + // test overloaded operator+= + r1.setZero(); + r2.setZero(); + r1.template triangularView() += m1; + r2 += m1up; + VERIFY_IS_APPROX(r1,r2); + // test overloaded operator= + m1.setZero(); + m1.template triangularView() = 3 * m2; + m3 = 3 * m2; + VERIFY_IS_APPROX(m3.template triangularView().toDenseMatrix(), m1); + + m1.setZero(); + m1.template triangularView() = 3 * m2; + VERIFY_IS_APPROX(m3.template triangularView().toDenseMatrix(), m1); + + m1.setZero(); + m1.template triangularView() = 3 * m2; + VERIFY_IS_APPROX(m3.template triangularView().toDenseMatrix(), m1); + + m1.setZero(); + m1.template triangularView() = 3 * m2; + VERIFY_IS_APPROX(m3.template triangularView().toDenseMatrix(), m1); + + m1.setRandom(); + m2 = m1.template triangularView(); + VERIFY(m2.isUpperTriangular()); + VERIFY(!m2.isLowerTriangular()); + m2 = m1.template triangularView(); + VERIFY(m2.isUpperTriangular()); + VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); + m2 = m1.template triangularView(); + VERIFY(m2.isUpperTriangular()); + m2.diagonal().cwise() -= Scalar(1); + VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); + m2 = m1.template triangularView(); + VERIFY(m2.isLowerTriangular()); + VERIFY(!m2.isUpperTriangular()); + m2 = m1.template triangularView(); + VERIFY(m2.isLowerTriangular()); + VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); + m2 = m1.template triangularView(); + VERIFY(m2.isLowerTriangular()); + m2.diagonal().cwise() -= Scalar(1); + VERIFY(m2.diagonal().isMuchSmallerThan(RealScalar(1))); // test swap m1.setOnes(); m2.setZero(); - m2.template triangularView().swap(m1); + m2.template triangularView().swap(m1); m3.setZero(); - m3.template triangularView().setOnes(); + m3.template triangularView().setOnes(); VERIFY_IS_APPROX(m2,m3); } @@ -137,12 +233,19 @@ template void triangular(const MatrixType& m) void test_triangular() { for(int i = 0; i < g_repeat ; i++) { - CALL_SUBTEST_1( triangular(Matrix()) ); - CALL_SUBTEST_2( triangular(Matrix()) ); - CALL_SUBTEST_3( triangular(Matrix3d()) ); - CALL_SUBTEST_4( triangular(MatrixXcf(4, 4)) ); - CALL_SUBTEST_5( triangular(Matrix,8, 8>()) ); - CALL_SUBTEST_6( triangular(MatrixXcd(17,17)) ); - CALL_SUBTEST_7( triangular(Matrix(5, 5)) ); + CALL_SUBTEST_1( triangular_square(Matrix()) ); + CALL_SUBTEST_2( triangular_square(Matrix()) ); + CALL_SUBTEST_3( triangular_square(Matrix3d()) ); + CALL_SUBTEST_4( triangular_square(MatrixXcf(4, 4)) ); + CALL_SUBTEST_5( triangular_square(Matrix,8, 8>()) ); + CALL_SUBTEST_6( triangular_square(MatrixXcd(17,17)) ); + CALL_SUBTEST_7( triangular_square(Matrix(5, 5)) ); + + CALL_SUBTEST_8( triangular_rect(Matrix()) ); + CALL_SUBTEST_9( triangular_rect(Matrix()) ); + CALL_SUBTEST_4( triangular_rect(MatrixXcf(4, 10)) ); + CALL_SUBTEST_6( triangular_rect(MatrixXcd(11, 3)) ); + CALL_SUBTEST_7( triangular_rect(Matrix(7, 6)) ); + } } -- cgit v1.2.3 From 88b551e89bd34367c12f0ba540ee1134c6897ba1 Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Thu, 19 Nov 2009 19:20:19 -0500 Subject: * fix compilation of unit-tests (sorry, had tested only 1 channel) * remove buggy (superfluous?) specialization in the meta-unroller --- Eigen/src/Core/TriangularMatrix.h | 15 --------------- test/triangular.cpp | 28 +++++++++++++++++----------- 2 files changed, 17 insertions(+), 26 deletions(-) (limited to 'test/triangular.cpp') diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index ed26683ea..34b029421 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -343,21 +343,6 @@ struct ei_triangular_assignment_selector } }; -template -struct ei_triangular_assignment_selector -{ - inline static void run(Derived1 &dst, const Derived2 &src) - { - if(Mode&UnitDiagBit) - { - if(ClearOpposite) - dst.coeffRef(0, 0) = 1; - } - else if(!(Mode & ZeroDiagBit)) - dst.copyCoeff(0, 0, src); - } -}; - // prevent buggy user code from causing an infinite recursion template struct ei_triangular_assignment_selector diff --git a/test/triangular.cpp b/test/triangular.cpp index 3c89a8abb..a10b457a7 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -141,7 +141,10 @@ template void triangular_rect(const MatrixType& m) { typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef Matrix VectorType; + enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; + typedef Matrix VectorType; + typedef Matrix RMatrixType; + int rows = m.rows(); int cols = m.cols(); @@ -153,10 +156,10 @@ template void triangular_rect(const MatrixType& m) r1(rows, cols), r2(rows, cols), mzero = MatrixType::Zero(rows, cols), - mones = MatrixType::Ones(rows, cols), - identity = Matrix + mones = MatrixType::Ones(rows, cols); + RMatrixType identity = Matrix ::Identity(rows, rows), - square = Matrix + square = Matrix ::Random(rows, rows); VectorType v1 = VectorType::Random(rows), v2 = VectorType::Random(rows), @@ -172,8 +175,6 @@ template void triangular_rect(const MatrixType& m) VERIFY(!m2.isLowerTriangular()); } -// VERIFY_IS_APPROX(m1up.transpose() * m2, m1.upper().transpose().lower() * m2); - // test overloaded operator+= r1.setZero(); r2.setZero(); @@ -227,25 +228,30 @@ template void triangular_rect(const MatrixType& m) m3.setZero(); m3.template triangularView().setOnes(); VERIFY_IS_APPROX(m2,m3); - } void test_triangular() { - for(int i = 0; i < g_repeat ; i++) { + for(int i = 0; i < g_repeat ; i++) + { + +#ifdef EIGEN_TEST_PART_7 + int r = ei_random(2,20); + int c = ei_random(2,20); +#endif + CALL_SUBTEST_1( triangular_square(Matrix()) ); CALL_SUBTEST_2( triangular_square(Matrix()) ); CALL_SUBTEST_3( triangular_square(Matrix3d()) ); CALL_SUBTEST_4( triangular_square(MatrixXcf(4, 4)) ); CALL_SUBTEST_5( triangular_square(Matrix,8, 8>()) ); CALL_SUBTEST_6( triangular_square(MatrixXcd(17,17)) ); - CALL_SUBTEST_7( triangular_square(Matrix(5, 5)) ); + CALL_SUBTEST_7( triangular_square(Matrix(r, r)) ); CALL_SUBTEST_8( triangular_rect(Matrix()) ); CALL_SUBTEST_9( triangular_rect(Matrix()) ); CALL_SUBTEST_4( triangular_rect(MatrixXcf(4, 10)) ); CALL_SUBTEST_6( triangular_rect(MatrixXcd(11, 3)) ); - CALL_SUBTEST_7( triangular_rect(Matrix(7, 6)) ); - + CALL_SUBTEST_7( triangular_rect(Matrix(r, c)) ); } } -- cgit v1.2.3 From eac3232095baf7c69e8ff7a6b636ffcf03db156b Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Thu, 19 Nov 2009 20:50:50 -0500 Subject: minor improvements in triangular stuff --- Eigen/src/Core/TriangularMatrix.h | 2 +- test/triangular.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) (limited to 'test/triangular.cpp') diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 34b029421..1286bf6a6 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -238,7 +238,7 @@ template class TriangularView DenseMatrixType toDenseMatrix() const { DenseMatrixType res(rows(), cols()); - res = *this; + evalToLazy(res); return res; } diff --git a/test/triangular.cpp b/test/triangular.cpp index a10b457a7..379f586c4 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -85,7 +85,7 @@ template void triangular_square(const MatrixType& m) m1 = MatrixType::Random(rows, cols); for (int i=0; i(); + while (ei_abs2(m1(i,i))<1e-1) m1(i,i) = ei_random(); Transpose trm4(m4); // test back and forward subsitution with a vector as the rhs -- cgit v1.2.3 From 6cbf662f14d8fd0606d9e9c6f7052e5dcc52b2ae Mon Sep 17 00:00:00 2001 From: Benoit Jacob Date: Thu, 19 Nov 2009 22:01:13 -0500 Subject: * don't laugh, but these bugs took me forever to fix. * expand unit tests to make sure to catch them: they nearly escaped the existing tests as these memory violations were highly dependent on the numbers of rows and cols. --- Eigen/src/Core/TriangularMatrix.h | 6 +++--- test/triangular.cpp | 29 +++++++++++++---------------- 2 files changed, 16 insertions(+), 19 deletions(-) (limited to 'test/triangular.cpp') diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 1286bf6a6..1a8f6493d 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -410,7 +410,7 @@ struct ei_triangular_assignment_selector @@ -451,9 +451,9 @@ struct ei_triangular_assignment_selector void triangular_rect(const MatrixType& m) m3 = 3 * m2; VERIFY_IS_APPROX(m3.template triangularView().toDenseMatrix(), m1); + m1.setZero(); m1.template triangularView() = 3 * m2; VERIFY_IS_APPROX(m3.template triangularView().toDenseMatrix(), m1); @@ -196,10 +197,10 @@ template void triangular_rect(const MatrixType& m) m1.template triangularView() = 3 * m2; VERIFY_IS_APPROX(m3.template triangularView().toDenseMatrix(), m1); + m1.setZero(); m1.template triangularView() = 3 * m2; VERIFY_IS_APPROX(m3.template triangularView().toDenseMatrix(), m1); - m1.setRandom(); m2 = m1.template triangularView(); VERIFY(m2.isUpperTriangular()); @@ -234,24 +235,20 @@ void test_triangular() { for(int i = 0; i < g_repeat ; i++) { - -#ifdef EIGEN_TEST_PART_7 - int r = ei_random(2,20); - int c = ei_random(2,20); -#endif + EIGEN_UNUSED int r = ei_random(2,20); + EIGEN_UNUSED int c = ei_random(2,20); CALL_SUBTEST_1( triangular_square(Matrix()) ); CALL_SUBTEST_2( triangular_square(Matrix()) ); CALL_SUBTEST_3( triangular_square(Matrix3d()) ); - CALL_SUBTEST_4( triangular_square(MatrixXcf(4, 4)) ); - CALL_SUBTEST_5( triangular_square(Matrix,8, 8>()) ); - CALL_SUBTEST_6( triangular_square(MatrixXcd(17,17)) ); - CALL_SUBTEST_7( triangular_square(Matrix(r, r)) ); - - CALL_SUBTEST_8( triangular_rect(Matrix()) ); - CALL_SUBTEST_9( triangular_rect(Matrix()) ); - CALL_SUBTEST_4( triangular_rect(MatrixXcf(4, 10)) ); - CALL_SUBTEST_6( triangular_rect(MatrixXcd(11, 3)) ); - CALL_SUBTEST_7( triangular_rect(Matrix(r, c)) ); + CALL_SUBTEST_4( triangular_square(Matrix,8, 8>()) ); + CALL_SUBTEST_5( triangular_square(MatrixXcd(r,r)) ); + CALL_SUBTEST_6( triangular_square(Matrix(r, r)) ); + + CALL_SUBTEST_7( triangular_rect(Matrix()) ); + CALL_SUBTEST_8( triangular_rect(Matrix()) ); + CALL_SUBTEST_9( triangular_rect(MatrixXcf(r, c)) ); + CALL_SUBTEST_5( triangular_rect(MatrixXcd(r, c)) ); + CALL_SUBTEST_6( triangular_rect(Matrix(r, c)) ); } } -- cgit v1.2.3