Skip to content

Commit 6be6f6a

Browse files
j-fuj-fu
committed
First version working with VoronoiFVM.jl
1 parent 30fe88e commit 6be6f6a

File tree

5 files changed

+495
-176
lines changed

5 files changed

+495
-176
lines changed

‎README.md‎

Lines changed: 4 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,9 @@
11
# ExtendableSparse
22

3-
Sparse matrix class which allows cheaper assembly by using
4-
a different data structure for extension.
3+
Sparse matrix class with efficient assembly.
4+
5+
While SparseMatrixCSC is extendable as well, extending it is expensive for larger problems,
6+
due to the implementation of extension without an intermediate data structure.
57

6-
Still work in progress.
78

89

‎examples/ExtendableSparseMatrixTest.jl‎

Lines changed: 32 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,6 @@ function randmatx!(A::AbstractSparseMatrix{Tv,Ti}, m::Ti, n::Ti,xnnz::Ti) where
1313
a=1.0+rand(Float64)
1414
A[i,j]+=a
1515
end
16-
1716
end
1817

1918
function testmatx!(A::AbstractSparseMatrix{Tv,Ti}, m::Ti, n::Ti, xnnz::Ti) where {Tv,Ti}
@@ -28,6 +27,7 @@ function testmatx!(A::AbstractSparseMatrix{Tv,Ti}, m::Ti, n::Ti, xnnz::Ti) where
2827
@printf("nnz(S)=%d nnz(A)=%d\n",nnz(S),nnz(A))
2928
@assert(nnz(S)==nnz(A))
3029
end
30+
flush!(A)
3131
@printf("nnz(S)=%d nnz(A)=%d\n",nnz(S),nnz(A))
3232
@assert(nnz(S)==nnz(A))
3333
(I,J,V)=findnz(S)
@@ -37,23 +37,28 @@ function testmatx!(A::AbstractSparseMatrix{Tv,Ti}, m::Ti, n::Ti, xnnz::Ti) where
3737
end
3838

3939
function check(;m=10,n=10,nnz=3)
40-
extmat1=SparseMatrixExtension{Float64,Int64}(m,n)
40+
extmat1=SparseMatrixExtension(Float64,Int64,m,n)
4141
testmatx!(extmat1,m,n,nnz)
42+
43+
extmat2=ExtendableSparseMatrix(Float64,Int64,m,n)
44+
testmatx!(extmat2,m,n,nnz)
4245
end
4346

4447

4548
function benchmark(;n=100,m=100,nnz=500, pyplot=false)
4649

47-
mat=spzeros(Float64,Int64,n,n)
50+
mat=spzeros(Float64,Int64,m,n)
4851
@time randmatx!(mat,m,n,nnz)
4952

50-
extmat=SparseMatrixExtension{Float64,Int64}(n,n)
53+
extmat=SparseMatrixExtension(Float64,Int64,m,n)
5154
@time randmatx!(extmat,m,n,nnz)
5255

5356

54-
xextmat=ExtendableSparseMatrixCSC(n,n,spzeros(Float64,Int64,n,n),SparseMatrixExtension{Float64,Int64}(n,n))
55-
@time randmatx!(xextmat,m,n,nnz)
56-
57+
xextmat=ExtendableSparseMatrix(Float64,Int64,m,n)
58+
@time begin
59+
randmatx!(xextmat,m,n,nnz)
60+
@inbounds flush!(xextmat)
61+
end
5762

5863

5964

@@ -67,4 +72,24 @@ function benchmark(;n=100,m=100,nnz=500, pyplot=false)
6772
# end
6873
end
6974

75+
function tsplice(;m=10, n=10, nnz=20)
76+
xnnz=nnz
77+
78+
A=ExtendableSparseMatrix(Float64,Int64,m,n)
79+
for inz=1:xnnz
80+
i=rand((1:m))
81+
j=rand((1:n))
82+
a=1.0+rand(Float64)
83+
A[i,j]+=a
84+
end
85+
flush!(A)
86+
for inz=1:2
87+
i=rand((1:m))
88+
j=rand((1:n))
89+
a=1.0+rand(Float64)
90+
A[i,j]+=a
91+
end
92+
flush!(A)
93+
end
94+
7095
end

‎src/ExtendableSparse.jl‎

Lines changed: 6 additions & 166 deletions
Original file line numberDiff line numberDiff line change
@@ -1,172 +1,12 @@
11
module ExtendableSparse
2+
using DocStringExtensions
23
using SparseArrays
4+
using LinearAlgebra
35

4-
mutable struct SparseMatrixExtension{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
5-
m::Ti
6-
n::Ti
7-
nnz::Ti
8-
colptr::Vector{Ti}
9-
rowval::Vector{Ti}
10-
nzval::Vector{Tv}
11-
12-
function SparseMatrixExtension{Tv,Ti}(m,n) where {Tv,Ti<:Integer}
13-
colptr=zeros(Ti,n)
14-
rowval=zeros(Ti,n)
15-
nzval=zeros(Tv,n)
16-
new(m,n,0,zeros(Ti,n),zeros(Ti,n),zeros(Tv,n))
17-
end
18-
end
6+
include("extension.jl")
7+
include("extendable.jl")
198

20-
function Base.setindex!(E::SparseMatrixExtension{Tv,Ti}, _v, _i::Integer, _j::Integer) where {Tv,Ti<:Integer}
21-
v = convert(Tv, _v)
22-
i = convert(Ti, _i)
23-
j = convert(Ti, _j)
9+
export SparseMatrixExtension,ExtendableSparseMatrix,flush!,nnz
2410

25-
if !((1 <= i <= E.m) & (1 <= j <= E.n))
26-
throw(BoundsError(E, (i,j)))
27-
end
28-
29-
if iszero(v)
30-
return E
31-
end
32-
33-
if E.rowval[j]==0
34-
E.rowval[j]=i
35-
E.nzval[j]=v
36-
E.nnz+=1
37-
return E
38-
end
39-
40-
k=j
41-
k0=j
42-
while k>0
43-
if E.rowval[k]==i
44-
E.nzval[k]=v
45-
return E
46-
end
47-
k0=k
48-
k=E.colptr[k]
49-
end
50-
push!(E.nzval,v)
51-
push!(E.rowval,i)
52-
push!(E.colptr,-1)
53-
E.colptr[k0]=length(E.nzval)
54-
E.nnz+=1
55-
return E
56-
end
57-
58-
function Base.getindex(E::SparseMatrixExtension{Tv,Ti},i::Integer, j::Integer) where {Tv,Ti<:Integer}
59-
if !((1 <= i <= E.m) & (1 <= j <= E.n))
60-
throw(BoundsError(E, (i,j)))
61-
end
62-
63-
k=j
64-
while k>0
65-
if E.rowval[k]==i
66-
return E.nzval[k]
67-
end
68-
k=E.colptr[k]
69-
end
70-
return zero(Tv)
71-
end
72-
73-
Base.size(E::SparseMatrixExtension) = (E.m, E.n)
74-
SparseArrays.nnz(E::SparseMatrixExtension)=E.nnz
75-
76-
77-
78-
#####################################################################################################
79-
mutable struct ExtendableSparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
80-
m::Int
81-
n::Int
82-
cscmatrix::SparseMatrixCSC{Tv,Ti}
83-
extmatrix::SparseMatrixExtension{Tv,Ti}
84-
end
85-
86-
function Base.setindex!(m::ExtendableSparseMatrixCSC, v, i::Integer, j::Integer)
87-
setindex!(m.cscmatrix,v,i,j)
88-
end
89-
90-
function Base.getindex(m::ExtendableSparseMatrixCSC,i::Integer, j::Integer)
91-
return getindex(m.cscmatrix,i,j)
92-
end
93-
94-
struct ColEntry{Tv,Ti<:Integer}
95-
j::Ti
96-
v::Tv
97-
end
98-
99-
isless(x::ColEntry,y::ColEntry)=(x.i<y.i)
100-
101-
function splice(E::ExtendableSparseMatrixCSC,S::SparseMatrixCSC)
102-
103-
104-
nnz=nnz(S)+nnz(E)
105-
colptr=Vector{Ti}(undef,S.m+1)
106-
rowval=Vector{Ti}(undef,nnz)
107-
nzval=Vector{Tv}(undef,nnz)
108-
109-
E_maxcol_ext=0
110-
S_maxcol=0
111-
for j=1:m
112-
lrow=0
113-
k=j
114-
while k>0
115-
lrow+=1
116-
k=E.colptr[k]
117-
end
118-
E_maxcol_ext=max(lrow,E_maxcol)
119-
S_maxcol=max(S.colptr[j+1]-S.colptr[j],S_maxcol)
120-
end
121-
122-
col=Vector{ColEntry}(undef,E_maxcol+S_Maxcol+10)
123-
maxcol=0
124-
125-
i=1
126-
for j=1:m
127-
# put extension entries into row and sort them
128-
k=j
129-
while k>0
130-
if colptr[k]>0
131-
lxcol+=1
132-
col[lxcol].j=colptr[k]
133-
col[lxcol].v=nzval[k]
134-
k=colptr[k]
135-
end
136-
end
137-
138-
sort!(col,lt=isless)
139-
# jointly sort old and mew entries into colptr
140-
141-
i0=i
142-
colptr[j]=i
143-
jcol=0
144-
k=S.colptr[j]
145-
while true
146-
if k<S.colptr[j+1] && icol>lxcol || S.colptr[k]<col[jcol].j
147-
rowval[i]=S.rowval[k]
148-
nzval[i]=S.nzval[k]
149-
k+=1
150-
i+=1
151-
continue
152-
end
153-
if jcol<lxcol
154-
rowval[i]=col[jcol].j
155-
nzval[i]=col[jcol].v
156-
jcol+=1
157-
i++
158-
continue
159-
end
160-
break
161-
end
162-
end
163-
maxrow=max(maxrow,i-i0)
164-
colptr[j]=i
165-
return SparseMatrixCSC(m,n,colptr,rowval,nzval)
166-
end
167-
168-
169-
170-
171-
export SparseMatrixExtension,ExtendableSparseMatrixCSC
11+
export xcolptrs,colptrs
17212
end # module

0 commit comments

Comments
 (0)