.NET Interop and the Math Kernel Library LAPACK Routines

bwhazel · ‎01-27-2011

I am building a program in C# and would like to use the DSYEV LAPACK routine within it. I have built a wrapper DLL in Fortran (code is below) linked against the MKL sequentially. This DLL is then linked into the C# via the System.Runtime.InteropServices.DllImport attribute.

Additionally, due to the differences between multi-dimensional arrays in C# and Fortran, I have added two additional methods to convert a 2-dimensional array into a 1-dimensional one and vice versa.

The Fortran DLL compiles without a problem, but the DSYEV function within C# is not functioning properly. The INFO parameter returns as -1, so a link seems to have been made, but the arrays A and W are returning with only 1 element when A should have N*N elements (to be reconverted back into 2D array) and W should have N elements.

More information on the DSYEV subroutine can be found here.

Basically, my question is: what am I doing wrong?

INTEROP DECLARATION:

[csharp]public static extern void DSYEV_INTEROP(string JOBZ, string UPLO, int N, ref double[] A, int LDA, ref double[] W, ref double[] WORK, int LWORK, ref int INFO);[/csharp]

FORTRAN DLL:

[fortran]SUBROUTINE DSYEV_INTEROP(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO)

    !DEC$ ATTRIBUTES DLLEXPORT::DSYEV_INTEROP
    !DEC$ ATTRIBUTES STDCALL::DSYEV_INTEROP
    !DEC$ ATTRIBUTES ALIAS:"DSYEV_INTEROP"::DSYEV_INTEROP

    !DEC$ ATTRIBUTES VALUE::JOBZ
    !DEC$ ATTRIBUTES VALUE::UPLO
    !DEC$ ATTRIBUTES VALUE::N
    !DEC$ ATTRIBUTES REFERENCE::A
    !DEC$ ATTRIBUTES VALUE::LDA
    !DEC$ ATTRIBUTES REFERENCE::W
    !DEC$ ATTRIBUTES REFERENCE::WORK
    !DEC$ ATTRIBUTES VALUE::LWORK
    !DEC$ ATTRIBUTES REFERENCE::INFO

    CHARACTER(LEN=1), INTENT(IN) :: JOBZ, UPLO
    INTEGER, INTENT(IN) :: N, LDA, LWORK
    INTEGER, INTENT(OUT) :: INFO
    !REAL*8, DIMENSION(:,:), INTENT(INOUT) :: A
    REAL*8, POINTER, DIMENSION(:,:) :: A
    REAL*8, DIMENSION(:), INTENT(INOUT) :: W, WORK

    CALL dsyev(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO)

END SUBROUTINE DSYEV_INTEROP[/fortran]

Gennady_F_Intel · ‎01-27-2011

Please check the first parameter first. the info == -1 points thatthe 1-th parameter had an illegal value.

Vladimir_Koldakov__I · ‎01-28-2011

Hello,
You do not need to write wrappers. Please see examples of using MKL in C#:
http://software.intel.com/en-us/articles/using-intel-mkl-in-your-c-program/
You need not also to build custom dll if you use MKL 10.3.
Thanks,
Vladimir

bwhazel · ‎02-16-2011

Many thanks for the responses, and apologies for my delayed reply! I am looking into this option for .NET. Is there a simple way of knowing which data types should be passed by reference or value?

Gennady_F_Intel · ‎02-16-2011

by reference.

as an example, the Pardiso calling will be looks like

phase = 11;

Pardiso.pardiso (pt, ref maxfct, ref mnum, ref mtype, ref phase,

ref n, a, ia, ja, idum, ref nrhs,

iparm, ref msglvl, ddum, ddum, ref error);

bwhazel · ‎02-17-2011

I've tried the example in the download: the C# code is written. I am, however, experiencing some difficulty building the custom mkl.dll for the DSYEV routine.

I am not using the makefile, but am using it as an example to follow. I first compiled the _fseeki64.c file as in the makefile using the Microsoft C compiler (I don't have the Intel one). I then tried a variation on the link command, having set additional environment variables for the MKL libraries and compiler libraries:

MKLROOT=C:\Program Files (x86)\Intel\ComposerXE-2011\mkl\lib\intel64
OMPLIB=C:\Program Files (x86)\Intel\ComposerXE-2011\compiler\lib\intel64

where each file was preceeded with the necessary environment variable:

link /DLL /MACHINE:X64 /def:user.def _fseeki64.obj %MKLROOT%\mkl_intel_lp64_dll.lib %MKLROOT%\mkl_intel_thread_dll.lib %MKLROOT%\mkl_core_dll.lib %OMPLIB%\libiomp5md.lib msvcrt.lib user32.lib /out:mkl.dll

All files are found, but then I get the error:

LINK : error LNK2001: unresolved external symbol _DllMainCRTStartup
mkl.dll : fatal error LNK1120: 1 unresolved externals

Basically, what am I doing wrong?

Gennady_F_Intel · ‎02-17-2011

may be You didn't set the MKLROOT?

I checked how it works.

nmake ia32 MKLROOT=C:\APPS\Intel\MKL\10.2.7

Microsoft Program Maintenance Utility Version 8.00.50727.42

Add path of the MKL libraries to the lib environment variable

set lib=%MKLROOT%\ia32\lib;%lib%

MKL entries for custom dll

Workaround for pardiso

Microsoft 32-bit C/C++ Optimizing Compiler Version 14.00.50727.42 for 80x86

_fseeki64.c

Build MKL custom dll

nmake mkl.dll MACHINE=IX86 MKL_LIB="mkl_intel_c_dll.lib mkl_intel_threa

d_dll.lib mkl_core_dll.lib" MSLIB=user32.lib

Microsoft Program Maintenance Utility Version 8.00.50727.42

link /DLL /MACHINE:IX86 /def:user.def _fseeki64.obj mkl_intel_c_dll.lib

mkl_intel_thread_dll.lib mkl_core_dll.lib libiomp5md.lib user32.lib /out:mkl.dl

l

Microsoft Incremental Linker Version 8.00.50727.42

Creating library mkl.lib and object mkl.exp

Add path of the mkl.dll and path of the MKL binaries to the path environment var

iable

set path=c:\!!;%MKLROOT%\ia32\bin;%path%

Build and run examples

nmake /a dgemm.exe dgeev.exe dfti_d1.exe pardiso.exe vddiv.exe

Microsoft Program Maintenance Utility Version 8.00.50727.42

Compile dgemm.cs

csc .\dgemm.cs

Microsoft Visual C# 2005 Compiler version 8.00.50727.4927

for Microsoft Windows 2005 Framework version 2.0.50727

Run dgemm example

dgemm.exe

....................

bwhazel · ‎03-13-2011

Right - I've got it working! I didn't realise I could just link the provided DLLs rather than build a custom one from the static libraries!

Many thanks for all the assistance!

trombif · ‎04-12-2011

Hi,
I am using eigen values/vectors decomposition DSYEV in mkl for a university project and I need results to be reproductible.
I wrote a little wrapper (according to examples) to call DSYEV from C# :


namespace Eigen
{
 class Solve
 {
 public double[,] decompo(double[,] A)
 {
 int lwork = -1;
 int n = A.GetLength(0);
 int lda = n;
 double[] A_bis = new double[n * n];
 double[] w = new double;
 for (int i = 0; i < n; i++)
 {
 int cte = i * n;
 for (int j = 0; j < n; j++)
 {
 A_bis[cte + j] = A[j, i];
 }
 }
 lwork = -1;
 double[] work_1 = new double[1];
 int info = LAPACK_mkl.dsyev('V', 'U', n, A_bis, lda, w, work_1, lwork);
 lwork = (int)work_1[0];
 double[] work = new double[lwork];
 info = LAPACK_mkl.dsyev('V', 'U', n, A_bis, lda, w, work, lwork);
 double[,] res2 = new double[n, n + 1];
 for (int i = 0; i < n; i++)
 {
 res2[i, 0] = w;
 for (int j = 0; j < n; j++)
 {
 res2[i, j + 1] = A_bis[i + j * n];
 }
 }
 GC.Collect(0);
 return res2;
 }
 }
}
namespace mkl
{
 public sealed class LAPACK_mkl
 {
 private LAPACK_mkl() { }
 public static int dsyev(char jobz, char uplo, int N, double[] A, int LDA, double[] w, double[] work, int lwork)
 {
 LAPACKNative.kmp_set_warnings_off();
 int num = Environment.ProcessorCount;
 LAPACKNative.omp_set_num_threads(ref num);
 int info = -1;
 LAPACKNative.dsyev(ref jobz, ref uplo, ref N, A, ref LDA, w, work, ref lwork, ref info);
 return info;
 }
 }
 /** LAPACK native declarations */
 [SuppressUnmanagedCodeSecurity]
 internal sealed class LAPACKNative
 {
 [DllImport("libiomp5md", EntryPoint = "omp_set_num_threads")]
 internal static extern void omp_set_num_threads(ref int num);
 [DllImport("libiomp5md", EntryPoint = "kmp_set_warnings_off")]
 internal static extern void kmp_set_warnings_off();
 [DllImport("mkl", EntryPoint = "DSYEV", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
 internal static extern void dsyev(ref char jobz, ref char uplo, ref int n, double[] A, ref int lda, double[] w, double[] work, ref int lwork, ref int info);
 private LAPACKNative()
 {
 kmp_set_warnings_off();
 int num = Environment.ProcessorCount; 
 omp_set_num_threads(ref num); 
 }
 }
}



I red something about memory alignment and 16-byte boundaries... does this apply here ?

Thank you for any answer or comment on my code

TimP · ‎04-12-2011

Yes, _aligned_malloc(size,16) rather than the default new for the arrays passed to MKL should improve numerical consistency and performance, at least on 32-bit Windows. On 64-bit Windows, perhaps new would default to 16-byte alignment.
I've seen an assertion that Intel C++ new would observe alignments specified by _declspec(align(16)). The Microsoft page states clearly that you must use _declspec, with _declspec(align(32)) preferred for improved caching and compatibility with AVX.
You omitted needed header definitions, which I wouldn't consider obvious.

trombif · ‎04-12-2011

First, thank you for the quick reply !

As I mentionned in my previous post, this is C# code... and_aligned_malloc you are talking about is C++ (am I wrong?). Is there such a command, or a trick to work around and create arrays which are16-byte aligned, for C# ?

And here is a list of definitions :

using System;

using System.Collections.Generic;

using System.Linq;

using System.Text;

using System.Security;

using System.Runtime.InteropServices;

using mkl;

Thank you

TimP · ‎04-12-2011

_aligned_malloc is probably simply a wrapper around malloc. The usual technique is to allocate a region, say 12 bytes longer than the data size (as you should be able to depend on 4-byte aligned malloc). Then you use the sub-region beginning at the first 16-byte aligned address in that region. If you free the region, you must free the entire region, including the bytes you skip. You shouldn't need a loop to find the aligned address, such as some of the web search references on c# data alignment show. You simply add 15 to the base pointer and mask off the low order bits, or find the remainder %16 of the base pointer and figure out how far up the region to start. For AVX, you probably want the extra 28 bytes, add 31 and mask to 32-byte alignment. I'm not fluent in C#, so should shut up now.

TimP · ‎04-12-2011

I discussed with a colleague who uses C#, and we suggest you use the mkl_malloc() function rather than c# new.

trombif · ‎04-13-2011

I am not fluent in C# either: the main resaon I have chosen C# instead of C++ is memory management, and it occurs to me that it becomes indeed more complex ...

Do you have a code example on how to use mkl_malloc() in C# ?

Thank you very much

TimP · ‎04-13-2011

No, I could get from my colleague an example of ipp malloc in C#. It should work like a plain malloc() with the addition of the alignment argument (by value).

trombif · ‎04-13-2011

That would be great !

TimP · ‎04-13-2011

[bash] [DllImport("ippcore-7.0.dll",EntryPoint="ippMalloc")] private static extern
 IntPtr Malloc_x32 
 (int length );
 [DllImport("ippcoreem64t-7.0.dll",EntryPoint="ippMalloc")] private static extern
 IntPtr Malloc_x64 
 (int length );
 public IntPtr Malloc (int length ) {
  if (IntPtr.Size == 8)  
   return Malloc_x64(length);
  else
   return Malloc_x32(length);
 }
[/bash]

He checks whether to use the 32-bit or X64 ipp .dll (according to IntPtr.Size), then uses the ippMalloc, which is said to return a 64-byte aligned region. You could use that one, if you don't mind involving another big .dll, or use mkl_malloc(length, alignment)

For that matter, I would think you could still use the _aligned_malloc provided by Visual Studio.

We just noticed that VS2010 doc says the standard default malloc() provides 16-byte alignment. No idea if that would carry over to c# new. Prior to VS2010, malloc was 16-byte aligned only in the X64 version.

trombif · ‎04-14-2011

Okay, thank you.

Using dependency walker on mkl_rt.dll, I foundMKL_MALLOC function, so it is ok.

In mkl documentation, it seems that syntax is :

a_ptr=mkl_malloc(alloc_size,alignment);

So i defined this :

[csharp]        public static unsafe IntPtr mkl_malloc(int length)
        {
            LAPACKNative.kmp_set_warnings_off();
            int align = 32;
            IntPtr pointer = LAPACKNative.Malloc(length, align);
            return pointer;
        }
        [DllImport("mkl_rt.dll", EntryPoint = "MKL_MALLOC", CallingConvention = CallingConvention.Cdecl), SuppressUnmanagedCodeSecurity]
        internal static extern IntPtr Malloc(int length, int align);[/csharp]

Problem is, when I call this function this way :

IntPtr pointer_A_bis = LAPACK_mkl.mkl_malloc(n * n * sizeof(double));

I get an error concerning "an attempt to to read or write in protected memory. It often indicates another memory is damaged" ...

Any idea ?

trombif · ‎05-10-2011

Finally, here is what I did (as I could not make mkl_malloc work...)

[csharp]public unsafe double[] linsolve(ref double[,] A, ref double[] B)
        {
            int m = A.GetLength(0);
            int n = A.GetLength(1);
            int lwork = -1;
            int lda = m;
            int ldb = m;
            int nrhs = 1;
            int rank = 0;
            int align = 128;
            IntPtr pointer_A_bis_ini = Marshal.AllocHGlobal(m * n * sizeof(double) + align);
            IntPtr pointer_A_bis = pointer_A_bis_ini;
            IntPtr pointer_B_bis_ini = Marshal.AllocHGlobal(n * nrhs * sizeof(double) + align);
            IntPtr pointer_B_bis = pointer_B_bis_ini;
            Int64 adresse_A_bis = pointer_A_bis.ToInt64() % align;
            Int64 adresse_B_bis = pointer_B_bis.ToInt64() % align;
            int offset = 0;
            while (adresse_A_bis != 0)
            {
                offset++;
                pointer_A_bis = new IntPtr(pointer_A_bis_ini.ToInt64() + offset);
                adresse_A_bis = pointer_A_bis.ToInt64() % align;
            }
            offset = 0;
            while (adresse_B_bis != 0)
            {
                offset++;
                pointer_B_bis = new IntPtr(pointer_B_bis_ini.ToInt64() + offset);
                adresse_B_bis = pointer_B_bis.ToInt64() % align;
            }
            double* A_bis = (double*)pointer_A_bis.ToPointer();
            double* B_bis = (double*)pointer_B_bis.ToPointer();
            for (int i = 0; i < nrhs; i++)
            {
                int cte = i * m;
                for (int j = 0; j < m; j++)
                {
                    B_bis[cte + j] = B;
                }
            }
            for (int i = 0; i < n; i++)
            {
                int cte = i * m;
                for (int j = 0; j < m; j++)
                {
                    A_bis[cte + j] = A[j, i];
                }
            }
            if (m == n)
            {
                IntPtr pointer_IPIV_ini = Marshal.AllocHGlobal(m * sizeof(int) + align);
                IntPtr pointer_IPIV = pointer_IPIV_ini;
                Int64 adresse_IPIV = pointer_IPIV.ToInt64() % align;
                offset = 0;
                while (adresse_IPIV != 0)
                {
                    offset++;
                    pointer_IPIV = new IntPtr(pointer_IPIV_ini.ToInt64() + offset);
                    adresse_IPIV = pointer_IPIV.ToInt64() % align;
                }
                int* IPIV = (int*)pointer_IPIV.ToPointer();
                for (int i = 0; i < m; i++)
                {
                    IPIV = 0;
                }
                LAPACK_mkl.dgetrf(m, m, A_bis, m, IPIV);
                LAPACK_mkl.dgetrs('N', m, nrhs, A_bis, m, IPIV, B_bis, m);
                double[] res2 = new double;
                for (int i = 0; i < n; i++)
                {
                    res2 = B_bis;
                }
                Marshal.FreeHGlobal(pointer_A_bis_ini);
                Marshal.FreeHGlobal(pointer_B_bis_ini);
                Marshal.FreeHGlobal(pointer_IPIV_ini);
                return res2;
            }
            else
            {
                int maxMN;
                if (m > n) { maxMN = m; } else { maxMN = n; }
                IntPtr pointer_JPVT_ini = Marshal.AllocHGlobal(n * sizeof(int) + align);
                IntPtr pointer_JPVT = pointer_JPVT_ini;
                Int64 adresse_JPVT = pointer_JPVT.ToInt64() % align;
                offset = 0;
                while (adresse_JPVT != 0)
                {
                    offset++;
                    pointer_JPVT = new IntPtr(pointer_JPVT_ini.ToInt64() + offset);
                    adresse_JPVT = pointer_JPVT.ToInt64() % align;
                }
                int* JPVT = (int*)pointer_JPVT.ToPointer();
                for (int i = 0; i < n; i++)
                {
                    JPVT = 0;
                }
                double RCOND = 2.0E-16;
                IntPtr pointer_work_1 = Marshal.AllocHGlobal(sizeof(double));
                double* work_1 = (double*)pointer_work_1.ToPointer();
                LAPACK_mkl.dgelsy(m, n, nrhs, A_bis, lda, B_bis, maxMN, JPVT, RCOND, rank, work_1, lwork);
                lwork = (int)work_1[0];
                IntPtr pointer_work_ini = Marshal.AllocHGlobal(lwork * sizeof(double) + align);
                IntPtr pointer_work = pointer_work_ini;
                Int64 adresse_work = pointer_work.ToInt64() % align;
                offset = 0;
                while (adresse_work != 0)
                {
                    offset++;
                    pointer_work = new IntPtr(pointer_work_ini.ToInt64() + offset);
                    adresse_work = pointer_work.ToInt64() % align;
                }
                double* work = (double*)pointer_work.ToPointer();
                LAPACK_mkl.dgelsy(m, n, nrhs, A_bis, m, B_bis, maxMN, JPVT, RCOND, rank, work, lwork);
                
                double[] res2 = new double;
                for (int i = 0; i < n; i++)
                {
                    res2 = B_bis;
                }
                Marshal.FreeHGlobal(pointer_work_1);
                Marshal.FreeHGlobal(pointer_work_ini);
                Marshal.FreeHGlobal(pointer_A_bis_ini);
                Marshal.FreeHGlobal(pointer_B_bis_ini);
                Marshal.FreeHGlobal(pointer_JPVT_ini);
                GC.Collect(0);
                return res2;
            }
        }[/csharp]

It seems it's working fine.

If anyone has a better idea, feel free to share !

Best regards, Eric