Skip to content

Commit f4ea67c

Browse files
committed
improve matrixproperties
1 parent 3924196 commit f4ea67c

1 file changed

Lines changed: 18 additions & 8 deletions

File tree

src/common/matrixproperties.jl

Lines changed: 18 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,10 @@ function isunitary(A; isapprox_kwargs...)
3131
return is_left_isometry(A; isapprox_kwargs...) &&
3232
is_right_isometry(A; isapprox_kwargs...)
3333
end
34-
# TODO: for matrices, isunitary corresponds to is_left_isometry + square.
34+
function isunitary(A::AbstractMatrix; isapprox_kwargs...)
35+
size(A, 1) == size(A, 2) || return false
36+
return is_left_isometry(A; isapprox_kwargs...)
37+
end
3538

3639
@doc """
3740
is_left_isometry(A; isapprox_kwargs...) -> Bool
@@ -42,8 +45,11 @@ The `isapprox_kwargs` can be used to control the tolerances of the equality.
4245
See also [`isisometry`](@ref) and [`is_right_isometry`](@ref).
4346
""" is_left_isometry
4447

45-
function is_left_isometry(A::AbstractMatrix; isapprox_kwargs...)
46-
return isapprox(A' * A, LinearAlgebra.I; isapprox_kwargs...)
48+
function is_left_isometry(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
49+
P = A' * A
50+
nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
51+
diagview(P) .-= 1
52+
return norm(P) <= max(atol, rtol * nP) # assume that the norm of I is `sqrt(n)`
4753
end
4854

4955
@doc """
@@ -55,8 +61,11 @@ The `isapprox_kwargs` can be used to control the tolerances of the equality.
5561
See also [`isisometry`](@ref) and [`is_left_isometry`](@ref).
5662
""" is_right_isometry
5763

58-
function is_right_isometry(A::AbstractMatrix; isapprox_kwargs...)
59-
return isapprox(A * A', LinearAlgebra.I; isapprox_kwargs...)
64+
function is_right_isometry(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
65+
P = A * A'
66+
nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
67+
diagview(P) .-= 1
68+
return norm(P) <= max(atol, rtol * nP) # assume that the norm of I is `sqrt(n)`
6069
end
6170

6271
"""
@@ -75,7 +84,7 @@ end
7584
function ishermitian_exact(A)
7685
return A == A'
7786
end
78-
function ishermitian_exact(A::AbstractMatrix; kwargs...)
87+
function ishermitian_exact(A::StridedMatrix; kwargs...)
7988
return _ishermitian_exact(A, Val(false); kwargs...)
8089
end
8190

@@ -95,11 +104,12 @@ end
95104
function isantihermitian_exact(A)
96105
return A == -A'
97106
end
98-
function isantihermitian_exact(A::AbstractMatrix; kwargs...)
107+
function isantihermitian_exact(A::StridedMatrix; kwargs...)
99108
return _ishermitian_exact(A, Val(true); kwargs...)
100109
end
101110

102-
# block implementation of exact checks
111+
# blocked implementation of exact checks for strided matrices
112+
# -----------------------------------------------------------
103113
function _ishermitian_exact(A::AbstractMatrix, anti::Val; blocksize = 32)
104114
n = size(A, 1)
105115
for j in 1:blocksize:n

0 commit comments

Comments
 (0)