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