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.
mkl_zdiamm
can be faster for band matrices, then general algorithm ofDot
for sparse matrices. Orzheev
, which can solve a problem withEigenvectors
for V9 and earlier. $\endgroup$