diff options
-rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 21 | ||||
-rw-r--r-- | Eigen/src/Cholesky/LLT.h | 5 | ||||
-rw-r--r-- | Eigen/src/LU/PartialPivLU.h | 12 | ||||
-rw-r--r-- | 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<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); 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<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); 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<k+bs; ++i) { Index piv = (row_transpositions[i] += k); diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 4cc09ec05..666cee20f 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -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()); + } } } |