Skip to content

Commit 5a0c548

Browse files
ranjananranjanan
authored
Merge pull request #1 from ranjanan/level
Construct Interpolation Operator
2 parents 5e25911 + 9ac252c commit 5a0c548

File tree

6 files changed

+216
-13
lines changed

6 files changed

+216
-13
lines changed

‎src/AMG.jl

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

33
include("strength.jl")
4-
export classical
4+
export strength_of_connection, Classical
55

66
include("splitting.jl")
7-
export RS
7+
export split_nodes, RS
88

99
include("gallery.jl")
1010
export poisson
1111

12+
include("classical.jl")
13+
1214
end # module

‎src/classical.jl

Lines changed: 166 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,166 @@
1+
struct Solver{S,T,P,PS}
2+
strength::S
3+
CF::T
4+
presmoother::P
5+
postsmoother::PS
6+
max_levels::Int64
7+
max_coarse::Int64
8+
end
9+
10+
struct Level{T}
11+
A::T
12+
end
13+
14+
function ruge_stuben(A::SparseMatrixCSC;
15+
strength = Classical(),
16+
CF = RS(),
17+
presmoother = GaussSiedel(),
18+
postsmoother = GaussSiedel(),
19+
max_levels = 10,
20+
max_coarse = 500)
21+
22+
s = Solver(strength, CF, presmoother,
23+
postsmoother, max_levels, max_levels)
24+
25+
levels = [Level(A)]
26+
27+
while length(levels) < max_levels && size(levels[end].A, 1)
28+
extend_heirarchy!(levels, strength, CF, A)
29+
end
30+
end
31+
32+
function extend_heirarchy!(levels::Vector{Level}, strength, CF, A)
33+
S = strength_of_connection(strength, A)
34+
splitting = split_nodes(CF, S)
35+
P, R = direct_interpolation(A, S, splitting)
36+
end
37+
38+
function direct_interpolation{T,V}(A::T, S::T, splitting::Vector{V})
39+
40+
fill!(S.nzval, 1.)
41+
S = A .* S
42+
Pp = rs_direct_interpolation_pass1(S, A, splitting)
43+
Pp .= Pp .+ 1
44+
45+
Px, Pj = rs_direct_interpolation_pass2(A, S, splitting, Pp)
46+
47+
Px .= abs.(Px)
48+
Pj .= Pj .+ 1
49+
50+
R = SparseMatrixCSC(maximum(Pj), size(A, 1), Pp, Pj, Px)
51+
P = R'
52+
53+
P, R
54+
end
55+
56+
57+
function rs_direct_interpolation_pass1(S, A, splitting)
58+
59+
Bp = zeros(Int, size(A.colptr))
60+
Sp = S.colptr
61+
Sj = S.rowval
62+
n_nodes = size(A, 1)
63+
nnz = 0
64+
for i = 1:n_nodes
65+
if splitting[i] == C_NODE
66+
nnz += 1
67+
else
68+
for jj = Sp[i]:Sp[i+1]-1
69+
if splitting[Sj[jj]] == C_NODE && Sj[jj] != i
70+
nnz += 1
71+
end
72+
end
73+
end
74+
Bp[i+1] = nnz
75+
end
76+
Bp
77+
end
78+
79+
80+
function rs_direct_interpolation_pass2{Tv, Ti}(A::SparseMatrixCSC{Tv,Ti},
81+
S::SparseMatrixCSC{Tv,Ti},
82+
splitting::Vector{Ti},
83+
Bp::Vector{Ti})
84+
85+
86+
Ap = A.colptr
87+
Aj = A.rowval
88+
Ax = A.nzval
89+
Sp = S.colptr
90+
Sj = S.rowval
91+
Sx = S.nzval
92+
Bj = zeros(Ti, Bp[end])
93+
Bx = zeros(Float64, Bp[end])
94+
n_nodes = size(A, 1)
95+
#Bp += 1
96+
97+
for i = 1:n_nodes
98+
if splitting[i] == C_NODE
99+
Bj[Bp[i]] = i
100+
Bx[Bp[i]] = 1
101+
else
102+
sum_strong_pos = 0
103+
sum_strong_neg = 0
104+
for jj = Sp[i]: Sp[i+1]-1
105+
if splitting[Sj[jj]] == C_NODE && Sj[jj] != i
106+
if Sx[jj] < 0
107+
sum_strong_neg += Sx[jj]
108+
else
109+
sum_strong_pos += Sx[jj]
110+
end
111+
end
112+
end
113+
114+
sum_all_pos = 0
115+
sum_all_neg = 0
116+
diag = 0;
117+
for jj = Ap[i]:Ap[i+1]
118+
if Aj[jj] == i
119+
diag += Ax[jj]
120+
else
121+
if Ax[jj] < 0
122+
sum_all_neg += Ax[jj];
123+
else
124+
sum_all_pos += Ax[jj];
125+
end
126+
end
127+
end
128+
129+
alpha = sum_all_neg / sum_strong_neg
130+
beta = sum_all_pos / sum_strong_pos
131+
132+
if sum_strong_pos == 0
133+
diag += sum_all_pos
134+
beta = 0
135+
end
136+
137+
neg_coeff = -alpha/diag;
138+
pos_coeff = -beta/diag;
139+
140+
nnz = Bp[i]
141+
for jj = Sp[i]:Sp[i+1]
142+
if splitting[Sj[jj]] == C_NODE && Sj[jj] != i
143+
Bj[nnz] = Sj[jj]
144+
if Sx[jj] < 0
145+
Bx[nnz] = neg_coeff * Sx[jj]
146+
else
147+
Bx[nnz] = pos_coeff * Sx[jj]
148+
nnz += 1
149+
end
150+
end
151+
end
152+
end
153+
end
154+
155+
m = zeros(Ti, n_nodes)
156+
sum = 0
157+
for i = 1:n_nodes
158+
m[i] = sum
159+
sum += splitting[i]
160+
end
161+
for i = 1:Bp[n_nodes]
162+
Bj[i] = m[Bj[i]]
163+
end
164+
165+
Bx, Bj, Bp
166+
end

‎src/multilevel.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
struct Level{Ti,Tv}
2+
A::SparseMatrixCSC{Ti,Tv}
3+
end

‎src/splitting.jl

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,11 @@ const F_NODE = 0
22
const C_NODE = 1
33
const U_NODE = 2
44

5-
function RS(S::SparseMatrixCSC)
5+
struct RS
6+
end
7+
8+
split_nodes(::RS, S::SparseMatrixCSC) = RS_CF_splitting(S - spdiagm(diag(S)))
9+
function RS_CF_splitting(S::SparseMatrixCSC)
610

711
m,n = size(S)
812

‎src/strength.jl

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,11 @@
1-
function classical(A::SparseMatrixCSC, θ::Float64)
1+
abstract type Strength end
2+
struct Classical{T} <: Strength
3+
θ::T
4+
end
5+
6+
function strength_of_connection{T}(c::Classical{T}, A::SparseMatrixCSC)
27

8+
θ = c.θ
39
I = Int[]
410
J = Int[]
511
V = Float64[]

‎test/runtests.jl

Lines changed: 31 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,20 +1,42 @@
11
using AMG
22
using Base.Test
33

4+
@testset "Strength of connection" begin
5+
46
# classical strength of connection
57
A = poisson(5)
6-
@test full(classical(A, 0.2)) == [ 1.0 0.5 0.0 0.0 0.0
7-
0.5 1.0 0.5 0.0 0.0
8-
0.0 0.5 1.0 0.5 0.0
9-
0.0 0.0 0.5 1.0 0.5
10-
0.0 0.0 0.0 0.5 1.0 ]
8+
X = strength_of_connection(Classical(0.2), A)
9+
@test full(X) == [ 1.0 0.5 0.0 0.0 0.0
10+
0.5 1.0 0.5 0.0 0.0
11+
0.0 0.5 1.0 0.5 0.0
12+
0.0 0.0 0.5 1.0 0.5
13+
0.0 0.0 0.0 0.5 1.0 ]
14+
15+
end
16+
17+
@testset "Splitting" begin
1118

1219
# Ruge-Stuben splitting
1320
S = poisson(7)
14-
S = S - spdiagm(diag(S)) #FIXME
15-
@test RS(S) == [0, 1, 0, 1, 0, 1, 0]
21+
@test split_nodes(RS(), S) == [0, 1, 0, 1, 0, 1, 0]
1622

1723
srand(0)
1824
S = sprand(10,10,0.1); S = S + S'
19-
S = S - spdiagm(diag(S))
20-
@test RS(S) == [0, 1, 1, 0, 0, 0, 0, 0, 1, 1]
25+
@test split_nodes(RS(), S) == [0, 1, 1, 0, 0, 0, 0, 0, 1, 1]
26+
27+
end
28+
29+
@testset "Interpolation" begin
30+
31+
# Direct Interpolation
32+
using AMG
33+
A = poisson(5)
34+
splitting = [1,0,1,0,1]
35+
P, R = AMG.direct_interpolation(A, A, splitting)
36+
@test P == [ 1.0 0.0 0.0
37+
0.5 0.5 0.0
38+
0.0 1.0 0.0
39+
0.0 0.5 0.5
40+
0.0 0.0 1.0 ]
41+
42+
end

0 commit comments

Comments
 (0)