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 *`

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.

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.

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.

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
```

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.