added linear algebra tests and example

Author: j-fu

.travis.yml

@@ -3,6 +3,7 @@ language: julia
 os:
   - osx
   - linux
+  - windows
 
 julia:
   - 1.0

Project.toml

@@ -4,7 +4,9 @@ authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>"]
 version = "0.1.3"
 
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

README.md

@@ -1,6 +1,6 @@
 # ExtendableSparse.jl
 
-[![Build Status](https://img.shields.io/travis/j-fu/ExtendableSparse.jl/master.svg?label=Linux+MacOSX)](https://travis-ci.org/j-fu/ExtendableSparse.jl)
+[![Build Status](https://img.shields.io/travis/j-fu/ExtendableSparse.jl/master.svg?label=Linux+MacOSX+Windows)](https://travis-ci.org/j-fu/ExtendableSparse.jl)
 [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://j-fu.github.io/ExtendableSparse.jl/stable)
 [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://j-fu.github.io/ExtendableSparse.jl/dev)
 
@@ -12,5 +12,4 @@ Without an intermediate data structure, efficient successive insertion/update of
 
 `ExtendableSparseMatrix` is aimed to work as a drop-in replacement to `SparseMatrixCSC` in finite element and finite volume codes.
 
-Currently it has the methods required for `AbstractSparseMatrix` (`getindex`, `setindex!`,`size`,`nnz`), as well as `lufact` and `mul!`, which is already sufficient for a number of interesting applications.
 

docs/make.jl

@@ -9,6 +9,7 @@ makedocs(sitename="ExtendableSparse.jl",
          repo="https://github.com/j-fu/ExtendableSparse.jl",
          pages=[
              "Home"=>"index.md",
+             "example.md",
              "changes.md",
              "api.md"
          ])

docs/src/changes.md

@@ -1,4 +1,8 @@
 # Changes
+## dev
+- more interface methods delegating to csc, in particular mul! and ldiv!
+- lazy creation of extendable part: don't create idle memory
+- nicer constructors
 
 ## V0.1, July 2019
 

docs/src/example.md

@@ -0,0 +1,43 @@
+# Examples
+
+An `ExtendableSparseMatrix` can serve as a drop-in replacement for
+`SparseMatrixCSC`, albeit with faster assembly.
+
+The code below provides a small benchmark example. 
+```julia
+using ExtendableSparse
+using SparseArrays
+
+function ext_create(n, row_nz)
+    A=ExtendableSparseMatrix(n,n)
+    sprand_sdd!(A,row_nz);
+end
+
+function csc_create(n, row_nz)
+    A=spzeros(n,n)
+    sprand_sdd!(A,row_nz);
+end
+
+csc_create(10,2);
+ext_create(10,2);
+n=65536
+row_nz=5
+@time Acsc=csc_create(n,row_nz);
+nnz(Acsc)
+@time Aext=ext_create(n,row_nz);
+nnz(Aext)
+b=rand(n);
+@time Acsc*(Acsc\b)≈b
+@time Aext*(Aext\b)≈b
+```
+
+The function [`sprand_sdd!`](@ref) fills a 
+sparse matrix with random entries such that it becomes strictly diagonally
+dominant and thus invertible and has a fixed number of nonzeros in
+its rows. Its  bandwidth is bounded by 2*sqrt(n), therefore it 
+resembles a typical matrix of a 2D piecewise linear FEM discretization.
+For filling a matrix a, the method  conveniently albeit naively
+uses just `a[i,j]=value`. This approach is considerably faster with 
+the [`ExtendableSparseMatrix`](@ref).
+
+

src/ExtendableSparse.jl

@@ -7,7 +7,7 @@ include("extension.jl")
 include("extendable.jl")
 include("sprand.jl")
 
-export SparseMatrixExtension,ExtendableSparseMatrix,flush!,nnz,sprand!
+export SparseMatrixExtension,ExtendableSparseMatrix,flush!,nnz, sprand!,sprand_sdd!
 
-export xcolptrs,colptrs
+export colptrs
 end # module

src/extendable.jl

@@ -39,7 +39,7 @@ function ExtendableSparseMatrix(::Type{Tv},::Type{Ti},m::Integer, n::Integer) wh
 end
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Create empty ExtendablSparseMatrix.
 This is a pendant to spzeros.
@@ -48,7 +48,7 @@ ExtendableSparseMatrix(::Type{Tv},m::Integer, n::Integer) where{Tv}=ExtendableSp
 
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Create empty ExtendablSparseMatrix.
 This is a pendant to spzeros.
@@ -69,7 +69,7 @@ end
 
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Return index corresponding to entry (i,j) in the array of nonzeros,
 if the entry exists, otherwise, return 0.
@@ -132,15 +132,15 @@ function Base.getindex(M::ExtendableSparseMatrix{Tv,Ti},i::Integer, j::Integer)
 end
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Size of ExtendableSparseMatrix.
 """
 Base.size(E::ExtendableSparseMatrix) = (E.cscmatrix.m, E.cscmatrix.n)
 
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Number of nonzeros of ExtendableSparseMatrix.
 """
@@ -165,13 +165,12 @@ end
 Base.isless(x::ColEntry{Tv, Ti},y::ColEntry{Tv, Ti}) where {Tv,Ti<:Integer} = (x.i<y.i)
 
 
+# Create new CSC matrix with sorted entries from  CSC matrix S and matrix extension E.
+#
+# This method assumes that there are no entries with the same
+# indices in E and S, therefore  it appears too dangerous for general use and
+# so we don't export it. Generalizations appear to be possible, though.
 function _splice(E::SparseMatrixExtension{Tv,Ti},S::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti<:Integer}
-    # Create new CSC matrix with sorted entries from  CSC matrix S and matrix extension E.
-    #
-    # This method assumes that there are no entries with the same
-    # indices in E and S, therefore  it appears too dangerous for general use and
-    # so we don't export it. Generalizations appear to be possible, though.
-
     @assert(S.m==E.m)
     @assert(S.n==E.n)
 
@@ -191,8 +190,8 @@ function _splice(E::SparseMatrixExtension{Tv,Ti},S::SparseMatrixCSC{Tv,Ti}) wher
         end
         E_maxcol=max(lcol,E_maxcol)
     end
-
-    # pre-allocate column 
+    
+    # pre-allocate column  data
     col=[ColEntry{Tv,Ti}(0,0) for i=1:E_maxcol]
 
 
@@ -265,7 +264,7 @@ end
 
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Flush and delegate to cscmatrix.
 """
@@ -276,9 +275,8 @@ end
 
 
 
-
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Flush and delegate to cscmatrix.
 """
@@ -289,29 +287,18 @@ end
 
 
 """
-$(TYPEDSIGNATURES)
-
-Flush and delegate to cscmatrix.
-"""
-function xcolptrs(E::ExtendableSparseMatrix)
-    @inbounds flush!(E)
-    return E.cscmatrix.colptr
-end
-
-
-"""
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Flush and delegate to cscmatrix.
 """
 function colptrs(E::ExtendableSparseMatrix)
-    flush!(E)
+    @inbounds flush!(E)
     return E.cscmatrix.colptr
 end
 
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Flush and delegate to cscmatrix.
 """
@@ -322,7 +309,7 @@ end
 
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Delegating LU factorization.
 """
@@ -332,7 +319,7 @@ function LinearAlgebra.lu(E::ExtendableSparseMatrix)
 end
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Delegating Matrix multiplication
 """
@@ -343,7 +330,7 @@ end
 
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Delegating Matrix ldiv
 """
@@ -353,7 +340,7 @@ function  LinearAlgebra.ldiv!(r::AbstractArray{T,1} where T, E::ExtendableSparse
 end
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Delegating Matrix multiplication
 """
@@ -364,7 +351,7 @@ end
 
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Delegating Matrix ldiv
 """

src/extension.jl

@@ -58,7 +58,7 @@ end
 
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
     
 Constructor of empty extension
 """
@@ -142,23 +142,23 @@ end
 
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Return tuple containing size of the matrix.
 """
 Base.size(E::SparseMatrixExtension) = (E.m, E.n)
 
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Return number of nonzero entries.
 """
 SparseArrays.nnz(E::SparseMatrixExtension)=E.nnz
 
 
 """
-$(TYPEDSIGNATURES)
+$(SIGNATURES)
 
 Dummy flush! method for Sparse matrix extension. Just
 used in thest methods

src/sprand.jl

@@ -13,5 +13,37 @@ function sprand!(A::AbstractSparseMatrix{Tv,Ti},xnnz::Int) where {Tv,Ti}
         a=1.0+rand(Tv)
         A[i,j]+=a
     end
+    A
+end
+
+"""
+$(SIGNATURES)
+
+Fill sparse matrix with random entries such that it becomes strictly diagonally
+dominant and thus invertible and has a fixed number of nonzeros in
+its rows. The matrix bandwidth is bounded by sqrt(n), therefore it 
+resembles a typical matrix of a 2D piecewise linear FEM discretization
+
+"""
+function sprand_sdd!(A::AbstractSparseMatrix{Tv,Ti},nnzrow::Int; dim=2) where {Tv,Ti}
+    m,n=size(A)
+    @assert m==n
+    nnzrow=min(n,nnzrow)
+    bandwidth=convert(Int,ceil(sqrt(n)))
+    for i=1:n
+        aii=0
+        for k=1:nnzrow
+            jmin=max(1,i-bandwidth)
+            jmax=min(n,i+bandwidth)
+            j=rand((jmin:jmax))
+            if i!=j
+                aij=rand()
+                A[i,j]=aij
+                aii+=abs(aij)
+            end
+        end
+        A[i,i]=aii+rand() # make it strictly diagonally dominant
+    end
+    A
 end