fix sparse LDLT with complexes

Eigen/src/Sparse/SparseLDLT.h

@@ -71,6 +71,8 @@ LDL License:
   *
   * \param MatrixType the type of the matrix of which we are computing the LDLT Cholesky decomposition
   *
+  * \warning the upper triangular part has to be specified. The rest of the matrix is not used. The input matrix must be column major.
+  *
   * \sa class LDLT, class LDLT
   */
 template<typename MatrixType, int Backend = DefaultBackend>
@@ -280,7 +282,7 @@ bool SparseLDLT<MatrixType,Backend>::_numeric(const MatrixType& a)
       Index i = Pinv ? Pinv[Ai[p]] : Ai[p]; /* get A(i,k) */
       if (i <= k)
       {
-        y[i] += Ax[p];            /* scatter A(i,k) into Y (sum duplicates) */
+        y[i] += ei_conj(Ax[p]);            /* scatter A(i,k) into Y (sum duplicates) */
         Index len;
         for (len = 0; tags[i] != k; i = m_parent[i])
         {
@@ -291,22 +293,23 @@ bool SparseLDLT<MatrixType,Backend>::_numeric(const MatrixType& a)
           pattern[--top] = pattern[--len];
       }
     }
+
     /* compute numerical values kth row of L (a sparse triangular solve) */
     m_diag[k] = y[k];                       /* get D(k,k) and clear Y(k) */
     y[k] = 0.0;
     for (; top < size; ++top)
     {
       Index i = pattern[top];       /* pattern[top:n-1] is pattern of L(:,k) */
-      Scalar yi = y[i];          /* get and clear Y(i) */
+      Scalar yi = (y[i]);             /* get and clear Y(i) */
       y[i] = 0.0;
       Index p2 = Lp[i] + m_nonZerosPerCol[i];
       Index p;
       for (p = Lp[i]; p < p2; ++p)
-        y[Li[p]] -= Lx[p] * yi;
+        y[Li[p]] -= ei_conj(Lx[p]) * (yi);
       Scalar l_ki = yi / m_diag[i];       /* the nonzero entry L(k,i) */
-      m_diag[k] -= l_ki * yi;
+      m_diag[k] -= l_ki * ei_conj(yi);
       Li[p] = k;                          /* store L(k,i) in column form of L */
-      Lx[p] = l_ki;
+      Lx[p] = (l_ki);
       ++m_nonZerosPerCol[i];              /* increment count of nonzeros in col i */
     }
     if (m_diag[k] == 0.0)
@@ -323,7 +326,7 @@ bool SparseLDLT<MatrixType,Backend>::_numeric(const MatrixType& a)
   return ok;  /* success, diagonal of D is all nonzero */
 }
 
-/** Computes b = L^-T L^-1 b */
+/** Computes b = L^-T D^-1 L^-1 b */
 template<typename MatrixType, int Backend>
 template<typename Derived>
 bool SparseLDLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
@@ -336,8 +339,6 @@ bool SparseLDLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
   if (m_matrix.nonZeros()>0) // otherwise L==I
     m_matrix.template triangularView<UnitLower>().solveInPlace(b);
   b = b.cwiseQuotient(m_diag);
-  // FIXME should be .adjoint() but it fails to compile...
-
   if (m_matrix.nonZeros()>0) // otherwise L==I
     m_matrix.adjoint().template triangularView<UnitUpper>().solveInPlace(b);
 

Eigen/src/Sparse/SparseLLT.h

@@ -195,14 +195,7 @@ bool SparseLLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
   ei_assert(size==b.rows());
 
   m_matrix.template triangularView<Lower>().solveInPlace(b);
-  // FIXME should be simply .adjoint() but it fails to compile...
+  m_matrix.adjoint().template triangularView<Upper>().solveInPlace(b);
-  if (NumTraits<Scalar>::IsComplex)
-  {
-    CholMatrixType aux = m_matrix.conjugate();
-    aux.transpose().template triangularView<Upper>().solveInPlace(b);
-  }
-  else
-    m_matrix.transpose().template triangularView<Upper>().solveInPlace(b);
 
   return true;
 }

test/sparse_solvers.cpp

@@ -149,26 +149,27 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
   }
 
   // test LDLT
-  if (!NumTraits<Scalar>::IsComplex)
   {
-    // TODO fix the issue with complex (see SparseLDLT::solveInPlace)
     SparseMatrix<Scalar> m2(rows, cols);
     DenseMatrix refMat2(rows, cols);
 
     DenseVector b = DenseVector::Random(cols);
     DenseVector refX(cols), x(cols);
 
-//     initSPD(density, refMat2, m2);
     initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0);
-    refMat2 += (refMat2.adjoint()).eval();
+    for(int i=0; i<rows; ++i)
-    refMat2.diagonal() *= 0.5;
+      m2.coeffRef(i,i) = refMat2(i,i) = ei_abs(ei_real(refMat2(i,i)));
 
-    refX = refMat2.llt().solve(b); // FIXME use LLT to compute the reference because LDLT seems to fail with large matrices
+    refX = refMat2.template selfadjointView<Upper>().llt().solve(b);
+    // FIXME use LLT to compute the reference because LDLT seems to fail with large matrices
     typedef SparseMatrix<Scalar,Upper|SelfAdjoint> SparseSelfAdjointMatrix;
     x = b;
     SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
     if (ldlt.succeeded())
       ldlt.solveInPlace(x);
+    else
+      std::cerr << "warning LDLT failed\n";
+
     VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default");
   }