refactor: fold the eigenvalue domain scan into a single helper call

Add a `_clamp_domain_eigenvalues!(D::Diagonal, domain_atol)` convenience
method that derives the default tolerance itself, so the real-input
branch of the `MatrixFunctionViaEig` kernels reduces to a one-line
preamble before delegating to the `DiagonalAlgorithm` kernel.

Author: lkdvos

src/implementations/logarithm.jl

@@ -58,8 +58,7 @@ function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEig)
     D, V = eig_full!(A, alg.eig_alg)
     diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol)
     if eltype(A) <: Real
-        atol = something(alg.domain_atol, default_domain_atol(diagview(D)))
-        _clamp_domain_eigenvalues!(diagview(D), atol)
+        _clamp_domain_eigenvalues!(D, alg.domain_atol)
         VlogD = V * logarithm!(D, D, diag_alg)
         logAc = rdiv!(VlogD, LinearAlgebra.lu!(V))
         return logA .= real.(logAc)

src/implementations/matrixfunctions.jl

@@ -22,6 +22,14 @@ function _clamp_domain_eigenvalues!(λ::AbstractVector{<:Real}, atol::Real)
     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)

src/implementations/power.jl

@@ -71,10 +71,7 @@ function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEig)
     D, V = eig_full!(A, alg.eig_alg)
     diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol)
     if eltype(A) <: Real
-        if !isinteger(p)
-            atol = something(alg.domain_atol, default_domain_atol(diagview(D)))
-            _clamp_domain_eigenvalues!(diagview(D), atol)
-        end
+        isinteger(p) || _clamp_domain_eigenvalues!(D, alg.domain_atol)
         VpD = V * power!(D, p, D, diag_alg)
         powAc = rdiv!(VpD, LinearAlgebra.lu!(V))
         return powA .= real.(powAc)

src/implementations/squareroot.jl

@@ -60,8 +60,7 @@ function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEig)
     D, V = eig_full!(A, alg.eig_alg)
     diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol)
     if eltype(A) <: Real
-        atol = something(alg.domain_atol, default_domain_atol(diagview(D)))
-        _clamp_domain_eigenvalues!(diagview(D), atol)
+        _clamp_domain_eigenvalues!(D, alg.domain_atol)
         VsqrtD = V * squareroot!(D, D, diag_alg)
         sqrtAc = rdiv!(VsqrtD, LinearAlgebra.lu!(V))
         return sqrtA .= real.(sqrtAc)