C Матрица умножения динамически распределенных матриц

Question

Я работаю с определенным ограничением выделения памяти матрицы, которое создает матрицы как таковые:

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

Я храню эту матрицу внутри массива структур, таких как:

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

Учитывая, что это так, и я не создаю матрицы с помощью двумерного массива, такого как MATRIX[SIZE][SIZE]: как правильно умножить две матрицы, созданные таким образом?

В этой текущей реализации, если я хочу сделать что-то вроде вычитания, я делаю это следующим образом:

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

Это достаточно просто, потому что размеры matrix_index1 и matrix_index2 идентичны, а возвращаемая матрица также имеет одинаковые размеры.

Как можно добиться умножения матриц с помощью этого метода построения матриц?

Answer 1

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

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

Проверено на onlinegdb , пример вывода:

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 для представления размера объекта и количества массивов.

Answer 2

Есть несколько проблем с этим кодом: он нечитабелен и очень недружелюбен к кэшу, что означает медленный.

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

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

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

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

Затем при доступе к «искаженному массиву» мы можем привести указатель к типу двумерного массива и использовать его каждый раз, когда получаем доступ к данным.:

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