Replace _slice by Base.:+, now fully fledged addition

Author: j-fu

README.md

@@ -12,7 +12,7 @@ Without an intermediate data structure, efficient successive insertion/update of
 
 The later is modeled after the linked list sparse matrix format described in the [whitepaper](https://www-users.cs.umn.edu/~saad/software/SPARSKIT/paper.ps) by Y. Saad. See also exercise P.3-16  in his [book](https://www-users.cs.umn.edu/~saad/IterMethBook_2ndEd.pdf).
 
-Any linear algebra method on `ExtendableSparseMatrix` starts with a `flush!` method which splices the LNK entries and the existing CSC entries into a new CSC struct and resets the LNK struct.
+Any linear algebra method on `ExtendableSparseMatrix` starts with a `flush!` method which adds the LNK entries and the existing CSC entries into a new CSC struct and resets the LNK struct.
 
 `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.
 

docs/make.jl

@@ -3,7 +3,7 @@ using Documenter, ExtendableSparse
 
 makedocs(sitename="ExtendableSparse.jl",
          modules = [ExtendableSparse],
-         doctest = false,
+         doctest = true,
          clean = true,
          authors = "J. Fuhrmann",
          repo="https://github.com/j-fu/ExtendableSparse.jl",

docs/src/changes.md

@@ -1,4 +1,10 @@
 # Changes
+
+## dev
+- What used to be `_splice`  is now `+` and allows now real addition (resulting in a CSC matrix)
+- Added constructors of LNK matrix from CSC matrix and vice versa
+- reorganized tests
+
 ## v0.2.1 Dec 22, 2019
 - Tried to track down the source from which I learned the linked list based struct in order
   to document this. Ended up with SPARSEKIT of Y.Saad, however I believe this 

docs/src/example.md

@@ -1,8 +1,35 @@
 # Examples
 
+
+## Usage idea
 An `ExtendableSparseMatrix` can serve as a drop-in replacement for
 `SparseMatrixCSC`, albeit with faster assembly.
 
+    
+```jldoctest
+using ExtendableSparse
+function example(n,bandwidth)
+    A=ExtendableSparseMatrix(n,n);
+    aii=0.0
+    for i=1:n
+        for j=max(1,i-bandwidth):5:min(n,i+bandwidth)
+            aij=-rand()
+            A[i,j]=aij
+            aii+=abs(aij)
+        end
+        A[i,i]=aii+0.1
+    end
+    b=rand(n)
+    x=A\b
+    A*x≈ b
+end
+example(1000,100)
+# output
+true
+```
+
+
+## Benchmark
 The code below provides a small benchmark example.
 ```julia
 using ExtendableSparse
@@ -35,7 +62,7 @@ b=rand(n);
 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 
+its rows (similar to the method in the first example). 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 

src/extendable.jl

@@ -138,132 +138,30 @@ Size of ExtendableSparseMatrix.
 Base.size(ext::ExtendableSparseMatrix) = (ext.cscmatrix.m, ext.cscmatrix.n)
 
 
-"""
-$(SIGNATURES)
-
-Number of nonzeros of ExtendableSparseMatrix.
-"""
-function SparseArrays.nnz(ext::ExtendableSparseMatrix)
-    ennz=0
-    if ext.lnkmatrix!=nothing
-        ennz=nnz(ext.lnkmatrix)
-    end
-    return nnz(ext.cscmatrix)+ennz
-end
-
-
-
-# Struct holding pair of value and row
-# number, for sorting
-struct ColEntry{Tv,Ti<:Integer}
-    rowval::Ti
-    nzval::Tv
-end
-
-# Comparison method for sorting
-Base.isless(x::ColEntry{Tv, Ti},y::ColEntry{Tv, Ti}) where {Tv,Ti<:Integer} = (x.rowval<y.rowval)
-
-"""
-$(SIGNATURES)
-
-Create new CSC matrix with sorted entries from a SparseMatrixCSC  csc and 
-[`SparseMatrixLNK`](@ref)  lnk.
-
-This method assumes that there are no entries with the same
-indices in lnk and csc, therefore  it appears too dangerous for general use and
-so we don't export it.  A generalization appears to be possible, though.
-"""
-function _splice(lnk::SparseMatrixLNK{Tv,Ti},csc::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti<:Integer}
-    @assert(csc.m==lnk.m)
-    @assert(csc.n==lnk.n)
-
-    xnnz=nnz(csc)+nnz(lnk)
-    colptr=Vector{Ti}(undef,csc.n+1)
-    rowval=Vector{Ti}(undef,xnnz)
-    nzval=Vector{Tv}(undef,xnnz)
-
-    # Detect the maximum column length of lnk
-    lnk_maxcol=0
-    for j=1:csc.n
-        lcol=0
-        k=j
-        while k>0
-            lcol+=1
-            k=lnk.colptr[k]
-        end
-        lnk_maxcol=max(lcol,lnk_maxcol)
-    end
-    
-    # pre-allocate column  data
-    col=[ColEntry{Tv,Ti}(0,0) for i=1:lnk_maxcol]
-
-
-    
-    inz=1 # counts the nonzero entries in the new matrix
-
-    # loop over all columns
-    for j=1:csc.n
-        # put extension entries into col and sort them
-        k=j
-        l_lnk_col=0
-        while k>0
-            if lnk.rowval[k]>0
-                l_lnk_col+=1
-                col[l_lnk_col]=ColEntry(lnk.rowval[k],lnk.nzval[k])
-            end
-            k=lnk.colptr[k]
-        end
-        sort!(col,1,l_lnk_col, Base.QuickSort, Base.Forward)
-
-
-        # jointly sort lnk and csc entries  into new matrix data
-        colptr[j]=inz
-        jcol=1 # counts the entries in col
-        jcsc=csc.colptr[j]  # counts entries in csc
-        while true
-            
-            # Insert entries from csc into new structure
-            # if the row number is before col[jcol]
-            if ( nnz(csc)>0 &&  (jcsc<csc.colptr[j+1]) ) && 
-                 (
-                     (jcol<=l_lnk_col && csc.rowval[jcsc]<col[jcol].rowval) || 
-                     (jcol>l_lnk_col)
-                 )
-                rowval[inz]=csc.rowval[jcsc]
-                nzval[inz]=csc.nzval[jcsc]
-                jcsc+=1
-                inz+=1
-
-            # Otherwise, insert next entry from col    
-            elseif jcol<=l_lnk_col
-                rowval[inz]=col[jcol].rowval
-                nzval[inz]=col[jcol].nzval
-                jcol+=1
-                inz+=1
-            else
-                break
-            end
-        end
-    end
-    colptr[csc.n+1]=inz
-    SparseMatrixCSC{Tv,Ti}(csc.m,csc.n,colptr,rowval,nzval)
-end
-
 
 """
 $(SIGNATURES)
 
-If there are new entries in extension, create new CSC matrix using [`_splice`](@ref)
-and reset linked list based extension.
+If there are new entries in extension, create new CSC matrix by adding the
+cscmatrix and linked list matrix and reset the linked list based extension.
 """
 function flush!(ext::ExtendableSparseMatrix{Tv,Ti}) where {Tv, Ti<:Integer}
     if ext.lnkmatrix!=nothing && nnz(ext.lnkmatrix)>0
-        ext.cscmatrix=_splice(ext.lnkmatrix,ext.cscmatrix)
+        ext.cscmatrix=ext.lnkmatrix+ext.cscmatrix
         ext.lnkmatrix=nothing
     end
     return ext
 end
 
+"""
+$(SIGNATURES)
+
+[`flush!`](@ref) and return number of nonzeros in ext.cscmatrix.
+"""
+function SparseArrays.nnz(ext::ExtendableSparseMatrix)
+    @inbounds flush!(ext)
+    return nnz(ext.cscmatrix)
+end
 
 
 """