diff options
Diffstat (limited to 'Eigen/src/Cholesky/LDLT.h')
-rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 37 |
1 files changed, 19 insertions, 18 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 9bd60d648..05b300c56 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -259,7 +259,7 @@ template<> struct ldlt_inplace<Lower> { transpositions.setIdentity(); if(sign) - *sign = real(mat.coeff(0,0))>0 ? 1:-1; + *sign = numext::real(mat.coeff(0,0))>0 ? 1:-1; return true; } @@ -278,22 +278,13 @@ template<> struct ldlt_inplace<Lower> // are compared; if any diagonal is negligible compared // to the largest overall, the algorithm bails. cutoff = abs(NumTraits<Scalar>::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; } @@ -309,11 +300,11 @@ template<> struct ldlt_inplace<Lower> for(int i=k+1;i<index_of_biggest_in_corner;++i) { Scalar tmp = mat.coeffRef(i,k); - mat.coeffRef(i,k) = conj(mat.coeffRef(index_of_biggest_in_corner,i)); - mat.coeffRef(index_of_biggest_in_corner,i) = conj(tmp); + mat.coeffRef(i,k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner,i)); + mat.coeffRef(index_of_biggest_in_corner,i) = numext::conj(tmp); } if(NumTraits<Scalar>::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) = numext::conj(mat.coeff(index_of_biggest_in_corner,k)); } // partition the matrix: @@ -334,6 +325,16 @@ template<> struct ldlt_inplace<Lower> } 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 = numext::real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0; + if(k == 0) + *sign = newSign; + else if(*sign != newSign) + *sign = 0; + } } return true; @@ -349,7 +350,7 @@ template<> struct ldlt_inplace<Lower> template<typename MatrixType, typename WDerived> static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w, const typename MatrixType::RealScalar& sigma=1) { - using internal::isfinite; + using numext::isfinite; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::Index Index; @@ -367,9 +368,9 @@ template<> struct ldlt_inplace<Lower> break; // Update the diagonal terms - RealScalar dj = real(mat.coeff(j,j)); + RealScalar dj = numext::real(mat.coeff(j,j)); Scalar wj = w.coeff(j); - RealScalar swj2 = sigma*abs2(wj); + RealScalar swj2 = sigma*numext::abs2(wj); RealScalar gamma = dj*alpha + swj2; mat.coeffRef(j,j) += swj2/alpha; @@ -380,7 +381,7 @@ template<> struct ldlt_inplace<Lower> 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*numext::conj(wj)/gamma)*w.tail(rs); } return true; } |