Skip to content

Commit 49df6b0

Browse files
j-fuj-fu
committed
added updateindex! method
provide fdrand and fdrand! matrix constructors automatic benchmarks in examples
1 parent 57a23f5 commit 49df6b0

File tree

9 files changed

+348
-85
lines changed

9 files changed

+348
-85
lines changed

‎README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@ Any linear algebra method on `ExtendableSparseMatrix` starts with a `flush!` met
1616

1717
`ExtendableSparseMatrix` is aimed to work as a drop-in replacement to `SparseMatrixCSC` in finite element and finite volume codes especally in those cases where the sparsity structure is hard to detect a priori and where working with an intermediadte COO representation appears to be not convenient.
1818

19+
In addition, the packages provides a method `updateindex!(A,op,v,i,j)' for both `SparseMatrixCSC` and for `ExtendableSparse` which allows to update a matrix element with one index search instead of two. It allows to replace e.g. `A[i,j]+=v` by `updateindex!(A,+,i,j,v)`. The former operation is lowered to `%1 = Base.getindex(A, 1, 2); %2 = %1 + 3; Base.setindex!(A, %2, 1, 2)` triggering two index searches, one for `getindex!` and another one for `setindex!`.
1920

2021

2122

‎docs/Project.toml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
[deps]
2+
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
23
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
34
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
45
DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"

‎docs/src/changes.md‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,20 @@
11
# Changes
2+
## v0.2.5, Jan 26, 2020
3+
- fixed allocations in Base.+
4+
- added updateindex! method
5+
- provide fdrand and fdrand! matrix constructors
6+
- automatic benchmarks in examples
7+
8+
## v0.2.4, Jan 19, 2020
9+
- Allow preconditioner creation directly from CSC Matrix
10+
- Rename AbstractPreconditioner to AbstractExtendablePreconditioner
11+
12+
## v0.2.3, Jan 15, 2020
13+
- Started to introduce preconditioners (undocumented)
14+
15+
## v0.2.3, Jan 8, 2020
16+
- added norm, cond, opnorm methods
17+
- resize! instead of push! when adding entries should trigger less allocation operations
218

319
## v0.2.2. Dec 23, 2019
420
- What used to be `_splice` is now `+` and allows now real addition (resulting in a CSC matrix)

‎docs/src/example.md‎

Lines changed: 89 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -1,73 +1,106 @@
11
# Examples
22

33

4-
## Usage idea
4+
## Matrix creation
55
An `ExtendableSparseMatrix` can serve as a drop-in replacement for
6-
`SparseMatrixCSC`, albeit with faster assembly.
7-
6+
`SparseMatrixCSC`, albeit with faster assembly during the buildup
7+
phase when using index based access. That means that code similar
8+
to the following example should be fast enough to avoid the assembly
9+
steps using the coordinate format:
810

9-
```jldoctest
10-
using ExtendableSparse
11-
function example(n,bandwidth)
12-
A=ExtendableSparseMatrix(n,n);
13-
aii=0.0
14-
for i=1:n
15-
for j=max(1,i-bandwidth):5:min(n,i+bandwidth)
16-
aij=-rand()
17-
A[i,j]=aij
18-
aii+=abs(aij)
19-
end
20-
A[i,i]=aii+0.1
21-
end
22-
b=rand(n)
23-
x=A\b
24-
A*x≈ b
11+
```@example
12+
using ExtendableSparse # hide
13+
n=3
14+
A=ExtendableSparseMatrix(n,n)
15+
for i=1:n-1
16+
A[i,i+1]=i
2517
end
26-
example(1000,100)
27-
# output
28-
true
18+
A
2919
```
3020

3121

32-
## Benchmark
33-
The code below provides a small benchmark example.
34-
```julia
35-
using ExtendableSparse
36-
using SparseArrays
22+
### Benchmark
3723

38-
function ext_create(n)
39-
A=ExtendableSparseMatrix(n,n)
40-
sprand_sdd!(A);
41-
end
24+
The method [`fdrand`](@ref) creates a matrix similar to the discetization
25+
matrix of a Poisson equation on a d-dimensional cube. The code uses the index
26+
access API for the creation of the matrix.
27+
This approach is considerably faster with
28+
the [`ExtendableSparseMatrix`](@ref) which uses a linked list based
29+
structure [`SparseMatrixLNK`](@ref) to grab new entries.
4230

43-
function csc_create(n)
44-
A=spzeros(n,n)
45-
sprand_sdd!(A);
46-
end
4731

48-
# Trigger JIT compilation before timing
49-
csc_create(10);
50-
ext_create(10);
51-
52-
n=90000
53-
@time Acsc=csc_create(n);
54-
nnz(Acsc)
55-
@time Aext=ext_create(n);
56-
nnz(Aext)
57-
b=rand(n);
58-
@time Acsc*(Acsc\b)≈b
59-
@time Aext*(Aext\b)≈b
32+
```@example
33+
using ExtendableSparse # hide
34+
using SparseArrays # hide
35+
using BenchmarkTools # hide
36+
37+
@benchmark fdrand(30,30,30, matrixtype=ExtendableSparseMatrix);
38+
```
39+
40+
```@example
41+
using ExtendableSparse # hide
42+
using SparseArrays # hide
43+
using BenchmarkTools # hide
44+
45+
@benchmark fdrand(30,30,30, matrixtype=SparseMatrixCSC);
46+
```
47+
48+
49+
## Matrix update
50+
For repeated calculations on the same sparsity structure (e.g. for time dependent
51+
problems or Newton iterations) it is convenient to skip all but the first creation steps
52+
and to just replace the values in the matrix after setting then elements of the `nzval`
53+
vector to zero. Typically in finite element and finite volume methods this step updates
54+
matrix entries (most of them several times) by adding values. In this case, the current indexing
55+
interface of Julia requires to access the matrix twice:
56+
57+
```@example
58+
using SparseArrays # hide
59+
A=spzeros(3,3)
60+
Meta.@lower A[1,2]+=3
61+
```
62+
For sparse matrices this requires to the index search in the structure twice.
63+
The packages provides the method [`updateindex!`](@ref) for both `SparseMatrixCSC` and
64+
for `ExtendableSparse` which allows to update a matrix element with one index search.
65+
66+
67+
### Benchmark for `SparseMatrixCSC`
68+
```@example
69+
using ExtendableSparse # hide
70+
using SparseArrays # hide
71+
using BenchmarkTools # hide
72+
73+
A=fdrand(30,30,30, matrixtype=SparseMatrixCSC);
74+
@benchmark fdrand!(A,30,30,30, update=(A,v,i,j)-> A[i,j]+=v);
75+
```
76+
```@example
77+
using ExtendableSparse # hide
78+
using SparseArrays # hide
79+
using BenchmarkTools # hide
80+
81+
A=fdrand(30,30,30, matrixtype=SparseMatrixCSC);
82+
@benchmark fdrand!(A,30,30,30, update=(A,v,i,j)-> updateindex!(A,+,v,i,j));
83+
```
84+
85+
### Benchmark for `ExtendableSparseMatrix`
86+
```@example
87+
using ExtendableSparse # hide
88+
using BenchmarkTools # hide
89+
90+
A=fdrand(30,30,30, matrixtype=ExtendableSparseMatrix);
91+
@benchmark fdrand!(A,30,30,30, update=(A,v,i,j)-> A[i,j]+=v);
92+
```
93+
```@example
94+
using ExtendableSparse # hide
95+
using BenchmarkTools # hide
96+
97+
A=fdrand(30,30,30, matrixtype=ExtendableSparseMatrix);
98+
@benchmark fdrand!(A,30,30,30, update=(A,v,i,j)-> updateindex!(A,+,v,i,j));
6099
```
61100

62-
The function [`sprand_sdd!`](@ref) fills a
63-
sparse matrix with random entries such that it becomes strictly diagonally
64-
dominant and thus invertible and has a fixed number of nonzeros in
65-
its rows (similar to the method in the first example). Its bandwidth is bounded by 2*sqrt(n), therefore it
66-
resembles a typical matrix of a 2D piecewise linear FEM discretization.
67-
For filling a matrix a, the method conveniently albeit naively
68-
uses just `a[i,j]=value`. This approach is considerably faster with
69-
the [`ExtendableSparseMatrix`](@ref) which uses a linked list based
70-
structure [`SparseMatrixLNK`](@ref) to grab new entries.
71101

102+
Note that the update process for `ExtendableSparse` is slightly slower
103+
than for `SparseMatrixCSC` due to the overhead which comes from checking
104+
the presence of new entries.
72105

73106

‎src/ExtendableSparse.jl‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -3,13 +3,14 @@ using DocStringExtensions
33
using SparseArrays
44
using LinearAlgebra
55

6+
include("sparsematrixcsc.jl")
67
include("sparsematrixlnk.jl")
78
include("extendable.jl")
89
include("preconditioners.jl")
910
include("sprand.jl")
1011

11-
export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, sprand!,sprand_sdd!
12-
export JacobiPreconditioner, ILU0Preconditioner, update!
12+
export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, sprand!,sprand_sdd!, fdrand,fdrand!
13+
export JacobiPreconditioner, ILU0Preconditioner, updateindex!
1314

1415
export colptrs
1516
end # module

‎src/extendable.jl‎

Lines changed: 30 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -75,27 +75,39 @@ end
7575
"""
7676
$(SIGNATURES)
7777
78-
Return index corresponding to entry [i,j] in the array of nonzeros,
79-
if the entry exists, otherwise, return 0.
80-
"""
81-
function findindex(csc::SparseMatrixCSC{T}, i::Integer, j::Integer) where T
82-
if !(1 <= i <= csc.m && 1 <= j <= csc.n); throw(BoundsError()); end
83-
r1 = Int(csc.colptr[j])
84-
r2 = Int(csc.colptr[j+1]-1)
85-
if r1>r2
86-
return 0
87-
end
88-
89-
# See sparsematrix.jl
90-
r1 = searchsortedfirst(csc.rowval, i, r1, r2, Base.Forward)
91-
if (r1>r2 ||csc.rowval[r1] != i)
92-
return 0
78+
Update element of the matrix with operation `op`.
79+
This can replace the following code and save one index
80+
search during acces:
81+
82+
```@example
83+
using ExtendableSparse # hide
84+
A=ExtendableSparseMatrix(3,3)
85+
A[1,2]+=0.1
86+
A
87+
```
88+
89+
```@example
90+
using ExtendableSparse # hide
91+
92+
A=ExtendableSparseMatrix(3,3)
93+
updateindex!(A,+,0.1,1,2)
94+
A
95+
```
96+
"""
97+
function updateindex!(ext::ExtendableSparseMatrix{Tv,Ti}, op,v, i,j) where{Tv,Ti<:Integer}
98+
k=findindex(ext.cscmatrix,i,j)
99+
if k>0
100+
ext.cscmatrix.nzval[k]=op(ext.cscmatrix.nzval[k],v)
101+
else
102+
if ext.lnkmatrix==nothing
103+
ext.lnkmatrix=SparseMatrixLNK{Tv, Ti}(ext.cscmatrix.m, ext.cscmatrix.n)
104+
end
105+
updateindex!(ext.lnkmatrix,op,v,i,j)
93106
end
94-
return r1
107+
ext
95108
end
96109

97110

98-
99111
"""
100112
$(SIGNATURES)
101113
@@ -158,6 +170,7 @@ function flush!(ext::ExtendableSparseMatrix{Tv,Ti}) where {Tv, Ti<:Integer}
158170
return ext
159171
end
160172

173+
161174
"""
162175
$(SIGNATURES)
163176

‎src/sparsematrixcsc.jl‎

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
"""
2+
$(SIGNATURES)
3+
4+
Return index corresponding to entry [i,j] in the array of nonzeros,
5+
if the entry exists, otherwise, return 0.
6+
"""
7+
function findindex(csc::SparseMatrixCSC{T}, i::Integer, j::Integer) where T
8+
if !(1 <= i <= csc.m && 1 <= j <= csc.n); throw(BoundsError()); end
9+
r1 = Int(csc.colptr[j])
10+
r2 = Int(csc.colptr[j+1]-1)
11+
if r1>r2
12+
return 0
13+
end
14+
15+
# See sparsematrix.jl
16+
r1 = searchsortedfirst(csc.rowval, i, r1, r2, Base.Forward)
17+
if (r1>r2 ||csc.rowval[r1] != i)
18+
return 0
19+
end
20+
return r1
21+
end
22+
23+
"""
24+
$(SIGNATURES)
25+
26+
Update element of the matrix with operation `op`.
27+
This can replace the following code and save one index
28+
search during acces:
29+
30+
```@example
31+
using ExtendableSparse # hide
32+
using SparseArrays # hide
33+
A=spzeros(3,3)
34+
A[1,2]+=0.1
35+
A
36+
```
37+
38+
```@example
39+
using ExtendableSparse # hide
40+
using SparseArrays # hide
41+
A=spzeros(3,3)
42+
updateindex!(A,+,0.1,1,2)
43+
A
44+
```
45+
"""
46+
function updateindex!(csc::SparseMatrixCSC{Tv,Ti}, op, v,i::Integer, j::Integer) where{Tv,Ti<:Integer}
47+
k=findindex(csc,i,j)
48+
if k>0 # update existing value
49+
csc.nzval[k]=op(csc.nzval[k],v)
50+
else # insert new value
51+
csc[i,j]= op(zero(Tv),v)
52+
end
53+
csc
54+
end
55+
56+
"""
57+
$(SIGNATURES)
58+
59+
Trival flush! method for allowing to run the code with both `ExtendableSparseMatrix` and
60+
`SparseMatrixCSC`.
61+
"""
62+
flush!(csc::SparseMatrixCSC)=csc

0 commit comments

Comments
 (0)