Cblas_sgemm produces wrong results

226 Views Asked by At

I have a problem using the sgemm function of cblas.

Here is the code:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cblas.h>

#define MATRIX_DIM 5

int main(){

    float *a_mat = calloc(MATRIX_DIM*MATRIX_DIM, sizeof(float));
    float *b_mat = calloc(MATRIX_DIM, sizeof(float));
    float *c_mat = calloc(MATRIX_DIM, sizeof(float));
    int i,j;

    for(i=0; i<MATRIX_DIM*MATRIX_DIM; i++) {
        a_mat[i] = 1.0f;
        b_mat[i] = 1.0f;
        c_mat[i] = 0.0f;
    }

    cblas_sgemm(CblasRowMajor, CblasNoTrans,
                 CblasNoTrans, MATRIX_DIM, MATRIX_DIM,
                 MATRIX_DIM, 1.0, a_mat,
                 MATRIX_DIM, b_mat, MATRIX_DIM,
                 1.0, c_mat, MATRIX_DIM);

    //RESULT PRINTING
    printf("Printing A MATRIX: \n");
    for(i=0; i<MATRIX_DIM; i++) {
        for(j=0; j<MATRIX_DIM; j++){
                printf("%0.1f ", a_mat[i*MATRIX_DIM+j]);
        }
        printf("\n");
    }
    printf("Printing B MATRIX: \n");
    for(i=0; i<MATRIX_DIM; i++) {
            for(j=0; j<MATRIX_DIM; j++){
                printf("%0.1f ", b_mat[i*MATRIX_DIM+j]);
        }
        printf("\n");
    }

    printf("\nPrinting the Results: \n");
    for(i=0; i<MATRIX_DIM;i++){
        for(j=0; j<MATRIX_DIM; j++){
                printf("%0.1f ", c_mat[i*MATRIX_DIM+j]);
        }
        printf("\n");
    }

    free(a_mat);
    free(b_mat);
    free(c_mat);

    return 0;
}

I am pretty some of the arguments I put in are wrong, but I really don't know which. The results should be a 5x5 matrix filled with 5.0. Instead the program responds with this:

6.0 6.0 6.0 16.0 86.0 
6.0 6.0 6.0 16.0 86.0 
16.0 36.0 6.0 46.0 86.0 
16.0 36.0 5.0 45.0 85.0 
20.0 80.0 5.0 45.0 85.0 

I know the rowmajor order or the transpose arguments might be wrong and I will figure those out later, but in this particular multiplication the answer should be 5.0 either way.

1

There are 1 best solutions below

4
On BEST ANSWER

Thanks to @AndrasDeak in the comments, all I needed was to allocate more space on the two matrices, which I previously overlooked.

So basically change:

float *b_mat = calloc(MATRIX_DIM, sizeof(float));
float *c_mat = calloc(MATRIX_DIM, sizeof(float));

To:

float *b_mat = calloc(MATRIX_DIM*MATRIX_DIM, sizeof(float));
float *c_mat = calloc(MATRIX_DIM*MATRIX_DIM, sizeof(float));

Since those are supposed to be 2-dimensional matrices and not vectors.