I am trying to implement mkl two-stage algorithm in Julia. “mkl_sparse_?_create_csc” function creates a handle for a matrix A in CSC format. In c language, the type of the handle is sparse_matrix_t. How to wrap this kind of type in Julia ccall?
Have you looked at the MKLSparse package, https://github.com/JuliaSparse/MKLSparse.jl, to see how things are done there? This was primarily wrapping an older version of the MKL sparse matrix calls that didn’t use the two-stage algorithms but it may give you some clue as how to start.
I believe that most of the two-stage algorithms return opaque pointers from functions like mkl_sparse_?_create_csc and you probably will need to manipulate them as the Julia equivalent of void *
2 Likes
According to your prescription, I have created a pointer of type Cvoid but my call giving me an error of no method matching. I am attaching my sample code here
I = [1, 4, 4, 5]; J = [1, 1, 2, 3]; V = [2.0, 1.0, 4.0, 6.0];
b=sparse(I,J,V);
p=Ptr{Cvoid}
ccall(("mkl_sparse_d_create_csc", librt),Cint (Ptr{Cvoid},Ptr{UInt8},
Ref{Cint},Ref{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Float64}),
p,"SPARSE_INDEX_BASE_ONE",b.m,b.n,
b.colptr,pointer(b.colptr,2),b.rowval,b.nzval)
MethodError: no method matching unsafe_convert(::Type{Ptr{Nothing}}, ::Type{Ptr{Nothing}})
Closest candidates are:
unsafe_convert(::Type{Ptr{Nothing}}, ::Base.RefValue{T}) where T at refvalue.jl:30
unsafe_convert(::Type{Ptr{Nothing}}, ::Base.RefArray{T,A,R} where R where A<:(AbstractArray{T,N} where N)) where T at refpointer.jl:90
unsafe_convert(::Type{Ptr{Nothing}}, ::Base.CFunction) at c.jl:36
I also looked into the python wrapper sparse_dot_mkl. I have tried to implement their logic for creating the matrix handle for the MKL routine, mkl_sparse_?_create_csc gives me the error -“The input parameters contain an invalid value”. I am showing a part of my code.
mutable struct sparse_matrix_t
end
I = [1, 4, 4, 5]; J = [1, 1, 2, 3]; V = [2.0, 1.0, 4.0, 6.0];
b=sparse(I,J,V);
d=sparse_matrix_t();
p=pointer_from_objref(d)
ccall(("mkl_sparse_d_create_csc", librt),Cint (Ptr{Cvoid},Ptr{UInt8},
Ref{Cint},Ref{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Float64}),
p,"SPARSE_INDEX_BASE_ONE",b.m,b.n,
b.colptr,pointer(b.colptr,2),b.rowval,b.nzval);
the value of the pointer p does not change during the operation.
types like sparse_index_base_t are enums
/* sparse matrix indexing: C-style or Fortran-style */
typedef enum
{
SPARSE_INDEX_BASE_ZERO = 0, /* C-style */
SPARSE_INDEX_BASE_ONE = 1 /* Fortran-style */
} sparse_index_base_t;
and would not be passed as character strings. I think that p will need to be passed as a Ref{Ptr{Cvoid}} so that the pointer itself can be changed. I will experiment and see if I can get something working.
1 Like
I did successfully change the value of the pointer by declaring an array of Ptr{Cvoid} and passing that
julia> using Base.Enums, MKL_jll, SparseArrays
julia> @enum sparse_index_base_t SPARSE_INDEX_BASE_ZERO=0 SPARSE_INDEX_BASE_ONE=1
julia> m = sparse(collect(1:10), [5, 2, 6, 3, 9, 1, 8, 5, 7, 10], ones(10))
julia> p = [Ptr{Cvoid}(0)]
1-element Array{Ptr{Nothing},1}:
Ptr{Nothing} @0x0000000000000000
julia> ccall((:mkl_sparse_d_create_csc, libmkl_rt), Cint, (Ref{Ptr{Cvoid}}, sparse_index_base_t, Int64, Int64, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Cdouble}),
p, SPARSE_INDEX_BASE_ONE, m.m, m.n, m.colptr, pointer(m.colptr, 2), m.rowval, m.nzval)
0
julia> p
1-element Array{Ptr{Nothing},1}:
Ptr{Nothing} @0x000000000443f000
Creating an array of one element may be more complicated than necessary. There should be a way of using Ref and I still need to check that pointer points to a valid sparse_matrix for MKL.
1 Like
Thank you.
But why are you creating the array of Ptr{Cvoid} in the 1st place?
There needs to be a location in memory that Julia allocates and can access which is then passed to the compiled code. One way of doing that is to allocate an array. As I said, it is probably possible to do that with a Ref{Ptr{Cvoid}} but I haven’t worked out exactly how.
It isn’t important exactly how the memory address is declared and passed around because it is an opaque pointer. All that the mkl_spblas.h file says about it is
struct sparse_matrix;
typedef struct sparse_matrix *sparse_matrix_t;
However, you do need to pass a pointer to a memory location, not the value in that memory location.
1 Like
To follow up on this, you can create p as below using Ref and avoid creation of an array.
julia> p = Ref(Ptr{Cvoid}(0))
Base.RefValue{Ptr{Nothing}}(Ptr{Nothing} @0x0000000000000000)
julia> p[]
Ptr{Nothing} @0x0000000000000000
julia> ccall((:mkl_sparse_d_create_csc, libmkl_rt), Cint, (Ref{Ptr{Cvoid}}, sparse_index_base_t, Int64, Int64, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Cdouble}), p, SPARSE_INDEX_BASE_ONE, m.m, m.n, m.colptr, pointer(m.colptr, 2), m.rowval, m.nzval)0
julia> p[]
Ptr{Nothing} @0x0000000004965000
julia> B = Matrix(Float64.(I(10)));
julia> C = zeros(Float64, 10, 10);
julia> ccall((:mkl_sparse_d_mm, libmkl_rt), sparse_status_t, (sparse_operation_t, Cdouble, Ptr{Cvoid}, matrix_descr, sparse_layout_t, Ptr{Cdouble}, BlasInt, BlasInt, Cdouble, Ptr{Cdouble}, BlasInt, BlasInt), SPARSE_OPERATION_NON_TRANSPOSE, 1.0, p[], matrix_descr(SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_FULL, SPARSE_DIAG_NON_UNIT), SPARSE_LAYOUT_COLUMN_MAJOR, B, 10, 10, 0.0, C, 10, 10)
SPARSE_STATUS_SUCCESS::sparse_status_t = 0
julia> C
10×10 Array{Float64,2}:
0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
The various enums are lightly modified code from mkl_spblas.h
using LinearAlgebra
using MKL_jll
using Base.Enums
using LinearAlgebra: BlasInt
function __init__()
ccall((:MKL_Set_Interface_Layer, libmkl_rt), Cint, (Cint,), Base.USE_BLAS64 ? 1 : 0)
end
# Basic types and constants for inspector-executor SpBLAS API
@enum(
sparse_status_t,
SPARSE_STATUS_SUCCESS = 0,
SPARSE_STATUS_NOT_INITIALIZED = 1,
SPARSE_STATUS_ALLOC_FAILED = 2,
SPARSE_STATUS_INVALID_VALUE = 3,
SPARSE_STATUS_EXECUTION_FAILED = 4,
SPARSE_STATUS_INTERNAL_ERROR = 5,
SPARSE_STATUS_NOT_SUPPORTED = 6
)
@enum(
sparse_operation_t,
SPARSE_OPERATION_NON_TRANSPOSE = 10,
SPARSE_OPERATION_TRANSPOSE = 11,
SPARSE_OPERATION_CONJUGATE_TRANSPOSE = 12
)
@enum(
sparse_matrix_type_t,
SPARSE_MATRIX_TYPE_GENERAL = 20,
SPARSE_MATRIX_TYPE_SYMMETRIC = 21,
SPARSE_MATRIX_TYPE_HERMITIAN = 22,
SPARSE_MATRIX_TYPE_TRIANGULAR = 23,
SPARSE_MATRIX_TYPE_DIAGONAL = 24,
SPARSE_MATRIX_TYPE_BLOCK_TRIANGULAR = 25,
SPARSE_MATRIX_TYPE_BLOCK_DIAGONAL = 26
)
@enum(
sparse_index_base_t,
SPARSE_INDEX_BASE_ZERO = 0,
SPARSE_INDEX_BASE_ONE = 1
)
@enum(
sparse_fill_mode_t,
SPARSE_FILL_MODE_LOWER = 40,
SPARSE_FILL_MODE_UPPER = 41,
SPARSE_FILL_MODE_FULL = 42
)
@enum(
sparse_diag_type_t,
SPARSE_DIAG_NON_UNIT = 50,
SPARSE_DIAG_UNIT = 51
)
@enum(
sparse_layout_t,
SPARSE_LAYOUT_ROW_MAJOR = 101,
SPARSE_LAYOUT_COLUMN_MAJOR = 102
)
@enum(
verbose_mode_t,
SPARSE_VERBOSE_OFF = 70,
SPARSE_VERBOSE_BASIC = 71,
SPARSE_VERBOSE_EXTENDED = 72
)
@enum(
sparse_memory_usage_t,
SPARSE_MEMORY_NONE = 80,
SPARSE_MEMORY_AGGRESSIVE = 81
)
@enum(
sparse_request_t,
SPARSE_STAGE_FULL_MULT = 90,
SPARSE_STAGE_NNZ_COUNT = 91,
SPARSE_STAGE_FINALIZE_MULT = 92,
SPARSE_STAGE_FULL_MULT_NO_VAL = 93,
SPARSE_STAGE_FINALIZE_MULT_NO_VAL = 94
)
struct matrix_descr
type::sparse_matrix_type_t
mode::sparse_fill_mode_t
diag::sparse_diag_type_t
end
1 Like
Previously I used p = Ref(Ptr{Cvoid}(0)) but I missed the actual value of sparse_index_base_t instead I used it as a character string. Thank you for your kind help.
can you please explain what are the difference of Ref(Ptr{Cvoid}), Ref(Ptr{Cvoid}(0)) and
pointer_from_objref(Ref(Ptr{Cvoid})) ? I can guess that Ref(Ptr{Cvoid}) referencing to the memory address of the *void pointer, while Ptr may or may not give the actual memory address.
In the mkl_spblas.h there are 11 input arguments for mkl_sparse_d_mm routine but in your example, you are giving 12 parameters as input arguments. Also for the non-symmetric matrix multiplication, I am getting the wrong answer from the mkl_sparse_d_mm routine.
julia> using Base.Enums, MKL_jll, SparseArrays
julia> I = [1, 4, 4, 5]; J = [1, 1, 2, 3]; V = [2.0, 1.0, 4.0, 6.0];
julia> m=sparse(I,J,V);
julia> p = Ref(Ptr{Cvoid}(0));
julia> ccall((:mkl_sparse_d_create_csc, libmkl_rt), Cint, (Ref{Ptr{Cvoid}}, sparse_index_base_t, Int64, Int64, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Cdouble}), p, SPARSE_INDEX_BASE_ONE, m.m, m.n, m.colptr, pointer(m.colptr, 2), m.rowval, m.nzval);
julia> B=rand(3,2)
3×2 Array{Float64,2}:
0.116917 0.206817
0.407721 0.265923
0.539494 0.77058
julia> C=zeros(5,2);
julia> ccall((:mkl_sparse_d_mm, libmkl_rt), sparse_status_t, (sparse_operation_t, Cdouble, Ptr{Cvoid}, matrix_descr, sparse_layout_t, Ptr{Cdouble}, Int64, Int64, Cdouble, Ptr{Cdouble}, Int64), SPARSE_OPERATION_NON_TRANSPOSE, 1.0, p[], matrix_descr(SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_FULL, SPARSE_DIAG_NON_UNIT), SPARSE_LAYOUT_COLUMN_MAJOR, B, 2, 3, 0.0, C, 5)
SPARSE_STATUS_SUCCESS::sparse_status_t = 0
julia> C
5×2 Array{Float64,2}:
0.233834 0.413634
0.0 0.0
0.0 0.0
2.15798 3.08232
0.206817 0.0
julia> m*B
5×2 Array{Float64,2}:
0.233834 0.413634
0.0 0.0
0.0 0.0
1.7478 1.27051
3.23696 4.62348
In c language, I am getting the correct answer for the same parameter values.
Ref{Ptr{Cvoid}} is a type, as is Ptr{Cvoid}. The expression Ptr{Cvoid}(0) is a call to the default constructor function for this type. It creates what is generally called a NULL pointer in C. That is, it is an address of an unspecified type that has been initialized to the 0 location in memory.
If I take a Ref to such an object it stores it in memory. So the sequence
julia> p = Ref(Ptr{Cvoid}(0))
Base.RefValue{Ptr{Nothing}}(Ptr{Nothing} @0x0000000000000000)
stores a zero address in memory and returns, in a Julia-accessible way, the (non-zero) address of this zero address. When that is passed in the call to :mkl_spars_d_create_csc the routine being called can allocate memory for the sparse_matrix_t struct and overwrite this location with the address of the new struct. That’s why we have to pass a pointer to a pointer or **void in C. It happens to be declared as *sparse_matrix_t but sparse_matrix_t is itself a pointer type so this is a pointer to a pointer.
Before calling the C function to initialize the storage the address in that location was
julia> p[]
Ptr{Nothing} @0x0000000000000000
because that is how it was initialized in Julia. It’s a convention that 0 is always an invalid address so trying to dereference a NULL pointer in C throws a memory error. After the call we have
julia> p[]
Ptr{Nothing} @0x0000000004965000
which is the location of the memory that was allocated to hold the sparse_matrix struct. Unfortunately, we don’t know what the internal structure of a sparse_matrix is so the only thing we can do with it is to pass around a pointer to the struct, which is what a sparse_matrix_t is.
For an opaque pointer like this, on the Julia side it is irrelevant how it is described as long as it is the right size to be an address. We are calling it a Ptr{Cvoid} just to remind ourselves that we are using it as an address to an opaque C struct.
You’re correct. I misread the declaration of the C function in the header file and appended an MKL_INT argument. Luckily for me, putting an extra argument on the stack doesn’t cause problems.
Regarding the results, can you create an example that does not use a random matrix like B?
I am just copping your example of mkl_sparse_d_mm with the same matrix input. But my result is quite different,
julia> m = sparse(collect(1:10), [5, 2, 6, 3, 9, 1, 8, 5, 7, 10], ones(10));
julia> p = Ref(Ptr{Cvoid}(0));
julia> ccall(("mkl_sparse_d_create_csc", librt), Cint, (Ref{Ptr{Cvoid}}, sparse_index_base_t, Int64, Int64, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Cdouble}), p, SPARSE_INDEX_BASE_ONE, m.m, m.n, m.colptr, pointer(m.colptr, 2), m.rowval, m.nzval);
julia> B = Matrix(Float64.(I(10)));
julia> C = zeros(Float64, 10, 10);
julia> ccall(("mkl_sparse_d_mm", librt), sparse_status_t, (sparse_operation_t, Cdouble, Ptr{Cvoid}, matrix_descr, sparse_layout_t, Ptr{Cdouble}, BlasInt, BlasInt, Cdouble, Ptr{Cdouble}, BlasInt), SPARSE_OPERATION_NON_TRANSPOSE, 1.0, p[], matrix_descr(SPARSE_MATRIX_TYPE_GENERAL, SPARSE_FILL_MODE_FULL, SPARSE_DIAG_NON_UNIT), SPARSE_LAYOUT_COLUMN_MAJOR, B, 10, 10, 0.0, C, 10)
SPARSE_STATUS_SUCCESS::sparse_status_t = 0
julia> C
10×10 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
Sorry but I can’t reproduce your results. I get
julia> C
10×10 Array{Float64,2}:
0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
julia> C == Matrix(m)
true
Can you check that the BLAS is using 64-bit ints?
julia> LinearAlgebra.BlasInt
Int64
and that you have initialized the interface layer with, e.g.
ccall((:MKL_Set_Interface_Layer, libmkl_rt), Cint, (Cint,), Base.USE_BLAS64 ? 1 : 0)
Now it is working fine after initializing the MKL interface level. I thought interface initialization is not mandatory. Does it replace the default openblas with mkl? Can you please explain the purpose of the ccall((:MKL_Set_Interface_Layer, libmkl_rt), Cint, (Cint,), Base.USE_BLAS64 ? 1 : 0) routine. I have seen it in the MKLSparse.jl but I can’t understand it.
Regarding replacing the default OpenBLAS with MKL, see the (as yet unregistered) package https://github.com/JuliaComputing/MKL.jl I have found it works well with Julia-1.5.2 on my Ubuntu system but fails in the precompilation on Julia-1.6.0-DEV.
I did just copy the ccall from https://github.com/JuliaSparse/MKLSparse.jl assuming, from the way the argument is constructed, that it would distinguish between a 32-bit BlasInt and a 64-bit BlasInt. However, when I look up the documentation for MKL_Set_Interface_Layer (something you could do too) that doesn’t seem to be what its purpose is and I haven’t been able to reconcile the values 0 and 1 for the arguments with the documentation. In the documentation the call seems to set up linkage conventions for different compilers.
1 Like