added more tests, version bump

- added simple iteration scheme
- fixed fdrand (boundray terms)
- added tests for fdrand, 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.3.6"
+version = "0.3.7"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"

src/ExtendableSparse.jl

@@ -6,11 +6,14 @@ using LinearAlgebra
 include("sparsematrixcsc.jl")
 include("sparsematrixlnk.jl")
 include("extendable.jl")
+export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, updateindex, colptrs
+
 include("preconditioners.jl")
-include("sprand.jl")
+export JacobiPreconditioner, ILU0Preconditioner, ParallelJacobiPreconditioner
 
-export SparseMatrixLNK,ExtendableSparseMatrix,flush!,nnz, sprand!,sprand_sdd!, fdrand,fdrand!
-export JacobiPreconditioner, ILU0Preconditioner, ParallelJacobiPreconditioner, updateindex!
+include("simple_iteration.jl")
+export simple,simple!
 
-export colptrs
+include("sprand.jl")
+export sprand!,sprand_sdd!, fdrand,fdrand!
 end # module

src/preconditioners.jl

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

src/simple_iteration.jl

@@ -0,0 +1,44 @@
+"""
+````        
+simple!(u,A,b;
+                 abstol::Real = zero(real(eltype(b))),
+                 reltol::Real = sqrt(eps(real(eltype(b)))),
+                 log=false,
+                 maxiter=100,
+                 P=nothing
+                 ) -> soution, [history]
+````
+
+Simple iteration scheme ``u_{i+1}= u_i - P^{-1} (A u_i -b)`` with similar API as the methods in IterativeSolvers.jl.
+
+"""
+function simple!(u,A,b;
+                 abstol::Real = zero(real(eltype(b))),
+                 reltol::Real = sqrt(eps(real(eltype(b)))),
+                 log=false,
+                 maxiter=100,
+                 Pl=nothing
+                 )
+    res=A*u-b # initial residual
+    upd=similar(res)
+    r0=norm(res) # residual norm
+    history=[r0]  # intialize history recording
+    for i=1:maxiter
+        ldiv!(upd,Pl,res)
+	u.-=upd # solve preconditioning system and update solution
+	mul!(res,A,u)         # calculate residual
+        res.-=b
+	r=norm(res)       # residual norm
+	push!(history,r)  # record in history 
+	if (r/r0)<reltol || r<abstol    # check for tolerance
+            break
+	end
+    end
+    if log 
+	return u,Dict( :resnorm =>history )
+    else
+        return u
+    end
+end
+
+simple(A,b;kwargs...)=simple!(zeros(eltype(b),length(b)),A,b;kwargs...)

src/sprand.jl

@@ -101,19 +101,19 @@ function fdrand!(A::AbstractMatrix,
                     update_pair(A,rand()*hy*hz/hx,l,l+1)
                 end
                 if i==1|| i==nx
-                    update(A,rand(),l,l)
+                    update(A,rand()*hy*hz,l,l)
                 end
                 if j<ny
                     update_pair(A,rand()*hx*hz/hy,l,l+nx)
                 end
                 if ny>2&&(j==1|| j==ny)
-                    update(A,rand(),l,l)
+                    update(A,rand()*hx*hz,l,l)
                 end
                 if k<nz
                     update_pair(A,rand()*hx*hy/hz,l,l+nxy)
                 end
                 if nz>2&&(k==1|| k==nz)
-                    update(A,rand(),l,l)
+                    update(A,rand()*hx*hy,l,l)
                 end
                 l=l+1
             end
@@ -162,7 +162,7 @@ Tested for Matrix, SparseMatrixCSC and ExtendableSparseMatrix.
 function fdrand(nx,ny,nz;
                 matrixtype::Type{Tv}=SparseMatrixCSC,
                 update = (A,v,i,j)-> A[i,j]+=v,
-                rand=()-> rand()) where Tv
+                rand= ()-> 0.1+rand()) where Tv
     N=nx*ny*nz
     if matrixtype==ExtendableSparseMatrix
         A=ExtendableSparseMatrix(N,N)