# What is wrong with this matrix multiplication?

I am attempting to write a matrix multiplication routine because I need to do some analysis in CUDA and I want to validate it with CPU code. I am trying to use dgemm to do this, but I seem to be doing it wrong.

My matrices are constructed as 1D arrays since the same thing is done in CUDA kernels, so I am attempting to call dgemm with 1D matrices.

When I use random numbers as my input matrices I get:

A = [ ...
0.840187717154710        0.394382926819093        0.783099223758606
0.798440033476073        0.911647357936784        0.197551369293384
0.335222755714889        0.768229594811904        0.277774710803188
];
B = [ ...
0.553969955795431        0.477397051862160        0.628870924761924
0.364784472791843        0.513400910195616        0.952229725174713
0.916195068003701        0.635711727959901        0.717296929432683
];
myC = [ ...
1.057423514989923        1.136811309272234        0.702808322919262
1.035616345918396        1.343436407318801        0.651590826816703
1.517807789358121        1.491925339118482        1.042304316032681
];


as my output, which, as I found out, is actually B*A. This happens consistently - when I attempt to multiply square matrices, A and B, I actually getB*A. I am testing this because when I do non-square matrices, it's just a mess - I can't even trace what is happening!

My C code is as follows:

void test_host_matmult() {
fprintf(stderr, "Inside test_host_matmult\n");
int nrA = 3;
int ncA = 3;

int nrB = 3;
int ncB = 3;

double* A;
double* B;
double* C;

A = (double*)malloc(nrA*ncA*sizeof(double));
B = (double*)malloc(nrB*ncB*sizeof(double));
C = (double*)malloc(nrA*ncB*sizeof(double));

//int icount = 0;
for(int ir = 0; ir < nrA; ir++) {
for(int ic = 0; ic < ncA; ic++) {
//A[ir*ncA + ic] = 10*(ir+1) + (ic+1);
A[ir*ncA + ic] = rand_double(0.0, 1.0);
//fprintf(stderr, "%i ... %i\n", icount++, 10*(ir+1) + (ic+1));
}
}
for(int ir = 0; ir < nrB; ir++) {
for(int ic = 0; ic < ncB; ic++) {
//B[ir*ncB + ic] = 10*(ir+1) + (ic+1);
B[ir*ncB + ic] = rand_double(0.0, 1.0);
}
}

matmult(A, nrA, ncA, B, nrB, ncB, C, nrA, ncB);

fprintf(stderr, "A = [ ...\n");
for(int ir = 0; ir < nrA; ir++) {
for(int ic = 0; ic < ncA; ic++) {
fprintf(stderr, "%25.15f", A[ir*ncA + ic]);
}
fprintf(stderr, "\n");
}
fprintf(stderr, "];\n");

fprintf(stderr, "B = [ ...\n");
for(int ir = 0; ir < nrB; ir++) {
for(int ic = 0; ic < ncB; ic++) {
fprintf(stderr, "%25.15f", B[ir*ncB + ic]);
}
fprintf(stderr, "\n");
}
fprintf(stderr, "];\n");

fprintf(stderr, "myC = [ ...\n");
for(int ir = 0; ir < nrA; ir++) {
for(int ic = 0; ic < ncB; ic++) {
fprintf(stderr, "%25.15f", C[ir*ncB + ic]);
}
fprintf(stderr, "\n");
}
fprintf(stderr, "];\n");
}


and here is the matmult function:

void matmult( double* A, int nrA, int ncA, double* B, int nrB, int ncB, double* C, int nrC, int ncC) {
// (nrA x ncA) . (nrB x ncB ) = (nrA x ncB)
if( !(nrA == nrC && ncB == ncC) ) {
fprintf(stderr, "ERROR: incorrect matrix sizes!!\n");
}
char TRANS = 'N';
int M = nrA;
int N = ncB;
int K = ncA;
double alpha = 1.0;
double beta = 0.0;
int lda = nrA;
int ldb = nrB;
int ldc = nrC;
dgemm_(&TRANS, &TRANS, &M, &N, &K, &alpha, A, &lda, B, &ldb, &beta, C, &ldc);

}


Note: I haven't run your code.

Perhaps this is a problem with the row-major/column-major conventions at play here: http://docs.nvidia.com/cuda/cublas/index.html#data-layout.

You seem to use row-major matrices, but BLAS uses the column-major convention. When you pass row-major $A$ and $B$ to dgemm, it implicitly interprets them as $A^t$ and $B^t$, giving the output matrix $A^tB^t$. When you read it back, you treat it as row-major, thus implicitly applying another transpose: $(A^tB^t)^t = BA$.

I suspect you also might have a problem in your code with non-square matrices, because then the leading dimension parameters and matrix sizes that you pass to dgemm would be wrong.

• I have suspected it was a row major/column major issue since the beginning but am unable to figure out what I need to change to get the expected results. Could you possibly identify where/what I need to change here? – drjrm3 Jul 6 '15 at 23:23
• The way your matrices are stored should be modified - A[irncA + ic]+ becomes A[ir + icnrA] for example (writes to column major format). – Jesse Chan Jul 6 '15 at 23:24
• @Laurbert515 If you insist on row-major order, can't you just pass the matrices in the "wrong" order ($B\times A$ instead of $A\times B$) to dgemm? Just make sure ld and matrix sizes are correct for all the implicit matrix transposes. – Kirill Jul 6 '15 at 23:25
• Thanks for the help. I wasn't 'insisting' on anything, I just didn't realize I needed to print them out in column major order as well. Silly, I know, but once I converted everything to column major, it works perfectly! – drjrm3 Jul 6 '15 at 23:35
• In case you need or want to keep them row-major, the easiest and most performing solution is simply toggling the TRANS parameter. – Federico Poloni Jul 7 '15 at 21:56