refactor: make domain-clamp helpers GPU compatible

Replace the scalar-indexing loops in `_clamp_domain_eigenvalues!` and
`_check_nonzero_eigenvalues` with reductions and broadcasts, so the
helpers work on GPU arrays. The `DomainError` now reports the worst
offending eigenvalue instead of the first one encountered.

Author: lkdvos

src/implementations/matrixfunctions.jl

@@ -5,20 +5,17 @@
 # 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)
-    for i in eachindex(λ)
-        x = λ[i]
-        if x < -atol
-            throw(
-                DomainError(
-                    x,
-                    "The matrix has a negative real eigenvalue beyond `domain_atol = $atol` and the result of this matrix function is complex. " *
-                        "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact."
-                )
+    λmin = minimum(λ; init = zero(eltype(λ)))
+    if λmin < -atol
+        throw(
+            DomainError(
+                λmin,
+                "The matrix has a negative real eigenvalue beyond `domain_atol = $atol` and the result of this matrix function is complex. " *
+                    "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact."
             )
-        elseif x < 0
-            λ[i] = zero(x)
-        end
+        )
     end
+    λ .= max.(λ, zero(eltype(λ)))
     return λ
 end
 
@@ -33,37 +30,32 @@ 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)
-    for i in eachindex(λ)
-        x = λ[i]
-        if abs(imag(x)) <= atol && real(x) < 0
-            if real(x) < -atol
-                throw(
-                    DomainError(
-                        x,
-                        "The matrix has an eigenvalue on the negative real axis beyond `domain_atol = $atol` and the result of this matrix function is complex. " *
-                            "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact."
-                    )
-                )
-            else
-                λ[i] = zero(x)
-            end
-        end
+    onaxis = x -> abs(imag(x)) <= atol && real(x) < 0
+    λmin = mapreduce(x -> onaxis(x) ? real(x) : zero(real(x)), min, λ; init = zero(real(eltype(λ))))
+    if λmin < -atol
+        throw(
+            DomainError(
+                λmin,
+                "The matrix has an eigenvalue on the negative real axis beyond `domain_atol = $atol` and the result of this matrix function is complex. " *
+                    "Pass a complex matrix to obtain the principal value, or increase `domain_atol` if the eigenvalue is a rounding artifact."
+            )
+        )
     end
+    λ .= ifelse.(onaxis.(λ), zero(eltype(λ)), λ)
     return λ
 end
 
 # Reject (numerically) zero eigenvalues for functions that are undefined there,
 # e.g. `logarithm` and `power` with a negative fractional power.
 function _check_nonzero_eigenvalues(λ, atol::Real)
-    for x in λ
-        if abs(x) <= atol
-            throw(
-                DomainError(
-                    x,
-                    "The matrix has a (numerically) zero eigenvalue within `domain_atol = $atol`, for which this matrix function is not defined."
-                )
+    amin = minimum(abs, λ; init = typemax(real(eltype(λ))))
+    if amin <= atol
+        throw(
+            DomainError(
+                amin,
+                "The matrix has a (numerically) zero eigenvalue within `domain_atol = $atol`, for which this matrix function is not defined."
             )
-        end
+        )
     end
     return λ
 end