devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

reduce float warnings (comparisons and implicit conversions)

Browse Source
This commit is contained in:
Erik Schultheis
2022-01-26 18:16:19 +00:00
committed by Rasmus Munk Larsen
parent 51311ec651
commit d271a7d545
41 changed files with 152 additions and 133 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Eigen/src/Cholesky/LDLT.h
View File

@@ -440,7 +440,7 @@ template<> struct ldlt_inplace<Lower>
// Update the terms of L
Index rs = size-j-1;
w.tail(rs) -= wj * mat.col(j).tail(rs);
if(gamma != 0)
if(!numext::is_exactly_zero(gamma))
mat.col(j).tail(rs) += (sigma*numext::conj(wj)/gamma)*w.tail(rs);
}
return true;

2
Eigen/src/Cholesky/LLT.h
View File

@@ -297,7 +297,7 @@ static Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const
if(rs)
{
temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs);
if(gamma != 0)
if(!numext::is_exactly_zero(gamma))
mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs);
}
}

10
Eigen/src/Core/ConditionEstimator.h
View File

@@ -162,12 +162,12 @@ rcond_estimate_helper(typename Decomposition::RealScalar matrix_norm, const Deco
{
typedef typename Decomposition::RealScalar RealScalar;
eigen_assert(dec.rows() == dec.cols());
if (dec.rows() == 0) return NumTraits<RealScalar>::infinity();
if (matrix_norm == RealScalar(0)) return RealScalar(0);
if (dec.rows() == 1) return RealScalar(1);
if (dec.rows() == 0) return NumTraits<RealScalar>::infinity();
if (numext::is_exactly_zero(matrix_norm)) return RealScalar(0);
if (dec.rows() == 1) return RealScalar(1);
const RealScalar inverse_matrix_norm = rcond_invmatrix_L1_norm_estimate(dec);
return (inverse_matrix_norm == RealScalar(0) ? RealScalar(0)
: (RealScalar(1) / inverse_matrix_norm) / matrix_norm);
return (numext::is_exactly_zero(inverse_matrix_norm) ? RealScalar(0)
: (RealScalar(1) / inverse_matrix_norm) / matrix_norm);
}
} // namespace internal

3
Eigen/src/Core/GeneralProduct.h
View File

@@ -219,7 +219,6 @@ template<> struct gemv_dense_selector<OnTheRight,ColMajor,true>
typedef typename Lhs::Scalar LhsScalar;
typedef typename Rhs::Scalar RhsScalar;
typedef typename Dest::Scalar ResScalar;
typedef typename Dest::RealScalar RealScalar;
typedef internal::blas_traits<Lhs> LhsBlasTraits;
typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType;
@@ -264,7 +263,7 @@ template<> struct gemv_dense_selector<OnTheRight,ColMajor,true>
{
gemv_static_vector_if<ResScalar,ActualDest::SizeAtCompileTime,ActualDest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
const bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0));
const bool alphaIsCompatible = (!ComplexByReal) || (numext::is_exactly_zero(numext::imag(actualAlpha)));
const bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),

14
Eigen/src/Core/MathFunctionsImpl.h
View File

@@ -119,7 +119,7 @@ RealScalar positive_real_hypot(const RealScalar& x, const RealScalar& y)
EIGEN_USING_STD(sqrt);
RealScalar p, qp;
p = numext::maxi(x,y);
if(p==RealScalar(0)) return RealScalar(0);
if(numext::is_exactly_zero(p)) return RealScalar(0);
qp = numext::mini(y,x) / p;
return p * sqrt(RealScalar(1) + qp*qp);
}
@@ -169,8 +169,8 @@ EIGEN_DEVICE_FUNC std::complex<T> complex_sqrt(const std::complex<T>& z) {
return
(numext::isinf)(y) ? std::complex<T>(NumTraits<T>::infinity(), y)
: x == zero ? std::complex<T>(w, y < zero ? -w : w)
: x > zero ? std::complex<T>(w, y / (2 * w))
: numext::is_exactly_zero(x) ? std::complex<T>(w, y < zero ? -w : w)
: x > zero ? std::complex<T>(w, y / (2 * w))
: std::complex<T>(numext::abs(y) / (2 * w), y < zero ? -w : w );
}
@@ -208,10 +208,10 @@ EIGEN_DEVICE_FUNC std::complex<T> complex_rsqrt(const std::complex<T>& z) {
const T woz = w / abs_z;
// Corner cases consistent with 1/sqrt(z) on gcc/clang.
return
abs_z == zero ? std::complex<T>(NumTraits<T>::infinity(), NumTraits<T>::quiet_NaN())
: ((numext::isinf)(x) || (numext::isinf)(y)) ? std::complex<T>(zero, zero)
: x == zero ? std::complex<T>(woz, y < zero ? woz : -woz)
: x > zero ? std::complex<T>(woz, -y / (2 * w * abs_z))
numext::is_exactly_zero(abs_z) ? std::complex<T>(NumTraits<T>::infinity(), NumTraits<T>::quiet_NaN())
: ((numext::isinf)(x) || (numext::isinf)(y)) ? std::complex<T>(zero, zero)
: numext::is_exactly_zero(x) ? std::complex<T>(woz, y < zero ? woz : -woz)
: x > zero ? std::complex<T>(woz, -y / (2 * w * abs_z))
: std::complex<T>(numext::abs(y) / (2 * w * abs_z), y < zero ? woz : -woz );
}

2
Eigen/src/Core/ProductEvaluators.h
View File

@@ -460,7 +460,7 @@ protected:
void eval_dynamic_impl(Dst& dst, const LhsT& lhs, const RhsT& rhs, const Func &func, const Scalar& s /* == 1 */, false_type)
{
EIGEN_UNUSED_VARIABLE(s);
eigen_internal_assert(s==Scalar(1));
eigen_internal_assert(numext::is_exactly_one(s));
call_restricted_packet_assignment_no_alias(dst, lhs.lazyProduct(rhs), func);
}

4
Eigen/src/Core/products/TriangularMatrixMatrix.h
View File

@@ -453,12 +453,12 @@ struct triangular_product_impl<Mode,LhsIsTriangular,Lhs,false,Rhs,false>
// Apply correction if the diagonal is unit and a scalar factor was nested:
if ((Mode&UnitDiag)==UnitDiag)
{
if (LhsIsTriangular && lhs_alpha!=LhsScalar(1))
if (LhsIsTriangular && !numext::is_exactly_one(lhs_alpha))
{
Index diagSize = (std::min)(lhs.rows(),lhs.cols());
dst.topRows(diagSize) -= ((lhs_alpha-LhsScalar(1))*a_rhs).topRows(diagSize);
}
else if ((!LhsIsTriangular) && rhs_alpha!=RhsScalar(1))
else if ((!LhsIsTriangular) && !numext::is_exactly_one(rhs_alpha))
{
Index diagSize = (std::min)(rhs.rows(),rhs.cols());
dst.leftCols(diagSize) -= (rhs_alpha-RhsScalar(1))*a_lhs.leftCols(diagSize);

7
Eigen/src/Core/products/TriangularMatrixVector.h
View File

@@ -211,7 +211,6 @@ template<int Mode> struct trmv_selector<Mode,ColMajor>
typedef typename Lhs::Scalar LhsScalar;
typedef typename Rhs::Scalar RhsScalar;
typedef typename Dest::Scalar ResScalar;
typedef typename Dest::RealScalar RealScalar;
typedef internal::blas_traits<Lhs> LhsBlasTraits;
typedef typename LhsBlasTraits::DirectLinearAccessType ActualLhsType;
@@ -237,7 +236,7 @@ template<int Mode> struct trmv_selector<Mode,ColMajor>
gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
bool alphaIsCompatible = (!ComplexByReal) || (numext::imag(actualAlpha)==RealScalar(0));
bool alphaIsCompatible = (!ComplexByReal) || numext::is_exactly_zero(numext::imag(actualAlpha));
bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
@@ -278,7 +277,7 @@ template<int Mode> struct trmv_selector<Mode,ColMajor>
dest = MappedDest(actualDestPtr, dest.size());
}
if ( ((Mode&UnitDiag)==UnitDiag) && (lhs_alpha!=LhsScalar(1)) )
if ( ((Mode&UnitDiag)==UnitDiag) && !numext::is_exactly_one(lhs_alpha) )
{
Index diagSize = (std::min)(lhs.rows(),lhs.cols());
dest.head(diagSize) -= (lhs_alpha-LhsScalar(1))*rhs.head(diagSize);
@@ -337,7 +336,7 @@ template<int Mode> struct trmv_selector<Mode,RowMajor>
dest.data(),dest.innerStride(),
actualAlpha);
if ( ((Mode&UnitDiag)==UnitDiag) && (lhs_alpha!=LhsScalar(1)) )
if ( ((Mode&UnitDiag)==UnitDiag) && !numext::is_exactly_one(lhs_alpha) )
{
Index diagSize = (std::min)(lhs.rows(),lhs.cols());
dest.head(diagSize) -= (lhs_alpha-LhsScalar(1))*rhs.head(diagSize);

16
Eigen/src/Core/util/Meta.h
View File

@@ -397,6 +397,8 @@ struct aligned_storage {
} // end namespace internal
template<typename T> struct NumTraits;
namespace numext {
#if defined(EIGEN_GPU_COMPILE_PHASE)
@@ -429,6 +431,20 @@ template<> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC
bool equal_strict(const double& x,const double& y) { return std::equal_to<double>()(x,y); }
#endif
/**
* \internal Performs an exact comparison of x to zero, e.g. to decide whether a term can be ignored.
* Use this to to bypass -Wfloat-equal warnings when exact zero is what needs to be tested.
*/
template<typename X> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC
bool is_exactly_zero(const X& x) { return equal_strict(x, typename NumTraits<X>::Literal{0}); }
/**
* \internal Performs an exact comparison of x to one, e.g. to decide whether a factor needs to be multiplied.
* Use this to to bypass -Wfloat-equal warnings when exact one is what needs to be tested.
*/
template<typename X> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC
bool is_exactly_one(const X& x) { return equal_strict(x, typename NumTraits<X>::Literal{1}); }
template<typename X, typename Y> EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC
bool not_equal_strict(const X& x,const Y& y) { return x != y; }

2
Eigen/src/Eigenvalues/ComplexSchur.h
View File

@@ -308,7 +308,7 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
// In this case, det==0, and all we have to do is checking that eival2_norm!=0
if(eival1_norm > eival2_norm)
eival2 = det / eival1;
else if(eival2_norm!=RealScalar(0))
else if(!numext::is_exactly_zero(eival2_norm))
eival1 = det / eival2;
// choose the eigenvalue closest to the bottom entry of the diagonal

10
Eigen/src/Eigenvalues/RealQZ.h
View File

@@ -239,7 +239,7 @@ namespace Eigen {
for (Index i=dim-1; i>=j+2; i--) {
JRs G;
// kill S(i,j)
if(m_S.coeff(i,j) != 0)
if(!numext::is_exactly_zero(m_S.coeff(i, j)))
{
G.makeGivens(m_S.coeff(i-1,j), m_S.coeff(i,j), &m_S.coeffRef(i-1, j));
m_S.coeffRef(i,j) = Scalar(0.0);
@@ -250,7 +250,7 @@ namespace Eigen {
m_Q.applyOnTheRight(i-1,i,G);
}
// kill T(i,i-1)
if(m_T.coeff(i,i-1)!=Scalar(0))
if(!numext::is_exactly_zero(m_T.coeff(i, i - 1)))
{
G.makeGivens(m_T.coeff(i,i), m_T.coeff(i,i-1), &m_T.coeffRef(i,i));
m_T.coeffRef(i,i-1) = Scalar(0.0);
@@ -288,7 +288,7 @@ namespace Eigen {
while (res > 0)
{
Scalar s = abs(m_S.coeff(res-1,res-1)) + abs(m_S.coeff(res,res));
if (s == Scalar(0.0))
if (numext::is_exactly_zero(s))
s = m_normOfS;
if (abs(m_S.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
break;
@@ -318,7 +318,7 @@ namespace Eigen {
using std::abs;
using std::sqrt;
const Index dim=m_S.cols();
if (abs(m_S.coeff(i+1,i))==Scalar(0))
if (numext::is_exactly_zero(abs(m_S.coeff(i + 1, i))))
return;
Index j = findSmallDiagEntry(i,i+1);
if (j==i-1)
@@ -629,7 +629,7 @@ namespace Eigen {
{
for(Index i=0; i<dim-1; ++i)
{
if(m_S.coeff(i+1, i) != Scalar(0))
if(!numext::is_exactly_zero(m_S.coeff(i + 1, i)))
{
JacobiRotation<Scalar> j_left, j_right;
internal::real_2x2_jacobi_svd(m_T, i, i+1, &j_left, &j_right);

6
Eigen/src/Eigenvalues/RealSchur.h
View File

@@ -314,7 +314,7 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMa
Scalar considerAsZero = numext::maxi<Scalar>( norm * numext::abs2(NumTraits<Scalar>::epsilon()),
(std::numeric_limits<Scalar>::min)() );
if(norm!=Scalar(0))
if(!numext::is_exactly_zero(norm))
{
while (iu >= 0)
{
@@ -517,7 +517,7 @@ inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Inde
Matrix<Scalar, 2, 1> ess;
v.makeHouseholder(ess, tau, beta);
if (beta != Scalar(0)) // if v is not zero
if (!numext::is_exactly_zero(beta)) // if v is not zero
{
if (firstIteration && k > il)
m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
@@ -537,7 +537,7 @@ inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Inde
Matrix<Scalar, 1, 1> ess;
v.makeHouseholder(ess, tau, beta);
if (beta != Scalar(0)) // if v is not zero
if (!numext::is_exactly_zero(beta)) // if v is not zero
{
m_matT.coeffRef(iu-1, iu-2) = beta;
m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);

14
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@@ -447,7 +447,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
// map the matrix coefficients to [-1:1] to avoid over- and underflow.
mat = matrix.template triangularView<Lower>();
RealScalar scale = mat.cwiseAbs().maxCoeff();
if(scale==RealScalar(0)) scale = RealScalar(1);
if(numext::is_exactly_zero(scale)) scale = RealScalar(1);
mat.template triangularView<Lower>() /= scale;
m_subdiag.resize(n-1);
m_hcoeffs.resize(n-1);
@@ -526,7 +526,7 @@ ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag
}
// find the largest unreduced block at the end of the matrix.
while (end>0 && subdiag[end-1]==RealScalar(0))
while (end>0 && numext::is_exactly_zero(subdiag[end - 1]))
{
end--;
}
@@ -538,7 +538,7 @@ ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag
if(iter > maxIterations * n) break;
start = end - 1;
while (start>0 && subdiag[start-1]!=0)
while (start>0 && !numext::is_exactly_zero(subdiag[start - 1]))
start--;
internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n);
@@ -843,12 +843,12 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta
// RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
// This explain the following, somewhat more complicated, version:
RealScalar mu = diag[end];
if(td==RealScalar(0)) {
if(numext::is_exactly_zero(td)) {
mu -= numext::abs(e);
} else if (e != RealScalar(0)) {
} else if (!numext::is_exactly_zero(e)) {
const RealScalar e2 = numext::abs2(e);
const RealScalar h = numext::hypot(td,e);
if(e2 == RealScalar(0)) {
if(numext::is_exactly_zero(e2)) {
mu -= e / ((td + (td>RealScalar(0) ? h : -h)) / e);
} else {
mu -= e2 / (td + (td>RealScalar(0) ? h : -h));
@@ -859,7 +859,7 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta
RealScalar z = subdiag[start];
// If z ever becomes zero, the Givens rotation will be the identity and
// z will stay zero for all future iterations.
for (Index k = start; k < end && z != RealScalar(0); ++k)
for (Index k = start; k < end && !numext::is_exactly_zero(z); ++k)
{
JacobiRotation<RealScalar> rot;
rot.makeGivens(x, z);

4
Eigen/src/Householder/Householder.h
View File

@@ -124,7 +124,7 @@ void MatrixBase<Derived>::applyHouseholderOnTheLeft(
{
*this *= Scalar(1)-tau;
}
else if(tau!=Scalar(0))
else if(!numext::is_exactly_zero(tau))
{
Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace,cols());
Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols());
@@ -162,7 +162,7 @@ void MatrixBase<Derived>::applyHouseholderOnTheRight(
{
*this *= Scalar(1)-tau;
}
else if(tau!=Scalar(0))
else if(!numext::is_exactly_zero(tau))
{
Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows());
Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1);

6
Eigen/src/Jacobi/Jacobi.h
View File

@@ -234,13 +234,13 @@ void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar
{
using std::sqrt;
using std::abs;
if(q==Scalar(0))
if(numext::is_exactly_zero(q))
{
m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1);
m_s = Scalar(0);
if(r) *r = abs(p);
}
else if(p==Scalar(0))
else if(numext::is_exactly_zero(p))
{
m_c = Scalar(0);
m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1);
@@ -468,7 +468,7 @@ void /*EIGEN_DONT_INLINE*/ apply_rotation_in_the_plane(DenseBase<VectorX>& xpr_x
OtherScalar c = j.c();
OtherScalar s = j.s();
if (c==OtherScalar(1) && s==OtherScalar(0))
if (numext::is_exactly_one(c) && numext::is_exactly_zero(s))
return;
apply_rotation_in_the_plane_selector<

2
Eigen/src/LU/FullPivLU.h
View File

@@ -519,7 +519,7 @@ void FullPivLU<MatrixType>::computeInPlace()
row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
col_of_biggest_in_corner += k; // need to add k to them.
if(biggest_in_corner==Score(0))
if(numext::is_exactly_zero(biggest_in_corner))
{
// before exiting, make sure to initialize the still uninitialized transpositions
// in a sane state without destroying what we already have.

4
Eigen/src/LU/PartialPivLU.h
View File

@@ -378,7 +378,7 @@ struct partial_lu_impl
row_transpositions[k] = PivIndex(row_of_biggest_in_col);
if(biggest_in_corner != Score(0))
if(!numext::is_exactly_zero(biggest_in_corner))
{
if(k != row_of_biggest_in_col)
{
@@ -404,7 +404,7 @@ struct partial_lu_impl
{
Index k = endk;
row_transpositions[k] = PivIndex(k);
if (Scoring()(lu(k, k)) == Score(0) && first_zero_pivot == -1)
if (numext::is_exactly_zero(Scoring()(lu(k, k))) && first_zero_pivot == -1)
first_zero_pivot = k;
}

2
Eigen/src/QR/ColPivHouseholderQR.h
View File

@@ -552,7 +552,7 @@ void ColPivHouseholderQR<MatrixType>::computeInPlace()
// http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
// and used in LAPACK routines xGEQPF and xGEQP3.
// See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html
if (m_colNormsUpdated.coeffRef(j) != RealScalar(0)) {
if (!numext::is_exactly_zero(m_colNormsUpdated.coeffRef(j))) {
RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j);
temp = (RealScalar(1) + temp) * (RealScalar(1) - temp);
temp = temp < RealScalar(0) ? RealScalar(0) : temp;

2
Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
View File

@@ -139,7 +139,7 @@ class SPQR : public SparseSolverBase<SPQR<MatrixType_> >
{
RealScalar max2Norm = 0.0;
for (int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm());
if(max2Norm==RealScalar(0))
if(numext::is_exactly_zero(max2Norm))
max2Norm = RealScalar(1);
pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon();
}

40
Eigen/src/SVD/BDCSVD.h
View File

@@ -282,7 +282,7 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
return *this;
}
if(scale==Literal(0)) scale = Literal(1);
if(numext::is_exactly_zero(scale)) scale = Literal(1);
MatrixX copy;
if (m_isTranspose) copy = matrix.adjoint()/scale;
else copy = matrix/scale;
@@ -621,7 +621,10 @@ void BDCSVD<MatrixType>::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, Ma
// but others are interleaved and we must ignore them at this stage.
// To this end, let's compute a permutation skipping them:
Index actual_n = n;
while(actual_n>1 && diag(actual_n-1)==Literal(0)) {--actual_n; eigen_internal_assert(col0(actual_n)==Literal(0)); }
while(actual_n>1 && numext::is_exactly_zero(diag(actual_n - 1))) {
--actual_n;
eigen_internal_assert(numext::is_exactly_zero(col0(actual_n)));
}
Index m = 0; // size of the deflated problem
for(Index k=0;k<actual_n;++k)
if(abs(col0(k))>considerZero)
@@ -753,11 +756,11 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
Index actual_n = n;
// Note that here actual_n is computed based on col0(i)==0 instead of diag(i)==0 as above
// because 1) we have diag(i)==0 => col0(i)==0 and 2) if col0(i)==0, then diag(i) is already a singular value.
while(actual_n>1 && col0(actual_n-1)==Literal(0)) --actual_n;
while(actual_n>1 && numext::is_exactly_zero(col0(actual_n - 1))) --actual_n;
for (Index k = 0; k < n; ++k)
{
if (col0(k) == Literal(0) || actual_n==1)
if (numext::is_exactly_zero(col0(k)) || actual_n == 1)
{
// if col0(k) == 0, then entry is deflated, so singular value is on diagonal
// if actual_n==1, then the deflated problem is already diagonalized
@@ -778,7 +781,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
// recall that at this stage we assume that z[j]!=0 and all entries for which z[j]==0 have been put aside.
// This should be equivalent to using perm[]
Index l = k+1;
while(col0(l)==Literal(0)) { ++l; eigen_internal_assert(l<actual_n); }
while(numext::is_exactly_zero(col0(l))) { ++l; eigen_internal_assert(l < actual_n); }
right = diag(l);
}
@@ -813,7 +816,8 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
{
// check that after the shift, f(mid) is still negative:
RealScalar midShifted = (right - left) / RealScalar(2);
if(shift==right)
// we can test exact equality here, because shift comes from `... ? left : right`
if(numext::equal_strict(shift, right))
midShifted = -midShifted;
RealScalar fMidShifted = secularEq(midShifted, col0, diag, perm, diagShifted, shift);
if(fMidShifted>0)
@@ -826,7 +830,8 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
// initial guess
RealScalar muPrev, muCur;
if (shift == left)
// we can test exact equality here, because shift comes from `... ? left : right`
if (numext::equal_strict(shift, left))
{
muPrev = (right - left) * RealScalar(0.1);
if (k == actual_n-1) muCur = right - left;
@@ -849,7 +854,7 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
// rational interpolation: fit a function of the form a / mu + b through the two previous
// iterates and use its zero to compute the next iterate
bool useBisection = fPrev*fCur>Literal(0);
while (fCur!=Literal(0) && abs(muCur - muPrev) > Literal(8) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(muCur), abs(muPrev)) && abs(fCur - fPrev)>NumTraits<RealScalar>::epsilon() && !useBisection)
while (!numext::is_exactly_zero(fCur) && abs(muCur - muPrev) > Literal(8) * NumTraits<RealScalar>::epsilon() * numext::maxi<RealScalar>(abs(muCur), abs(muPrev)) && abs(fCur - fPrev) > NumTraits<RealScalar>::epsilon() && !useBisection)
{
++m_numIters;
@@ -869,8 +874,9 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
muCur = muZero;
fCur = fZero;
if (shift == left && (muCur < Literal(0) || muCur > right - left)) useBisection = true;
if (shift == right && (muCur < -(right - left) || muCur > Literal(0))) useBisection = true;
// we can test exact equality here, because shift comes from `... ? left : right`
if (numext::equal_strict(shift, left) && (muCur < Literal(0) || muCur > right - left)) useBisection = true;
if (numext::equal_strict(shift, right) && (muCur < -(right - left) || muCur > Literal(0))) useBisection = true;
if (abs(fCur)>abs(fPrev)) useBisection = true;
}
@@ -881,7 +887,8 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
std::cout << "useBisection for k = " << k << ", actual_n = " << actual_n << "\n";
#endif
RealScalar leftShifted, rightShifted;
if (shift == left)
// we can test exact equality here, because shift comes from `... ? left : right`
if (numext::equal_strict(shift, left))
{
// to avoid overflow, we must have mu > max(real_min, |z(k)|/sqrt(real_max)),
// the factor 2 is to be more conservative
@@ -959,7 +966,8 @@ void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& d
// Instead fo abbording or entering an infinite loop,
// let's just use the middle as the estimated zero-crossing:
muCur = (right - left) * RealScalar(0.5);
if(shift == right)
// we can test exact equality here, because shift comes from `... ? left : right`
if(numext::equal_strict(shift, right))
muCur = -muCur;
}
}
@@ -1004,7 +1012,7 @@ void BDCSVD<MatrixType>::perturbCol0
// The offset permits to skip deflated entries while computing zhat
for (Index k = 0; k < n; ++k)
{
if (col0(k) == Literal(0)) // deflated
if (numext::is_exactly_zero(col0(k))) // deflated
zhat(k) = Literal(0);
else
{
@@ -1077,7 +1085,7 @@ void BDCSVD<MatrixType>::computeSingVecs
for (Index k = 0; k < n; ++k)
{
if (zhat(k) == Literal(0))
if (numext::is_exactly_zero(zhat(k)))
{
U.col(k) = VectorType::Unit(n+1, k);
if (m_compV) V.col(k) = VectorType::Unit(n, k);
@@ -1123,7 +1131,7 @@ void BDCSVD<MatrixType>::deflation43(Eigen::Index firstCol, Eigen::Index shift,
RealScalar c = m_computed(start, start);
RealScalar s = m_computed(start+i, start);
RealScalar r = numext::hypot(c,s);
if (r == Literal(0))
if (numext::is_exactly_zero(r))
{
m_computed(start+i, start+i) = Literal(0);
return;
@@ -1163,7 +1171,7 @@ void BDCSVD<MatrixType>::deflation44(Eigen::Index firstColu , Eigen::Index first
<< m_computed(firstColm + i+1, firstColm+i+1) << " "
<< m_computed(firstColm + i+2, firstColm+i+2) << "\n";
#endif
if (r==Literal(0))
if (numext::is_exactly_zero(r))
{
m_computed(firstColm + i, firstColm + i) = m_computed(firstColm + j, firstColm + j);
return;

6
Eigen/src/SVD/JacobiSVD.h
View File

@@ -377,7 +377,7 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
const RealScalar precision = NumTraits<Scalar>::epsilon();
if(n==0)
if(numext::is_exactly_zero(n))
{
// make sure first column is zero
work_matrix.coeffRef(p,p) = work_matrix.coeffRef(q,p) = Scalar(0);
@@ -684,7 +684,7 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
m_info = InvalidInput;
return *this;
}
if(scale==RealScalar(0)) scale = RealScalar(1);
if(numext::is_exactly_zero(scale)) scale = RealScalar(1);
/*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
@@ -777,7 +777,7 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
{
Index pos;
RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
if(maxRemainingSingularValue == RealScalar(0))
if(numext::is_exactly_zero(maxRemainingSingularValue))
{
m_nonzeroSingularValues = i;
break;

6
Eigen/src/SparseCore/TriangularSolver.h
View File

@@ -116,7 +116,7 @@ struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Lower,ColMajor>
for(Index i=0; i<lhs.cols(); ++i)
{
Scalar& tmp = other.coeffRef(i,col);
if (tmp!=Scalar(0)) // optimization when other is actually sparse
if (!numext::is_exactly_zero(tmp)) // optimization when other is actually sparse
{
LhsIterator it(lhsEval, i);
while(it && it.index()<i)
@@ -151,7 +151,7 @@ struct sparse_solve_triangular_selector<Lhs,Rhs,Mode,Upper,ColMajor>
for(Index i=lhs.cols()-1; i>=0; --i)
{
Scalar& tmp = other.coeffRef(i,col);
if (tmp!=Scalar(0)) // optimization when other is actually sparse
if (!numext::is_exactly_zero(tmp)) // optimization when other is actually sparse
{
if(!(Mode & UnitDiag))
{
@@ -241,7 +241,7 @@ struct sparse_solve_triangular_sparse_selector<Lhs,Rhs,Mode,UpLo,ColMajor>
{
tempVector.restart();
Scalar& ci = tempVector.coeffRef(i);
if (ci!=Scalar(0))
if (!numext::is_exactly_zero(ci))
{
// find
typename Lhs::InnerIterator it(lhs, i);