non-allocating ishermitian_approx

Author: lkdvos

src/common/matrixproperties.jl

@@ -69,31 +69,32 @@ is_right_isometric(A; kwargs...) = is_left_isometric(A'; kwargs...)
 Test whether a linear map is Hermitian, i.e. `A = A'`.
 The `isapprox_kwargs` can be used to control the tolerances of the equality.
 """
-function ishermitian(A; atol::Real = 0, rtol::Real = 0, norm = LinearAlgebra.norm, kwargs...)
+function ishermitian(A; atol::Real = 0, rtol::Real = 0, kwargs...)
     if iszero(atol) && iszero(rtol)
         return ishermitian_exact(A; kwargs...)
     else
-        return 2 * norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
+        return ishermitian_approx(A; atol, rtol, kwargs...)
     end
 end
-function ishermitian_exact(A)
-    return A == A'
-end
-function ishermitian_exact(A::StridedMatrix; kwargs...)
-    return strided_ishermitian_exact(A, Val(false); kwargs...)
+
+ishermitian_exact(A) = A == A'
+ishermitian_exact(A::StridedMatrix; kwargs...) = strided_ishermitian_exact(A, Val(false); kwargs...)
+function ishermitian_approx(A; atol, rtol, kwargs...)
+    return 2 * norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
 end
+ishermitian_approx(A::StridedMatrix; kwargs...) = strided_ishermitian_approx(A, Val(false); kwargs...)
 
 """
     isantihermitian(A; isapprox_kwargs...)
 
 Test whether a linear map is anti-Hermitian, i.e. `A = -A'`.
 The `isapprox_kwargs` can be used to control the tolerances of the equality.
 """
-function isantihermitian(A; atol::Real = 0, rtol::Real = 0, norm = LinearAlgebra.norm, kwargs...)
+function isantihermitian(A; atol::Real = 0, rtol::Real = 0, kwargs...)
     if iszero(atol) && iszero(rtol)
         return isantihermitian_exact(A; kwargs...)
     else
-        return 2 * norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
+        return isantihermitian_approx(A; atol, rtol, kwargs...)
     end
 end
 function isantihermitian_exact(A)
@@ -102,6 +103,10 @@ end
 function isantihermitian_exact(A::StridedMatrix; kwargs...)
     return strided_ishermitian_exact(A, Val(true); kwargs...)
 end
+function isantihermitian_approx(A; atol, rtol, kwargs...)
+    return 2 * norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
+end
+isantihermitian_approx(A::StridedMatrix; kwargs...) = strided_ishermitian_approx(A, Val(true); kwargs...)
 
 # blocked implementation of exact checks for strided matrices
 # -----------------------------------------------------------
@@ -139,3 +144,51 @@ function _ishermitian_exact_offdiag(Al, Au, ::Val{anti}) where {anti}
     end
     return true
 end
+
+
+function strided_ishermitian_approx(
+        A::AbstractMatrix, anti::Val;
+        blocksize = 32, atol::Real = default_hermitian_tol(A), rtol::Real = 0
+    )
+    n = size(A, 1)
+    ϵ = abs2(zero(eltype(A)))
+    ϵmax = oftype(ϵ, rtol > 0 ? max(atol, rtol * norm(A)) : atol)^2
+    for j in 1:blocksize:n
+        jb = min(blocksize, n - j + 1)
+        ϵ += _ishermitian_approx_diag(view(A, j:(j + jb - 1), j:(j + jb - 1)), anti)
+        ϵ < ϵmax || return false
+        for i in 1:blocksize:(j - 1)
+            ib = blocksize
+            ϵ += _ishermitian_approx_offdiag(
+                view(A, i:(i + ib - 1), j:(j + jb - 1)),
+                view(A, j:(j + jb - 1), i:(i + ib - 1)),
+                anti
+            )
+            ϵ < ϵmax || return false
+        end
+    end
+    return true
+end
+
+function _ishermitian_approx_diag(A, ::Val{anti}) where {anti}
+    n = size(A, 1)
+    ϵ = abs2(zero(eltype(A)))
+    @inbounds for j in 1:n
+        @simd for i in 1:j
+            val = anti ? (A[i, j] + adjoint(A[j, i])) : (A[i, j] - adjoint(A[j, i]))
+            ϵ += abs2(val)
+        end
+    end
+    return ϵ
+end
+function _ishermitian_approx_offdiag(Al, Au, ::Val{anti}) where {anti}
+    m, n = size(Al) # == reverse(size(Al))
+    ϵ = abs2(zero(eltype(Al)))
+    @inbounds for j in 1:n
+        @simd for i in 1:m
+            val = anti ? (Al[i, j] + adjoint(Au[j, i])) : (Al[i, j] - adjoint(Au[j, i]))
+            ϵ += abs2(val)
+        end
+    end
+    return ϵ
+end