Merged the latest version of the code from eigen/eigen

This commit is contained in:
Benoit Steiner
2014-02-18 18:51:24 -08:00
37 changed files with 504 additions and 223 deletions

View File

@@ -263,8 +263,8 @@ template<typename MatrixType> void cholesky_bug241(const MatrixType& m)
// LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal.
// This test checks that LDLT reports correctly that matrix is indefinite.
// See http://forum.kde.org/viewtopic.php?f=74&t=106942
template<typename MatrixType> void cholesky_indefinite(const MatrixType& m)
// See http://forum.kde.org/viewtopic.php?f=74&t=106942 and bug 736
template<typename MatrixType> void cholesky_definiteness(const MatrixType& m)
{
eigen_assert(m.rows() == 2 && m.cols() == 2);
MatrixType mat;
@@ -280,6 +280,24 @@ template<typename MatrixType> void cholesky_indefinite(const MatrixType& m)
VERIFY(!ldlt.isNegative());
VERIFY(!ldlt.isPositive());
}
{
mat << 0, 0, 0, 0;
LDLT<MatrixType> ldlt(mat);
VERIFY(ldlt.isNegative());
VERIFY(ldlt.isPositive());
}
{
mat << 0, 0, 0, 1;
LDLT<MatrixType> ldlt(mat);
VERIFY(!ldlt.isNegative());
VERIFY(ldlt.isPositive());
}
{
mat << -1, 0, 0, 0;
LDLT<MatrixType> ldlt(mat);
VERIFY(ldlt.isNegative());
VERIFY(!ldlt.isPositive());
}
}
template<typename MatrixType> void cholesky_verify_assert()
@@ -309,7 +327,7 @@ void test_cholesky()
CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) );
CALL_SUBTEST_3( cholesky(Matrix2d()) );
CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) );
CALL_SUBTEST_3( cholesky_indefinite(Matrix2d()) );
CALL_SUBTEST_3( cholesky_definiteness(Matrix2d()) );
CALL_SUBTEST_4( cholesky(Matrix3f()) );
CALL_SUBTEST_5( cholesky(Matrix4d()) );
s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE);

View File

@@ -279,6 +279,13 @@ template<typename Scalar, int Mode, int Options> void transformations()
t1 = Eigen::Scaling(s0,s0,s0) * t1;
VERIFY_IS_APPROX(t0.matrix(), t1.matrix());
t0 = t3;
t0.scale(s0);
t1 = t3 * Eigen::Scaling(s0);
VERIFY_IS_APPROX(t0.matrix(), t1.matrix());
t0.prescale(s0);
t1 = Eigen::Scaling(s0) * t1;
VERIFY_IS_APPROX(t0.matrix(), t1.matrix());
t0.setIdentity();
t0.prerotate(q1).prescale(v0).pretranslate(v0);

View File

@@ -54,6 +54,8 @@ template<typename MatrixType> void qr_invertible()
{
using std::log;
using std::abs;
using std::pow;
using std::max;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef typename MatrixType::Scalar Scalar;
@@ -65,7 +67,7 @@ template<typename MatrixType> void qr_invertible()
if (internal::is_same<RealScalar,float>::value)
{
// let's build a matrix more stable to inverse
MatrixType a = MatrixType::Random(size,size*2);
MatrixType a = MatrixType::Random(size,size*4);
m1 += a * a.adjoint();
}
@@ -81,8 +83,11 @@ template<typename MatrixType> void qr_invertible()
m3 = qr.householderQ(); // get a unitary
m1 = m3 * m1 * m3;
qr.compute(m1);
VERIFY_IS_APPROX(absdet, qr.absDeterminant());
VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant());
// This test is tricky if the determinant becomes too small.
// Since we generate random numbers with magnitude rrange [0,1], the average determinant is 0.5^size
VERIFY_IS_MUCH_SMALLER_THAN( abs(absdet-qr.absDeterminant()), (max)(RealScalar(pow(0.5,size)),(max)(abs(absdet),abs(qr.absDeterminant()))) );
}
template<typename MatrixType> void qr_verify_assert()

View File

@@ -154,16 +154,16 @@ initSparse(double density,
sparseMat.finalize();
}
template<typename Scalar> void
template<typename Scalar,int Options,typename Index> void
initSparse(double density,
Matrix<Scalar,Dynamic,1>& refVec,
SparseVector<Scalar>& sparseVec,
SparseVector<Scalar,Options,Index>& sparseVec,
std::vector<int>* zeroCoords = 0,
std::vector<int>* nonzeroCoords = 0)
{
sparseVec.reserve(int(refVec.size()*density));
sparseVec.setZero();
for(int i=0; i<refVec.size(); i++)
for(Index i=0; i<refVec.size(); i++)
{
Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0);
if (v!=Scalar(0))
@@ -178,10 +178,10 @@ initSparse(double density,
}
}
template<typename Scalar> void
template<typename Scalar,int Options,typename Index> void
initSparse(double density,
Matrix<Scalar,1,Dynamic>& refVec,
SparseVector<Scalar,RowMajor>& sparseVec,
SparseVector<Scalar,Options,Index>& sparseVec,
std::vector<int>* zeroCoords = 0,
std::vector<int>* nonzeroCoords = 0)
{

View File

@@ -270,6 +270,14 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.row(0).dot(refM2.row(0)));
else
VERIFY_IS_APPROX(m1.innerVector(0).dot(refM2.row(0)), refM1.col(0).dot(refM2.row(0)));
DenseVector rv = DenseVector::Random(m1.cols());
DenseVector cv = DenseVector::Random(m1.rows());
Index r = internal::random<Index>(0,m1.rows()-2);
Index c = internal::random<Index>(0,m1.cols()-1);
VERIFY_IS_APPROX(( m1.template block<1,Dynamic>(r,0,1,m1.cols()).dot(rv)) , refM1.row(r).dot(rv));
VERIFY_IS_APPROX(m1.row(r).dot(rv), refM1.row(r).dot(rv));
VERIFY_IS_APPROX(m1.col(c).dot(cv), refM1.col(c).dot(cv));
VERIFY_IS_APPROX(m1.conjugate(), refM1.conjugate());
VERIFY_IS_APPROX(m1.real(), refM1.real());

View File

@@ -13,8 +13,9 @@ template<typename SparseMatrixType, typename DenseMatrix, bool IsRowMajor=Sparse
template<typename SparseMatrixType, typename DenseMatrix> struct test_outer<SparseMatrixType,DenseMatrix,false> {
static void run(SparseMatrixType& m2, SparseMatrixType& m4, DenseMatrix& refMat2, DenseMatrix& refMat4) {
int c = internal::random(0,m2.cols()-1);
int c1 = internal::random(0,m2.cols()-1);
typedef typename SparseMatrixType::Index Index;
Index c = internal::random<Index>(0,m2.cols()-1);
Index c1 = internal::random<Index>(0,m2.cols()-1);
VERIFY_IS_APPROX(m4=m2.col(c)*refMat2.col(c1).transpose(), refMat4=refMat2.col(c)*refMat2.col(c1).transpose());
VERIFY_IS_APPROX(m4=refMat2.col(c1)*m2.col(c).transpose(), refMat4=refMat2.col(c1)*refMat2.col(c).transpose());
}
@@ -22,8 +23,9 @@ template<typename SparseMatrixType, typename DenseMatrix> struct test_outer<Spar
template<typename SparseMatrixType, typename DenseMatrix> struct test_outer<SparseMatrixType,DenseMatrix,true> {
static void run(SparseMatrixType& m2, SparseMatrixType& m4, DenseMatrix& refMat2, DenseMatrix& refMat4) {
int r = internal::random(0,m2.rows()-1);
int c1 = internal::random(0,m2.cols()-1);
typedef typename SparseMatrixType::Index Index;
Index r = internal::random<Index>(0,m2.rows()-1);
Index c1 = internal::random<Index>(0,m2.cols()-1);
VERIFY_IS_APPROX(m4=m2.row(r).transpose()*refMat2.col(c1).transpose(), refMat4=refMat2.row(r).transpose()*refMat2.col(c1).transpose());
VERIFY_IS_APPROX(m4=refMat2.col(c1)*m2.row(r), refMat4=refMat2.col(c1)*refMat2.row(r));
}
@@ -37,9 +39,9 @@ template<typename SparseMatrixType> void sparse_product()
{
typedef typename SparseMatrixType::Index Index;
Index n = 100;
const Index rows = internal::random<int>(1,n);
const Index cols = internal::random<int>(1,n);
const Index depth = internal::random<int>(1,n);
const Index rows = internal::random<Index>(1,n);
const Index cols = internal::random<Index>(1,n);
const Index depth = internal::random<Index>(1,n);
typedef typename SparseMatrixType::Scalar Scalar;
enum { Flags = SparseMatrixType::Flags };
@@ -244,6 +246,7 @@ void test_sparse_product()
CALL_SUBTEST_1( (sparse_product<SparseMatrix<double,RowMajor> >()) );
CALL_SUBTEST_2( (sparse_product<SparseMatrix<std::complex<double>, ColMajor > >()) );
CALL_SUBTEST_2( (sparse_product<SparseMatrix<std::complex<double>, RowMajor > >()) );
CALL_SUBTEST_3( (sparse_product<SparseMatrix<float,ColMajor,long int> >()) );
CALL_SUBTEST_4( (sparse_product_regression_test<SparseMatrix<double,RowMajor>, Matrix<double, Dynamic, Dynamic, RowMajor> >()) );
}
}

View File

@@ -9,14 +9,14 @@
#include "sparse.h"
template<typename Scalar> void sparse_vector(int rows, int cols)
template<typename Scalar,typename Index> void sparse_vector(int rows, int cols)
{
double densityMat = (std::max)(8./(rows*cols), 0.01);
double densityVec = (std::max)(8./float(rows), 0.1);
typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
typedef Matrix<Scalar,Dynamic,1> DenseVector;
typedef SparseVector<Scalar> SparseVectorType;
typedef SparseMatrix<Scalar> SparseMatrixType;
typedef SparseVector<Scalar,0,Index> SparseVectorType;
typedef SparseMatrix<Scalar,0,Index> SparseMatrixType;
Scalar eps = 1e-6;
SparseMatrixType m1(rows,rows);
@@ -101,9 +101,10 @@ template<typename Scalar> void sparse_vector(int rows, int cols)
void test_sparse_vector()
{
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_1( sparse_vector<double>(8, 8) );
CALL_SUBTEST_2( sparse_vector<std::complex<double> >(16, 16) );
CALL_SUBTEST_1( sparse_vector<double>(299, 535) );
CALL_SUBTEST_1(( sparse_vector<double,int>(8, 8) ));
CALL_SUBTEST_2(( sparse_vector<std::complex<double>, int>(16, 16) ));
CALL_SUBTEST_1(( sparse_vector<double,long int>(299, 535) ));
CALL_SUBTEST_1(( sparse_vector<double,short>(299, 535) ));
}
}

View File

@@ -55,8 +55,16 @@ template<typename MatrixType> void stable_norm(const MatrixType& m)
Index rows = m.rows();
Index cols = m.cols();
Scalar big = internal::random<Scalar>() * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4));
Scalar small = internal::random<Scalar>() * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4));
// get a non-zero random factor
Scalar factor = internal::random<Scalar>();
while(numext::abs2(factor)<RealScalar(1e-4))
factor = internal::random<Scalar>();
Scalar big = factor * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4));
factor = internal::random<Scalar>();
while(numext::abs2(factor)<RealScalar(1e-4))
factor = internal::random<Scalar>();
Scalar small = factor * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4));
MatrixType vzero = MatrixType::Zero(rows, cols),
vrand = MatrixType::Random(rows, cols),
@@ -91,7 +99,7 @@ template<typename MatrixType> void stable_norm(const MatrixType& m)
VERIFY_IS_APPROX(vsmall.blueNorm(), sqrt(size)*abs(small));
VERIFY_IS_APPROX(vsmall.hypotNorm(), sqrt(size)*abs(small));
// Test compilation of cwise() version
// Test compilation of cwise() version
VERIFY_IS_APPROX(vrand.colwise().stableNorm(), vrand.colwise().norm());
VERIFY_IS_APPROX(vrand.colwise().blueNorm(), vrand.colwise().norm());
VERIFY_IS_APPROX(vrand.colwise().hypotNorm(), vrand.colwise().norm());

View File

@@ -130,10 +130,11 @@ void run_fixed_size_test(int num_elements)
// MUST be positive because in any other case det(cR_t) may become negative for
// odd dimensions!
const Scalar c = abs(internal::random<Scalar>());
// Also if c is to small compared to t.norm(), problem is ill-posed (cf. Bug 744)
const Scalar c = internal::random<Scalar>(0.5, 2.0);
FixedMatrix R = randMatrixSpecialUnitary<Scalar>(dim);
FixedVector t = Scalar(50)*FixedVector::Random(dim,1);
FixedVector t = Scalar(32)*FixedVector::Random(dim,1);
HomMatrix cR_t = HomMatrix::Identity(dim+1,dim+1);
cR_t.block(0,0,dim,dim) = c*R;
@@ -149,9 +150,9 @@ void run_fixed_size_test(int num_elements)
HomMatrix cR_t_umeyama = umeyama(src_block, dst_block);
const Scalar error = ( cR_t_umeyama*src - dst ).array().square().sum();
const Scalar error = ( cR_t_umeyama*src - dst ).squaredNorm();
VERIFY(error < Scalar(10)*std::numeric_limits<Scalar>::epsilon());
VERIFY(error < Scalar(16)*std::numeric_limits<Scalar>::epsilon());
}
void test_umeyama()