Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

I am currently trying to link a program against the Intel MKL 11.0 library instead of using NetLIB or OpenBLAS. Doing this I recognized the following error which I can not explain to my self at the moment. Consider the following C code example computing a complex scalar product using zdotc:

#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
double complex zdotc_(int *n, double complex *X, int *incx, double complex *Y, int *INCY ); 
int main ( ) {
    int n = 5; 
    int incx = 1, incy = 1; 
    double complex x[5] = {1,I,2,2+I,3}; 
    double complex y[5] = {I,3,I*3, 2+2*I, 9}; 
    double complex ret; 
    ret = zdotc_(&n,x,&incx,y,&incy); 
    printf("n   = %d\n", n); 
    printf("ret = %lg + %lgi\n", creal(ret), cimag(ret)); 
    return 0; 
}

I compiled this example using the command line flags given by the MKL Advisor. I select "GNU C/C++, 32 Bit Integer, Dynamic Linking, GNU OpenMP". The resulting command line is:

gcc zdotc_test.c -o zdot_mkl_gcc -O2  -L$MKLROOT/lib/intel64 -lmkl_intel_lp64 -lmkl_gnu_thread -lmkl_core -ldl -lpthread -lm  -fopenmp -m64 -I$MKLROOT/include

The output of this program is:

n   = 0
ret = 0 + 1.07933e+21i

which is obviously wrong and especially why is n altered?

If I select GNU Fortran instead of GNU C/C++, I have to replace -lmkl_intel_lp64 by -lmkl_gf_lp64 and then the correct output

n   = 5
ret = 33 + 6i

is produced.

So my question is: where are the detailed differences between those to interfaces and why does the first one produced this error?

share|improve this question
    
You probably want to use the C interface instead of the Fortran one (it comes with MKL) as it has far fewer portability issues. For the Fortran interface you would have to deal with the Fortran calling convention, which varies with platform and compiler. –  Rufflewind Jul 2 at 0:38

1 Answer 1

up vote 3 down vote accepted

The difference has to do with the calling convention (ABI) differences. ZDOTC is a problematic function because it returns a double complex, which is usually considered a struct rather than a "simple" data type. Thus, it can either be returned on the stack as a return value, or by reference as an implicit first argument to the "C" style function. I am guessing that what you are seeing here is that in the _intel_ case, it expects the return value to be an implicit first argument, so it returns the value and ends up modifying n, and the return value is then just completely wrong. In that case the correct declaration is

void zdotc_(double complex *retval, int *n, double complex *X, int *incx, double complex *Y, int *INCY );

It's stupid stuff like this that makes me implement everything in C++ if possible. In particular, I tend to implement BLAS level 1 in C++, where performance is limited by memory bandwidth anyways. For BLAS levels 2 and 3, I will call actual BLAS, which does not have functions with this problem.

I believe I've seen benchmarks before comparing naive dot product implementations with optimized BLAS, and if you use __restrict type keywords, most optimizing compilers will produce code that's as fast as the optimized BLAS.

My own version is located here, where you would need the corresponding source file to provide forwarding of higher level BLAS, located here. I would refer you to Eigen, but my code is slightly easier to port to C99 since it's less heavily templated.

share|improve this answer
    
Ok this small ABI change seems to explain the effect I have seen in my software, that all functions relying on zdotc and zdotu are crashing. So I have to write a configure check which variant of zdotc is working in the given BLAS implementation. But in this way the MKL is small kind of incompatible whit the Quasi-Standard BLAS ABI. Because the calling convention of this functions changed. –  Grisu Feb 26 '13 at 10:45

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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