Matrix Multiplication problem: Wrong results in a 1st row of a resulting matrix C ( C[8][8] = A[8][8] x B[8][8] )

[SergeyKostrov]

I detected some problem with SGEMM and DGEMM functions. In essence, there are wrong numbers in a 1st row of a resulting matrix C ( C[8][8] = A[8][8] x B[8][8] ).

It is verified with 32-bit and 64-bit versions of MKL:
...
#define __INTEL_MKL__ 10
#define __INTEL_MKL_MINOR__ 3
#define __INTEL_MKL_UPDATE__ 12
...
and
...
#define __INTEL_MKL__ 11
#define __INTEL_MKL_MINOR__ 0
#define __INTEL_MKL_UPDATE__ 2
...
using Intel, Microsoft, MinGW and Borland C++ compilers.

[SergeyKostrov]

Here are matricies: A[8][8]= 5 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 x B[8][8]= 7 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 Here is a result calculated usimg an online matrix-calculator: [ www.bluebit.gr/matrix-calculator/matrix_multiplication.aspx ] C[8][8]= 70.0 80.0 120.0 160.0 200.0 240.0 280.0 320.0 42.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 42.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 42.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 42.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 42.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 42.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 42.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 Here is a result calculated by: [ MKL SGEMM and DGEMM functions ] C[8][8]= 70.0 84.0 126.0 168.0 210.0 252.0 294.0 336.0 40.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 40.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 40.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 40.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 40.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 40.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 40.0 72.0 108.0 144.0 180.0 216.0 252.0 288.0 The C/C++ test case is based on a well known MKL example for these functions ( nothing is special / alpha=1 and beta=0 ). Please do your own verifications. Thanks in advance.

[SergeyKostrov]

Also, if I use 'T' and 'T' for first two input parameters in both functions then a resulting transposed matrix C is correct. In overall, wrong numbers look like some rounding error ( I could be wrong here ) and it needs to be investigated.

[mecej4]

You have, apparently, overlooked the fact that C stores matrices row-by-row, whereas Fortran stores them column-by-column. You did not show any code, but your results are consistent with passing C-native 2-D arrays to routines (such as SGEMM/DGEMM) that expect Fortran-native 2-D arrays.

There are routines in MKL with the prefix cblas_  that you may consider for using in combination with C-native 2-D arrays.

[SergeyKostrov]

Mecej4, Just run the test with values for matricies I gave. If I don't put values 5 and 7 repectively ( 1st values ) it calculates the C properly.

[SergeyKostrov]

>>...You did not show any code... Just run the test with values for matricies I gave. If I don't put values 5 and 7 repectively ( 1st values ) it calculates the C properly.

[SergeyKostrov]

A message to a Moderator: Many of my posts are delayed and please take a look. Thank you in advance.

[SergeyKostrov]

Note: Corrected a 1st value in the 2nd B matrix. Just run the test with: A[8][8]= 5 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 x B[8][8]= 7 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 and after that with: A[8][8]= 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 x B[8][8]= 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8

[SergeyKostrov]

>>...I detected some problem with SGEMM and DGEMM functions. In essence, there are wrong numbers in a 1st row of >>a resulting matrix C ( C[8][8] = A[8][8] x B[8][8] )... A manual verification was also completed for all values of the 1st row ( here are calculations for 1st and 8th values ): ( 5*7 ) + ( 2*1 ) + ( 3*1 ) + ( 4*1 ) + ( 5*1 ) + ( 6*1 ) + ( 7*1 ) + ( 8*1 ) = 70 ... ( 5*8 ) + ( 2*8 ) + ( 3*8 ) + ( 4*8 ) + ( 5*8 ) + ( 6*8 ) + ( 7*8 ) + ( 8*8 ) = 320 Note: Not 336 and I expect that Intel software developer needs to verify that simple test case with MKL functions for matrix multiplication in its own test application ( that is, independent ).

[SergeyKostrov]

>>...I detected some problem with SGEMM and DGEMM functions. In essence, there are wrong numbers in >>a 1st row of a resulting matrix C ( C[8][8] = A[8][8] x B[8][8] )... I wonder if Intel software engineers had a chance to verify ( independently ) these two functions with 8x8 matricies and data I've provided? Thanks in advance. PS: I'm very sorry but my test case can not be provided.

[Murat_G_Intel]

Hi Sergey,

As mecej4 already hinted, C programming language stores the 2D arrays in row-major fashion. If you would like to use 2D C arrays for dgemm, you need to call cblas_dgemm in row-major ordering. Namely, cblas_dgemm(order=CblasRowMajor, transa, transb, m, n, k ...). If you just call dgemm, it assumes column-major matrix ordering.

I created a simple program which calls cblas_dgemm and dgemm for your input values. As you can see from the output, cblas_dgemm gives expected values, whereas, dgemm call gives the "wrong" results mentioned in your post. However, the results are not actually wrong, it's just the dgemm uses what you think as rows of your matrix as columns. I'm attaching the simple program.

Here is the output from it:

Results from CBLAS row-major ordering
70, 80, 120, 160, 200, 240, 280, 320,
42, 72, 108, 144, 180, 216, 252, 288,
42, 72, 108, 144, 180, 216, 252, 288,
42, 72, 108, 144, 180, 216, 252, 288,
42, 72, 108, 144, 180, 216, 252, 288,
42, 72, 108, 144, 180, 216, 252, 288,
42, 72, 108, 144, 180, 216, 252, 288,
42, 72, 108, 144, 180, 216, 252, 288,

Results from Fortran-style col-major ordering
70, 84, 126, 168, 210, 252, 294, 336,
40, 72, 108, 144, 180, 216, 252, 288,
40, 72, 108, 144, 180, 216, 252, 288,
40, 72, 108, 144, 180, 216, 252, 288,
40, 72, 108, 144, 180, 216, 252, 288,
40, 72, 108, 144, 180, 216, 252, 288,
40, 72, 108, 144, 180, 216, 252, 288,
40, 72, 108, 144, 180, 216, 252, 288,

Best Wishes,

Efe

[SergeyKostrov]

>>...However, the results are not actually wrong, it's just the dgemm uses what you think as rows of your matrix as columns... Hi Efe, Thank you very much and this is exactly what I expected to confirm. Best regards, Sergey

[Abdul__J_]

Hi,

can you tell me that what would be the values for leading dimensions lda, ldb and ldc in case of non-square matrix-multiplication? If A is m x n B is n x k then C will be m x k so what would be lda, ldb and ldc for cblas_gemm?

[mecej4]

Abdul J wrote:
can you tell me that what would be the values for leading dimensions lda, ldb and ldc 

There is no single answer to this question, as these dimensions depend on whether the matrices are stored by rows or by columns. Please see the example .../mkl/examples/cblas/source/cblas_sgemmx.c.

If the matrix is stored by rows, as is more common in C, lda should be equal to the max. possible row length = max. number of columns of A. If you allocate the matrices to be the exact size needed and not larger, lda = the the number of columns of A.

And, it goes similarly for the other matrices.