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，这样你得到相邻的记忆。如果您让编译器计算数组索引而不是手动计算，这通常也有助于提高性能。

我们可以通过在结构的末尾使用一个灵活的数组成员来显着简化这一点，然后在我们访问它时将其用作老式的“错位数组”。 “重整数组”在 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]) );

然后在访问“重整数组”时，我们可以转换为指向二维数组类型的指针，并在每次访问数据时使用该指针：

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