re-organized+added tests

bugfix in getindex

Author: j-fu

examples/ExtendableSparseMatrixTest.jl

@@ -5,8 +5,9 @@ using LinearAlgebra
 
 include("extension.jl")
 include("extendable.jl")
+include("sprand.jl")
 
-export SparseMatrixExtension,ExtendableSparseMatrix,flush!,nnz
+export SparseMatrixExtension,ExtendableSparseMatrix,flush!,nnz,sprand!
 
 export xcolptrs,colptrs
 end # module

src/ExtendableSparse.jl

@@ -32,6 +32,9 @@ ExtendableSparseMatrix(::Type{Tv}, ::Type{Ti}, m::Integer, n::Integer) where {Tv
 
 
 
+
+
+
 """
 $(TYPEDSIGNATURES)
 
@@ -43,18 +46,21 @@ function findindex(S::SparseMatrixCSC{T}, i::Integer, j::Integer) where T
     r1 = Int(S.colptr[j])
     r2 = Int(S.colptr[j+1]-1)
     if r1>r2
-        return 0
+        return zero(T)
     end
 
     # See sparsematrix.jl
     r1 = searchsortedfirst(S.rowval, i, r1, r2, Base.Forward)
-    if (r1>length(S.rowval) ||S.rowval[r1] != i)
-        return 0
+    if (r1>r2 ||S.rowval[r1] != i)
+        return zero(T)
     end
     return r1
 end
 
 
+
+
+
 """
 $(TYPEDSIGNATURES)
 
@@ -70,6 +76,10 @@ function Base.setindex!(M::ExtendableSparseMatrix, v, i::Integer, j::Integer)
     end
 end
 
+
+
+
+
 """
 $(TYPEDSIGNATURES)
 
@@ -86,6 +96,10 @@ function Base.getindex(M::ExtendableSparseMatrix,i::Integer, j::Integer)
 end
 
 
+
+
+
+
 """
 $(TYPEDSIGNATURES)
 
@@ -94,6 +108,10 @@ Matrix size.
 Base.size(E::ExtendableSparseMatrix) = (E.cscmatrix.m, E.cscmatrix.n)
 
 
+
+
+
+
 """
 $(TYPEDSIGNATURES)
 
@@ -103,6 +121,11 @@ SparseArrays.nnz(E::ExtendableSparseMatrix)=(nnz(E.cscmatrix)+nnz(E.extmatrix))
 
 
 
+
+
+
+
+
 # Struct holding pair of value and row
 # number, for sorting
 mutable struct ColEntry{Tv,Ti<:Integer}
@@ -114,7 +137,6 @@ end
 Base.isless(x::ColEntry{Tv, Ti},y::ColEntry{Tv, Ti}) where {Tv,Ti<:Integer} = (x.i<y.i)
 
 
-
 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.
     #
@@ -130,7 +152,7 @@ function _splice(E::SparseMatrixExtension{Tv,Ti},S::SparseMatrixCSC{Tv,Ti}) wher
     rowval=Vector{Ti}(undef,xnnz)
     nzval=Vector{Tv}(undef,xnnz)
 
-    # Detect the maximum column lenght of E
+    # Detect the maximum column length of E
     E_maxcol=0
     for j=1:S.n
         lcol=0
@@ -143,7 +165,7 @@ function _splice(E::SparseMatrixExtension{Tv,Ti},S::SparseMatrixCSC{Tv,Ti}) wher
     end
 
     # pre-allocate column 
-    col=[ColEntry{Tv,Ti}(0,0) for i=1:E_maxcol+1]
+    col=[ColEntry{Tv,Ti}(0,0) for i=1:E_maxcol]
 
 
 
@@ -199,6 +221,9 @@ function _splice(E::SparseMatrixExtension{Tv,Ti},S::SparseMatrixCSC{Tv,Ti}) wher
 end
 
 
+
+
+
 """
 $(TYPEDSIGNATURES)
 
@@ -215,6 +240,9 @@ end
 
 
 
+
+
+
 """
 $(TYPEDSIGNATURES)
 
@@ -225,6 +253,9 @@ function SparseArrays.nonzeros(E::ExtendableSparseMatrix)
     return nonzeros(E.cscmatrix)
 end
 
+
+
+
 """
 $(TYPEDSIGNATURES)
 
@@ -236,6 +267,9 @@ function SparseArrays.rowvals(E::ExtendableSparseMatrix)
 end
 
 
+
+
+
 """
 $(TYPEDSIGNATURES)
 
@@ -246,6 +280,9 @@ function xcolptrs(E::ExtendableSparseMatrix)
     return E.cscmatrix.colptr
 end
 
+
+
+
 """
 $(TYPEDSIGNATURES)
 
@@ -257,6 +294,16 @@ function colptrs(E::ExtendableSparseMatrix)
 end
 
 
+"""
+$(TYPEDSIGNATURES)
+
+Flush and delegate to cscmatrix.
+"""
+function SparseArrays.findnz(E::ExtendableSparseMatrix)
+    flush!(E)
+    return findnz(E.cscmatrix)
+end
+
 
 """
 $(TYPEDSIGNATURES)
@@ -270,6 +317,8 @@ function LinearAlgebra.lu(E::ExtendableSparseMatrix)
 end
 
 
+
+
 """
 $(TYPEDSIGNATURES)
 

src/extendable.jl

@@ -0,0 +1,17 @@
+"""
+$(TYPEDSIGNATURES)
+
+Fill empty sparse matrix A with random nonzero elements from interval [1,2]
+using incremental assembly.
+
+"""
+function sprand!(A::AbstractSparseMatrix{Tv,Ti},xnnz::Int) where {Tv,Ti}
+    m,n=size(A)
+    for i=1:xnnz
+        i=rand((1:m))
+        j=rand((1:n))
+        a=1.0+rand(Tv)
+        A[i,j]+=a
+    end
+end
+

src/sprand.jl

@@ -0,0 +1,70 @@
+module ExtendableSparseTest
+using ExtendableSparse
+using SparseArrays
+using Printf
+
+
+
+function randtest!(A::AbstractSparseMatrix{Tv,Ti},xnnz::Int,nsplice::Int) where {Tv,Ti}
+    m,n=size(A)
+    S=spzeros(m,n)
+    for isplice=1:nsplice
+        for inz=1:xnnz
+            i=rand((1:m))
+            j=rand((1:n))
+            a=1.0+rand(Float64)
+            S[i,j]+=a
+            A[i,j]+=a
+            @assert(nnz(S)==nnz(A))
+        end
+        flush!(A)
+        for j=1:n
+            @assert(issorted(A.cscmatrix.rowval[A.cscmatrix.colptr[j]:A.cscmatrix.colptr[j+1]-1]))
+        end
+        @assert(nnz(S)==nnz(A))
+
+        (I,J,V)=findnz(S)
+        for inz=1:nnz(S)
+            @assert(A[I[inz],J[inz]]==V[inz])
+        end
+
+        (I,J,V)=findnz(A)
+        for inz=1:nnz(A)
+            @assert(S[I[inz],J[inz]]==V[inz])
+        end
+
+    end
+    return true
+end
+
+
+
+function check(;m=1000,n=1000,nnz=5000,nsplice=1)
+    mat=ExtendableSparseMatrix(Float64,Int64,m,n)
+    return randtest!(mat,nnz,nsplice)
+end
+
+
+
+function benchmark(;n=10000,m=10000,nnz=50000)
+    
+    println("SparseMatrixCSC:")
+    mat=spzeros(Float64,Int64,m,n)
+    @time sprand!(mat,nnz)
+    
+    println("SparseMatrixExtension:")
+    extmat=SparseMatrixExtension(Float64,Int64,m,n)
+    @time sprand!(extmat,nnz)
+    
+    println("ExtendableSparseMatrix:")
+    xextmat=ExtendableSparseMatrix(Float64,Int64,m,n)
+    @time begin
+        sprand!(xextmat,nnz)
+        @inbounds flush!(xextmat)
+    end
+    return
+    b
+end
+
+
+end