Merge pull request #10 from jkrch/master

Added parallel jacobi

Author: j-fu

src/ExtendableSparse.jl

@@ -10,7 +10,7 @@ include("preconditioners.jl")
 include("sprand.jl")
 
 export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, sprand!,sprand_sdd!, fdrand,fdrand!
-export JacobiPreconditioner, ILU0Preconditioner, updateindex!
+export JacobiPreconditioner, ILU0Preconditioner, ParallelJacobiPreconditioner, updateindex!
 
 export colptrs
 end # module

src/parallel_jacobi.jl

@@ -0,0 +1,37 @@
+struct ParallelJacobiPreconditioner{Tv, Ti} <: AbstractExtendablePreconditioner{Tv,Ti}
+    extmatrix::ExtendableSparseMatrix{Tv,Ti}
+    invdiag::Array{Tv,1}
+end
+
+function update!(precon::ParallelJacobiPreconditioner)
+    cscmatrix=precon.extmatrix.cscmatrix
+    invdiag=precon.invdiag
+    n=cscmatrix.n
+    Threads.@threads for i=1:n
+        @inbounds  invdiag[i]=1.0/cscmatrix[i,i]
+    end
+    precon
+end
+
+function ParallelJacobiPreconditioner(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
+
+ParallelJacobiPreconditioner(cscmatrix::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}=ParallelJacobiPreconditioner(ExtendableSparseMatrix(cscmatrix))
+
+
+function  LinearAlgebra.ldiv!(u::AbstractArray{T,1} where T, precon::ParallelJacobiPreconditioner, v::AbstractArray{T,1} where T)
+    invdiag=precon.invdiag
+    n=length(invdiag)
+    Threads.@threads for i=1:n
+        @inbounds u[i]=invdiag[i]*v[i]
+    end
+end
+
+function LinearAlgebra.ldiv!(precon::ParallelJacobiPreconditioner, v::AbstractArray{T,1} where T)
+    ldiv!(v, precon, v)
+end

src/preconditioners.jl

@@ -25,4 +25,6 @@ timestamp!(precon::AbstractExtendablePreconditioner)= precon.pattern_timestamp=t
 
 include("jacobi.jl")
 include("ilu0.jl")
+
+include("parallel_jacobi.jl")