mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-02-17 18:09:55 +08:00
Code Issues Packages Projects Releases Wiki Activity

fix LDLT, now it really only uses a given triangular part!

Browse Source
This commit is contained in:
Gael Guennebaud 2010-06-09 14:01:06 +02:00
parent 201bd253ad
commit e242ac9345
4 changed files with 52 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
Eigen/src/Cholesky/LDLT.h
View File

@ -69,7 +69,6 @@ template<typename _MatrixType, int _UpLo> class LDLT
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef typename MatrixType::Index Index;
// typedef typename ei_plain_col_type<MatrixType, Index>::type IntColVectorType;
typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> TmpMatrixType;
typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
@ -251,10 +250,26 @@ template<> struct ei_ldlt_inplace<Lower>
transpositions.coeffRef(k) = index_of_biggest_in_corner;
if(k != index_of_biggest_in_corner)
{
mat.row(k).swap(mat.row(index_of_biggest_in_corner));
mat.col(k).swap(mat.col(index_of_biggest_in_corner));
// apply the transposition while taking care to consider only
// the lower triangular part
Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element
mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k));
mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s));
std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner));
for(int i=k+1;i<index_of_biggest_in_corner;++i)
{
Scalar tmp = mat.coeffRef(i,k);
mat.coeffRef(i,k) = ei_conj(mat.coeffRef(index_of_biggest_in_corner,i));
mat.coeffRef(index_of_biggest_in_corner,i) = ei_conj(tmp);
}
if(NumTraits<Scalar>::IsComplex)
mat.coeffRef(index_of_biggest_in_corner,k) = ei_conj(mat.coeff(index_of_biggest_in_corner,k));
}
// partition the matrix:
// A00 | - | -
// lu = A10 | A11 | -
// A20 | A21 | A22
Index rs = size - k - 1;
Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1);
Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);

5
Eigen/src/Cholesky/LLT.h
View File

@ -208,9 +208,12 @@ template<> struct ei_llt_inplace<Lower>
for (Index k=0; k<size; k+=blockSize)
{
// partition the matrix:
// A00 | - | -
// lu = A10 | A11 | -
// A20 | A21 | A22
Index bs = std::min(blockSize, size-k);
Index rs = size - k - bs;
Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs);
Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs);
Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);

12
Eigen/src/LU/PartialPivLU.h
View File

@ -339,9 +339,9 @@ struct ei_partial_lu_impl
Index tsize = size - k - bs; // trailing size
// partition the matrix:
// A00 | A01 | A02
// lu = A10 | A11 | A12
// A20 | A21 | A22
// A00 | A01 | A02
// lu = A_0 | A_1 | A_2 = A10 | A11 | A12
// A20 | A21 | A22
BlockType A_0(lu,0,0,rows,k);
BlockType A_2(lu,0,k+bs,rows,tsize);
BlockType A11(lu,k,k,bs,bs);
@ -350,8 +350,8 @@ struct ei_partial_lu_impl
BlockType A22(lu,k+bs,k+bs,trows,tsize);
Index nb_transpositions_in_panel;
// recursively calls the blocked LU algorithm with a very small
// blocking size:
// recursively call the blocked LU algorithm on [A11^T A21^T]^T
// with a very small blocking size:
if(!blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
row_transpositions+k, nb_transpositions_in_panel, 16))
{
@ -364,7 +364,7 @@ struct ei_partial_lu_impl
}
nb_transpositions += nb_transpositions_in_panel;
// update permutations and apply them to A10
// update permutations and apply them to A_0
for(Index i=k; i<k+bs; ++i)
{
Index piv = (row_transpositions[i] += k);

44
test/cholesky.cpp
View File

@ -121,14 +121,18 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
VERIFY_IS_APPROX(symm * matX, matB);
}
int sign = ei_random<int>()%2 ? 1 : -1;
if(sign == -1)
// LDLT
{
symm = -symm; // test a negative matrix
}
int sign = ei_random<int>()%2 ? 1 : -1;
if(sign == -1)
{
symm = -symm; // test a negative matrix
}
SquareMatrixType symmUp = symm.template triangularView<Upper>();
SquareMatrixType symmLo = symm.template triangularView<Lower>();
{
LDLT<SquareMatrixType,Lower> ldltlo(symmLo);
VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix());
vecX = ldltlo.solve(vecB);
@ -143,20 +147,20 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
matX = ldltup.solve(matB);
VERIFY_IS_APPROX(symm * matX, matB);
// if(MatrixType::RowsAtCompileTime==Dynamic)
// {
// // note : each inplace permutation requires a small temporary vector (mask)
//
// // check inplace solve
// matX = matB;
// VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0);
// VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval());
//
//
// matX = matB;
// VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0);
// VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval());
// }
if(MatrixType::RowsAtCompileTime==Dynamic)
{
// note : each inplace permutation requires a small temporary vector (mask)
// check inplace solve
matX = matB;
VERIFY_EVALUATION_COUNT(matX = ldltlo.solve(matX), 0);
VERIFY_IS_APPROX(matX, ldltlo.solve(matB).eval());
matX = matB;
VERIFY_EVALUATION_COUNT(matX = ldltup.solve(matX), 0);
VERIFY_IS_APPROX(matX, ldltup.solve(matB).eval());
}
}
}