aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Eigen/src/LU/FullPivLU.h7
-rw-r--r--test/inverse.cpp2
-rw-r--r--test/lu.cpp65
-rw-r--r--test/main.h12
-rw-r--r--test/qr_colpivoting.cpp4
-rw-r--r--test/qr_fullpivoting.cpp2
6 files changed, 42 insertions, 50 deletions
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index ec551645b..0a305d52b 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -422,8 +422,11 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
// when k==0, biggest_in_corner is the biggest coeff absolute value in the original matrix
if(k == 0) cutoff = biggest_in_corner * NumTraits<Scalar>::epsilon();
- // if the pivot (hence the corner) is exactly zero, terminate to avoid generating nan/inf values
- if(ei_abs(biggest_in_corner) < cutoff)
+ // if the pivot (hence the corner) is "zero", terminate to avoid generating nan/inf values.
+ // Notice that using an exact comparison (biggest_in_corner==0) here, as Golub-van Loan do in
+ // their pseudo-code, results in numerical instability! The cutoff here has been validated
+ // by running the unit test 'lu' with many repetitions.
+ if(biggest_in_corner < cutoff)
{
// before exiting, make sure to initialize the still uninitialized transpositions
// in a sane state without destroying what we already have.
diff --git a/test/inverse.cpp b/test/inverse.cpp
index 3f6138e0c..1e567ad14 100644
--- a/test/inverse.cpp
+++ b/test/inverse.cpp
@@ -42,7 +42,7 @@ template<typename MatrixType> void inverse(const MatrixType& m)
m2(rows, cols),
mzero = MatrixType::Zero(rows, cols),
identity = MatrixType::Identity(rows, rows);
- createRandomProjectionOfRank(rows,rows,rows,m1);
+ createRandomPIMatrixOfRank(rows,rows,rows,m1);
m2 = m1.inverse();
VERIFY_IS_APPROX(m1, m2.inverse() );
diff --git a/test/lu.cpp b/test/lu.cpp
index 02f6ec805..442202a33 100644
--- a/test/lu.cpp
+++ b/test/lu.cpp
@@ -28,9 +28,6 @@ using namespace std;
template<typename MatrixType> void lu_non_invertible()
{
- static int times_called = 0;
- times_called++;
-
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
/* this test covers the following files:
@@ -55,11 +52,16 @@ template<typename MatrixType> void lu_non_invertible()
cols2 = cols = MatrixType::ColsAtCompileTime;
}
+ enum {
+ RowsAtCompileTime = MatrixType::RowsAtCompileTime,
+ ColsAtCompileTime = MatrixType::ColsAtCompileTime
+ };
typedef typename ei_kernel_retval_base<FullPivLU<MatrixType> >::ReturnType KernelMatrixType;
typedef typename ei_image_retval_base<FullPivLU<MatrixType> >::ReturnType ImageMatrixType;
- typedef Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> DynamicMatrixType;
- typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime>
+ typedef Matrix<typename MatrixType::Scalar, ColsAtCompileTime, ColsAtCompileTime>
CMatrixType;
+ typedef Matrix<typename MatrixType::Scalar, RowsAtCompileTime, RowsAtCompileTime>
+ RMatrixType;
int rank = ei_random<int>(1, std::min(rows, cols)-1);
@@ -68,26 +70,21 @@ template<typename MatrixType> void lu_non_invertible()
MatrixType m1(rows, cols), m3(rows, cols2);
CMatrixType m2(cols, cols2);
- createRandomProjectionOfRank(rank, rows, cols, m1);
+ createRandomPIMatrixOfRank(rank, rows, cols, m1);
FullPivLU<MatrixType> lu;
- // The special value 0.01 below works well in tests. Keep in mind that we're only computing the rank of projections.
- // So it's not clear at all the epsilon should play any role there.
+ // The special value 0.01 below works well in tests. Keep in mind that we're only computing the rank
+ // of singular values are either 0 or 1.
+ // So it's not clear at all that the epsilon should play any role there.
lu.setThreshold(RealScalar(0.01));
lu.compute(m1);
- // FIXME need better way to construct trapezoid matrices. extend triangularView to support rectangular.
- DynamicMatrixType u(rows,cols);
- for(int i = 0; i < rows; i++)
- for(int j = 0; j < cols; j++)
- u(i,j) = i>j ? Scalar(0) : lu.matrixLU()(i,j);
- DynamicMatrixType l(rows,rows);
- for(int i = 0; i < rows; i++)
- for(int j = 0; j < rows; j++)
- l(i,j) = (i<j || j>=cols) ? Scalar(0)
- : i==j ? Scalar(1)
- : lu.matrixLU()(i,j);
+ MatrixType u(rows,cols);
+ u = lu.matrixLU().template triangularView<Upper>();
+ RMatrixType l = RMatrixType::Identity(rows,rows);
+ l.block(0,0,rows,std::min(rows,cols)).template triangularView<StrictlyLower>()
+ = lu.matrixLU().block(0,0,rows,std::min(rows,cols));
VERIFY_IS_APPROX(lu.permutationP() * m1 * lu.permutationQ(), l*u);
@@ -101,20 +98,8 @@ template<typename MatrixType> void lu_non_invertible()
VERIFY(!lu.isSurjective());
VERIFY((m1 * m1kernel).isMuchSmallerThan(m1));
VERIFY(m1image.fullPivLu().rank() == rank);
+ VERIFY_IS_APPROX(m1 * m1.adjoint() * m1image, m1image);
- // The following test is damn hard to get to succeed over a large number of repetitions.
- // We're checking that the image is indeed the image, i.e. adding it as new columns doesn't increase the rank.
- // Since we've already tested rank() above, the point here is not to test rank(), it is to test image().
- // Since image() is implemented in a very simple way that doesn't leave much room for choice, the occasional
- // errors that we get here (one in 1e+4 repetitions roughly) are probably just a sign that it's a really
- // hard test, so we just limit how many times it's run.
- if(times_called < 100)
- {
- DynamicMatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols());
- sidebyside << m1, m1image;
- VERIFY(sidebyside.fullPivLu().rank() == rank);
- }
-
m2 = CMatrixType::Random(cols,cols2);
m3 = m1*m2;
m2 = CMatrixType::Random(cols,cols2);
@@ -128,20 +113,18 @@ template<typename MatrixType> void lu_invertible()
/* this test covers the following files:
LU.h
*/
+ typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
int size = ei_random<int>(1,200);
MatrixType m1(size, size), m2(size, size), m3(size, size);
- m1 = MatrixType::Random(size,size);
-
- if (ei_is_same_type<RealScalar,float>::ret)
- {
- // let's build a matrix more stable to inverse
- MatrixType a = MatrixType::Random(size,size*2);
- m1 += a * a.adjoint();
- }
+ FullPivLU<MatrixType> lu;
+ lu.setThreshold(0.01);
+ do {
+ m1 = MatrixType::Random(size,size);
+ lu.compute(m1);
+ } while(!lu.isInvertible());
- FullPivLU<MatrixType> lu(m1);
VERIFY(0 == lu.dimensionOfKernel());
VERIFY(lu.kernel().cols() == 1); // the kernel() should consist of a single (zero) column vector
VERIFY(size == lu.rank());
diff --git a/test/main.h b/test/main.h
index 6d296b2e3..96324de33 100644
--- a/test/main.h
+++ b/test/main.h
@@ -148,7 +148,7 @@ namespace Eigen
#define EIGEN_INTERNAL_DEBUGGING
#define EIGEN_NICE_RANDOM
-#include <Eigen/QR> // required for createRandomProjectionOfRank
+#include <Eigen/QR> // required for createRandomPIMatrixOfRank
#define VERIFY(a) do { if (!(a)) { \
@@ -342,8 +342,13 @@ inline bool test_isUnitary(const MatrixBase<Derived>& m)
return m.isUnitary(test_precision<typename ei_traits<Derived>::Scalar>());
}
+/** Creates a random Partial Isometry matrix of given rank.
+ *
+ * A partial isometry is a matrix all of whose singular values are either 0 or 1.
+ * This is very useful to test rank-revealing algorithms.
+ */
template<typename MatrixType>
-void createRandomProjectionOfRank(int desired_rank, int rows, int cols, MatrixType& m)
+void createRandomPIMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m)
{
typedef typename ei_traits<MatrixType>::Scalar Scalar;
enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
@@ -360,7 +365,8 @@ void createRandomProjectionOfRank(int desired_rank, int rows, int cols, MatrixTy
if(desired_rank == 1)
{
- m = VectorType::Random(rows) * VectorType::Random(cols).transpose();
+ // here we normalize the vectors to get a partial isometry
+ m = VectorType::Random(rows).normalized() * VectorType::Random(cols).normalized().transpose();
return;
}
diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp
index abee32184..96cc66316 100644
--- a/test/qr_colpivoting.cpp
+++ b/test/qr_colpivoting.cpp
@@ -36,7 +36,7 @@ template<typename MatrixType> void qr()
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
MatrixType m1;
- createRandomProjectionOfRank(rank,rows,cols,m1);
+ createRandomPIMatrixOfRank(rank,rows,cols,m1);
ColPivHouseholderQR<MatrixType> qr(m1);
VERIFY_IS_APPROX(rank, qr.rank());
VERIFY(cols - qr.rank() == qr.dimensionOfKernel());
@@ -64,7 +64,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
typedef typename MatrixType::Scalar Scalar;
int rank = ei_random<int>(1, std::min(int(Rows), int(Cols))-1);
Matrix<Scalar,Rows,Cols> m1;
- createRandomProjectionOfRank(rank,Rows,Cols,m1);
+ createRandomPIMatrixOfRank(rank,Rows,Cols,m1);
ColPivHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);
VERIFY_IS_APPROX(rank, qr.rank());
VERIFY(Cols - qr.rank() == qr.dimensionOfKernel());
diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp
index 60255f94c..7ad3af1fe 100644
--- a/test/qr_fullpivoting.cpp
+++ b/test/qr_fullpivoting.cpp
@@ -35,7 +35,7 @@ template<typename MatrixType> void qr()
typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
MatrixType m1;
- createRandomProjectionOfRank(rank,rows,cols,m1);
+ createRandomPIMatrixOfRank(rank,rows,cols,m1);
FullPivHouseholderQR<MatrixType> qr(m1);
VERIFY_IS_APPROX(rank, qr.rank());
VERIFY(cols - qr.rank() == qr.dimensionOfKernel());
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-04 16:55:52 +0000