-
Notifications
You must be signed in to change notification settings - Fork 4
/
matrices.jl
85 lines (70 loc) · 2.4 KB
/
matrices.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
const TransposeOrAdjoint{T,M} = Union{Transpose{T,M},Adjoint{T,M}}
"""
matrix_versions(A::AbstractMatrix)
Return various versions of the same matrix:
- dense and sparse
- transpose and adjoint
Used for internal testing.
"""
function matrix_versions(A)
A_dense = Matrix(A)
A_sparse = sparse(A)
versions = [
A_dense,
transpose(Matrix(transpose(A_dense))),
adjoint(Matrix(adjoint(A_dense))),
A_sparse,
transpose(sparse(transpose(A_sparse))),
adjoint(sparse(adjoint(A_sparse))),
]
if issymmetric(A)
lower_triangles = [
Matrix(LowerTriangular(A_dense)), sparse(LowerTriangular(A_sparse))
]
upper_triangles = [
Matrix(UpperTriangular(A_dense)), sparse(UpperTriangular(A_sparse))
]
symmetric_versions = vcat(
Symmetric.(versions),
Hermitian.(versions),
Symmetric.(lower_triangles, :L),
Symmetric.(upper_triangles, :U),
)
append!(versions, symmetric_versions)
end
return versions
end
"""
respectful_similar(A::AbstractMatrix)
respectful_similar(A::AbstractMatrix, ::Type{T})
Like `Base.similar` but returns a transpose or adjoint when `A` is a transpose or adjoint.
"""
respectful_similar(A::AbstractMatrix) = respectful_similar(A, eltype(A))
respectful_similar(A::AbstractMatrix, ::Type{T}) where {T} = similar(A, T)
function respectful_similar(A::Transpose, ::Type{T}) where {T}
return transpose(respectful_similar(parent(A), T))
end
function respectful_similar(A::Adjoint, ::Type{T}) where {T}
return adjoint(respectful_similar(parent(A), T))
end
function respectful_similar(A::Union{Symmetric,Hermitian}, ::Type{T}) where {T}
return respectful_similar(sparse(A), T)
end
"""
same_pattern(A, B)
Perform a partial equality check on the sparsity patterns of `A` and `B`:
- if the return is `true`, they might have the same sparsity pattern but we're not sure
- if the return is `false`, they definitely don't have the same sparsity pattern
"""
same_pattern(A, B) = size(A) == size(B)
function same_pattern(
A::Union{SparseMatrixCSC,SparsityPatternCSC},
B::Union{SparseMatrixCSC,SparsityPatternCSC},
)
return size(A) == size(B) && nnz(A) == nnz(B)
end
function check_same_pattern(A, S)
if !same_pattern(A, S)
throw(DimensionMismatch("`A` and `S` must have the same sparsity pattern."))
end
end