C Matrix Multiplication Dynamically Allocated Matrices

Question

I am working within a specific matrix memory allocation constraint that creates matrices as such:

float * matrix_data = (float *) malloc(rows * cols * sizeof(float));

I am storing this matrix inside an array of structs like this:

#define MAX_MATRICES 100

struct matrix{
    char matrixName[50];
    int rows;
    int columns;
    float* matrix_data;
};
typedef struct matrix matrix_t;

matrix_t our_matrix[MAX_MATRICES];

Given that this is the case, and I am not creating matrices by way of a 2d array like MATRIX[SIZE][SIZE]: what is the correct way to multiply two matrices created in this way?

With this current implementation if I want to do something like subtract, I do it as follows:

int max_col = our_matrix[matrix_index1].columns;
      free(our_matrix[number_of_matrices].matrix_data);
      our_matrix[number_of_matrices].data = (float *) malloc(our_matrix[matrix_index1].rows * our_matrix[matrix_index1].columns * sizeof(float)); 
      float *data1 = our_matrix[matrix_index1].matrix_data;
      float *data2 = our_matrix[matrix_index2].matrix_data;

      int col, row;
      for(col = 1; col <= our_matrix[matrix_index2].columns; col++){
        for(row = 1; row <= our_matrix[matrix_index2].rows; row++){
          our_matrix[number_of_matrices].matrix_data[(col-1) + (row-1) * max_col] =
            (data1[(col-1) + (row-1) * (max_col)]) - (data2[(col-1) + (row-1) * (max_col)]);  
        }
      }

This is simple enough because the dimensions of matrix_index1 and matrix_index2 are identical and the matrix that they return is also of identical dimensions.

How might I achieve matrix multiplication with this method of matrix construction?

Answer 1

There's several problems with this code: it is unreadable and it is very cache-unfriendly, meaning slow.

Regarding cache, you should always iterate over the outer-most dimension of the 2D array (doesn't matter if you call that one row or column), and you should only have a single malloc call in the code, so that you get adjacent memory. If you let the compiler calculate array indices instead of doing that manually, that usually helps performance as well.

We can simplify this significantly by using a flexible array member at the end of the struct, then using that as an old school "mangled array" whenever we access it. "Mangled array" meaning the array type is single dimension in the C syntax, but we access it as if it is a 2D array.

The struct type would be:

typedef struct 
{
  char   name[50];
  size_t columns;
  size_t rows;
  float  data[];
} matrix_t;

And we alloc memory for it once, with a single call:

matrix_t* matrix = malloc( sizeof *matrix + sizeof(float[c][r]) );

Then when accessing the "mangled array", we can cast to a pointer to 2D array type and use that one every time we access the data:

float (*data)[r] = (float(*)[r]) matrix->data;

Full example:

#include <stdlib.h>
#include <stdio.h>

typedef struct 
{
  char   name[50];
  size_t columns;
  size_t rows;
  float  data[];
} matrix_t;

int main (void)
{
  size_t c = 3;
  size_t r = 5;
  matrix_t* matrix = malloc(sizeof *matrix + sizeof(float[c][r]));

  float (*data)[r] = (float(*)[r]) matrix->data;
  for(size_t i=0; i<c; i++)
  {
    for(size_t j=0; j<r; j++)
    {
      data[i][j] = (float)i+j; // assign some random value
      printf("%.2f ", data[i][j]);
    }
    printf("\n");
  }

  free(matrix);
}

Answer 2

Write proper abstraction, then work your way up. It will be way easier:

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

struct matrix_s {
    char matrixName[50];
    size_t columns;
    size_t rows;
    float* data;
};

typedef struct matrix_s matrix_t;

void m_init(matrix_t *t, size_t columns, size_t rows) {
    t->rows = rows;
    t->columns = columns;
    t->data = calloc(rows * columns, sizeof(*t->data));
    if (t->data == NULL) abort();
}

size_t m_columns(const matrix_t *t) {
    return t->columns;
}

size_t m_rows(const matrix_t *t) {
    return t->rows;
}

// matrix_get 
// (x,y) = (col,row) always in that order
float *m_get(const matrix_t *t, size_t x, size_t y) {
    assert(x < m_columns(t));
    assert(y < m_rows(t));
    // __UNCONST
    // see for example `char *strstr(const char *haystack, ...` 
    // it takes `const char*` but returns `char*` nonetheless.
    return (float*)&t->data[t->rows * x + y];
}

// fill matrix with a fancy patterns just so it's semi-unique
void m_init_seq(matrix_t *t, size_t columns, size_t rows) {
    m_init(t, columns, rows);
    for (size_t i = 0; i < t->columns; ++i) {
        for (size_t j = 0; j < t->rows; ++j) {
            *m_get(t, i, j) = i + 100 * j;
        }
    }
}

void m_print(const matrix_t *t) {
    printf("matrix %p\n", (void*)t->data);
    for (size_t i = 0; i < t->columns; ++i) {
        for (size_t j = 0; j < t->rows; ++j) {
            printf("%5g\t", *m_get(t, i, j));
        }
        printf("\n");
    }
    printf("\n");
}

void m_multiply(matrix_t *out, const matrix_t *a, const matrix_t *b) {
    assert(m_columns(b) == m_rows(a));
    assert(m_columns(out) == m_columns(a));
    assert(m_rows(out) == m_rows(b));
    // Index from 0, not from 1
    // don't do `(col-1) + (row-1)` strange things
    for (size_t col = 0; col < m_columns(out); ++col) {
        for (size_t row = 0; row < m_rows(out); ++row) {
            float sum = 0;
            for (size_t i = 0; i < m_rows(a); ++i) {
                sum += *m_get(a, col, i) * *m_get(b, i, row);
            }
            *m_get(out, col, row) = sum;
        }
    }
}

int main()
{
    matrix_t our_matrix[100];

    m_init_seq(&our_matrix[0], 4, 2);
    m_init_seq(&our_matrix[1], 2, 3);

    m_print(&our_matrix[0]);
    m_print(&our_matrix[1]);

    m_init(&our_matrix[2], 4, 3);
    m_multiply(&our_matrix[2], &our_matrix[0], &our_matrix[1]);

    m_print(&our_matrix[2]);

    return 0;
}

Tested on onlinegdb, example output:

matrix 0xf9d010
    0     100   
    1     101   
    2     102   
    3     103   

matrix 0xf9d040
    0     100     200   
    1     101     201   

matrix 0xf9d060
  100   10100   20100   
  101   10301   20501   
  102   10502   20902   
  103   10703   21303   

Without the abstraction it's just a big mess. That would be something along:

  int col, row;
  for(col = 0; col < our_matrix[number_of_matrices].columns; col++){
    for(row = 0; row < our_matrix[number_of_matrices].rows; row++){
        for (size_t i = 0; i < our_matrix[matrix_index1].rows; ++i) {
            our_matrix[number_of_matrices].data[col * our_matrix[number_of_matrices].columns + row] = 
                our_matrix[matrix_index1].data[col * our_matrix[matrix_index1].columns + i] +
                our_matrix[matrix_index2].data[i * our_matrix[matrix_index2].columns + row];
        }  
    }
  }

Notes:

• Iterating from 0 up to < is way easier to read then all the (col-1) * ... + (row-1).
• Remember to check if the indexes are our of bounds. It's easy to do even with a simple assertion, ex. assert(row < matrix->rows && col < matrix->cols);
• Use size_t type to represent object size and array count.