Merge pull request #13 from gridap/adding_lu_to_CSR

Implemented support for lu, lu! + SparseMatrixCSR ...

Author: fverdugo

Project.toml

@@ -1,11 +1,12 @@
 name = "SparseMatricesCSR"
 uuid = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
 authors = ["Víctor Sande <[email protected]>", "Francesc Verdugo <[email protected]>"]
-version = "0.6.2"
+version = "0.6.3"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [compat]
 julia = "1"

src/SparseMatricesCSR.jl

@@ -2,9 +2,10 @@ module SparseMatricesCSR
 
 using SparseArrays
 using LinearAlgebra
+using SuiteSparse
 
 import Base: convert, copy, size, getindex, show, count, *, IndexStyle
-import LinearAlgebra: mul!
+import LinearAlgebra: mul!, lu, lu!
 import SparseArrays: nnz, getnzval, nonzeros, nzrange
 import SparseArrays: findnz, rowvals, getnzval, issparse
 

src/SparseMatrixCSR.jl

@@ -118,6 +118,32 @@ function Base.copy(a::SparseMatrixCSR{Bi}) where Bi
   SparseMatrixCSR{Bi}(a.m,a.n,copy(a.rowptr),copy(a.colval),copy(a.nzval))
 end
 
+_copy_and_increment(x) = copy(x) .+ 1
+
+function LinearAlgebra.lu(a::SparseMatrixCSR{0})
+  rowptr = _copy_and_increment(a.rowptr)
+  colval = _copy_and_increment(a.colval)
+  Transpose(lu(SparseMatrixCSC(a.m,a.n,rowptr,colval,a.nzval)))
+end
+
+function LinearAlgebra.lu(a::SparseMatrixCSR{1})
+  Transpose(lu(SparseMatrixCSC(a.m,a.n,a.rowptr,a.colval,a.nzval)))
+end
+
+function LinearAlgebra.lu!(
+    translu::Transpose{T,<:SuiteSparse.UMFPACK.UmfpackLU{T}},
+    a::SparseMatrixCSR{1}) where {T}
+  Transpose(lu!(translu.parent,SparseMatrixCSC(a.m,a.n,a.rowptr,a.colval,a.nzval)))
+end
+
+function LinearAlgebra.lu!(
+  translu::Transpose{T,<:SuiteSparse.UMFPACK.UmfpackLU{T}},
+  a::SparseMatrixCSR{0}) where {T}
+  rowptr = _copy_and_increment(a.rowptr)
+  colval = _copy_and_increment(a.colval)
+  Transpose(lu!(translu.parent,SparseMatrixCSC(a.m,a.n,rowptr,colval,a.nzval)))
+end
+
 size(S::SparseMatrixCSR) = (S.m, S.n)
 IndexStyle(::Type{<:SparseMatrixCSR}) = IndexCartesian()
 function getindex(A::SparseMatrixCSR{Bi,T}, i0::Integer, i1::Integer) where {Bi,T}

test/SparseMatrixCSR.jl

@@ -82,9 +82,21 @@ function test_csr(Bi,Tv,Ti)
   @test y ≈ z
 
   @test CSR*x ≈ CSC*x
+end
 
+function test_lu(Bi,I,J,V)
+  CSR=sparsecsr(Val(Bi),I,J,V)
+  CSC=sparse(I,J,V)
+  x=rand(3)
+  @test norm(CSR\x-CSC\x) < 1.0e-14
+  fact=lu(CSR)
+  lu!(fact,CSR)
+  y=similar(x)
+  ldiv!(y,fact,x)
+  @test norm(y-CSC\x) < 1.0e-14
 end
 
+
 for Bi in (0,1)
   for Tv in (Float32,Float64)
     for Ti in (Int32,Int64)
@@ -93,4 +105,10 @@ for Bi in (0,1)
   end
 end
 
+I = [1,1,2,2,2,3,3]
+J = [1,2,1,2,3,2,3]
+V = [4.0,1.0,-1.0,4.0,1.0,-1.0,4.0]
+test_lu(0,I,J,V)
+test_lu(1,I,J,V)
+
 end # module