search

JuliaSparse / SuiteSparseGraphBLAS.jl

Sparse, General Linear Algebra for Graphs!
MIT License
103 stars 17 forks source link

Extension: GxB Matrix import/export CSC #13

Closed fcdimitr closed 3 years ago

fcdimitr commented 3 years ago

I am trying to extend your wrapper to provide additional functionality, specifically, to be able to pass CSC matrices seamlessly between Julia and C (without going through the extractTuples and build functions).

The GxB_Matrix_import_CSC (and export respectively) allows us to work with pointers to the CSC vectors. See https://github.com/DrTimothyAldenDavis/GraphBLAS/blob/stable/Doc/GraphBLAS_UserGuide.pdf page 98, for example.

My initial try is the following:

function GxB_Matrix_import_CSC(Aj)

  A     = GrB_Matrix{Float64}()
  A_ptr = pointer_from_objref(A)

  type = GrB_FP64;

  nr = convert( UInt64, size(Aj,1) )
  nc = convert( UInt64, size(Aj,2) )
  nv = convert( UInt64, SparseArrays.nnz( Aj ) )

  ne = convert( Int64, -1 )

  Ap = Ref( Ptr{Cvoid}(pointer( ZeroBasedIndices( Aj.colptr.-1 ) ) ) )
  Ai = Ref( Ptr{Cvoid}(pointer( ZeroBasedIndices( Aj.rowval.-1 ) ) ) )
  X  = Ref( Ptr{Cvoid}(pointer( Aj.nzval                         ) ) )

  desc = GrB_NULL

  return GrB_Info(
    ccall(
      dlsym(SuiteSparseGraphBLAS.graphblas_lib, "GxB_Matrix_import_CSC"),
      Cint,
      # ----- INPUT SPECIFICATION
      (Ptr{Cvoid},              # GrB_Matrix *A
       Ptr{Cvoid},              # GrB_Type   type
       UInt64,                  # GrB_Index  nrows
       UInt64,                  # GrB_Index  ncols
       UInt64,                  # GrB_Index  nvals
       Int64,                   # int64_t    nonempty
       Ptr{Ptr{Cvoid}},         # GrB_Index  **Ap
       Ptr{Ptr{Cvoid}},         # GrB_Index  **Ai
       Ptr{Ptr{Cvoid}},         # void       **Ax
       Ptr{Cvoid}),             # (descriptor, unused)
      # ----- INPUT PARAMS
      A_ptr,                    # GrB_Matrix *A
      type.p,                   # GrB_Type   type
      nr,                       # GrB_Index  nrows
      nc,                       # GrB_Index  ncols
      nv,                       # GrB_Index  nvals
      ne,                       # int64_t    nonempty
      Ap,                       # GrB_Index  **Ap
      Ai,                       # GrB_Index  **Ai
      X,                        # void       **Ax
      desc.p                    # (descriptor, unused)
    )
  )

end

I am working with Double pointers in this example. When running the following code

Aj = sprand(Float64,10,10,0.1);
GxB_Matrix_import_CSC(Aj);

I get a segmentation violation:


signal (11): Segmentation fault: 11
in expression starting at REPL[10]:1
GxB_Matrix_import_CSC at ...04f00df9359a326cadce4d8a1c49ccbecc391293/lib/libgraphblas.3.1.1.dylib (unknown line)
GxB_Matrix_import_CSC at ...graphblastest.jl:86
unknown function (ip: 0x107052745)
jl_apply at /Users/julia/buildbot/worker/package_macos64/build/src/./julia.h:1690 [inlined]
do_call at /Users/julia/buildbot/worker/package_macos64/build/src/interpreter.c:117
eval_body at /Users/julia/buildbot/worker/package_macos64/build/src/interpreter.c:0
jl_interpret_toplevel_thunk at /Users/julia/buildbot/worker/package_macos64/build/src/interpreter.c:660
jl_toplevel_eval_flex at /Users/julia/buildbot/worker/package_macos64/build/src/toplevel.c:840
jl_toplevel_eval_flex at /Users/julia/buildbot/worker/package_macos64/build/src/toplevel.c:790
jl_toplevel_eval at /Users/julia/buildbot/worker/package_macos64/build/src/toplevel.c:849 [inlined]
jl_toplevel_eval_in at /Users/julia/buildbot/worker/package_macos64/build/src/toplevel.c:883
eval at ./boot.jl:331
eval_user_input at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:134
repl_backend_loop at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:195
start_repl_backend at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:180
#run_repl#37 at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:292
run_repl at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288
#807 at ./client.jl:399
jfptr_YY.807_70113.clone_1 at /Applications/Julia-1.5.app/Contents/Resources/julia/lib/julia/sys.dylib (unknown line)
jl_apply at /Users/julia/buildbot/worker/package_macos64/build/src/./julia.h:1690 [inlined]
do_apply at /Users/julia/buildbot/worker/package_macos64/build/src/builtins.c:655
jl_f__apply at /Users/julia/buildbot/worker/package_macos64/build/src/builtins.c:669 [inlined]
jl_f__apply_latest at /Users/julia/buildbot/worker/package_macos64/build/src/builtins.c:705
#invokelatest#1 at ./essentials.jl:710 [inlined]
invokelatest at ./essentials.jl:709 [inlined]
run_main_repl at ./client.jl:383
exec_options at ./client.jl:313
_start at ./client.jl:506
jfptr__start_57081.clone_1 at /Applications/Julia-1.5.app/Contents/Resources/julia/lib/julia/sys.dylib (unknown line)
true_main at /usr/local/bin/julia (unknown line)
main at /usr/local/bin/julia (unknown line)
Allocations: 2317213 (Pool: 2316562; Big: 651); GC: 3

Any ideas what I am doing wrong?

abhinavmehndiratta commented 3 years ago

Hi @fcdimitr you need to wrap the values in a struct here you go -

using SparseArrays

struct foo
    x::Float64
end

function GxB_Matrix_import_CSC(Aj)
  A = GrB_Matrix{Float64}()
  A_ptr = pointer_from_objref(A)
  type = GrB_FP64
  nr = convert(UInt64, size(Aj, 1))
  nc = convert(UInt64, size(Aj, 2))
  nv = convert(UInt64, SparseArrays.nnz(Aj))
  ne = convert(Int64, -1)

  Ap = Ref(pointer(ZeroBasedIndices(Aj.colptr.-1)))
  Ai = Ref(pointer(ZeroBasedIndices(Aj.rowval.-1 ))) 
  X  = Ref(pointer(foo.(Aj.nzval)))

  desc = GrB_NULL

  g = GrB_Info(
    ccall(
      dlsym(graphblas_lib, "GxB_Matrix_import_CSC"),
      Cint,
      # ----- INPUT SPECIFICATION
      (Ptr{Cvoid},              # GrB_Matrix *A
       Ptr{Cvoid},              # GrB_Type   type
       UInt64,                  # GrB_Index  nrows
       UInt64,                  # GrB_Index  ncols
       UInt64,                  # GrB_Index  nvals
       Int64,                   # int64_t    nonempty
       Ptr{Cvoid},              # GrB_Index  **Ap
       Ptr{Cvoid},              # GrB_Index  **Ai
       Ptr{Cvoid},              # void       **Ax
       Ptr{Cvoid}),             # (descriptor, unused)
      # ----- INPUT PARAMS
      A_ptr,                    # GrB_Matrix *A
      type.p,                   # GrB_Type   type
      nr,                       # GrB_Index  nrows
      nc,                       # GrB_Index  ncols
      nv,                       # GrB_Index  nvals
      ne,                       # int64_t    nonempty
      Ap,                       # GrB_Index  **Ap
      Ai,                       # GrB_Index  **Ai
      X,                        # void       **Ax
      desc.p                    # (descriptor, unused)
    )
  )
  g != GrB_SUCCESS && return g
  return A
end

Since the C function expects a pointer to a pointer, this will work.

julia>  using GraphBLASInterface, SuiteSparseGraphBLAS

julia> GrB_init(GrB_NONBLOCKING)
GrB_SUCCESS::GrB_Info = 0

julia> using SparseArrays

julia> Aj = sprand(Int64,10,10,0.1);

julia> f = GxB_Matrix_import_CSC(Aj);

julia> @GxB_fprint(f, GxB_COMPLETE)

  10x10 GraphBLAS double matrix, sparse by col:
  f, 12 entries

    (6,0)    -6.51981090010754e+18
    (2,1)    2.38012960675141e+18
    (2,2)    7.15625156581955e+18
    (6,2)    -8.4509572234965e+17
    (0,5)    -5.14419131400837e+18
    (7,6)    -8.76743228265509e+18
    (6,7)    -4.87354216758125e+18
    (8,7)    -5.25682341711608e+18
    (5,8)    -8.96788900640334e+18
    (9,8)    6.97693687580616e+18
    (4,9)    -1.12098026882762e+18
    (8,9)    -9.16149892069805e+18

julia> Aj
10×10 SparseMatrixCSC{Int64,Int64} with 12 stored entries:
  [7 ,  1]  =  -6519810900107535143
  [3 ,  2]  =  2380129606751411208
  [3 ,  3]  =  7156251565819545931
  [7 ,  3]  =  -845095722349649729
  [1 ,  6]  =  -5144191314008365360
  [8 ,  7]  =  -8767432282655090949
  [7 ,  8]  =  -4873542167581250904
  [9 ,  8]  =  -5256823417116077675
  [6 ,  9]  =  -8967889006403336730
  [10,  9]  =  6976936875806161528
  [5 , 10]  =  -1120980268827621848
  [9 , 10]  =  -9161498920698047516
abhinavmehndiratta commented 3 years ago

If you want, you can check out this fork https://github.com/cvdlab/SuiteSparseGraphBLAS.jl, they've built an interface on top of this library.

fcdimitr commented 3 years ago

Thank you very much @abhinavmehndiratta, both for the solution, and for the link to the updated fork. I am closing the issue.