tolerance in Float64

Author: lkdvos

src/implementations/matrixfunctions.jl

@@ -71,29 +71,26 @@ end
     )
 end
 
+_clamp_domain_eigenvalues!(D::Diagonal, atol::Real) =
+    _clamp_domain_eigenvalues!(diagview(D), atol)
+
 # Clamp real eigenvalues that are negative within `atol` (rounding artifacts) to zero,
 # and throw a `DomainError` for eigenvalues that are genuinely negative, since then the
 # result cannot be expressed with the same (real) scalar type.
 function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Real}, atol::Real)
     λmin = minimum(λ; init = zero(eltype(λ)))
+    atol = atol < 0 ? default_domain_atol(λ) : oftype(λmin, atol)
     λmin < -atol && throw_negative_eigenvalue(λmin, atol, "a negative real eigenvalue")
     λ .= max.(λ, zero(eltype(λ)))
     return λ
 end
 
-# Convenience method for the eigenvalues of a decomposition, deriving the default
-# tolerance from the eigenvalues themselves when `domain_atol` is `nothing`.
-function _clamp_domain_eigenvalues!(D::Diagonal, domain_atol::Union{Nothing, Real})
-    λ = diagview(D)
-    atol = something(domain_atol, default_domain_atol(λ))
-    return _clamp_domain_eigenvalues!(λ, atol)
-end
-
 # Complex eigenvalues of a real matrix: only eigenvalues (numerically) on the negative
 # real axis obstruct a real result; complex-conjugate pairs do not.
 function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Complex}, atol::Real)
     onaxis = x -> abs(imag(x)) <= atol && real(x) < 0
     λmin = mapreduce(x -> onaxis(x) ? real(x) : zero(real(x)), min, λ; init = zero(real(eltype(λ))))
+    atol = atol < 0 ? default_domain_atol(λ) : oftype(λmin, atol)
     λmin < -atol && throw_negative_eigenvalue(λmin, atol, "an eigenvalue on the negative real axis")
     λ .= ifelse.(onaxis.(λ), zero(eltype(λ)), λ)
     return λ

src/interface/matrixfunctions.jl

@@ -38,17 +38,16 @@ For matrix functions with a restricted domain (e.g. [`squareroot`](@ref) and [`l
 as rounding artifacts and clamped to the domain boundary, with `nothing` denoting the default
 tolerance [`default_domain_atol`](@ref).
 """
-struct MatrixFunctionViaEigh{A <: AbstractAlgorithm, T <: Union{Nothing, Real}} <: AbstractAlgorithm
+struct MatrixFunctionViaEigh{A <: AbstractAlgorithm} <: AbstractAlgorithm
     eigh_alg::A
-    domain_atol::T
-end
-function MatrixFunctionViaEigh(eigh_alg::AbstractAlgorithm; domain_atol::Union{Nothing, Real} = nothing)
-    return MatrixFunctionViaEigh(eigh_alg, domain_atol)
+    domain_atol::Float64 # negative value for runtime defaults
 end
+MatrixFunctionViaEigh(eigh_alg::AbstractAlgorithm; domain_atol::Real = -1.0) =
+    MatrixFunctionViaEigh(eigh_alg, Float64(domain_atol))
 function Base.show(io::IO, alg::MatrixFunctionViaEigh)
     print(io, "MatrixFunctionViaEigh(")
     _show_alg(io, alg.eigh_alg)
-    isnothing(alg.domain_atol) || print(io, "; domain_atol=", alg.domain_atol)
+    alg.domain_atol < 0 || print(io, "; domain_atol=", alg.domain_atol)
     return print(io, ")")
 end
 
@@ -62,16 +61,15 @@ For matrix functions with a restricted domain (e.g. [`squareroot`](@ref) and [`l
 as rounding artifacts and clamped to the domain boundary, with `nothing` denoting the default
 tolerance [`default_domain_atol`](@ref).
 """
-struct MatrixFunctionViaEig{A <: AbstractAlgorithm, T <: Union{Nothing, Real}} <: AbstractAlgorithm
+struct MatrixFunctionViaEig{A <: AbstractAlgorithm} <: AbstractAlgorithm
     eig_alg::A
-    domain_atol::T
-end
-function MatrixFunctionViaEig(eig_alg::AbstractAlgorithm; domain_atol::Union{Nothing, Real} = nothing)
-    return MatrixFunctionViaEig(eig_alg, domain_atol)
+    domain_atol::Float64 # negative value for runtime defaults
 end
+MatrixFunctionViaEig(eig_alg::AbstractAlgorithm; domain_atol::Real = -1.0) =
+    MatrixFunctionViaEig(eig_alg, Float64(domain_atol))
 function Base.show(io::IO, alg::MatrixFunctionViaEig)
     print(io, "MatrixFunctionViaEig(")
     _show_alg(io, alg.eig_alg)
-    isnothing(alg.domain_atol) || print(io, "; domain_atol=", alg.domain_atol)
+    alg.domain_atol < 0 || print(io, "; domain_atol=", alg.domain_atol)
     return print(io, ")")
 end