asArray and misc changes

Author: Wimmerer

Project.toml

@@ -4,7 +4,6 @@ authors = ["Will Kimmerer <[email protected]", "Abhinav Mehndiratta <abhina
 version = "0.4.0"
 
 [deps]
-CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
@@ -15,9 +14,9 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SSGraphBLAS_jll = "7ed9a814-9cab-54e9-8e9e-d9e95b4d61b1"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 
 [compat]
-CEnum = "0.4"
 MacroTools = "0.5"
 SSGraphBLAS_jll = "5.1"
 julia = "1.6"

README.md

@@ -3,12 +3,6 @@
 # SuiteSparseGraphBLAS.jl
 A fast, general sparse linear algebra and graph computation package, based on SuiteSparse:GraphBLAS.
 
-## v0.5
-v0.5 is planned to release in August, after the entire GraphBLAS interface is available and some memory pressure issues have been sorted out. Check back then for more benchmarks, a more Julian interface, automatic differentiation support with ChainRules.jl and better integration with the wider ecosystem!
-
-If you're fine with sharp edges then give v0.4 a try, and let me know of any issues you find.
-The docs should provide enough information to run virtually the entire set of GraphBLAS functions.
-
 ### Installation
 ```julia
 using Pkg

benchmarks/Manifest.toml

@@ -107,6 +107,12 @@ version = "0.4.2"
 deps = ["Random"]
 uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
 
+[[IntelOpenMP_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c"
+uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
+version = "2018.0.3+2"
+
 [[InteractiveUtils]]
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
@@ -150,6 +156,10 @@ git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 version = "0.21.2"
 
+[[LazyArtifacts]]
+deps = ["Artifacts", "Pkg"]
+uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
+
 [[LibCURL]]
 deps = ["LibCURL_jll", "MozillaCACerts_jll"]
 uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
@@ -182,6 +192,18 @@ version = "0.3.4"
 [[Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
+[[MKLSparse]]
+deps = ["LinearAlgebra", "MKL_jll", "SparseArrays"]
+git-tree-sha1 = "c5fd60128db7d4e0ad2358d1017d293c81a49335"
+uuid = "0c723cd3-b8cd-5d40-b370-ba682dde9aae"
+version = "1.1.0"
+
+[[MKL_jll]]
+deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
+git-tree-sha1 = "5455aef09b40e5020e1520f551fa3135040d4ed0"
+uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
+version = "2021.1.1+2"
+
 [[MacroTools]]
 deps = ["Markdown", "Random"]
 git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"

benchmarks/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+MKLSparse = "0c723cd3-b8cd-5d40-b370-ba682dde9aae"
 MatrixMarket = "4d4711f2-db25-561a-b6b3-d35e7d4047d3"
 SuiteSparseGraphBLAS = "c2e53296-7b14-11e9-1210-bddfa8111e1d"
 SuiteSparseMatrixCollection = "ac199af8-68bc-55b8-82c4-7abd6f96ed98"

benchmarks/benchmarks.jl

@@ -26,12 +26,17 @@ using BenchmarkTools
 using SparseArrays
 using LinearAlgebra
 
+
 #OPTIONS SET 1:
 # Maximum number of samples taken for each benchmark
 BenchmarkTools.DEFAULT_PARAMETERS.samples = 10
 # Total amount of time allowed for each benchmark, minimum of 1 sample taken.
 BenchmarkTools.DEFAULT_PARAMETERS.seconds = 60
 
+# Comment or uncomment this line to disable or enable MKLSparse respectively.
+# This will only work for SpMM and SpMV and only operates on CSC.
+using MKLSparse
+
 # Change this to change the size of the dense RHS of csrtimesfull and csctimesfull
 const sizefullrhs = [1,2,4]
 
@@ -288,6 +293,8 @@ if length(ARGS) != 0
         if splitext(ARGS[1])[2] == ".mtx"
             singlebench(ARGS[1])
         else
+            lines = readlines(ARGS[1])
+            filter!((x) -> !occursin("#", x), lines)
             singlebench.(readlines(ARGS[1]))
         end
     elseif tryparse(Int64, ARGS[1]) !== nothing

benchmarks/do10

@@ -3,13 +3,4 @@
 # FIXME: the paths would need to point somewhere else.
 # better yet, the sparse matrix IDs should be used so
 # they can be directly downloaded from sparse.tamu.edu.
-julia benchmarks.jl /home/faculty/d/davis/Freescale2.mtx
-julia benchmarks.jl /home/faculty/d/davis/indochina-2004.mtx
-julia benchmarks.jl /home/faculty/d/davis/nd24k.mtx
-julia benchmarks.jl /home/faculty/d/davis/nd6k.mtx
-julia benchmarks.jl /home/faculty/d/davis/roadNet-CA.mtx
-julia benchmarks.jl /raid/GAP/GAP-twitter/GAP-twitter.mtx
-julia benchmarks.jl /raid/GAP/GAP-road/GAP-road.mtx
-julia benchmarks.jl /raid/GAP/GAP-web/GAP-web.mtx
-julia benchmarks.jl /raid/GAP/GAP-urand/GAP-urand.mtx
-julia benchmarks.jl /raid/GAP/GAP-kron/GAP-kron.mtx
+julia benchmarks.jl mtx.txt

benchmarks/mtx.txt

@@ -0,0 +1,12 @@
+# You may provide paths or SuiteSparse ID numbers here.
+# A hashtag denotes a commented line.
+/home/faculty/d/davis/Freescale2.mtx
+/home/faculty/d/davis/indochina-2004.mtx
+/home/faculty/d/davis/nd24k.mtx
+/home/faculty/d/davis/nd6k.mtx
+/home/faculty/d/davis/roadNet-CA.mtx
+/raid/GAP/GAP-twitter/GAP-twitter.mtx
+/raid/GAP/GAP-road/GAP-road.mtx
+/raid/GAP/GAP-web/GAP-web.mtx
+/raid/GAP/GAP-urand/GAP-urand.mtx
+/raid/GAP/GAP-kron/GAP-kron.mtx

src/SuiteSparseGraphBLAS.jl

@@ -19,7 +19,6 @@ using SparseArrays: nonzeroinds
 using MacroTools
 using LinearAlgebra
 using Random: randsubseq, default_rng, AbstractRNG, GLOBAL_RNG
-using CEnum
 using SpecialFunctions: lgamma, gamma, erf, erfc
 using Base.Broadcast
 include("abstracts.jl")
@@ -76,6 +75,8 @@ include("operations/resize.jl")
 include("print.jl")
 include("import.jl")
 include("export.jl")
+include("pack.jl")
+include("unpack.jl")
 include("options.jl")
 #EXPERIMENTAL
 include("operations/argminmax.jl")
@@ -89,6 +90,7 @@ include("chainrules/selectrules.jl")
 include("chainrules/constructorrules.jl")
 #include("random.jl")
 include("misc.jl")
+include("asjulia.jl")
 export libgb
 export UnaryOps, BinaryOps, Monoids, Semirings #Submodules
 export UnaryOp, BinaryOp, Monoid, Semiring #UDFs
@@ -130,6 +132,7 @@ function __init__()
     # In the future this should hopefully allow us to do no-copy passing of arrays between Julia and SS:GrB.
     # In the meantime it helps Julia respond to memory pressure from SS:GrB and finalize things in a timely fashion.
     libgb.GxB_init(libgb.GrB_NONBLOCKING, cglobal(:jl_malloc), cglobal(:jl_calloc), cglobal(:jl_realloc), cglobal(:jl_free), true)
+    gbset(:nthreads, Sys.CPU_THREADS ÷ 2)
     # Eagerly load selectops constants.
     _loadselectops()
     # Set printing done by SuiteSparse:GraphBLAS to base-1 rather than base-0.

src/asjulia.jl

@@ -0,0 +1,100 @@
+function asArray(f::Function, A::GBVecOrMat{T}; dropzeros=false, freeunpacked=false) where {T}
+    if gbget(A, SPARSITY_STATUS) != GBDENSE
+        X = similar(A)
+        if X isa GBVector
+            X[:] = zero(T)
+        else
+            X[:,:] = zero(T)
+        end
+        #I don't like this, it defeats the purpose of this method, which is to make no copies.
+        # But somehow maintaining the input A in its original form is key to the to_vec implementation
+        # for ChainRules. Temporarily it's fine, it's no worse than it originally was.
+        # TODO: fix this issue with the ChainRules code.
+        A = eadd(X, A)
+    end
+    array = _unpackdensematrix!(A)
+    result = try
+        f(array, A)
+    finally
+        if freeunpacked
+            ccall(:jl_free, Cvoid, (Ptr{T},), pointer(array))
+        else
+            _packdensematrix!(A, array)
+            if dropzeros
+                select!(nonzeros, A)
+            end
+        end
+    end
+    return result
+end
+
+function asSparseMatrixCSC(f::Function, A::GBMatrix{T}; freeunpacked=false) where {T}
+    colptr, rowidx, values =  _unpackcscmatrix!(A)
+    array = SparseMatrixCSC{T, libgb.GrB_Index}(size(A, 1), size(A, 2), colptr, rowidx, values)
+    result = try
+        f(array, A)
+    finally
+        if freeunpacked
+            ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), pointer(colptr))
+            ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), pointer(rowidx))
+            ccall(:jl_free, Cvoid, (Ptr{T},), pointer(values))
+        else
+            _packcscmatrix!(A, colptr, rowidx, values)
+        end
+    end
+    return result
+end
+
+function asSparseVector(f::Function, A::GBVector{T}; freeunpacked=false) where {T}
+    colptr, rowidx, values =  _unpackcscmatrix!(A)
+    vector = SparseVector{T, libgb.GrB_Index}(size(A, 1), rowidx, values)
+    result = try
+        f(vector, A)
+    finally
+        if freeunpacked
+            ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), pointer(colptr))
+            ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), pointer(rowidx))
+            ccall(:jl_free, Cvoid, (Ptr{T},), pointer(values))
+        else
+            _packcscmatrix!(A, colptr, rowidx, values)
+        end
+    end
+    return result
+end
+
+
+function Base.Matrix(A::GBMatrix)
+    return asArray(A) do arr, _
+        return copy(arr)
+    end
+end
+
+function Matrix!(A::GBMatrix)
+    return asArray(A; freeunpacked=true) do arr, _
+        return copy(arr)
+    end
+end
+
+function Base.Vector(v::GBVector)
+    return asArray(v) do vec, _
+        return copy(vec)
+    end
+end
+
+function Vector!(v::GBVector)
+    return asArray(v; freeunpacked=true) do vec, _
+        return copy(vec)
+    end
+end
+
+function SparseArrays.SparseMatrixCSC(A::GBMatrix)
+    return asArray(A) do arr, _
+        return copy(arr)
+    end
+end
+
+function SparseArrays.SparseVector(v::GBVector)
+    return asArray(v) do arr, _
+        return copy(arr)
+    end
+end

src/chainrules/chainruleutils.jl

@@ -25,24 +25,3 @@ end
 
 # LinearAlgebra.norm doesn't like the nothings.
 LinearAlgebra.norm(A::GBArray, p::Real=2) = norm(nonzeros(A), p)
-
-# Broadcast b into the rows of A. WARNING: THIS DOES NOT MATCH JULIA.
-function broadcast_emul!(C, A, b, op; mask = nothing, accum = nothing, desc = nothing)
-    B = diagm(b)
-    mul!(C, A, B, (any, op); mask, accum, desc)
-    return C
-end
-
-function broadcast_emul(A, b, op; mask = nothing, accum = nothing, desc = nothing)
-    B = diagm(b)
-    mul(A, B, (any, op); mask, accum, desc)
-end
-
-function find_k(A, B::GBArray, op, minmax; mask = nothing, accum = nothing, desc = nothing)
-    K = []
-    for col ∈ axes(B, 2)
-        intermediate = broadcast_emul(A, B[:, col], op; mask, accum, desc)
-        push!(K, argminmax(intermediate, minmax, 2))
-    end
-    return hcat(K...)
-end