Wrong results given by zgemm_() in clapack

Open discussion regarding features, bugs, issues, vendors, etc.

Wrong results given by zgemm_() in clapack

Postby robinchm » Wed Feb 19, 2014 5:15 pm

I am testing zgemm_() subroutine and write the following code to compute C=A*B, of which all are 2*2 matrices.
--------------------------------------------------------------------------

#include<stdio.h>
#include<f2c.h>
#include<clapack.h>
#include<blaswrap.h>

int main(){
char transa = 'N';
char transb = 'N';
integer m = 2;
integer n = 2;
integer k = 2;
doublecomplex alpha;
doublecomplex a[2][2];
integer lda = 2;
doublecomplex b[2][2];
integer ldb = 2;
doublecomplex beta;
doublecomplex c[2][2];
integer ldc = 2;

int i, j;

a[0][0].r = 1;
a[0][1].r = 1;
a[1][0].r = 1;
a[1][1].r = 1;
b[0][0].r = 1;
b[0][1].r = 0;
b[1][0].r = 1;
b[1][1].r = 2;

alpha.r = 1; alpha.i = 0;
beta.r = 0; beta.i = 0;

// Print A, B, C before multiplication
for(i=0;i<2;i++){
for(j=0;j<2;j++) printf("%+f i %+f ", a[i][j].r, a[i][j].i);
putchar('\n');
}
for(i=0;i<2;i++){
for(j=0;j<2;j++) printf("%+f i %+f ", b[i][j].r, b[i][j].i);
putchar('\n');
}
for(i=0;i<2;i++){
for(j=0;j<2;j++) printf("%+f i %+f ", c[i][j].r, c[i][j].i);
putchar('\n');
}

zgemm_(&transa, &transb, &m, &n, &k, &alpha, a, &lda, b, &ldb, &beta, c, &ldc);

// Print A, B, C after multiplication
for(i=0;i<2;i++){
for(j=0;j<2;j++) printf("%+f i %+f ", a[i][j].r, a[i][j].i);
putchar('\n');
}
for(i=0;i<2;i++){
for(j=0;j<2;j++) printf("%+f i %+f ", b[i][j].r, b[i][j].i);
putchar('\n');
}
for(i=0;i<2;i++){
for(j=0;j<2;j++) printf("%+f i %+f ", c[i][j].r, c[i][j].i);
putchar('\n');
}

return 0;
}

--------------------------------------------------------------------------
The screen output is:
--------------------------------------------------------------------------
+1.000000 i +0.000000 +1.000000 i +0.000000 // A before
+1.000000 i +0.000000 +1.000000 i +0.000000
+1.000000 i +0.000000 +0.000000 i +0.000000 // B before
+1.000000 i +0.000000 +2.000000 i +0.000000
+0.000000 i +0.000000 +0.000000 i +0.000000 // C before
+0.000000 i +0.000000 +0.000000 i +0.000000
+1.000000 i +0.000000 +1.000000 i +0.000000 // A after
+1.000000 i +0.000000 +1.000000 i +0.000000
+1.000000 i +0.000000 +0.000000 i +0.000000 // B after
+1.000000 i +0.000000 +2.000000 i +0.000000
+1.000000 i +0.000000 +1.000000 i +0.000000 // C after
+3.000000 i +0.000000 +3.000000 i +0.000000
--------------------------------------------------------------------------

As can be seen, the returned C is clearly wrong. Which point am I missing?
robinchm
 
Posts: 3
Joined: Wed Feb 19, 2014 5:08 pm

Re: Wrong results given by zgemm_() in clapack

Postby Julien Langou » Wed Feb 19, 2014 6:13 pm

Everything is transposed. LAPACK is essentially working in column major format, not row major format. You entered your data in row major format, and you print it in row major format. So you entered A' (as far as LAPACK sees it) and B', then LAPACK performed ( A' * B' ), and you print in row major so in the LAPACK words that's (A'* B' )'. So all is good. Well anyway. LAPACK works with column major format. This is your problem as far as I read your code.

You may want to use the CBLAS interface. This is a C interface for the BLAS and it supports row major format. See:
http://www.netlib.org/blas/blast-forum/cblas.tgz

Cheers,
Julien.
Julien Langou
 
Posts: 745
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Wrong results given by zgemm_() in clapack

Postby robinchm » Wed Feb 19, 2014 9:31 pm

Thanks! I thought clapack has converted from fortran's col-major to c's row-major. Never thought about this!

I can transpose A and B with "transa" and "transb" easily, and in this case is C the correct row-major matrix I need?

Regarding CBLAS, it has to be compiled by "gfortran" as shown in the example, while clapack is compiled with C. (I do not have the root access so flags like -lcblas do not work quite well) Anyway since clapack covers blas I am happy with it, once I got this col-major issue done.

Thanks again!
robinchm
 
Posts: 3
Joined: Wed Feb 19, 2014 5:08 pm

Re: Wrong results given by zgemm_() in clapack

Postby Julien Langou » Wed Feb 19, 2014 9:51 pm

I can transpose A and B with "transa" and "transb" easily, and in this case is C the correct row-major matrix I need?

Yup. Just that then you want to do B^T * A^T so do not forget to switch the operands A and B in the function call. If you do this, you will pretty much do exactly what the CBLAS would do for you. (But transpose and switch operands.) But yes you can do it yourself. Julien.
Julien Langou
 
Posts: 745
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA


Return to User Discussion

Who is online

Users browsing this forum: No registered users and 1 guest