From e242ac9345e728d8e347bb0d48497891e3891276 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 9 Jun 2010 14:01:06 +0200 Subject: [PATCH] fix LDLT, now it really only uses a given triangular part! --- Eigen/src/Cholesky/LDLT.h | 21 +++++++++++++++--- Eigen/src/Cholesky/LLT.h | 5 ++++- Eigen/src/LU/PartialPivLU.h | 12 +++++----- test/cholesky.cpp | 44 ++++++++++++++++++++----------------- 4 files changed, 52 insertions(+), 30 deletions(-) diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 79657f086..e9fa4a5d2 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -69,7 +69,6 @@ template class LDLT typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; typedef typename MatrixType::Index Index; -// typedef typename ei_plain_col_type::type IntColVectorType; typedef Matrix TmpMatrixType; typedef Transpositions TranspositionType; @@ -251,10 +250,26 @@ template<> struct ei_ldlt_inplace 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::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 A21(mat,k+1,k,rs,1); Block A10(mat,k,0,1,k); diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 6e853436c..6e7660ddf 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -208,9 +208,12 @@ template<> struct ei_llt_inplace for (Index k=0; k A11(m,k, k, bs,bs); Block A21(m,k+bs,k, rs,bs); Block A22(m,k+bs,k+bs,rs,rs); diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index a9172289c..7e7516a5c 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -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 void cholesky(const MatrixType& m) VERIFY_IS_APPROX(symm * matX, matB); } - int sign = ei_random()%2 ? 1 : -1; - - if(sign == -1) + // LDLT { - symm = -symm; // test a negative matrix - } + int sign = ei_random()%2 ? 1 : -1; + + if(sign == -1) + { + symm = -symm; // test a negative matrix + } + + SquareMatrixType symmUp = symm.template triangularView(); + SquareMatrixType symmLo = symm.template triangularView(); - { LDLT ldltlo(symmLo); VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix()); vecX = ldltlo.solve(vecB); @@ -143,20 +147,20 @@ template 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()); + } } }