Merge pull request #39 from JuliaDiffEq/myb/fact

Add `lu_instance`

Author: ChrisRackauckas

Project.toml

@@ -1,6 +1,6 @@
 name = "ArrayInterface"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "2.4.1"
+version = "2.5.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -16,7 +16,8 @@ BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "LabelledArrays", "StaticArrays", "BandedMatrices", "BlockBandedMatrices"]
+test = ["Test", "LabelledArrays", "StaticArrays", "BandedMatrices", "BlockBandedMatrices", "SuiteSparse"]

README.md

@@ -60,6 +60,11 @@ A `getindex` which is always allowed.
 
 A `setindex!` which is always allowed.
 
+## lu_instance(A)
+
+Return an instance of the LU factorization object with the correct type
+cheaply.
+
 ## List of things to add
 
 - https://github.com/JuliaLang/julia/issues/22216

src/ArrayInterface.jl

@@ -346,8 +346,8 @@ abstract type ColoringAlgorithm end
 """
     fast_matrix_colors(A)
 
-    Query whether a matrix has a fast algorithm for getting the structural
-    colors of the matrix.
+Query whether a matrix has a fast algorithm for getting the structural
+colors of the matrix.
 """
 fast_matrix_colors(A) = false
 fast_matrix_colors(A::AbstractArray) = fast_matrix_colors(typeof(A))
@@ -356,8 +356,8 @@ fast_matrix_colors(A::Type{<:Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagona
 """
     matrix_colors(A::Union{Array,UpperTriangular,LowerTriangular})
 
-    The color vector for dense matrix and triangular matrix is simply
-    `[1,2,3,...,size(A,2)]`
+The color vector for dense matrix and triangular matrix is simply
+`[1,2,3,...,size(A,2)]`
 """
 function matrix_colors(A::Union{Array,UpperTriangular,LowerTriangular})
     eachindex(1:size(A,2)) # Vector size matches number of rows
@@ -379,12 +379,48 @@ function matrix_colors(A::Union{Tridiagonal,SymTridiagonal})
     _cycle(1:3,size(A,2))
 end
 
+"""
+  lu_instance(A) -> lu_factorization_instance
+
+Return an instance of the LU factorization object with the correct type
+cheaply.
+"""
+function lu_instance(A::Matrix{T}) where T
+  noUnitT = typeof(zero(T))
+  luT = LinearAlgebra.lutype(noUnitT)
+  ipiv = Vector{LinearAlgebra.BlasInt}(undef, 0)
+  info = zero(LinearAlgebra.BlasInt)
+  return LU{luT}(similar(A, 0, 0), ipiv, info)
+end
+
+"""
+  lu_instance(a::Number) -> a
+
+Return the number.
+"""
+lu_instance(a::Number) = a
+
 function __init__()
 
+  @require SuiteSparse="4607b0f0-06f3-5cda-b6b1-a6196a1729e9" begin
+    lu_instance(jac_prototype::SparseMatrixCSC) = SuiteSparse.UMFPACK.UmfpackLU(Ptr{Cvoid}(), Ptr{Cvoid}(), 1, 1,
+                                                                                      jac_prototype.colptr[1:1],
+                                                                                      jac_prototype.rowval[1:1],
+                                                                                      jac_prototype.nzval[1:1],
+                                                                                      0)
+  end
+
   @require StaticArrays="90137ffa-7385-5640-81b9-e52037218182" begin
     ismutable(::Type{<:StaticArrays.StaticArray}) = false
     can_setindex(::Type{<:StaticArrays.StaticArray}) = false
     ismutable(::Type{<:StaticArrays.MArray}) = true
+    function lu_instance(_A::StaticArrays.StaticMatrix{N,N}) where {N}
+      A = StaticArrays.SArray(_A)
+      L = LowerTriangular(A)
+      U = UpperTriangular(A)
+      p = StaticArrays.SVector{N,Int}(1:N)
+      return StaticArrays.LU(L, U, p)
+    end
   end
 
   @require LabelledArrays="2ee39098-c373-598a-b85f-a56591580800" begin

test/runtests.jl

@@ -1,6 +1,6 @@
 using ArrayInterface, Test
 using Base: setindex
-import ArrayInterface: has_sparsestruct, findstructralnz, fast_scalar_indexing
+import ArrayInterface: has_sparsestruct, findstructralnz, fast_scalar_indexing, lu_instance
 @test ArrayInterface.ismutable(rand(3))
 
 using StaticArrays
@@ -127,3 +127,16 @@ rowind,colind=findstructralnz(BBB)
         @test all(isone, y2)
     end
 end
+
+using SuiteSparse
+@testset "lu_instance" begin
+  for A in [
+    randn(5, 5),
+    @SMatrix(randn(5, 5)),
+    @MMatrix(randn(5, 5)),
+    sprand(50, 50, 0.5)
+  ]
+    @test lu_instance(A) isa typeof(lu(A))
+  end
+  @test lu_instance(1) === 1
+end