diff options
author | 2008-08-09 19:26:14 +0000 | |
---|---|---|
committer | 2008-08-09 19:26:14 +0000 | |
commit | d6e88f81551d5d9c1b66f4fc13d2f2211cb689ff (patch) | |
tree | e29d73a7fe043c84818e2b8f8b5cc573bc7cf13d | |
parent | 4fa40367e9bf55ea8b2ad1040b3fb73f94e4f12f (diff) |
* add LU unit-test. Seems like we have very good numerical stability!
* some cleanup, and grant me a copyright line on the determinant test.
-rw-r--r-- | Eigen/src/LU/Inverse.h | 2 | ||||
-rw-r--r-- | Eigen/src/LU/LU.h | 38 | ||||
-rw-r--r-- | test/CMakeLists.txt | 1 | ||||
-rw-r--r-- | test/determinant.cpp | 1 | ||||
-rw-r--r-- | test/lu.cpp | 119 |
5 files changed, 146 insertions, 15 deletions
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h index 6f6256d92..f2bf98f8a 100644 --- a/Eigen/src/LU/Inverse.h +++ b/Eigen/src/LU/Inverse.h @@ -161,7 +161,7 @@ struct ei_compute_inverse static inline void run(const MatrixType& matrix, MatrixType* result) { LU<MatrixType> lu(matrix); - lu.solve(MatrixType::Identity(matrix.rows(), matrix.cols()), result); + lu.computeInverse(result); } }; diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index e603d09dc..a7d5cd2e8 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -38,12 +38,9 @@ * are permutation matrices. * * This decomposition provides the generic approach to solving systems of linear equations, computing - * the rank, invertibility, inverse, and determinant. However for the case when invertibility is - * assumed, we have a specialized variant (see MatrixBase::inverse()) achieving better performance. + * the rank, invertibility, inverse, kernel, and determinant. * - * \sa MatrixBase::lu(), MatrixBase::determinant(), MatrixBase::rank(), MatrixBase::kernelDim(), - * MatrixBase::kernelBasis(), MatrixBase::solve(), MatrixBase::isInvertible(), - * MatrixBase::inverse(), MatrixBase::computeInverse() + * \sa MatrixBase::lu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse() */ template<typename MatrixType> class LU { @@ -141,6 +138,18 @@ template<typename MatrixType> class LU return isInjective() && isSurjective(); } + inline void computeInverse(MatrixType *result) const + { + solve(MatrixType::Identity(m_lu.rows(), m_lu.cols()), result); + } + + inline MatrixType inverse() const + { + MatrixType result; + computeInverse(&result); + return result; + } + protected: MatrixType m_lu; IntColVectorType m_p; @@ -163,7 +172,7 @@ LU<MatrixType>::LU(const MatrixType& matrix) IntRowVectorType cols_transpositions(matrix.cols()); int number_of_transpositions = 0; - RealScalar biggest; + RealScalar biggest = RealScalar(0); for(int k = 0; k < size; k++) { int row_of_biggest_in_corner, col_of_biggest_in_corner; @@ -224,7 +233,7 @@ void LU<MatrixType>::computeKernel(Matrix<typename MatrixType::Scalar, > *result) const { ei_assert(!isInvertible()); - const int dimker = dimensionOfKernel(), rows = m_lu.rows(), cols = m_lu.cols(); + const int dimker = dimensionOfKernel(), cols = m_lu.cols(); result->resize(cols, dimker); /* Let us use the following lemma: @@ -284,21 +293,22 @@ bool LU<MatrixType>::solve( * Step 4: result = Qd; */ - ei_assert(b.rows() == m_lu.rows()); - const int smalldim = std::min(m_lu.rows(), m_lu.cols()); + const int rows = m_lu.rows(); + ei_assert(b.rows() == rows); + const int smalldim = std::min(rows, m_lu.cols()); typename OtherDerived::Eval c(b.rows(), b.cols()); // Step 1 - for(int i = 0; i < m_lu.rows(); i++) c.row(m_p.coeff(i)) = b.row(i); + for(int i = 0; i < rows; i++) c.row(m_p.coeff(i)) = b.row(i); // Step 2 Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime, MatrixType::MaxRowsAtCompileTime, - MatrixType::MaxRowsAtCompileTime> l(m_lu.rows(), m_lu.rows()); + MatrixType::MaxRowsAtCompileTime> l(rows, rows); l.setZero(); - l.corner(Eigen::TopLeft,m_lu.rows(),smalldim) - = m_lu.corner(Eigen::TopLeft,m_lu.rows(),smalldim); + l.corner(Eigen::TopLeft,rows,smalldim) + = m_lu.corner(Eigen::TopLeft,rows,smalldim); l.template marked<UnitLower>().inverseProductInPlace(c); // Step 3 @@ -330,7 +340,7 @@ bool LU<MatrixType>::solve( * \sa class LU */ template<typename Derived> -const LU<typename MatrixBase<Derived>::EvalType> +inline const LU<typename MatrixBase<Derived>::EvalType> MatrixBase<Derived>::lu() const { return eval(); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index aa3860575..cb52ba150 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -99,6 +99,7 @@ EI_ADD_TEST(map) EI_ADD_TEST(array) EI_ADD_TEST(triangular) EI_ADD_TEST(cholesky) +EI_ADD_TEST(lu "-O2") EI_ADD_TEST(determinant) EI_ADD_TEST(inverse) EI_ADD_TEST(qr) diff --git a/test/determinant.cpp b/test/determinant.cpp index dc9f5eb54..ab9460b15 100644 --- a/test/determinant.cpp +++ b/test/determinant.cpp @@ -1,6 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. Eigen itself is part of the KDE project. // +// Copyright (C) 2008 Benoit Jacob <jacob@math.jussieu.fr> // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> // // Eigen is free software; you can redistribute it and/or diff --git a/test/lu.cpp b/test/lu.cpp new file mode 100644 index 000000000..5b0795d08 --- /dev/null +++ b/test/lu.cpp @@ -0,0 +1,119 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008 Benoit Jacob <jacob@math.jussieu.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" +#include <Eigen/LU> + +template<typename Derived> +void doSomeRankPreservingOperations(Eigen::MatrixBase<Derived>& m) +{ + for(int a = 0; a < 3*(m.rows()+m.cols()); a++) + { + double d = Eigen::ei_random<double>(-1,1); + int i = Eigen::ei_random<int>(0,m.rows()-1); // i is a random row number + int j; + do { + j = Eigen::ei_random<int>(0,m.rows()-1); + } while (i==j); // j is another one (must be different) + m.row(i) += d * m.row(j); + + i = Eigen::ei_random<int>(0,m.cols()-1); // i is a random column number + do { + j = Eigen::ei_random<int>(0,m.cols()-1); + } while (i==j); // j is another one (must be different) + m.col(i) += d * m.col(j); + } +} + +template<typename MatrixType> void lu_non_invertible() +{ + /* this test covers the following files: + LU.h + */ + int rows = ei_random<int>(10,200), cols = ei_random<int>(10,200), cols2 = ei_random<int>(10,200); + int rank = ei_random<int>(1, std::min(rows, cols)-1); + + MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1); + m1.setRandom(); + if(rows <= cols) + for(int i = rank; i < rows; i++) m1.row(i).setZero(); + else + for(int i = rank; i < cols; i++) m1.col(i).setZero(); + doSomeRankPreservingOperations(m1); + + LU<MatrixType> lu(m1); + VERIFY(cols - rank == lu.dimensionOfKernel()); + VERIFY(rank == lu.rank()); + VERIFY(!lu.isInjective()); + VERIFY(!lu.isInvertible()); + VERIFY(lu.isSurjective() == (lu.rank() == rows)); + VERIFY((m1 * lu.kernel()).isMuchSmallerThan(m1)); + lu.computeKernel(&k); + VERIFY((m1 * k).isMuchSmallerThan(m1)); + m2.setRandom(); + m3 = m1*m2; + m2.setRandom(); + lu.solve(m3, &m2); + VERIFY_IS_APPROX(m3, m1*m2); + m3.setRandom(); + VERIFY(!lu.solve(m3, &m2)); +} + +template<typename MatrixType> void lu_invertible() +{ + /* this test covers the following files: + LU.h + */ + int size = ei_random<int>(10,200); + + MatrixType m1(size, size), m2(size, size), m3(size, size); + m1.setRandom(); + + LU<MatrixType> lu(m1); + VERIFY(0 == lu.dimensionOfKernel()); + VERIFY(size == lu.rank()); + VERIFY(lu.isInjective()); + VERIFY(lu.isSurjective()); + VERIFY(lu.isInvertible()); + m3.setRandom(); + lu.solve(m3, &m2); + VERIFY_IS_APPROX(m3, m1*m2); + VERIFY_IS_APPROX(m2, lu.inverse()*m3); + m3.setRandom(); + VERIFY(lu.solve(m3, &m2)); +} + +void test_lu() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( lu_non_invertible<MatrixXf>() ); + CALL_SUBTEST( lu_non_invertible<MatrixXd>() ); + CALL_SUBTEST( lu_non_invertible<MatrixXcf>() ); + CALL_SUBTEST( lu_non_invertible<MatrixXcd>() ); + CALL_SUBTEST( lu_invertible<MatrixXf>() ); + CALL_SUBTEST( lu_invertible<MatrixXd>() ); + CALL_SUBTEST( lu_invertible<MatrixXcf>() ); + CALL_SUBTEST( lu_invertible<MatrixXcd>() ); + } +} |