What is wrong with this matrix multiplication?

Question

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);

}

Answer 1

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.