rework tests

Author: lkdvos

src/common/matrixproperties.jl

@@ -80,7 +80,7 @@ end
 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))
+    return norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
 end
 ishermitian_approx(A::StridedMatrix; kwargs...) = strided_ishermitian_approx(A, Val(false); kwargs...)
 
@@ -104,7 +104,7 @@ 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))
+    return norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
 end
 isantihermitian_approx(A::StridedMatrix; kwargs...) = strided_ishermitian_approx(A, Val(true); kwargs...)
 
@@ -159,7 +159,7 @@ function strided_ishermitian_approx(
         ϵ² < ϵ²max || return false
         for i in 1:blocksize:(j - 1)
             ib = blocksize
-            ϵ² += _ishermitian_approx_offdiag(
+            ϵ² += 2 * _ishermitian_approx_offdiag(
                 view(A, i:(i + ib - 1), j:(j + jb - 1)),
                 view(A, j:(j + jb - 1), i:(i + ib - 1)),
                 anti
@@ -175,7 +175,7 @@ function _ishermitian_approx_diag(A, ::Val{anti}) where {anti}
     ϵ² = 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]))
+            val = _project_hermitian(A[i, j], A[j, i], !anti)
             ϵ² += abs2(val) * (1 + Int(i < j))
         end
     end
@@ -186,9 +186,9 @@ function _ishermitian_approx_offdiag(Al, Au, ::Val{anti}) where {anti}
     ϵ² = 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]))
+            val = _project_hermitian(Al[i, j], Au[j, i], !anti)
             ϵ² += abs2(val)
         end
     end
-    return 2ϵ²
+    return ϵ²
 end

src/implementations/projections.jl

@@ -85,14 +85,16 @@ function project_hermitian_native!(A::AbstractMatrix, B::AbstractMatrix, anti::V
     return B
 end
 
+@inline function _project_hermitian(Aij::Number, Aji::Number, anti::Bool)
+    return anti ? (Aij - Aji') / 2 : (Aij + Aji') / 2
+end
 function _project_hermitian_offdiag!(
         Au::AbstractMatrix, Al::AbstractMatrix, Bu::AbstractMatrix, Bl::AbstractMatrix, ::Val{anti}
     ) where {anti}
-
     m, n = size(Au) # == reverse(size(Au))
     return @inbounds for j in 1:n
         @simd for i in 1:m
-            val = anti ? (Au[i, j] - adjoint(Al[j, i])) / 2 : (Au[i, j] + adjoint(Al[j, i])) / 2
+            val = _project_hermitian(Au[i, j], Al[j, i], anti)
             Bu[i, j] = val
             aval = adjoint(val)
             Bl[j, i] = anti ? -aval : aval
@@ -104,7 +106,7 @@ function _project_hermitian_diag!(A::AbstractMatrix, B::AbstractMatrix, ::Val{an
     n = size(A, 1)
     @inbounds for j in 1:n
         @simd for i in 1:(j - 1)
-            val = anti ? (A[i, j] - adjoint(A[j, i])) / 2 : (A[i, j] + adjoint(A[j, i])) / 2
+            val = _project_hermitian(A[i, j], A[j, i], anti)
             B[i, j] = val
             aval = adjoint(val)
             B[j, i] = anti ? -aval : aval

test/projections.jl

@@ -2,7 +2,7 @@ using MatrixAlgebraKit
 using Test
 using TestExtras
 using StableRNGs
-using LinearAlgebra: LinearAlgebra, Diagonal, norm
+using LinearAlgebra: LinearAlgebra, Diagonal, norm, normalize!
 
 const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
 
@@ -43,6 +43,20 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
         @test isantihermitian(Ba)
         @test Ba ≈ Aa
     end
+
+    # test approximate error calculation
+    A = normalize!(randn(rng, T, m, m))
+    Ah = project_hermitian(A)
+    Ah_approx = Ah + noisefactor * Aa
+    ϵ = norm(project_antihermitian(Ah_approx))
+    @test !ishermitian(Ah_approx; atol = (999 // 1000) * ϵ)
+    @test ishermitian(Ah_approx; atol = (1001 // 1000) * ϵ)
+
+    Aa = project_antihermitian(A)
+    Aa_approx = Aa + noisefactor * Ah
+    ϵ = norm(project_hermitian(Aa_approx))
+    @test !isantihermitian(Aa_approx; atol = (999 // 1000) * ϵ)
+    @test isantihermitian(Aa_approx; atol = (1001 // 1000) * ϵ)
 end
 
 @testset "project_isometric! for T = $T" for T in BLASFloats