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

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());
+    }
   }
 
 }