Aside from that, here is my test code:

- Code: Select all
`//~returns random matrix`

double* matrix_random(double *A, int row, int col)

{

int i,j;

srand(time(NULL));

for(i=0;i<row;i++)

{

for(j=0;j<col;j++)

{

A[i*col + j] = rand() % 5;

}

}

return A;

}

dgemm_test(int rowsA, int colsB, int common)

{

int i,j,k;

char trans = 'T';

char norm = 'N';

double one = 1.0;

double zero = 0.0;

double neg = -1.0;

double *A, *B, *C, *D;

A = malloc(sizeof(double)*rowsA*common);

B = malloc(sizeof(double)*colsB*common);

C = calloc(rowsA*colsB,sizeof(double));

D = calloc(rowsA*colsB,sizeof(double));

A = matrix_random(A,rowsA,common);

B = matrix_random(B,common,colsB);

//~ Do a matrix multiply: A*B = C

for(i=0;i<rowsA;i++)

{

for(j=0;j<colsB;j++)

{

C[i*colsB+j]=0;

for(k=0;k<common;k++)

{

C[i*colsB+j]+=A[i*common+k]*B[k*colsB+j];

}

}

}

dgemm_(&trans, &trans, &rowsA, &colsB, &common, &one, A, &common, B, &colsB, &zero, C, &rowsA);

}

int main(void)

{

dgemm_test(4, 3, 2);

return 0;

}

Example output:

A, resulting 4 by 2 matrix

0, 0,

0, 4,

3, 1,

3, 2,

B, resulting 2 by 3 matrix

0, 0, 0,

4, 3, 1,

C, resulting 4 by 3 matrix

0, 0, 0,

16, 12, 4,

4, 3, 1,

8, 6, 2,

D, resulting 4 by 3 matrix

0, 16, 4,

8, 0, 12,

3, 6, 0,

4, 1, 2,

If I print C and D here, they are just the same except transposed. Using the 'n' function does not help. It only specifies whether or not the input (A and B) is transposed. Because the output of A*B is transposed, I have to perform a transpose function before I can add another matrix to it. Is there a way to get the C result to not be transposed?