GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-19 16:19:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug fix in inverse for 1x1 matrix,

Browse Source 
some compilation fixes in sparse_solvers
This commit is contained in:
Gael Guennebaud 2009-07-27 18:09:56 +02:00
parent 94cc30180e
commit 5f3606bce9
4 changed files with 14 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Eigen/src/LU/Inverse.h
View File

@ -309,9 +309,8 @@ struct ei_compute_inverse_with_check<MatrixType, 1>
{ {
static inline bool run(const MatrixType& matrix, MatrixType* result) static inline bool run(const MatrixType& matrix, MatrixType* result)
{ {
if( 0 == result->coeffRef(0,0) ) return false;
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
if( matrix.coeff(0,0) == Scalar(0) ) return false;
result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0); result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
return true; return true;
} }

4
Eigen/src/Sparse/SparseLDLT.h
View File

@ -333,12 +333,12 @@ bool SparseLDLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
return false; return false;
if (m_matrix.nonZeros()>0) // otherwise L==I if (m_matrix.nonZeros()>0) // otherwise L==I
m_matrix.template triangular<LowerTriangular|UnitDiagBit>().solveInPlace(b); m_matrix.template triangular<UnitLowerTriangular>().solveInPlace(b);
b = b.cwise() / m_diag; b = b.cwise() / m_diag;
// FIXME should be .adjoint() but it fails to compile... // FIXME should be .adjoint() but it fails to compile...
if (m_matrix.nonZeros()>0) // otherwise L==I if (m_matrix.nonZeros()>0) // otherwise L==I
m_matrix.transpose().template triangular<UpperTriangular|UnitDiagBit>().solveInPlace(b); m_matrix.transpose().template triangular<UnitUpperTriangular>().solveInPlace(b);
return true; return true;
} }

5
test/inverse.cpp
View File

@ -76,10 +76,7 @@ template<typename MatrixType> void inverse(const MatrixType& m)
VectorType v3 = VectorType::Random(rows); VectorType v3 = VectorType::Random(rows);
MatrixType m3 = v3*v3.transpose(), m4(rows,cols); MatrixType m3 = v3*v3.transpose(), m4(rows,cols);
invertible = m3.computeInverseWithCheck( &m4 ); invertible = m3.computeInverseWithCheck( &m4 );
if( 1 == rows ){ VERIFY( rows==1 ? invertible : !invertible );
VERIFY( invertible ); }
else{
VERIFY( !invertible ); }
} }
void test_inverse() void test_inverse()

19
test/sparse_solvers.cpp
View File

@ -65,12 +65,12 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
// lower - dense // lower - dense
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().solveTriangular(vec2), VERIFY_IS_APPROX(refMat2.template triangularView<LowerTriangular>().solve(vec2),
m2.template triangular<LowerTriangular>().solve(vec3)); m2.template triangular<LowerTriangular>().solve(vec3));
// upper - dense // upper - dense
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.template marked<UpperTriangular>().solveTriangular(vec2), VERIFY_IS_APPROX(refMat2.template triangularView<UpperTriangular>().solve(vec2),
m2.template triangular<UpperTriangular>().solve(vec3)); m2.template triangular<UpperTriangular>().solve(vec3));
// TODO test row major // TODO test row major
@ -81,21 +81,21 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
// lower - sparse // lower - sparse
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular);
initSparse<Scalar>(density, refMatB, matB); initSparse<Scalar>(density, refMatB, matB);
refMat2.template marked<LowerTriangular>().solveTriangularInPlace(refMatB); refMat2.template triangularView<LowerTriangular>().solveInPlace(refMatB);
m2.template triangular<LowerTriangular>().solveInPlace(matB); m2.template triangular<LowerTriangular>().solveInPlace(matB);
VERIFY_IS_APPROX(matB.toDense(), refMatB); VERIFY_IS_APPROX(matB.toDense(), refMatB);
// upper - sparse // upper - sparse
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular);
initSparse<Scalar>(density, refMatB, matB); initSparse<Scalar>(density, refMatB, matB);
refMat2.template marked<UpperTriangular>().solveTriangularInPlace(refMatB); refMat2.template triangularView<UpperTriangular>().solveInPlace(refMatB);
m2.template triangular<UpperTriangular>().solveInPlace(matB); m2.template triangular<UpperTriangular>().solveInPlace(matB);
VERIFY_IS_APPROX(matB, refMatB); VERIFY_IS_APPROX(matB, refMatB);
// test deprecated API // test deprecated API
initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().solveTriangular(vec2), VERIFY_IS_APPROX(refMat2.template triangularView<LowerTriangular>().solve(vec2),
m2.template marked<LowerTriangular>().solveTriangular(vec3)); m2.template triangular<LowerTriangular>().solve(vec3));
} }
// test LLT // test LLT
@ -127,6 +127,7 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
x = b; x = b;
SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x); SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x);
VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)"); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)");
// TODO fix TAUCS with complexes
x = b; x = b;
SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x); SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x);
VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)"); VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)");
@ -151,7 +152,7 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
refMat2 += refMat2.adjoint(); refMat2 += refMat2.adjoint();
refMat2.diagonal() *= 0.5; refMat2.diagonal() *= 0.5;
refMat2.ldlt().solve(b, &refX); refMat2.llt().solve(b, &refX); // FIXME use LLT to compute the reference because LDLT seems to fail with large matrices
typedef SparseMatrix<Scalar,UpperTriangular|SelfAdjoint> SparseSelfAdjointMatrix; typedef SparseMatrix<Scalar,UpperTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
x = b; x = b;
SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2); SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
@ -228,8 +229,8 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
void test_sparse_solvers() void test_sparse_solvers()
{ {
for(int i = 0; i < g_repeat; i++) { for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST( sparse_solvers<double>(8, 8) ); // CALL_SUBTEST( sparse_solvers<double>(8, 8) );
CALL_SUBTEST( sparse_solvers<std::complex<double> >(16, 16) ); CALL_SUBTEST( sparse_solvers<std::complex<double> >(16, 16) );
CALL_SUBTEST( sparse_solvers<double>(101, 101) ); // CALL_SUBTEST( sparse_solvers<double>(100, 100) );
} }
} }