Skip to content

Commit f25bc64

Browse files
ranjananranjanan
committed
Some reorg
1 parent 1a39841 commit f25bc64

File tree

5 files changed

+103
-126
lines changed

5 files changed

+103
-126
lines changed

‎src/AMG.jl

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,12 @@
11
module AMG
22

3+
include("strength.jl")
4+
export classical
5+
36
include("mg.jl")
47
export RS
58

9+
include("gallery.jl")
10+
export poisson
11+
612
end # module

‎src/gallery.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
poisson(n) = sparse(Tridiagonal(fill(-1, n-1), fill(2, n), fill(-1, n-1)))

‎src/mg.jl renamed to ‎src/splitting.jl

Lines changed: 26 additions & 126 deletions
Original file line numberDiff line numberDiff line change
@@ -1,104 +1,5 @@
1-
#=function mg(a::SparseMatrixCSC, b::Vector, x = zeros(size(b)))
2-
3-
r = a*x - b
4-
5-
# Smoothing steps
6-
smooth!(x, A, b)
7-
8-
# Coarsening steps
9-
while level < nlevels
10-
Ac[level] = coarsen(A[level-1])
11-
f[level] = coarse(f[level-1])
12-
end
13-
x = solve(A[end], f[end])
14-
15-
# Interpolation step
16-
interpolate!(x, nlevels)
17-
18-
if norm(r) < tol
19-
break
20-
end
21-
22-
x
23-
end
24-
25-
struct Level
26-
A::SparseMatrixCSC
27-
R::RestrictionMatrix
28-
P::ProlongationMatrix
29-
end
30-
31-
struct RestrictionMatrix
32-
r::SparseMatrixCSC
33-
end
34-
35-
struct ProlongationMatrix
36-
p::SparseMatrixCSC
37-
end=#
38-
39-
function classical(A::SparseMatrixCSC, θ::Float64)
40-
41-
I = Int[]
42-
J = Int[]
43-
V = Float64[]
44-
45-
m, n = size(A)
46-
47-
for i = 1:n
48-
neighbors = A[:,i]
49-
m = find_max_off_diag(neighbors, i)
50-
threshold = θ * m
51-
for j in nzrange(A, i)
52-
row = A.rowval[j]
53-
val = A.nzval[j]
54-
if abs(val) >= threshold
55-
push!(I, row)
56-
push!(J, i)
57-
push!(V, abs(val))
58-
end
59-
end
60-
end
61-
S = sparse(I, J, V)
62-
63-
scale_cols_by_largest_entry(S)
64-
end
65-
66-
function find_max_off_diag(neighbors, col)
67-
max_offdiag = 0
68-
for (i,v) in enumerate(neighbors)
69-
if col != i
70-
max_offdiag = max(max_offdiag, abs(v))
71-
end
72-
end
73-
max_offdiag
74-
end
75-
76-
function scale_cols_by_largest_entry(A::SparseMatrixCSC)
77-
78-
m,n = size(A)
79-
80-
I = zeros(Int, size(A.nzval))
81-
J = similar(I)
82-
V = zeros(size(A.nzval))
83-
84-
k = 1
85-
for i = 1:n
86-
m = maximum(A[:,i])
87-
for j in nzrange(A, i)
88-
row = A.rowval[j]
89-
val = A.nzval[j]
90-
I[k] = row
91-
J[k] = i
92-
V[k] = val / m
93-
k += 1
94-
end
95-
end
96-
97-
sparse(I,J,V)
98-
end
99-
1001
function RS(S::SparseMatrixCSC)
101-
2+
1023
m,n = size(S)
1034

1045
n_nodes = n
@@ -107,17 +8,17 @@ function RS(S::SparseMatrixCSC)
1078
Tj = S.rowval
1089
Sp = Tp
10910
Sj = Tj
110-
11+
11112
# compute lambdas
11213
for i = 1:n
11314
lambda[i] = Tp[i+1] - Tp[i]
11415
end
115-
16+
11617
interval_ptr = zeros(Int, n+1)
11718
interval_count = zeros(Int, n+1)
11819
index_to_node = zeros(Int,n)
11920
node_to_index = zeros(Int,n)
120-
21+
12122
for i = 1:n
12223
interval_count[lambda[i]+1] += 1
12324
end
@@ -135,96 +36,96 @@ function RS(S::SparseMatrixCSC)
13536
interval_count[lambda_i] += 1
13637
end
13738
splitting = fill(2, n)
138-
39+
13940
# all nodes with no neighbors become F nodes
14041
for i = 1:n
14142
# check if diagonal or check if no neighbors
14243
if lambda[i] == 0 || (lambda[i] == 1 && Tj[Tp[i]] == i)
14344
splitting[i] = F_NODE
14445
end
14546
end
146-
47+
14748
# Now add elements to C and F, in descending order of lambda
14849
for top_index = n_nodes:-1:1
14950
i = index_to_node[top_index]
15051
lambda_i = lambda[i] + 1
151-
52+
15253
# if (n_nodes == 4)
15354
# std::cout << "selecting node #" << i << " with lambda " << lambda[i] << std::endl;
154-
55+
15556
# remove i from its interval
15657
interval_count[lambda_i] -= 1
157-
58+
15859
if splitting[i] == F_NODE
15960
continue
16061
else
161-
62+
16263
@assert splitting[i] == U_NODE
163-
64+
16465
splitting[i] = C_NODE
165-
66+
16667
# For each j in S^T_i /\ U
16768
for jj = Tp[i]:Tp[i+1]-1
16869

16970
j = Tj[jj]
170-
71+
17172
if splitting[j] == U_NODE
17273
splitting[j] = F_NODE
173-
74+
17475
# For each k in S_j /\ U
17576
for kk = Sp[j]: Sp[j+1]-1
17677
k = Sj[kk]
177-
78+
17879
if splitting[k] == U_NODE
17980
# move k to the end of its current interval
18081
lambda[k] >= n_nodes - 1 && continue
181-
82+
18283
lambda_k = lambda[k] + 1
18384
old_pos = node_to_index[k]
18485
new_pos = interval_ptr[lambda_k] + interval_count[lambda_k]# - 1
185-
86+
18687
node_to_index[index_to_node[old_pos]] = new_pos
18788
node_to_index[index_to_node[new_pos]] = old_pos
18889
(index_to_node[old_pos], index_to_node[new_pos]) = (index_to_node[new_pos], index_to_node[old_pos])
189-
90+
19091
# update intervals
19192
interval_count[lambda_k] -= 1
19293
interval_count[lambda_k+1] += 1 # invalid write!
19394
interval_ptr[lambda_k+1] = new_pos - 1
194-
95+
19596
# increment lambda_k
19697
lambda[k] += 1
19798
end
19899
end
199100
end
200101
end
201-
102+
202103
# For each j in S_i /\ U
203104
for jj = Sp[i]: Sp[i+1]-1
204105

205106
j = Sj[jj]
206107

207108
if splitting[j] == U_NODE # decrement lambda for node j
208109

209-
lambda[j] == 0 && continue
210-
110+
lambda[j] == 0 && continue
111+
211112
# assert(lambda[j] > 0);//this would cause problems!
212-
113+
213114
# move j to the beginning of its current interval
214115
lambda_j = lambda[j] + 1
215116
old_pos = node_to_index[j]
216117
new_pos = interval_ptr[lambda_j]
217-
118+
218119
node_to_index[index_to_node[old_pos]] = new_pos
219120
node_to_index[index_to_node[new_pos]] = old_pos
220121
(index_to_node[old_pos],index_to_node[new_pos]) = (index_to_node[new_pos],index_to_node[old_pos])
221-
122+
222123
# update intervals
223124
interval_count[lambda_j] -= 1
224125
interval_count[lambda_j-1] += 1
225126
interval_ptr[lambda_j] += 1
226127
interval_ptr[lambda_j-1] = interval_ptr[lambda_j] - interval_count[lambda_j-1]
227-
128+
228129
# decrement lambda_j
229130
lambda[j] -= 1
230131
end
@@ -233,4 +134,3 @@ function RS(S::SparseMatrixCSC)
233134
end
234135
splitting
235136
end
236-
poisson(n) = sparse(Tridiagonal(fill(-1, n-1), fill(2, n), fill(-1, n-1)))

‎src/strength.jl

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
function classical(A::SparseMatrixCSC, θ::Float64)
2+
3+
I = Int[]
4+
J = Int[]
5+
V = Float64[]
6+
7+
m, n = size(A)
8+
9+
for i = 1:n
10+
neighbors = A[:,i]
11+
m = find_max_off_diag(neighbors, i)
12+
threshold = θ * m
13+
for j in nzrange(A, i)
14+
row = A.rowval[j]
15+
val = A.nzval[j]
16+
if abs(val) >= threshold
17+
push!(I, row)
18+
push!(J, i)
19+
push!(V, abs(val))
20+
end
21+
end
22+
end
23+
S = sparse(I, J, V)
24+
25+
scale_cols_by_largest_entry(S)
26+
end
27+
28+
function find_max_off_diag(neighbors, col)
29+
max_offdiag = 0
30+
for (i,v) in enumerate(neighbors)
31+
if col != i
32+
max_offdiag = max(max_offdiag, abs(v))
33+
end
34+
end
35+
max_offdiag
36+
end
37+
38+
function scale_cols_by_largest_entry(A::SparseMatrixCSC)
39+
40+
m,n = size(A)
41+
42+
I = zeros(Int, size(A.nzval))
43+
J = similar(I)
44+
V = zeros(size(A.nzval))
45+
46+
k = 1
47+
for i = 1:n
48+
m = maximum(A[:,i])
49+
for j in nzrange(A, i)
50+
row = A.rowval[j]
51+
val = A.nzval[j]
52+
I[k] = row
53+
J[k] = i
54+
V[k] = val / m
55+
k += 1
56+
end
57+
end
58+
59+
sparse(I,J,V)
60+
end

‎test/thing.txt

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
0 0 0 0 0 0 0 0 0.7670696322682211 0
2+
0 0 0 0 0.838117753907359 0 0 0.1196525672223625 0 0
3+
0 0 0 0 0.9147120238969264 0 0 0 0 0
4+
0 0 0 0.13574455851185352 0 0 0 0 0 0
5+
0 0 0 0 0 0 0 0 0.8019235854122897 0
6+
0 0 0 0 0 0 0 0 0.035344549147287685 0
7+
0 0 0 0 0 0 0 0 0 0
8+
0 0 0.6052967398293401 0 0 0 0 0 0 0
9+
0 0 0 0 0 0 0.7228497594213787 0 0.48466052213279887 0
10+
0 0 0 0 0 0.30007495800798534 0 0 0 0

0 commit comments

Comments
 (0)