Started to introduce preconditioners

Author: j-fu

Project.toml

@@ -1,7 +1,7 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>"]
-version = "0.2.3"
+version = "0.2.4"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

src/ExtendableSparse.jl

@@ -5,9 +5,11 @@ using LinearAlgebra
 
 include("sparsematrixlnk.jl")
 include("extendable.jl")
+include("preconditioners.jl")
 include("sprand.jl")
 
 export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, sprand!,sprand_sdd!
+export JacobiPreconditioner, ILU0Preconditioner
 
 export colptrs
 end # module

src/extendable.jl

@@ -17,6 +17,11 @@ mutable struct ExtendableSparseMatrix{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv
     Linked list structure holding data of extension
     """
     lnkmatrix::Union{SparseMatrixLNK{Tv,Ti},Nothing}
+
+    """
+    Time stamp of last pattern update
+    """
+    pattern_timestamp::Float64
 end
 
 
@@ -26,7 +31,7 @@ $(SIGNATURES)
 Create empty ExtendableSparseMatrix.
 """
 function ExtendableSparseMatrix{Tv,Ti}(m::Integer, n::Integer) where{Tv,Ti<:Integer}
-    ExtendableSparseMatrix{Tv,Ti}(spzeros(Tv,Ti,m,n),nothing)
+    ExtendableSparseMatrix{Tv,Ti}(spzeros(Tv,Ti,m,n),nothing,time())
 end
 
 """
@@ -62,7 +67,7 @@ $(SIGNATURES)
   Create ExtendableSparseMatrix from SparseMatrixCSC
 """
 function ExtendableSparseMatrix(csc::SparseMatrixCSC{Tv,Ti}) where{Tv,Ti<:Integer}
-    return ExtendableSparseMatrix{Tv,Ti}(csc, nothing)
+    return ExtendableSparseMatrix{Tv,Ti}(csc, nothing, time())
 end
 
 
@@ -148,6 +153,7 @@ function flush!(ext::ExtendableSparseMatrix{Tv,Ti}) where {Tv, Ti<:Integer}
     if ext.lnkmatrix!=nothing && nnz(ext.lnkmatrix)>0
         ext.cscmatrix=ext.lnkmatrix+ext.cscmatrix
         ext.lnkmatrix=nothing
+        ext.pattern_timestamp
     end
     return ext
 end

src/ilu0.jl

@@ -0,0 +1,85 @@
+mutable struct ILU0Preconditioner{Tv, Ti} <: AbstractPreconditioner{Tv,Ti}
+    extmatrix::ExtendableSparseMatrix{Tv,Ti}
+    xdiag::Array{Tv,1}
+    idiag::Array{Ti,1}
+    pattern_timestamp::Float64
+end
+
+
+function ILU0Preconditioner(extmatrix::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti}
+    @assert size(extmatrix,1)==size(extmatrix,2)
+    flush!(extmatrix)
+    n=size(extmatrix,1)
+    xdiag=Array{Tv,1}(undef,n)
+    idiag=Array{Ti,1}(undef,n)
+    precon=ILU0Preconditioner{Tv, Ti}(extmatrix,xdiag,idiag,0.0)
+    update!(precon)
+end
+
+function update!(precon::ILU0Preconditioner{Tv,Ti}) where {Tv,Ti}
+    cscmatrix=precon.extmatrix.cscmatrix
+    colptr=cscmatrix.colptr
+    rowval=cscmatrix.rowval
+    nzval=cscmatrix.nzval
+    n=cscmatrix.n
+    xdiag=precon.xdiag
+    idiag=precon.idiag
+
+    # Find main diagonal index and
+    # copy main diagonal values
+    if need_symbolic_update(precon)
+        @inbounds for j=1:n
+            @inbounds for k=colptr[j]:colptr[j+1]-1
+                i=rowval[k]
+                if i==j
+                    idiag[j]=k
+                    break
+                end
+            end
+        end
+        timestamp!(precon)
+    end
+    
+    @inbounds for j=1:n
+        xdiag[j]=one(Tv)/nzval[idiag[j]]
+        @inbounds for k=idiag[j]+1:colptr[j+1]-1
+            i=rowval[k]
+            for l=colptr[i]:colptr[i+1]-1
+                if rowval[l]==j
+                    xdiag[i]-=nzval[l]*xdiag[j]*nzval[k]
+                    break
+                end
+            end
+        end
+    end
+    precon
+end
+
+
+function  LinearAlgebra.ldiv!(u::AbstractArray{T,1}, precon::ILU0Preconditioner, v::AbstractArray{T,1}) where T
+    cscmatrix=precon.extmatrix.cscmatrix
+    colptr=cscmatrix.colptr
+    rowval=cscmatrix.rowval
+    n=cscmatrix.n
+    nzval=cscmatrix.nzval
+    xdiag=precon.xdiag
+    idiag=precon.idiag
+    
+    @inbounds for j=1:n
+        x=zero(T)
+        @inbounds for k=colptr[j]:idiag[j]-1
+            x+=nzval[k]*u[rowval[k]]
+        end
+        u[j]=xdiag[j]*(v[j]-x)
+    end
+    
+    @inbounds for j=n:-1:1
+        x=zero(T)
+        @inbounds for k=idiag[j]+1:colptr[j+1]-1
+            x+=u[rowval[k]]*nzval[k]
+        end
+        u[j]-=x*xdiag[j]
+    end
+end
+
+

src/jacobi.jl

@@ -0,0 +1,31 @@
+struct JacobiPreconditioner{Tv, Ti} <: AbstractPreconditioner{Tv,Ti}
+    extmatrix::ExtendableSparseMatrix{Tv,Ti}
+    invdiag::Array{Tv,1}
+end
+
+function update!(precon::JacobiPreconditioner)
+    cscmatrix=precon.extmatrix.cscmatrix
+    invdiag=precon.invdiag
+    n=cscmatrix.n
+    @inbounds for i=1:n
+        invdiag[i]=1.0/cscmatrix[i,i]
+    end
+    precon
+end
+
+function JacobiPreconditioner(extmatrix::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti}
+    @assert size(extmatrix,1)==size(extmatrix,2)
+    flush!(extmatrix)
+    invdiag=Array{Tv,1}(undef,extmatrix.cscmatrix.n)
+    precon=JacobiPreconditioner{Tv, Ti}(extmatrix,invdiag)
+    update!(precon)
+end
+
+
+function  LinearAlgebra.ldiv!(u::AbstractArray{T,1} where T, precon::JacobiPreconditioner, v::AbstractArray{T,1} where T)
+    invdiag=precon.invdiag
+    n=length(invdiag)
+    @inbounds @simd for i=1:n
+        u[i]=invdiag[i]*v[i]
+    end
+end

src/preconditioners.jl

@@ -0,0 +1,7 @@
+abstract type AbstractPreconditioner{Tv,Ti} end
+
+need_symbolic_update(precon::AbstractPreconditioner)=precon.extmatrix.pattern_timestamp>precon.pattern_timestamp
+timestamp!(precon::AbstractPreconditioner)= precon.pattern_timestamp=time()
+
+include("jacobi.jl")
+include("ilu0.jl")