refactor: delegate eigenvalue transforms to the DiagonalAlgorithm kernels

In the `MatrixFunctionViaEig(h)` implementations of `squareroot`,
`logarithm` and `power`, apply the scalar function to the eigenvalue
matrix by calling the function's own `DiagonalAlgorithm` kernel (which
already contains the domain checks and tolerance clamping), mirroring
how `exponential!` recurses into `exponential!((τ, D), D)`. The eig
kernels only retain the negative-real-axis scan of the complex
eigenvalues for real input, which has no `Diagonal` counterpart.

Author: lkdvos

src/implementations/logarithm.jl

@@ -48,31 +48,23 @@ end
 function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEigh)
     check_input(logarithm!, A, logA, alg)
     D, V = eigh_full!(A, alg.eigh_alg)
-    λ = diagview(D)
-    atol = something(alg.domain_atol, default_domain_atol(λ))
-    _check_nonzero_eigenvalues(λ, atol)
-    _clamp_domain_eigenvalues!(λ, atol)
-    λ .= log.(λ)
-    VD = V * D
+    VD = V * logarithm!(D, D, DiagonalAlgorithm(; domain_atol = alg.domain_atol))
     mul!(logA, VD, V')
     return project_hermitian!(logA)
 end
 
 function logarithm!(A::AbstractMatrix, logA, alg::MatrixFunctionViaEig)
     check_input(logarithm!, A, logA, alg)
     D, V = eig_full!(A, alg.eig_alg)
-    λ = diagview(D)
-    atol = something(alg.domain_atol, default_domain_atol(λ))
-    _check_nonzero_eigenvalues(λ, atol)
+    diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol)
     if eltype(A) <: Real
-        _clamp_domain_eigenvalues!(λ, atol)
-        λ .= log.(λ)
-        VD = V * D
-        logAc = rdiv!(VD, LinearAlgebra.lu!(V))
+        atol = something(alg.domain_atol, default_domain_atol(diagview(D)))
+        _clamp_domain_eigenvalues!(diagview(D), atol)
+        VlogD = V * logarithm!(D, D, diag_alg)
+        logAc = rdiv!(VlogD, LinearAlgebra.lu!(V))
         return logA .= real.(logAc)
     else
-        λ .= log.(λ)
-        logA .= V .* transpose(λ)
+        logA .= V .* transpose(diagview(logarithm!(D, D, diag_alg)))
         return rdiv!(logA, LinearAlgebra.lu!(V))
     end
 end

src/implementations/power.jl

@@ -54,43 +54,32 @@ end
 function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEigh)
     check_input(power!, A, p, powA, alg)
     D, V = eigh_full!(A, alg.eigh_alg)
-    λ = diagview(D)
+    diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol)
     if isinteger(p)
-        p < 0 && any(iszero, λ) && throw(LinearAlgebra.SingularException(0))
-        λ .= λ .^ p
-        VD = V * D
+        VD = V * power!(D, p, D, diag_alg)
         mul!(powA, VD, V')
         return project_hermitian!(powA)
     else
-        atol = something(alg.domain_atol, default_domain_atol(λ))
-        p < 0 && _check_nonzero_eigenvalues(λ, atol)
-        _clamp_domain_eigenvalues!(λ, atol)
         # `A^p = (V * D^(p/2)) * (V * D^(p/2))'` is hermitian by construction
-        λ .= λ .^ (p / 2)
-        Vs = rmul!(V, D)
+        Vs = rmul!(V, power!(D, p / 2, D, diag_alg))
         return _mul_herm!(powA, Vs)
     end
 end
 
 function power!(A::AbstractMatrix, p::Real, powA, alg::MatrixFunctionViaEig)
     check_input(power!, A, p, powA, alg)
     D, V = eig_full!(A, alg.eig_alg)
-    λ = diagview(D)
-    if isinteger(p)
-        p < 0 && any(iszero, λ) && throw(LinearAlgebra.SingularException(0))
-    else
-        atol = something(alg.domain_atol, default_domain_atol(λ))
-        p < 0 && _check_nonzero_eigenvalues(λ, atol)
-        eltype(A) <: Real && _clamp_domain_eigenvalues!(λ, atol)
-    end
+    diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol)
     if eltype(A) <: Real
-        λ .= λ .^ p
-        VD = V * D
-        powAc = rdiv!(VD, LinearAlgebra.lu!(V))
+        if !isinteger(p)
+            atol = something(alg.domain_atol, default_domain_atol(diagview(D)))
+            _clamp_domain_eigenvalues!(diagview(D), atol)
+        end
+        VpD = V * power!(D, p, D, diag_alg)
+        powAc = rdiv!(VpD, LinearAlgebra.lu!(V))
         return powA .= real.(powAc)
     else
-        λ .= λ .^ p
-        powA .= V .* transpose(λ)
+        powA .= V .* transpose(diagview(power!(D, p, D, diag_alg)))
         return rdiv!(powA, LinearAlgebra.lu!(V))
     end
 end

src/implementations/squareroot.jl

@@ -48,29 +48,25 @@ end
 function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEigh)
     check_input(squareroot!, A, sqrtA, alg)
     D, V = eigh_full!(A, alg.eigh_alg)
-    λ = diagview(D)
-    atol = something(alg.domain_atol, default_domain_atol(λ))
-    _clamp_domain_eigenvalues!(λ, atol)
+    diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol)
     # `sqrt(A) = (V * D^(1/4)) * (V * D^(1/4))'` is hermitian by construction
-    λ .= sqrt.(sqrt.(λ))
-    Vs = rmul!(V, D)
+    sqrtD = squareroot!(D, D, diag_alg)
+    Vs = rmul!(V, squareroot!(sqrtD, sqrtD, diag_alg))
     return _mul_herm!(sqrtA, Vs)
 end
 
 function squareroot!(A::AbstractMatrix, sqrtA, alg::MatrixFunctionViaEig)
     check_input(squareroot!, A, sqrtA, alg)
     D, V = eig_full!(A, alg.eig_alg)
-    λ = diagview(D)
+    diag_alg = DiagonalAlgorithm(; domain_atol = alg.domain_atol)
     if eltype(A) <: Real
-        atol = something(alg.domain_atol, default_domain_atol(λ))
-        _clamp_domain_eigenvalues!(λ, atol)
-        λ .= sqrt.(λ)
-        VD = V * D
-        sqrtAc = rdiv!(VD, LinearAlgebra.lu!(V))
+        atol = something(alg.domain_atol, default_domain_atol(diagview(D)))
+        _clamp_domain_eigenvalues!(diagview(D), atol)
+        VsqrtD = V * squareroot!(D, D, diag_alg)
+        sqrtAc = rdiv!(VsqrtD, LinearAlgebra.lu!(V))
         return sqrtA .= real.(sqrtAc)
     else
-        λ .= sqrt.(λ)
-        sqrtA .= V .* transpose(λ)
+        sqrtA .= V .* transpose(diagview(squareroot!(D, D, diag_alg)))
         return rdiv!(sqrtA, LinearAlgebra.lu!(V))
     end
 end