C Матричное умножение Динамически размещаемые матрицы

Напишите правильную абстракцию, а затем работайте над собой. Будет намного проще:

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

Протестировано на онлайнgdb, пример вывода:

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   

Без абстракции это просто большой беспорядок. Это было бы что-то вроде:

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

Примечания:

• Итерация от 0 до < читается намного легче, чем все (col-1) * ... + (row-1).
• Не забудьте проверить, соответствуют ли индексы нашим границам. Это легко сделать даже с помощью простого утверждения, например. assert(row < matrix->rows && col < matrix->cols);
• Используйте тип size_t для представления размера объекта и количества массивов.

С этим кодом есть несколько проблем: он нечитаем и очень неудобен для кэширования, то есть медленный.

Что касается кеша, вы всегда должны перебирать самое внешнее измерение 2D-массива (не имеет значения, вызываете ли вы эту одну строку или столбец), и у вас должен быть только один вызов malloc в коде, так что вы получаете смежные Память. Если вы позволите компилятору вычислять индексы массива вместо того, чтобы делать это вручную, это также обычно повышает производительность.

Мы можем значительно упростить это, используя элемент гибкого массива в конце структуры, а затем используя его как «искаженный массив» старой школы всякий раз, когда мы к нему обращаемся. «Искаженный массив» означает, что тип массива является одномерным в синтаксисе C, но мы обращаемся к нему, как если бы это был двумерный массив.

Тип структуры будет:

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

И выделяем память под него один раз, одним вызовом:

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

Затем при доступе к «искаженному массиву» мы можем привести указатель к типу 2D-массива и использовать его каждый раз, когда мы обращаемся к данным:

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

Полный пример:

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