16
$\begingroup$

MKL has a lot of functions. The question is: how to use these functions directly from Mathematica?

Furthermore, Mathematica shipped with some MKL libraries, which is located in

$InstallationDirectory/SystemFiles/Libraries/

Can I use these internal MKL libraries? Is there any restriction?

$\endgroup$
  • $\begingroup$ What functions do you need? Many of Mathematica's built in functions rely on the MKL. This includes I believe (not sure) most of machine-precision linear algebra, fft and random number generation (possible it's not the default RNG method but it can be selected). $\endgroup$ – Szabolcs Sep 7 '14 at 17:09
  • $\begingroup$ @Szabolcs For example, I thought that mkl_zdiamm can be faster for band matrices, then general algorithm of Dot for sparse matrices. Or zheev, which can solve a problem with Eigenvectors for V9 and earlier. $\endgroup$ – ybeltukov Sep 7 '14 at 17:15
  • $\begingroup$ Yes, it's definitely interesting to be able to do this. $\endgroup$ – Szabolcs Sep 7 '14 at 17:16
13
$\begingroup$

I can partially answer to my own questions. Amazingly, but it is easier to use the internal MKL. Let us consider my question about multiplication a band matrix by a dense matrix. The corresponding function is mkl_zdiamm.

I wrote the following code (diamm.c)

#include <stdio.h>
#include <stdlib.h>
#include <WolframLibrary.h>
#include <mkl.h>

DLLEXPORT mint WolframLibrary_getVersion( ) {
    return WolframLibraryVersion;
}

DLLEXPORT int WolframLibrary_initialize( WolframLibraryData libData) {
    return 0;
}

DLLEXPORT void WolframLibrary_uninitialize( WolframLibraryData libData) {
    return;
}

DLLEXPORT int version(WolframLibraryData libData, mint argc, MArgument *args, MArgument res) {
  char* buf = (char*)malloc(200*sizeof(char));
  mkl_get_version_string(buf, 200);
  MArgument_setUTF8String(res, buf);
  return LIBRARY_NO_ERROR; 
}


DLLEXPORT int zdiamm(WolframLibraryData libData, mint argc, MArgument *args, MArgument res) 
{

    freopen("/tmp/diamm.log", "w", stdout);
    setbuf(stdout, NULL); // print immideately

    MTensor A  = MArgument_getMTensor(args[0]);
    MTensor ID = MArgument_getMTensor(args[1]);
    MTensor B  = MArgument_getMTensor(args[2]);
    MTensor C;

    const mint *dimsA =  libData->MTensor_getDimensions(A);
    const mint *dimsID = libData->MTensor_getDimensions(ID);
    const mint *dimsB =  libData->MTensor_getDimensions(B);

    MKL_INT m = dimsB[1], n = dimsB[0], nd = dimsID[0];

    int err = LIBRARY_NO_ERROR;

    err = libData->MTensor_new(MType_Complex, 2, dimsB, &C);
    if (err) return err;    

    MKL_Complex16 *dataA = (MKL_Complex16*) libData->MTensor_getComplexData(A);
    MKL_INT *dataID      = (MKL_INT*)       libData->MTensor_getIntegerData(ID);
    MKL_Complex16 *dataB = (MKL_Complex16*) libData->MTensor_getComplexData(B);
    MKL_Complex16 *dataC = (MKL_Complex16*) libData->MTensor_getComplexData(C);

    MKL_Complex16 alpha = {1.0, 0.0}, beta = {0.0, 0.0}, zero = {0.0, 0.0};

    int i;

    for (i = 0; i < n*m; i++) {
        dataC[i] = zero;
    }

    char transa = 'N';
    char *matdescra = "G";

    mkl_zdiamm(&transa, &m, &n, &m, &alpha, matdescra, dataA, &m, dataID, &nd, dataB, &m, &beta, dataC, &m);

    MArgument_setMTensor(res, C);
    return LIBRARY_NO_ERROR;
}

and compile it with (I use Linux)

gcc -DMKL_ILP64 -shared -fPIC -o libdiamm.so -I/usr/local/Wolfram/Mathematica/10.0/SystemFiles/IncludeFiles/C/ -I /opt/intel/composerxe-2013.3.174/mkl/include/ diamm.c

Useful notes:

  • There is no explicit reference to MKL libraries. Internal MKL libraries are already linked by Mathematica. Unfortunately, it doesn't allow me to use external MKL libraries.
  • I use headers from external MKL, because Mathematica shipped without MKL headers.
  • It seems to me that it works without MKL license.
  • Be careful with -DMKL_ILP64. It means "64-bit integers" and it is necessary for 64-bit systems (at least for 64-bit Linux). See also: MKL Link Line Advisor.

Now let's go back to Mathematica

AppendTo[$LibraryPath, NotebookDirectory[]];
version = LibraryFunctionLoad["libdiamm", "version", {}, "UTF8String"];
zdiamm = LibraryFunctionLoad["libdiamm", "zdiamm", 
     {{_Complex, 2}, {_Integer, 1}, {_Complex, 2}}, {_Complex, 2}];

version[]

Intel(R) Math Kernel Library Version 11.1.2 Product Build 20140122 \ for Intel(R) 64 architecture applications

  • I assume that libdiamm.so is in the same directory.
  • It is really internal MKL since I have installed a bit different version (11.1.3).
  • Do not forget to LibraryUnload["libdiamm"] after recompiling the code.

It remains to check the zdiamm function (I omit the details, because it is not relevant to the question under consideration)

$HistoryLength = 0;
n = 10000;
b = 300;
k = 300;
a = SparseArray[
   Flatten[#, 1] &@Table[{i, Mod[i + j, n, 1]}, {i, n}, {j, -b, b}] ->
     RandomComplex[{-1 - I, 1 + I}, n (2 b + 1)]];
u = RandomComplex[{-1 - I, 1 + I}, {n, k}];
v = a.u;

id = Union[#2 - #] & @@ Transpose[a@"NonzeroPositions"];
val = Normal@SparseArray[
     Transpose@{#, Length@a + #2 - #} & @@ 
       Transpose[a@"NonzeroPositions"] -> a@"NonzeroValues"][[All, id + Length@a]];
v3 = Transpose@zdiamm[Transpose@val, id, Transpose@u];

Max@Abs[v3 - v]

1.31642*10^-13

It works and return the same result! Unfortunately, mkl_zdiamm is several times slower then Mathematica's Dot for sparse matrices.

$\endgroup$
  • $\begingroup$ Where do you get the MKL headers? I thought MKL was not free. (Haven't tried the code yet, but I plan to try it on OS X.) $\endgroup$ – Szabolcs Sep 7 '14 at 17:11
  • $\begingroup$ @Szabolcs I have free noncommercial license. It seems to me that it is free only for Linux. Ref (English version is under reconstruction). $\endgroup$ – ybeltukov Sep 7 '14 at 17:31
  • $\begingroup$ It appears that you are using functions (such as mkl_get_version_string), which are not accessible at link time. You don't seem to be linking to the MKL. How is it then that the linker doesn't complain? Is it due to the prototype of these functions in mkl.h? If yes, what does the prototype need to look like? I was trying to replicate this on OS X with no success. I do not have the MKL installed, so I was trying to figure out the correct prototype based on headers I found on the internet. $\endgroup$ – Szabolcs Sep 5 '16 at 9:05
  • $\begingroup$ @Szabolcs I did not provide any specific linker options. Is it possible to download MKL (the licence is not necessary) and then figure out the necessary headers? $\endgroup$ – ybeltukov Sep 6 '16 at 14:05
  • $\begingroup$ I tried it on Linux now. It appears that it is possible to create a shared library even if there are undefined symbols if compiling with -fPIC. I assume these symbols will be resolved when the shared library is loaded. On OS X this does not seem to work after a few naive tries: it always gives an error about undefined symbols. $\endgroup$ – Szabolcs Sep 6 '16 at 14:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.