Skip to content

Commit cb6d13f

Browse files
committed
Add gmg as the test case and fix bug in finding coarse levels
1 parent 0341389 commit cb6d13f

File tree

2 files changed

+51
-1
lines changed

2 files changed

+51
-1
lines changed

src/classical.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,7 @@ function ruge_stuben{Ti,Tv}(A::SparseMatrixCSC{Ti,Tv};
2020

2121
levels = Vector{Level{Ti,Tv}}()
2222

23-
while length(levels) + 1 < max_levels && size(A, 1) < max_coarse
23+
while length(levels) + 1 < max_levels && size(A, 1) > max_coarse
2424
A = extend_heirarchy!(levels, strength, CF, A)
2525
#if size(A, 1) <= max_coarse
2626
# break

test/gmg.jl

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
import AMG: Level, MultiLevel, GaussSeidel
2+
3+
function multigrid(A::SparseMatrixCSC{T,V}; max_levels = 10, max_coarse = 10,
4+
presmoother = GaussSeidel(), postsmoother = GaussSeidel()) where {T,V}
5+
6+
levels = Vector{Level{T,V}}()
7+
8+
while length(levels) + 1 < max_levels && size(A, 1) > max_coarse
9+
A = extend!(levels, A)
10+
#=if size(A, 1) <= max_coarse
11+
# push!(levels, Level(A,spzeros(T,0,0),spzeros(T,0,0)))
12+
break
13+
end=#
14+
end
15+
16+
MultiLevel(levels, A, presmoother, postsmoother)
17+
end
18+
19+
function extend!(levels, A::SparseMatrixCSC{Ti,Tv}) where {Ti,Tv}
20+
21+
size_F = size(A, 1)
22+
size_C = rem(size_F,2) == 0 ? div((size_F-1), 2) +1 : div((size_F-1), 2)
23+
total_len = size_C + 2 * (size_C - 1)
24+
25+
I = Vector{Tv}(total_len)
26+
J = Vector{Tv}(total_len)
27+
V = Vector{Ti}(total_len)
28+
29+
l = 1
30+
for k = 1:size_C
31+
I[l], J[l], V[l] = 2k, k, 1
32+
l += 1
33+
end
34+
for k = 1:size_C - 1
35+
I[l], J[l], V[l] = 2k+1, k, 0.5
36+
I[l+1], J[l+1], V[l+1] = 2k+1, k+1, 0.5
37+
l += 2
38+
end
39+
40+
P = sparse(I, J, V, size_F, size_C)
41+
42+
R = transpose(P)
43+
44+
push!(levels, Level(A, P, R))
45+
46+
R * A * P
47+
end
48+
49+
Base.show(io::IO, level::Level) = print(io, "Level")
50+

0 commit comments

Comments
 (0)