4
$\begingroup$

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

}
Improve this question
$\endgroup$
4
$\begingroup$

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.

Improve this answer
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – drjrm3 Jul 6 '15 at 23:23
  • 2
    $\begingroup$ The way your matrices are stored should be modified - A[irncA + ic]+ becomes A[ir + icnrA] for example (writes to column major format). $\endgroup$ – Jesse Chan Jul 6 '15 at 23:24
  • 1
    $\begingroup$ @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. $\endgroup$ – Kirill Jul 6 '15 at 23:25
  • $\begingroup$ 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! $\endgroup$ – drjrm3 Jul 6 '15 at 23:35
  • 2
    $\begingroup$ In case you need or want to keep them row-major, the easiest and most performing solution is simply toggling the TRANS parameter. $\endgroup$ – Federico Poloni Jul 7 '15 at 21:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.