Hello,
Does anybody know how to feed correctly std::complex matrices to zgemm function in mkl BLAS? I can compile the code below, but the result is wrong:
====================================================================================
#include <complex>
#include <iostream>
#include <iomanip>
#define MKL_Complex16 std::complex<double>
#include "mkl.h"
typedef std::complex<double> Complex;
const int Arows = 3, Acols = 2, Brows = 3, Bcols = 2, Crows = 2, Ccols = 2;
int main()
{
const char *notrans = "n";
const char *trans = "t";
const char *conj = "c";
Complex A[Arows][Acols], A1[Arows][Acols], B[Brows][Bcols], B1[Brows][Bcols], C[Crows][Ccols], C1[Crows][Ccols];
Complex alpha(1.0, 0.0);
Complex beta(0.0, 0.0);
for (int i = 0; i < Arows; i++){
for(int j = 0; j <Acols; j++){
A[i][j] = Complex(i+1,j+1);
A1[i][j] = Complex(i+1,j+1);
}
}
for (int i = 0; i < Brows; i++){
for(int j = 0; j < Bcols; j++){
B[i][j] = -Complex(j+1,i+1);
B1[i][j] = -Complex(j+1,i+1);
}
}
for (int i = 0; i < Crows; i++){
for(int j = 0; j < Ccols; j++){
C[i][j] = Complex(0.0, 0.0);
C1[i][j] = Complex(0.0, 0.0);
}
}
Complex* pma = &(A[0][0]);
Complex* pmb = &(B[0][0]);
Complex* pmc = &(C[0][0]);
std::cout << "Matrix A" << std::endl;
for (int i = 0; i < Arows; i++){
for(int j = 0; j < Acols; j++){
std::cout << std::setw(4) << A[i][j];
}
std::cout << std::endl;
}
std::cout << "Matrix B" << std::endl;
for (int i = 0; i < Brows; i++){
for(int j = 0; j < Bcols; j++){
std::cout << std::setw(4) << B[i][j];
}
std::cout << std::endl;
}
int iArows = Arows, iAcols = Acols, iBrows = Brows, iBcols = Bcols, iCrows = Crows, iCcols = Ccols;
zgemm( notrans, notrans, &iArows, &iBcols, &iAcols, &alpha, pma, &iArows, pmb, &iBrows, &beta, pmc, &iArows);
for (int i = 0; i < Arows; i++){
for (int j = 0; j < Bcols; j++){
for (int k = 0; k < Acols; k++){
C1[i][j] += A1[i][k] * B1[k][j];
};
};
};
std::cout << "From zgemm" << std::endl;
for (int i = 0; i < Crows; i++){
for(int j = 0; j < Ccols; j++){
std::cout <<std::setw(4) << C[i][j];
}
std::cout << std::endl;
}
std::cout << "From direct" << std::endl;
for (int i = 0; i < Crows; i++){
for(int j = 0; j < Ccols; j++){
std::cout << std::setw(4) << C1[i][j];
}
std::cout << std::endl;
}
return 0;
}
=============================================================================================================
Result:
dmitry@linux-4doj:~/tests> ./zgemm-clean
Matrix A
(1,1)(1,2)
(2,1)(2,2)
(3,1)(3,2)
Matrix B
(-1,-1)(-2,-1)
(-1,-2)(-2,-2)
(-1,-3)(-2,-3)
From zgemm
(-2,-8)(-4,-8)
(-5,-10)(4,-12)
From direct
(5,-22)(2,-26)
(1,-9)(-3,-12)
Any help is very appreciated.
Dmitry