MKL cgemm

[hotimhou]

Hi,
I would like to find out if anyone has any problems using the function cgemm in MKL. I can't seem to get it right, even though I've tried different combinations of the command. Attached here are the C code for a 4x4 matrix premultiplying a 4x6 matrix.

The results are here:

A = [
0.0+1.0*j 1.0+2.0*j 2.0+3.0*j 3.0+4.0*j ;
1.0+2.0*j 2.0+3.0*j 3.0+4.0*j 4.0+5.0*j ;
2.0+3.0*j 3.0+4.0*j 4.0+5.0*j 5.0+6.0*j ;
3.0+4.0*j 4.0+5.0*j 5.0+6.0*j 6.0+7.0*j ;
]

B = [
0.0+1.0*j 1.0+2.0*j 2.0+3.0*j 3.0+4.0*j 4.0+5.0*j 5.0+6.0*j ;
1.0+2.0*j 2.0+3.0*j 3.0+4.0*j 4.0+5.0*j 5.0+6.0*j 6.0+7.0*j ;
2.0+3.0*j 3.0+4.0*j 4.0+5.0*j 5.0+6.0*j 6.0+7.0*j 7.0+8.0*j ;
3.0+4.0*j 4.0+5.0*j 5.0+6.0*j 6.0+7.0*j 7.0+8.0*j 8.0+9.0*j ;
]

C = [
-16.0+40.0*j -20.0+56.0*j -24.0+72.0*j -28.0+88.0*j -22.0+44.0*j -26.0+72.0*j ;
-30.0+100.0*j -34.0+128.0*j -28.0+88.0*j -32.0+128.0*j -36.0+168.0*j -40.0+208.0*j ;
-24.0+72.0*j -28.0+104.0*j -32.0+136.0*j -36.0+168.0*j -30.0+76.0*j -34.0+120.0*j ;
-38.0+164.0*j -42.0+208.0*j -36.0+120.0*j -40.0+176.0*j -44.0+232.0*j -48.0+288.0*j ;
]

The results are different from those I got from Octave or Matlab.
Can anyone advise me what I did wrong? Thanks in advance!

TH

[Todd_R_Intel]

It looks like you are setting up a row major data structure and passing it to the CGEMM function which expects a column-major data structure. Either you'll need to initialize A and B so thatthey store by column rather than row and have pA, pB, and pC point to columns instead of rows -or- you could try out the cblas_cgemm routine described in an appendix of the reference manual.

-Todd

[hotimhou]

Hi Todd,
Thanks for your help. I tried cblas_cgemm as well. I could compile the code, but it produces a segmentation fault when cblas_cgemm was called. Basically, I just replaced the cgemm call in the previous code with cblas_cgemm, with an additional parameter specifying the matrices being row-major. The program was linked with:
lmkl_ia32 -lguide -lpthread -lm
the way it was documented. Any problem with such a scheme?

[hotimhou]

As a follow up to my previous post, I have corrected the code that posted earlier. For straight forward complex matrix multiplication of the form C := Alpha*A*B + Beta*C, let RowA, ColA be the number of rows and columns of A, RowB, ColB be the number of rows and columns of B respectively. C is assumed to be compliant to the dimensions of multiplication of A and B. Alpha and Beta are both complex values. The matrix is stored in a row-major type array.

Then using cblas_cgemm, the resultant C can be found in this manner:

cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, RowA, ColB, ColA, Α, A, ColA, B, ColB, Β, C, ColB);
(Refer to mkl_cblas.h for the prototype of cblas_cgemm.)

Here, lda, ldb, and ldc are the number of entries for each row. I hope this info is useful for anyone out there trying to perform complex matrix multiplication.

Message Edited by hotimhou@yahoo.com on 07-30-2005 02:38 AM