C# calling Intel MKL cblas_dgemm_batch

Question

I can call Intel MKL cblas_dgem from C#, see the following code:

[DllImport("custom_mkl", CallingConvention = CallingConvention.Cdecl, ExactSpelling = true, SetLastError = false)]
internal static extern void cblas_dgemm(
    int Order, int TransA, int TransB, MKL_INT M, MKL_INT N, MKL_INT K,
    double alpha, [In] double[,] A, MKL_INT lda, [In] double[,] B, MKL_INT ldb,
    double beta, [In, Out] double[,] C, MKL_INT ldc);

and

void cblas_dgemm (const CBLAS_LAYOUT Layout, const CBLAS_TRANSPOSE transa, const CBLAS_TRANSPOSE transb, const MKL_INT m, const MKL_INT n, const MKL_INT k, const double alpha, const double *a, const MKL_INT lda, const double *b, const MKL_INT ldb, const double beta, double *c, const MKL_INT ldc);

But I'm not able to call cblas_dgemm_batch from C#, see the following code:

[DllImport("custom_mkl", CallingConvention = CallingConvention.Cdecl, ExactSpelling = true, SetLastError = false)] // not working
internal static extern void cblas_dgemm_batch(
    int Layout, [In] int[] transa_array, [In] int[] transb_array, [In] MKL_INT[] m_array, [In] MKL_INT[] n_array, [In] MKL_INT[]  k_array, 
    [In] double[] alpha_array, [In] double[][,] a_array, [In] MKL_INT[] lda_array, [In] double[][,] b_array, [In] MKL_INT[] ldb_array,
    [In] double[] beta_array, [In, Out] double[][,] c_array, [In] MKL_INT[] ldc_array, MKL_INT group_count, [In] MKL_INT[] group_size);

and

void cblas_dgemm_batch (const CBLAS_LAYOUT Layout, const CBLAS_TRANSPOSE* transa_array, const CBLAS_TRANSPOSE* transb_array, const MKL_INT* m_array, const MKL_INT* n_array, const MKL_INT* k_array, const double* alpha_array, const double **a_array, const MKL_INT* lda_array, const double **b_array, const MKL_INT* ldb_array, const double* beta_array, double **c_array, const MKL_INT* ldc_array, const MKL_INT group_count, const MKL_INT* group_size);

I'm getting the following error message:

• System.Runtime.InteropServices.MarshalDirectiveException
• Cannot marshal 'parameter #8': There is no marshaling support for nested arrays.

I can understand that the problem are the nested array parameters. This parameter should be array of pointers to arrays. But how can I call cblas_dgemm_batch from C#?

Answer 1

Using the following custom marshaler for the jagged arrays is the solution:

class JaggedArrayMarshaler : ICustomMarshaler
{
    static ICustomMarshaler GetInstance(string cookie)
    {
        return new JaggedArrayMarshaler();
    }
    GCHandle[] handles;
    GCHandle buffer;
    Array[] array;
    public void CleanUpManagedData(object ManagedObj)
    {
    }
    public void CleanUpNativeData(IntPtr pNativeData)
    {
        buffer.Free();
        foreach (GCHandle handle in handles) handle.Free();
    }
    public int GetNativeDataSize()
    {
        return IntPtr.Size;
    }
    public IntPtr MarshalManagedToNative(object ManagedObj)
    {
        array = (Array[])ManagedObj;
        handles = new GCHandle[array.Length];
        for (int i = 0; i < array.Length; i++)
            handles[i] = GCHandle.Alloc(array[i], GCHandleType.Pinned);
        IntPtr[] pointers = new IntPtr[handles.Length];
        for (int i = 0; i < handles.Length; i++)
            pointers[i] = handles[i].AddrOfPinnedObject();
        buffer = GCHandle.Alloc(pointers, GCHandleType.Pinned);
        return buffer.AddrOfPinnedObject();
    }
    public object MarshalNativeToManaged(IntPtr pNativeData)
    {
        return array;
    }
}

and using the above marshaler:

[DllImport("custom_mkl", CallingConvention = CallingConvention.Cdecl, ExactSpelling = true, SetLastError = false)]
internal static extern void cblas_dgemm_batch(
    int Layout, [In] int[] transa_array, [In] int[] transb_array, [In] MKL_INT[] m_array, [In] MKL_INT[] n_array, [In] MKL_INT[] k_array,
    [In] double[] alpha_array, 
    [MarshalAs(UnmanagedType.CustomMarshaler, MarshalTypeRef = typeof(JaggedArrayMarshaler))][In] double[][,] a_array, [In] MKL_INT[] lda_array, 
    [MarshalAs(UnmanagedType.CustomMarshaler, MarshalTypeRef = typeof(JaggedArrayMarshaler))][In] double[][,] b_array, [In] MKL_INT[] ldb_array,
    [In] double[] beta_array, 
    [MarshalAs(UnmanagedType.CustomMarshaler, MarshalTypeRef = typeof(JaggedArrayMarshaler))][In, Out] double[][,] c_array, 
    [In] MKL_INT[] ldc_array, MKL_INT group_count, [In] MKL_INT[] group_size);

I'm using the following code to test it:

public static double[][,] Dot(double[][,] a, double[][,] b)
{
    int n0 = a.Length;
    if (b.Length != n0) throw new System.Exception("Group size must be the same");
    int Order = 101; // row-major arrays
    int[] TransA = new int[n0];
    int[] TransB = new int[n0];
    MKL_INT[] M = new MKL_INT[n0];
    MKL_INT[] N = new MKL_INT[n0];
    MKL_INT[] K = new MKL_INT[n0];
    double[] alpha = new double[n0];
    double[] beta = new double[n0];
    double[][,] c = new double[n0][,];
    MKL_INT GroupCount = n0;
    MKL_INT[] GroupSize = new MKL_INT[n0];
    for (int i0 = 0; i0 < n0; i0++)
    {
        int n1 = a[i0].GetLength(0);
        int n2 = a[i0].GetLength(1);
        int n3 = b[i0].GetLength(0);
        int n4 = b[i0].GetLength(1);
        if (n2 != n3) throw new System.Exception("Inner matrix dimensions must agree");
        TransA[i0] = 111; // trans='N'
        TransB[i0] = 111; // trans='N'
        M[i0] = n1; N[i0] = n4; K[i0] = n2;
        alpha[i0] = 1; beta[i0] = 0;
        c[i0] = new double[n1, n4];
        GroupSize[i0] = 1;
    }
    MKL_INT[] lda = K;
    MKL_INT[] ldb = N;
    MKL_INT[] ldc = N;
    _mkl.cblas_dgemm_batch(Order, TransA, TransB, M, N, K, alpha, a, lda, b, ldb, beta, c, ldc, GroupCount, GroupSize);
    return c;
}

and

double[,] A0 = new double[,] { { 1, 2 }, { 3, 4 } };
double[,] A1 = new double[,] { { 5, 6 }, { 7, 8 } };
double[,] B0 = new double[,] { { 9, 10 }, { 11, 12 } };
double[,] B1 = new double[,] { { 13, 14 }, { 15, 16 } };
double[][,] A = new double[][,] { A0, A1 };
double[][,] B = new double[][,] { B0, B1 };
double[][,] C = MKL.Dot(A, B);

When I run the code it works. I can see in the debugger that the method MarshalManagedToNative is called 3 times (as expected) before the cblas_dgemm_batch is called.