Skip to content

Commit 3ad76ba

Browse files
MatrixIdentification.jl
1 parent 42299dd commit 3ad76ba

File tree

1 file changed

+95
-0
lines changed

1 file changed

+95
-0
lines changed
Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
using LinearAlgebra
2+
3+
## compute the percentage banded for a matrix given a bandwidth
4+
function compute_bandedness(A, bandwidth)
5+
n = size(A, 1)
6+
total_band_positions = 0
7+
non_zero_in_band = 0
8+
bandwidth = bandwidth
9+
for r in 1:n
10+
for c in 1:n
11+
if abs(r - c) <= bandwidth
12+
total_band_positions += 1 # This position belongs to the band
13+
if A[r, c] != 0
14+
non_zero_in_band += 1 # This element is non-zero in the band
15+
end
16+
end
17+
end
18+
end
19+
20+
percentage_filled = non_zero_in_band / total_band_positions * 100
21+
return percentage_filled
22+
end
23+
24+
## use the compute banded function to compute the bandedness given different bandsizes
25+
function isbanded(A)
26+
n = size(A, 1)
27+
count = 5
28+
potential_band_sizes = []
29+
30+
while count < n - 5:
31+
push!(potential_band_sizes, count += 5)
32+
end
33+
34+
for bandsize in potential_band_sizes
35+
percentage_filled = compute_bandedness(A, bandsize)
36+
if percentage_filled >= 0.9
37+
return true
38+
end
39+
end
40+
41+
return false
42+
43+
## compute the sparsity for a given matrix
44+
function compute_sparsity(A)
45+
n = size(A, 1)
46+
percentage_sparsity = length(nonzeros(A)) / n^2
47+
return percentage_sparsity
48+
end
49+
50+
## get the percentage banded for a bandwidth of 5 and percentage sparsity
51+
function getstructure(A::Matrix)::Any
52+
percentage_banded = compute_bandedness(A, 5)
53+
percentage_sparsity = compute_sparsity(A)
54+
55+
return (percentage_banded, percentage_sparsity)
56+
end
57+
58+
## return the best type of matrix for a given sparse matrix
59+
function sparsestructure(A::SparseMatrixCSC)::Any
60+
sym = issymmetric(A)
61+
herm = ishermitian(A)
62+
banded = isbanded(A)
63+
posdef = isposdef(A)
64+
lower_triangular = istril(A)
65+
upper_triangular = istriu(A)
66+
67+
68+
n = size(A, 1)
69+
70+
if banded
71+
return BandedMatrix(A)
72+
end
73+
74+
if sym
75+
return Symmetric(A)
76+
end
77+
78+
if herm
79+
return Hermitian(A)
80+
end
81+
82+
if lower_triangular
83+
return LowerTriangular(A)
84+
end
85+
86+
if upper_triangular
87+
return UpperTriangular(A)
88+
end
89+
end
90+
91+
92+
93+
94+
95+

0 commit comments

Comments
 (0)