topic Vector-vector product, cblas_sgemm problem in Intel® oneAPI Math Kernel Library

Vector-vector product, cblas_sgemm problem

apolo74 — Thu, 21 Jul 2011 10:13:47 GMT

I have a couple of vectors ( ax=[5, 2, 6, 1, 7] and ay=[8, 1, 7] ). The result of doing ay' * ax should be

[ 40 16 48 8 56

5 2 6 1 7

35 14 42 7 49]

But when running this operation using MKL (CBLAS):

cblas_sgemm( CblasRowMajor, CblasTrans, CblasNoTrans, 3,

5, 5, 1.0f, ay, 3, ax, 5,0.0f, aa, 5 );

being 'aa' the matrix where I'm planning to save the result of size 15. I get ALMOST what I'm supposed to get, since every time I run the code the element in the middle (row 2, col 3, aa[7] = 6) gives noise values, but all the other values are as they should. If I plot 'aa' before calling cblas_sgemm(), then some times I get the right result and some others I get 6.0001 or some others it is aa[8]=1.0001. So it seems that I'm missing some initialization, but according to the examples I don't see any initizalization either. Can someone tell me what I'm doing wrong? Thanks for any help,

Boris

Vector-vector product, cblas_sgemm problem

mecej4 — Fri, 22 Jul 2011 03:32:37 GMT

Since CBlas is merely an interface to the Fortran-77 BLAS routines, it is easy to get some of the parameters wrong, especially to a user unfamiliar with Fortran-77 conventions regarding multidimensional arrays.

That is why there is a cautionary note at the beginning of Appendix D:

"Warning

Users of the CBLAS interface should be aware that the CBLAS are just a C interface to the BLAS, which is based on the FORTRAN standard and subject to the FORTRAN standard restrictions. In particular, the output parameters should not be referenced through more than one argument."

"In the descriptions of CBLAS interfaces, links provided for each function group lead to the descriptions of the respective Fortran-interface BLAS functions."

From the MKL documentation page for ?GEMM, we have

A, B and C are matrices:

op(A) is an m-by-k matrix,

op(B) is a k-by-n matrix,

C is an m-by-n matrix.

For your example, op(A)=ay' is 3 X 1, op(B)=ax is 1 X 5. Thus, arguments #4, #5, #6, which are m,n,k, should be 3,5,1, not 3,5,5 as you have.

One could probe the consequences of how the incorrect arguments that you provided caused the mystefying behavior in the results that you observed, but I do not think that it is worthwhile. Your conjecture as to initialization being the cause turns out to be incorrect since, when argument beta = 0, the matrix C need not be initialized, as the documentation clearly states.

The next time you post a similar query, it would make it easier to analyze the reported problem if you gave a complete example that could be compiled and run to reproduce the problem.

[cpp]#include 
#include 

main(){
float ax[]={5, 2, 6, 1, 7}, ay[]={8, 1, 7}; 
float aa[5*3]; 
int i,j;
cblas_sgemm( CblasRowMajor, CblasTrans, CblasNoTrans, 
             3, 5, 1, 
			 1.0f, 
			 ay, 3, ax, 5, 
			 0.0f, 
			 aa, 5 );
for(i=0; i<3; i++){
   for(j=0; j<5; j++)printf("%9.4f ",aa[5*i+j]);
   printf("n");
   }
}   
[/cpp]

Vector-vector product, cblas_sgemm problem

TimP — Fri, 22 Jul 2011 12:36:57 GMT

I'm a little hesitant to rely on the beta=0f. eliminating the need for initialization, although I do this myself. If the uninitialized data patterns happen to include NaN or Inf, you would be depending on the library code explicitly zeroing the data for that case. Also, it's conceivable that the result might depend on whether gradual underflow is enabled (no-ftz), as unitialized data typically includes de-normals.

Vector-vector product, cblas_sgemm problem

Murat_G_Intel — Fri, 22 Jul 2011 17:41:57 GMT

Tim, I understand your hesitation. However, when the beta=0.0, the library code should ignore the values of the C matrix. In other words, the values of C Matrix are overwritten. Therefore, as long as you set beta=0.0, Nan, Inf or subnormals inC matrix should not change the results.

Vector-vector product, cblas_sgemm problem

mecej4 — Fri, 22 Jul 2011 20:24:10 GMT

In other words, the code would do something similar to

if(beta != 0.0){
for(i...) for(j...)
C += ...
}
else{
for(i...) for(j...)
C = ...
}

I have seen code with this structure in the Fortran-77 sources of BLAS and Lapack.

Vector-vector product, cblas_sgemm problem

Murat_G_Intel — Fri, 22 Jul 2011 20:38:42 GMT

Precisely, mecej4. The code treats beta == 0.0 as a special case and does C = alpha*op(A)*op(B). It DOES NOT explicity multiply the entries of C with 0.0.

This behaviour is described in MKL Reference Manual, under the entry for C matrix.All Level3 BLAS should havethe same behaviour except from TRSM and TRMM. There is no C matrix for TRSM and TRMM.

Vector-vector product, cblas_sgemm problem

apolo74 — Mon, 25 Jul 2011 10:28:31 GMT

Dear mecej4 and all, thanks for your help and sorry for bothering you with this simple error (for you at least)... in my defense I started using MKL a couple of days ago and I'm not that familiar with the logic and documentation. Thanks again for the explanation and next time I'll build some simple code of my problem.

Boris