diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -254,17 +254,17 @@ template<> struct ldlt_inplace typedef typename MatrixType::Index Index; eigen_assert(mat.rows()==mat.cols()); const Index size = mat.rows(); if (size <= 1) { transpositions.setIdentity(); if(sign) - *sign = real(mat.coeff(0,0))>0 ? 1:-1; + *sign = math::real(mat.coeff(0,0))>0 ? 1:-1; return true; } RealScalar cutoff(0), biggest_in_corner; for (Index k = 0; k < size; ++k) { // Find largest diagonal element @@ -273,52 +273,43 @@ template<> struct ldlt_inplace index_of_biggest_in_corner += k; if(k == 0) { // The biggest overall is the point of reference to which further diagonals // are compared; if any diagonal is negligible compared // to the largest overall, the algorithm bails. cutoff = abs(NumTraits::epsilon() * biggest_in_corner); - - if(sign) - *sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1; - } - else if(sign) - { - // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right - int newSign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0; - if(newSign != *sign) - *sign = 0; } // Finish early if the matrix is not full rank. if(biggest_in_corner < cutoff) { for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i; + if(sign) *sign = 0; break; } transpositions.coeffRef(k) = index_of_biggest_in_corner; if(k != 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) = conj(mat.coeff(index_of_biggest_in_corner,k)); + mat.coeffRef(index_of_biggest_in_corner,k) = math::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); @@ -329,16 +320,26 @@ template<> struct ldlt_inplace { temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint(); mat.coeffRef(k,k) -= (A10 * temp.head(k)).value(); if(rs>0) A21.noalias() -= A20 * temp.head(k); } if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff)) A21 /= mat.coeffRef(k,k); + + if(sign) + { + // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right + int newSign = math::real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0; + if(k == 0) + *sign = newSign; + else if(*sign != newSign) + *sign = 0; + } } return true; } // Reference for the algorithm: Davis and Hager, "Multiple Rank // Modifications of a Sparse Cholesky Factorization" (Algorithm 1) // Trivial rearrangements of their computations (Timothy E. Holy) @@ -362,30 +363,30 @@ template<> struct ldlt_inplace // Apply the update for (Index j = 0; j < size; j++) { // Check for termination due to an original decomposition of low-rank if (!(isfinite)(alpha)) break; // Update the diagonal terms - RealScalar dj = real(mat.coeff(j,j)); + RealScalar dj = math::real(mat.coeff(j,j)); Scalar wj = w.coeff(j); RealScalar swj2 = sigma*abs2(wj); RealScalar gamma = dj*alpha + swj2; mat.coeffRef(j,j) += swj2/alpha; alpha += swj2/dj; // Update the terms of L Index rs = size-j-1; w.tail(rs) -= wj * mat.col(j).tail(rs); if(gamma != 0) - mat.col(j).tail(rs) += (sigma*conj(wj)/gamma)*w.tail(rs); + mat.col(j).tail(rs) += (sigma*math::conj(wj)/gamma)*w.tail(rs); } return true; } template static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1) { // Apply the permutation to the input w diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -227,17 +227,17 @@ static typename MatrixType::Index llt_ra } } else { temp = vec; RealScalar beta = 1; for(Index j=0; j struct llt_inplace { @@ -272,17 +272,17 @@ template struct llt_inp for(Index k = 0; k < size; ++k) { Index rs = size-k-1; // remaining size Block A21(mat,k+1,k,rs,1); Block A10(mat,k,0,1,k); Block A20(mat,k+1,0,rs,k); - RealScalar x = real(mat.coeff(k,k)); + RealScalar x = math::real(mat.coeff(k,k)); if (k>0) x -= A10.squaredNorm(); if (x<=RealScalar(0)) return k; mat.coeffRef(k,k) = x = sqrt(x); if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint(); if (rs>0) A21 *= RealScalar(1)/x; } return -1; diff --git a/test/cholesky.cpp b/test/cholesky.cpp --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -263,20 +263,28 @@ template void chole // LDLT is not guaranteed to work for indefinite matrices, but happens to work fine if matrix is diagonal. // This test checks that LDLT reports correctly that matrix is indefinite. // See http://forum.kde.org/viewtopic.php?f=74&t=106942 template void cholesky_indefinite(const MatrixType& m) { eigen_assert(m.rows() == 2 && m.cols() == 2); MatrixType mat; - mat << 1, 0, 0, -1; - LDLT ldlt(mat); - VERIFY(!ldlt.isNegative()); - VERIFY(!ldlt.isPositive()); + { + mat << 1, 0, 0, -1; + LDLT ldlt(mat); + VERIFY(!ldlt.isNegative()); + VERIFY(!ldlt.isPositive()); + } + { + mat << 1, 2, 2, 1; + LDLT ldlt(mat); + VERIFY(!ldlt.isNegative()); + VERIFY(!ldlt.isPositive()); + } } template void cholesky_verify_assert() { MatrixType tmp; LLT llt; VERIFY_RAISES_ASSERT(llt.matrixL())