aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
authorGravatar Benoit Jacob <jacob.benoit.1@gmail.com>2008-08-09 19:26:14 +0000
committerGravatar Benoit Jacob <jacob.benoit.1@gmail.com>2008-08-09 19:26:14 +0000
commitd6e88f81551d5d9c1b66f4fc13d2f2211cb689ff (patch)
treee29d73a7fe043c84818e2b8f8b5cc573bc7cf13d
parent4fa40367e9bf55ea8b2ad1040b3fb73f94e4f12f (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.h2
-rw-r--r--Eigen/src/LU/LU.h38
-rw-r--r--test/CMakeLists.txt1
-rw-r--r--test/determinant.cpp1
-rw-r--r--test/lu.cpp119
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>() );
+ }
+}