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 的维度相同，并且它们 return 的矩阵也具有相同的维度。

如何使用这种矩阵构造方法实现矩阵乘法？

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

这段代码有几个问题：它不可读而且非常 cache-unfriendly，意味着很慢。

关于缓存，您应该始终遍历二维数组的 outer-most 维度（调用一行或一列无关紧要），并且您应该只在代码，以便您获得相邻的内存。如果让编译器计算数组索引而不是手动计算，通常也有助于提高性能。

我们可以通过在结构的末尾使用一个灵活的数组成员来显着简化这一点，然后在我们访问它时将其用作旧学校 "mangled array"。 "Mangled array"表示数组类型在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]) );

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

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

C 矩阵乘法动态分配矩阵

C Matrix Multiplication Dynamically Allocated Matrices

c

matrix

ansi-c