diff options
author | Gael Guennebaud <g.gael@free.fr> | 2009-07-06 23:43:20 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2009-07-06 23:43:20 +0200 |
commit | 1aea45335ff60c6a5ed6dc06cd798d050eff661a (patch) | |
tree | d4dfb36ebe76b45f25dd984ac681ac9f92894515 /test/triangular.cpp | |
parent | 08e419dcb151da41f169304f751e5467cf0c7b4a (diff) |
* bybye Part, welcome TriangularView and SelfAdjointView.
* move solveTriangular*() to TriangularView::solve*()
* move .llt() to SelfAdjointView
* add a high level wrapper to the efficient selfadjoint * vector product
* improve LLT so that we can specify which triangular part is meaningless
=> there are still many things to do (doc, cleaning, improve the matrix products, etc.)
Diffstat (limited to 'test/triangular.cpp')
-rw-r--r-- | test/triangular.cpp | 57 |
1 files changed, 31 insertions, 26 deletions
diff --git a/test/triangular.cpp b/test/triangular.cpp index 850c3eee0..1829a2578 100644 --- a/test/triangular.cpp +++ b/test/triangular.cpp @@ -1,4 +1,4 @@ -// This file is part of Eigen, a lightweight C++ template library +// This file is triangularView of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@gmail.com> @@ -51,8 +51,8 @@ template<typename MatrixType> void triangular(const MatrixType& m) v2 = VectorType::Random(rows), vzero = VectorType::Zero(rows); - MatrixType m1up = m1.template part<Eigen::UpperTriangular>(); - MatrixType m2up = m2.template part<Eigen::UpperTriangular>(); + MatrixType m1up = m1.template triangularView<Eigen::UpperTriangular>(); + MatrixType m2up = m2.template triangularView<Eigen::UpperTriangular>(); if (rows*cols>1) { @@ -66,22 +66,22 @@ template<typename MatrixType> void triangular(const MatrixType& m) // test overloaded operator+= r1.setZero(); r2.setZero(); - r1.template part<Eigen::UpperTriangular>() += m1; + r1.template triangularView<Eigen::UpperTriangular>() += m1; r2 += m1up; VERIFY_IS_APPROX(r1,r2); // test overloaded operator= m1.setZero(); - m1.template part<Eigen::UpperTriangular>() = (m2.transpose() * m2).lazy(); + m1.template triangularView<Eigen::UpperTriangular>() = (m2.transpose() * m2).lazy(); m3 = m2.transpose() * m2; - VERIFY_IS_APPROX(m3.template part<Eigen::LowerTriangular>().transpose(), m1); + VERIFY_IS_APPROX(m3.template triangularView<Eigen::LowerTriangular>().transpose().toDense(), m1); // test overloaded operator= m1.setZero(); - m1.template part<Eigen::LowerTriangular>() = (m2.transpose() * m2).lazy(); - VERIFY_IS_APPROX(m3.template part<Eigen::LowerTriangular>(), m1); + m1.template triangularView<Eigen::LowerTriangular>() = (m2.transpose() * m2).lazy(); + VERIFY_IS_APPROX(m3.template triangularView<Eigen::LowerTriangular>().toDense(), m1); - VERIFY_IS_APPROX(m3.template part<DiagonalBits>(), m3.diagonal().asDiagonal()); + // VERIFY_IS_APPROX(m3.template triangularView<DiagonalBits>(), m3.diagonal().asDiagonal()); m1 = MatrixType::Random(rows, cols); for (int i=0; i<rows; ++i) @@ -89,37 +89,42 @@ template<typename MatrixType> void triangular(const MatrixType& m) Transpose<MatrixType> trm4(m4); // test back and forward subsitution - m3 = m1.template part<Eigen::LowerTriangular>(); - VERIFY(m3.template marked<Eigen::LowerTriangular>().solveTriangular(m3).cwise().abs().isIdentity(test_precision<RealScalar>())); - VERIFY(m3.transpose().template marked<Eigen::UpperTriangular>() - .solveTriangular(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>())); + m3 = m1.template triangularView<Eigen::LowerTriangular>(); + VERIFY(m3.template triangularView<Eigen::LowerTriangular>().solve(m3).cwise().abs().isIdentity(test_precision<RealScalar>())); + VERIFY(m3.transpose().template triangularView<Eigen::UpperTriangular>() + .solve(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>())); // check M * inv(L) using in place API m4 = m3; - m3.transpose().template marked<Eigen::UpperTriangular>().solveTriangularInPlace(trm4); + m3.transpose().template triangularView<Eigen::UpperTriangular>().solveInPlace(trm4); VERIFY(m4.cwise().abs().isIdentity(test_precision<RealScalar>())); - m3 = m1.template part<Eigen::UpperTriangular>(); - VERIFY(m3.template marked<Eigen::UpperTriangular>().solveTriangular(m3).cwise().abs().isIdentity(test_precision<RealScalar>())); - VERIFY(m3.transpose().template marked<Eigen::LowerTriangular>() - .solveTriangular(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>())); + m3 = m1.template triangularView<Eigen::UpperTriangular>(); + VERIFY(m3.template triangularView<Eigen::UpperTriangular>().solve(m3).cwise().abs().isIdentity(test_precision<RealScalar>())); + VERIFY(m3.transpose().template triangularView<Eigen::LowerTriangular>() + .solve(m3.transpose()).cwise().abs().isIdentity(test_precision<RealScalar>())); // check M * inv(U) using in place API m4 = m3; - m3.transpose().template marked<Eigen::LowerTriangular>().solveTriangularInPlace(trm4); + m3.transpose().template triangularView<Eigen::LowerTriangular>().solveInPlace(trm4); VERIFY(m4.cwise().abs().isIdentity(test_precision<RealScalar>())); - m3 = m1.template part<Eigen::UpperTriangular>(); - VERIFY(m2.isApprox(m3 * (m3.template marked<Eigen::UpperTriangular>().solveTriangular(m2)), largerEps)); - m3 = m1.template part<Eigen::LowerTriangular>(); - VERIFY(m2.isApprox(m3 * (m3.template marked<Eigen::LowerTriangular>().solveTriangular(m2)), largerEps)); + m3 = m1.template triangularView<Eigen::UpperTriangular>(); + VERIFY(m2.isApprox(m3 * (m3.template triangularView<Eigen::UpperTriangular>().solve(m2)), largerEps)); + m3 = m1.template triangularView<Eigen::LowerTriangular>(); + VERIFY(m2.isApprox(m3 * (m3.template triangularView<Eigen::LowerTriangular>().solve(m2)), largerEps)); - VERIFY((m1.template part<Eigen::UpperTriangular>() * m2.template part<Eigen::UpperTriangular>()).isUpperTriangular()); + // check solve with unit diagonal + m3 = m1.template triangularView<Eigen::UnitUpperTriangular>(); + VERIFY(m2.isApprox(m3 * (m1.template triangularView<Eigen::UnitUpperTriangular>().solve(m2)), largerEps)); + +// VERIFY(( m1.template triangularView<Eigen::UpperTriangular>() +// * m2.template triangularView<Eigen::UpperTriangular>()).isUpperTriangular()); // test swap m1.setOnes(); m2.setZero(); - m2.template part<Eigen::UpperTriangular>().swap(m1); + m2.template triangularView<Eigen::UpperTriangular>().swap(m1); m3.setZero(); - m3.template part<Eigen::UpperTriangular>().setOnes(); + m3.template triangularView<Eigen::UpperTriangular>().setOnes(); VERIFY_IS_APPROX(m2,m3); } |