From d92de9574a2aec4af51ad04c0dc5cd2eb39e45bf Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 3 Jun 2010 11:56:08 +0200 Subject: [PATCH] fix sparse LDLT with complexes --- Eigen/src/Sparse/SparseLDLT.h | 27 ++++++++++++++------------- Eigen/src/Sparse/SparseLLT.h | 9 +-------- test/sparse_solvers.cpp | 13 +++++++------ 3 files changed, 22 insertions(+), 27 deletions(-) diff --git a/Eigen/src/Sparse/SparseLDLT.h b/Eigen/src/Sparse/SparseLDLT.h index ae1a96b4f..a6785d0af 100644 --- a/Eigen/src/Sparse/SparseLDLT.h +++ b/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 @@ -213,7 +215,7 @@ void SparseLDLT::_symbolic(const MatrixType& a) m_parent[k] = -1; /* parent of k is not yet known */ tags[k] = k; /* mark node k as visited */ m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ - Index kk = P ? P[k] : k; /* kth original, or permuted, column */ + Index kk = P ? P[k] : k; /* kth original, or permuted, column */ Index p2 = Ap[kk+1]; for (Index p = Ap[kk]; p < p2; ++p) { @@ -269,10 +271,10 @@ bool SparseLDLT::_numeric(const MatrixType& a) for (Index k = 0; k < size; ++k) { /* compute nonzero pattern of kth row of L, in topological order */ - y[k] = 0.0; /* Y(0:k) is now all zero */ + y[k] = 0.0; /* Y(0:k) is now all zero */ Index top = size; /* stack for pattern is empty */ - tags[k] = k; /* mark node k as visited */ - m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ + tags[k] = k; /* mark node k as visited */ + m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ Index kk = (P) ? (P[k]) : (k); /* kth original, or permuted, column */ Index p2 = Ap[kk+1]; for (Index p = Ap[kk]; p < p2; ++p) @@ -280,7 +282,7 @@ bool SparseLDLT::_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::_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) */ + Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */ + 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::_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 template bool SparseLDLT::solveInPlace(MatrixBase &b) const @@ -336,8 +339,6 @@ bool SparseLDLT::solveInPlace(MatrixBase &b) const if (m_matrix.nonZeros()>0) // otherwise L==I m_matrix.template triangularView().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().solveInPlace(b); diff --git a/Eigen/src/Sparse/SparseLLT.h b/Eigen/src/Sparse/SparseLLT.h index 4ec3ee009..47d58f8e6 100644 --- a/Eigen/src/Sparse/SparseLLT.h +++ b/Eigen/src/Sparse/SparseLLT.h @@ -195,14 +195,7 @@ bool SparseLLT::solveInPlace(MatrixBase &b) const ei_assert(size==b.rows()); m_matrix.template triangularView().solveInPlace(b); - // FIXME should be simply .adjoint() but it fails to compile... - if (NumTraits::IsComplex) - { - CholMatrixType aux = m_matrix.conjugate(); - aux.transpose().template triangularView().solveInPlace(b); - } - else - m_matrix.transpose().template triangularView().solveInPlace(b); + m_matrix.adjoint().template triangularView().solveInPlace(b); return true; } diff --git a/test/sparse_solvers.cpp b/test/sparse_solvers.cpp index 00df1bffd..ea8aee718 100644 --- a/test/sparse_solvers.cpp +++ b/test/sparse_solvers.cpp @@ -149,26 +149,27 @@ template void sparse_solvers(int rows, int cols) } // test LDLT - if (!NumTraits::IsComplex) { - // TODO fix the issue with complex (see SparseLDLT::solveInPlace) SparseMatrix m2(rows, cols); DenseMatrix refMat2(rows, cols); DenseVector b = DenseVector::Random(cols); DenseVector refX(cols), x(cols); -// initSPD(density, refMat2, m2); initSparse(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0); - refMat2 += (refMat2.adjoint()).eval(); - refMat2.diagonal() *= 0.5; + for(int i=0; i().llt().solve(b); + // FIXME use LLT to compute the reference because LDLT seems to fail with large matrices typedef SparseMatrix SparseSelfAdjointMatrix; x = b; SparseLDLT ldlt(m2); if (ldlt.succeeded()) ldlt.solveInPlace(x); + else + std::cerr << "warning LDLT failed\n"; + VERIFY(refX.isApprox(x,test_precision()) && "LDLT: default"); }