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表示对象大小和数组计数的类型。