devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

* 4x4 inverse: revert to cofactors method

Browse Source 
* inverse tests: use createRandomMatrixOfRank, use more strict precision
* tests: createRandomMatrixOfRank: support 1x1 matrices
* determinant: nest the xpr
* Minor: add comment
This commit is contained in:
Benoit Jacob
2009-12-09 12:43:25 -05:00
parent f0315295e9
commit d2e44f2636
6 changed files with 43 additions and 121 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
test/inverse.cpp
View File

@@ -38,18 +38,11 @@ template<typename MatrixType> void inverse(const MatrixType& m)
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
MatrixType m1 = MatrixType::Random(rows, cols),
MatrixType m1(rows, cols),
m2(rows, cols),
mzero = MatrixType::Zero(rows, cols),
identity = MatrixType::Identity(rows, rows);
if (ei_is_same_type<RealScalar,float>::ret)
{
// let's build a more stable to inverse matrix
MatrixType a = MatrixType::Random(rows,cols);
m1 += m1 * m1.adjoint() + a * a.adjoint();
}
createRandomMatrixOfRank(rows,rows,rows,m1);
m2 = m1.inverse();
VERIFY_IS_APPROX(m1, m2.inverse() );

15
test/main.h
View File

@@ -353,13 +353,26 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType&
typedef Matrix<Scalar, Rows, Rows> MatrixAType;
typedef Matrix<Scalar, Cols, Cols> MatrixBType;
if(desired_rank == 0)
{
m.setZero(rows,cols);
return;
}
if(desired_rank == 1)
{
m = VectorType::Random(rows) * VectorType::Random(cols).transpose();
return;
}
MatrixAType a = MatrixAType::Random(rows,rows);
MatrixType d = MatrixType::Identity(rows,cols);
MatrixBType b = MatrixBType::Random(cols,cols);
// set the diagonal such that only desired_rank non-zero entries reamain
const int diag_size = std::min(d.rows(),d.cols());
d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank);
if(diag_size != desired_rank)
d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank);
HouseholderQR<MatrixAType> qra(a);
HouseholderQR<MatrixBType> qrb(b);

28
test/prec_inverse_4x4.cpp
View File

@@ -26,28 +26,6 @@
#include <Eigen/LU>
#include <algorithm>
Matrix4f inverse(const Matrix4f& m)
{
Matrix4f r;
r(0,0) = m.minor(0,0).determinant();
r(1,0) = -m.minor(0,1).determinant();
r(2,0) = m.minor(0,2).determinant();
r(3,0) = -m.minor(0,3).determinant();
r(0,2) = m.minor(2,0).determinant();
r(1,2) = -m.minor(2,1).determinant();
r(2,2) = m.minor(2,2).determinant();
r(3,2) = -m.minor(2,3).determinant();
r(0,1) = -m.minor(1,0).determinant();
r(1,1) = m.minor(1,1).determinant();
r(2,1) = -m.minor(1,2).determinant();
r(3,1) = m.minor(1,3).determinant();
r(0,3) = -m.minor(3,0).determinant();
r(1,3) = m.minor(3,1).determinant();
r(2,3) = -m.minor(3,2).determinant();
r(3,3) = m.minor(3,3).determinant();
return r / (m(0,0)*r(0,0) + m(1,0)*r(0,1) + m(2,0)*r(0,2) + m(3,0)*r(0,3));
}
template<typename MatrixType> void inverse_permutation_4x4()
{
typedef typename MatrixType::Scalar Scalar;
@@ -79,7 +57,7 @@ template<typename MatrixType> void inverse_general_4x4(int repeat)
do {
m = MatrixType::Random();
absdet = ei_abs(m.determinant());
} while(absdet < 10 * epsilon<Scalar>());
} while(absdet < epsilon<Scalar>());
MatrixType inv = m.inverse();
double error = double( (m*inv-MatrixType::Identity()).norm() * absdet / epsilon<Scalar>() );
error_sum += error;
@@ -89,8 +67,8 @@ template<typename MatrixType> void inverse_general_4x4(int repeat)
double error_avg = error_sum / repeat;
EIGEN_DEBUG_VAR(error_avg);
EIGEN_DEBUG_VAR(error_max);
VERIFY(error_avg < (NumTraits<Scalar>::IsComplex ? 8.4 : 1.4) );
VERIFY(error_max < (NumTraits<Scalar>::IsComplex ? 160.0 : 75.) );
VERIFY(error_avg < (NumTraits<Scalar>::IsComplex ? 8.0 : 1.0));
VERIFY(error_max < (NumTraits<Scalar>::IsComplex ? 64.0 : 16.0));
}
void test_prec_inverse_4x4()