Tension between C++'s complex and MKL (ZGEMM)

[barakb]

Hi there,

I am developing a code in C++ and need to multiply complex matrices together. I want to use BLAS's ZGEMM, but am getting strange results which, I think, are due to precisions issues:

My matrices are defined with complex and I use ZGEMM as an external function: I define it with either

extern "C" {void zgemm(const char*, const char*, const int*, const int*, const int*,const complex *, const complex *, const int*, const complex *, const int*, const complex*, complex *,const int*);}

or, due to possible inconsistencies between complex and MKL_Complex16, with

extern "C" {void zgemm_(const char*, const char*, const int*, const int*, const int*,const void*, void*
, const int*, void*, const int*, const void*, void*,const int*);}

and while the code compiles, I find that when I use ZGEMM iteratively (i.e. put as input matrices that were obtained as output before hand), then after many iterations (depending on the size of the matrices), results are wrong.

This leads me to suspect that the way I define zgemm above leads to incompatibilities in the way complex variables are defined in C++ and MKL.

To fix this I tried to follow the MKL's user manual suggestion and to add

#define MKL_Complex16 complex
#include "mkl.h"

to my code, but this does not compile! I get the compilation error

/local/opt/intel/mkl/10.0.2.018/include/mkl_types.h(43): error: expected a ";"
} MKL_Complex16;
^

So, does anyone know how to properly call from C++ MKL functions with complex variables ?

I don't think I do ;-)

Thanks,

Barak

[Gennady_F_Intel]

Barak,
please look at the post disccused earlier.
you can find zgemm_bug_small.cpp example there... will it help you?
--Gennady

[barakb]

Quoting - Gennady Fedorov (Intel)

Barak,
please look at the post disccused earlier.
you can find zgemm_bug_small.cpp example there... will it help you?
--Gennady

Hi Gennady,

I'm not sure it does. I am comparing the results of calling ZGEMM from C++ to those obtained with MTL, and am getting results which are _not_ exact at machine precision. Then, when I iterate the multiplications these small errors become significant.

I'm not sure how to debug this...

Do you know of any alternatives to ZGEMM? I actually only need to multiply ywo complex matrices together...

[Dmitry_B_Intel]

Hi barakb,

That ZGEMM from MKL and similar functionfrom MTL produce results that are not identical to the last bit is expected. Moreover, ZGEMM may produce non-identical results from run to run on the same machine if it is run in threaded mode (see for example a discussion onMKL numerical stability and threading). That the iterative algorithm diverges is more the property of that algorithm rather than deficiency of matrix-matrix multiplication. You may observedifferent resultsiterating even simpler logistics equation: x = lambda*x*(1 - x).

Regarding compilation error, the version of MKL you use (10.0.2) doesn't contain sentinels that check if MKL_Complex16 is already defined. They are introduced in a later MKL release. However, you may solve this compile problem with this trick:

#define MKL_Complex16 forget_this_name
#include "mkl_types.h"
#undef MKL_Complex16

#define MKL_Complex16 complex
#include "mkl.h"

Thanks
Dima