Skip to content

Commit 1b996ee

Browse files
j-fuj-fu
committed
Merge branch 'jkrch-master'
2 parents 0a8bff8 + 60d92df commit 1b996ee

File tree

5 files changed

+238
-1
lines changed

5 files changed

+238
-1
lines changed

‎src/ExtendableSparse.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@ include("extendable.jl")
1515
export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, updateindex!, rawupdateindex!, colptrs
1616

1717
include("factorizations.jl")
18-
export JacobiPreconditioner, ILU0Preconditioner, ParallelJacobiPreconditioner
18+
export JacobiPreconditioner, ILU0Preconditioner, ParallelJacobiPreconditioner, ParallelILU0Preconditioner, reorderlinsys
1919
export AbstractFactorization,LUFactorization, CholeskyFactorization
2020
export issolver
2121
export factorize!, update!

‎src/factorizations.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -127,6 +127,7 @@ end
127127
include("jacobi.jl")
128128
include("ilu0.jl")
129129
include("parallel_jacobi.jl")
130+
include("parallel_ilu0.jl")
130131
include("umfpack_lu.jl")
131132
include("cholmod_cholesky.jl")
132133

@@ -137,6 +138,7 @@ include("cholmod_cholesky.jl")
137138
@makefrommatrix ILU0Preconditioner
138139
@makefrommatrix JacobiPreconditioner
139140
@makefrommatrix ParallelJacobiPreconditioner
141+
@makefrommatrix ParallelILU0Preconditioner
140142

141143
end
142144

‎src/parallel_ilu0.jl‎

Lines changed: 210 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,210 @@
1+
mutable struct ParallelILU0Preconditioner{Tv, Ti} <: AbstractPreconditioner{Tv,Ti}
2+
A::ExtendableSparseMatrix{Tv,Ti}
3+
xdiag::Array{Tv,1}
4+
idiag::Array{Ti,1}
5+
phash::UInt64
6+
7+
coloring::Array{Array{Ti,1},1}
8+
coloring_index::Array{Array{Ti,1},1}
9+
coloring_index_reverse::Array{Array{Ti,1},1}
10+
11+
function ParallelILU0Preconditioner{Tv,Ti}() where {Tv,Ti}
12+
p=new()
13+
p.phash=0
14+
p
15+
end
16+
end
17+
18+
"""
19+
```
20+
ParallelILU0Preconditioner(;valuetype=Float64,indextype=Int64)
21+
ParallelILU0Preconditioner(matrix)
22+
```
23+
24+
Parallel ILU preconditioner with zero fill-in.
25+
"""
26+
ParallelILU0Preconditioner(;valuetype::Type=Float64, indextype::Type=Int64)=ParallelILU0Preconditioner{valuetype,indextype}()
27+
28+
29+
function update!(precon::ParallelILU0Preconditioner{Tv,Ti}) where {Tv,Ti}
30+
flush!(precon.A)
31+
32+
# Get coloring and reorder matrix
33+
precon.coloring=graphcol(precon.A.cscmatrix)
34+
precon.coloring_index, precon.coloring_index_reverse=coloringindex(precon.coloring)
35+
precon.A=ExtendableSparseMatrix(reordermatrix(precon.A.cscmatrix, precon.coloring))
36+
37+
cscmatrix=precon.A.cscmatrix
38+
colptr=cscmatrix.colptr
39+
rowval=cscmatrix.rowval
40+
nzval=cscmatrix.nzval
41+
n=cscmatrix.n
42+
43+
if precon.phash==0
44+
n=size(precon.A,1)
45+
precon.xdiag=Array{Tv,1}(undef,n)
46+
precon.idiag=Array{Ti,1}(undef,n)
47+
end
48+
49+
xdiag=precon.xdiag
50+
idiag=precon.idiag
51+
52+
53+
# Find main diagonal index and
54+
# copy main diagonal values
55+
if precon.phash != precon.A.phash
56+
@inbounds for j=1:n
57+
@inbounds for k=colptr[j]:colptr[j+1]-1
58+
i=rowval[k]
59+
if i==j
60+
idiag[j]=k
61+
break
62+
end
63+
end
64+
end
65+
precon.phash=precon.A.phash
66+
end
67+
68+
@inbounds for j=1:n
69+
xdiag[j]=one(Tv)/nzval[idiag[j]]
70+
@inbounds for k=idiag[j]+1:colptr[j+1]-1
71+
i=rowval[k]
72+
for l=colptr[i]:colptr[i+1]-1
73+
if rowval[l]==j
74+
xdiag[i]-=nzval[l]*xdiag[j]*nzval[k]
75+
break
76+
end
77+
end
78+
end
79+
end
80+
precon
81+
end
82+
83+
84+
function LinearAlgebra.ldiv!(u::AbstractArray{T,1}, precon::ParallelILU0Preconditioner, v::AbstractArray{T,1}) where T
85+
cscmatrix=precon.A.cscmatrix
86+
colptr=cscmatrix.colptr
87+
rowval=cscmatrix.rowval
88+
n=cscmatrix.n
89+
nzval=cscmatrix.nzval
90+
xdiag=precon.xdiag
91+
idiag=precon.idiag
92+
93+
coloring = precon.coloring
94+
coloring_index = precon.coloring_index
95+
coloring_index_reverse = precon.coloring_index_reverse
96+
97+
@inbounds for indset in coloring_index
98+
@inbounds Threads.@threads for j in indset
99+
x=zero(T)
100+
@inbounds for k=colptr[j]:idiag[j]-1
101+
x+=nzval[k]*u[rowval[k]]
102+
end
103+
u[j]=xdiag[j]*(v[j]-x)
104+
end
105+
end
106+
107+
@inbounds for indset in coloring_index_reverse
108+
@inbounds Threads.@threads for j in indset
109+
x=zero(T)
110+
@inbounds for k=idiag[j]+1:colptr[j+1]-1
111+
x+=u[rowval[k]]*nzval[k]
112+
end
113+
u[j]-=x*xdiag[j]
114+
end
115+
end
116+
end
117+
118+
119+
function LinearAlgebra.ldiv!(precon::ParallelILU0Preconditioner, v::AbstractArray{T,1} where T)
120+
ldiv!(v, precon, v)
121+
end
122+
123+
124+
# Returns an independent set of the graph of a matrix
125+
# Reference: https://research.nvidia.com/sites/default/files/pubs/2015-05_Parallel-Graph-Coloring/nvr-2015-001.pdf
126+
function indset(A::SparseMatrixCSC{Tv,Ti}, W::StridedVector{Ti}) where {Tv,Ti}
127+
# Random numbers for all vertices
128+
lenW = length(W)
129+
# r = sample(1:lenW, lenW, replace = false)
130+
r = rand(lenW)
131+
@inbounds for i = 1:lenW
132+
if W[i] == 0
133+
r[i] = 0
134+
end
135+
end
136+
# Empty independent set
137+
S = zeros(Int, lenW)
138+
# Get independent set by comparing random number of vertex with the random
139+
# numbers of all neighbor vertices
140+
@inbounds Threads.@threads for i in 1:lenW
141+
if W[i] != 0
142+
j = A.rowval[A.colptr[i]:A.colptr[i+1]-1]
143+
if all(x->x==1, r[i] .>= r[j])
144+
S[i] = W[i]
145+
end
146+
end
147+
end
148+
# Remove zero entries and return independent set
149+
return filter!(x->x≠0, S)
150+
end
151+
152+
# Returns coloring of the graph of a matrix
153+
# Reference: https://research.nvidia.com/sites/default/files/pubs/2015-05_Parallel-Graph-Coloring/nvr-2015-001.pdf
154+
function graphcol(A::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
155+
# Empty list for coloring
156+
C = []
157+
# Array of vertices
158+
W = [1:size(A)[1];]
159+
# Get all independent sets of the graph of the matrix
160+
while any(W .!= 0)
161+
# Get independent set
162+
S = indset(A + transpose(A), W)
163+
push!(C, S)
164+
# Remove entries in S from W
165+
@inbounds for s in S
166+
W[s] = 0
167+
end
168+
end
169+
# Return coloring
170+
return C
171+
end
172+
173+
174+
# Reorders a sparse matrix with provided coloring
175+
function reordermatrix(A::SparseMatrixCSC{Tv,Ti}, coloring::Array{Array{Int64,1},1}) where {Tv,Ti}
176+
c = collect(Iterators.flatten(coloring))
177+
return A[c,:][:,c]
178+
end
179+
180+
# Reorders a linear system with provided coloring
181+
function reorderlinsys(A::SparseMatrixCSC{Tv,Ti}, b::StridedVector{Tv}, coloring::Array{Array{Int64,1},1}) where {Tv,Ti}
182+
c = collect(Iterators.flatten(coloring))
183+
return A[c,:][:,c], b[c]
184+
end
185+
186+
187+
# Returns an array with the same structure of the input coloring and ordered
188+
# entries 1:length(coloring) and an array with the structure of
189+
# reverse(coloring) and ordered entries length(coloring):-1:1
190+
function coloringindex(coloring::Array{Array{Int64,1},1})
191+
# First array
192+
c = deepcopy(coloring)
193+
cnt = 1
194+
@inbounds for i in 1:length(c)
195+
@inbounds for j in 1:length(c[i])
196+
c[i][j] = cnt
197+
cnt += 1
198+
end
199+
end
200+
# Second array
201+
cc = deepcopy(reverse(coloring))
202+
@inbounds for i in 1:length(cc)
203+
@inbounds for j in 1:length(cc[i])
204+
cnt -= 1
205+
cc[i][j] = cnt
206+
end
207+
end
208+
# Return both
209+
return c, cc
210+
end

‎test/Project.toml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
AlgebraicMultigrid = "2169fc97-5a83-5252-b627-83903c6c433c"
33
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
44
IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
5+
IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
56
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
67
Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
78
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

‎test/runtests.jl‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ using BenchmarkTools
88
using Pardiso
99
using AlgebraicMultigrid
1010
using IncompleteLU
11+
using IterativeSolvers
1112

1213
##############################################################
1314
@testset "Constructors" begin
@@ -394,3 +395,26 @@ if Pardiso.PARDISO_LOADED[]
394395
end
395396
end
396397

398+
399+
400+
##############################################
401+
function test_parilu0(n)
402+
A = ExtendableSparseMatrix(n, n)
403+
sprand_sdd!(A)
404+
flush!(A)
405+
A = A.cscmatrix
406+
b = A * ones(n)
407+
P_par = ParallelILU0Preconditioner(A)
408+
A_reord, b_reord = reorderlinsys(A, b, P_par.coloring)
409+
P_ser = ILU0Preconditioner(A_reord)
410+
sol_ser, hist_ser = gmres(A_reord, b_reord, Pl=P_ser, log=true)
411+
sol_par, hist_par = gmres(A_reord, b_reord, Pl=P_par, log=true)
412+
sol_ser == sol_par && hist_ser.iters == hist_par.iters
413+
end
414+
415+
@testset "Parallel ILU0" begin
416+
@test test_parilu0(10)
417+
@test test_parilu0(100)
418+
@test test_parilu0(1000)
419+
end
420+

0 commit comments

Comments
 (0)