dgemm routine works incorrectly

[vich_ka]

Hi! Please, help me:) I don't know what's wrong... maybe my hands... maybe my programm:)))

So, elementary code how i used dgemm:

static double A[2][3],B[3][2],C[2][2];

int main(){
int i,j;

for(i=0;i<2;i++){
for(j=0;j<3;j++){
A=(double)(j+1);
}
}
for(i=0;i<3;i++){
for(j=0;j<2;j++){
B=(double)(i+1);
}
}
dgemm_64('N','N',2,2,3,1.0,&A[0][0],2,&B[0][0],3,0.0,&C[0][0],2);
for(i=0;i<2;i++){
for(j=0;j<2;j++){
printf("%le ",C);
}
}
}

Answer = [8 9;17 16], but it must be [14 14;14 14].

I don't know what to do:(

[Evgueni_P_Intel]

Quoting - vich_ka

Hi! Please, help me:) I don't know what's wrong... maybe my hands... maybe my programm:)))

So, elementary code how i used dgemm:

//begin

static double A[1][3],B[3][1],C[1][1];
int main(){
int i,j;
for(i=0;i<2;i++){
A[0]=(double) (i+1);
}

for(i=0;i<3;i++){
B[0]=(double) (i+1);
}
dgemm('N','N',1,1,3,1.0,&A[0][0],1,&B[0][0],3,0.0,&C[0][0],1);
printf("%le ",C[0][0]);
}

//end

Answer=5. But "in my opinion" it should be 14!

I don't know what to do:(

Most likely, increment 3 is incorrect for B.
The correct incrementwould be1.

[barragan_villanueva_]

Hi,

In MKL dgemm has FORTRAN-like interface,and all arguments should be passed by reference.
So, correct example is to be as follows:

#include
#include

static double A[1][3],B[3][1],C[1][1];
int main(){
int i,j;
int one=1;
int three=3;
double a=1.0;
double b=0.0;

for(i=0;i<3;i++){
A[0]=(double) (i+1);
}

for(i=0;i<3;i++){
B[0]=(double) (i+1);
}
dgemm("N","N",&one,&one,&three,&a,&A[0][0],&one,&B[0][0],&three,&b,&C[0][0],&one);
printf("%le ",C[0][0]);
}

which gives 1.400000e+01

Thanks
-- Victor

[vich_ka]

Quoting - Evgueni Petrov aka espetrov (Intel)

Most likely, increment 3 is incorrect for B.
The correct incrementwould be1.

I don't understand why it should be 1... but i tried to do it and have an error...

[vich_ka]

Quoting - Victor Pasko (Intel)

Hi,

In MKL dgemm has FORTRAN-like interface,and all arguments should be passed by reference.
So, correct example is to be as follows:

#include
#include

static double A[1][3],B[3][1],C[1][1];
int main(){
int i,j;
int one=1;
int three=3;
double a=1.0;
double b=0.0;

for(i=0;i<3;i++){
A[0]=(double) (i+1);
}

for(i=0;i<3;i++){
B[0]=(double) (i+1);
}
dgemm("N","N",&one,&one,&three,&a,&A[0][0],&one,&B[0][0],&three,&b,&C[0][0],&one);
printf("%le ",C[0][0]);
}

which gives 1.400000e+01

Thanks
-- Victor

thanks for your answer.

In my code was silly error

for(i=0;i<2;i++){
A[0]=(double) (i+1);
}

instead of

for(i=0;i<3;i++){
A[0]=(double) (i+1);
}

you've seen it. so, i set it right and now i have right answer = 14.

but my general problem isn't solved. i've corrected the first topic. please look it over again. there is new example, which incorrect.

i use "sunperf.h" library, which include "mkl_blas.h". i tried to pass arguments by reference, but have an error.

[Evgueni_P_Intel]

vich_ka, yes i was wrong about the increment.

[barragan_villanueva_]

Hi,

Your new code corresponds to the following operation:

A= {{1, 2, 3}, {1, 2, 3}}

B= {{1, 1}, {2, 2}, {3, 3}}

Therefore:

|1 3 2| |1 2| |8 17|
* |1 3| =
|2 1 3| |2 3| |9 16|

not as you want

|1 2 3||1 1| |14 14|
* |2 2| =
|1 2 3||3 3| |14 14|

So, the problem is in data layout

Thanks
-- Victor

[vich_ka]

Quoting - Victor Pasko (Intel)

Hi,

Your code corresponds to the following operation:

A= {{1, 2, 3}, {1, 2, 3}}

B= {{1, 1}, {2, 2}, {3, 3}}

Therefore:

|1 3 2| |1 2| |8 17|
* |1 3| =
|2 1 3| |2 3| |9 16|

not as you want

|1 2 3||1 1| |14 14|
* |2 2| =
|1 2 3||3 3| |14 14|

So, the problem is in data layout

Thanks
-- Victor

Thanks a lot, we don't know this feature. I thought that if i write "A[2][3]" it means that A is two-dimensional array instead of string:) Even in the help i read: " A (input) DOUBLE PRECISION array of DIMENSION ( LDA, ka )".

You really save us:) Thanks!

[Vladimir_Koldakov__I]

Hi,

It's really the MKL dgemm assumes column-major matrix ordering so the input data should be reorganized.

There is one more option: using cblas_dgemm interface that allows both column- and row-major matrix ordering:

void cblas_dgemm(const CBLAS_ORDER Order, const CBLAS_TRANSPOSE TransA,

const CBLAS_TRANSPOSE TransB, const MKL_INT M, const MKL_INT N,

const MKL_INT K, const double alpha, const double *A,

const MKL_INT lda, const double *B, const MKL_INT ldb,

const double beta, double *C, const MKL_INT ldc);

Thanks,

Vladimir

[barragan_villanueva_]

Quoting - vich_ka

i use "sunperf.h" library, which include "mkl_blas.h". i tried to pass arguments by reference, but have an error.

It looks like using "sunperf.h" calls directly MKL FORTRAN dgemm not cblas_dgemm where data layout can be C or FORTRAN.