Adding support for lu+SparseMatrixCSR{0}

Author: amartinhuertas

src/SparseMatrixCSR.jl

@@ -118,8 +118,12 @@ 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})
-  @assert false "Base.lu(a::SparseMatrixCSR{0}) not yet implemented"
+  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})
@@ -132,6 +136,14 @@ function LinearAlgebra.lu!(
   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

@@ -84,6 +84,19 @@ function test_csr(Bi,Tv,Ti)
   @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)
@@ -95,14 +108,7 @@ 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]
-CSR=sparsecsr(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
+test_lu(0,I,J,V)
+test_lu(1,I,J,V)
 
 end # module