矩阵分治内存优化

Question

我目前正在 C 上实现 Strassen 的矩阵乘法算法并让它正常工作。但是，我在每次递归时传递子矩阵的方式远非最佳，因为我创建了一个新数组并从原始矩阵中复制元素。这是我的实现

// this algorithm assumes that both matrices have the same size and are square matrices.
int **matrix_multiply_divide_conquer(const int **A, const int **B, int length, int threshold)
{
  int **result;
  int i, j, k;
  // allocate result space
  result = (int **)malloc(sizeof(int *) * length);
  for (i = 0; i < length; i++)
  {
    result[i] = (int *)malloc(sizeof(int) * length);
  }

  if (length == threshold)
  {
    // compute the dot product
    matrix_multiply((const int **)A, (const int **)B, result, length, length, length, false);
    return result;
  }
  else
  {
    // stored in an int given that the matrix length is guaranteed to be a power of 2
    int nLength = length / 2;
    int i;
    int **a11;
    int **a12;
    int **a21;
    int **a22;

    int **b11;
    int **b12;
    int **b21;
    int **b22;

    // allocate new sub matrices
    a11 = (int **)malloc(sizeof(int *) * nLength);
    a12 = (int **)malloc(sizeof(int *) * nLength);
    a21 = (int **)malloc(sizeof(int *) * nLength);
    a22 = (int **)malloc(sizeof(int *) * nLength);

    b11 = (int **)malloc(sizeof(int *) * nLength);
    b12 = (int **)malloc(sizeof(int *) * nLength);
    b21 = (int **)malloc(sizeof(int *) * nLength);
    b22 = (int **)malloc(sizeof(int *) * nLength);

    for (i = 0; i < nLength; i++)
    {
      a11[i] = (int *)malloc(sizeof(int) * nLength);
      a12[i] = (int *)malloc(sizeof(int) * nLength);
      a21[i] = (int *)malloc(sizeof(int) * nLength);
      a22[i] = (int *)malloc(sizeof(int) * nLength);

      b11[i] = (int *)malloc(sizeof(int) * nLength);
      b12[i] = (int *)malloc(sizeof(int) * nLength);
      b21[i] = (int *)malloc(sizeof(int) * nLength);
      b22[i] = (int *)malloc(sizeof(int) * nLength);
    }

    a11 = get_sub_matrix((const int **)A, 0, 0, nLength);
    a12 = get_sub_matrix((const int **)A, 0, nLength, nLength);
    a21 = get_sub_matrix((const int **)A, nLength, 0, nLength);
    a22 = get_sub_matrix((const int **)A, nLength, nLength, nLength);


    b11 = get_sub_matrix((const int **)B, 0, 0, nLength);
    b12 = get_sub_matrix((const int **)B, 0, nLength, nLength);
    b21 = get_sub_matrix((const int **)B, nLength, 0, nLength);
    b22 = get_sub_matrix((const int **)B, nLength, nLength, nLength);

    // combine the results
    int **c11 = matrix_addition((const int **)matrix_multiply_divide_conquer((const int **)a11, (const int **)b11, nLength, 4), (const int **)matrix_multiply_divide_conquer((const int **)a12, (const int **)b21, nLength, 4), nLength, nLength);
    int **c12 = matrix_addition((const int **)matrix_multiply_divide_conquer((const int **)a11, (const int **)b12, nLength, 4), (const int **)matrix_multiply_divide_conquer((const int **)a12, (const int **)b22, nLength, 4), nLength, nLength);
    int **c21 = matrix_addition((const int **)matrix_multiply_divide_conquer((const int **)a21, (const int **)b11, nLength, 4), (const int **)matrix_multiply_divide_conquer((const int **)a22, (const int **)b21, nLength, 4), nLength, nLength);
    int **c22 = matrix_addition((const int **)matrix_multiply_divide_conquer((const int **)a21, (const int **)b12, nLength, 4), (const int **)matrix_multiply_divide_conquer((const int **)a22, (const int **)b22, nLength, 4), nLength, nLength);

    // combine result quarters
    combine_sub_matrix((const int **)c11, result, nLength, 0, 0);
    combine_sub_matrix((const int **)c12, result, nLength, 0, nLength);
    combine_sub_matrix((const int **)c21, result, nLength, nLength, 0);
    combine_sub_matrix((const int **)c22, result, nLength, nLength, nLength);
    return result;
  }
}

这是 get_sub_matrix 方法：

int **get_sub_matrix(const int **A, int verticalOffset, int horizontalOffset, int size)
{
  int **sub_matrix;
  sub_matrix = (int **)malloc(sizeof(int *) * size);
  int i, j;
  for (i = 0; i < size; i++)
  {
    //each row will contain a number of columns size
    sub_matrix[i] = (int *)malloc(sizeof(int) * size);
  }
  for (i = 0; i < size; i++)
  {
    for (j = 0; j < size; j++)
    {
      sub_matrix[i][j] = A[i + verticalOffset][j + horizontalOffset];
    }
  }
  return sub_matrix;
}

我确信这不是传递子矩阵的最佳方式，因为它们只是被读取而从不被修改，所以我想问你是否可以使用指针来解决这个问题以及如何做到这一点（因为我还在纠结他们）。

我尝试了以下方法，但没有成功：

int **a11 = &A[0][0];
int **a12 = &A[0][nLength];
int **a21 = &A[nLength][0] ;
int **a22 = &A[nLength][nLength];

int **b11 = &B[0][0];
int **b12 = &B[0][nLength];
int **b21 = &B[nLength][0] ;
int **b22 = &B[nLength][nLength];

感谢您帮助我。

Answer 1

您现在正在处理一组行，

 int **x = malloc(sizeof(int *) * size);

然后为每一行分配行内存。

for (i = 0; i < length; i++)
{
    x[i] = malloc(sizeof(int) * length);
}

相反，您可以使用将单个内存块视为二维数组：

int *x = malloc(sizeof(int) * size * size);

这允许您指向子矩阵的开头并将其传递给您的其他函数。但请注意您必须如何进行矩阵索引。因为编译器现在不知道您将单个块视为 2 维数组，所以您必须自己显式地进行索引。因此，给定一个大小为 length 的矩阵 length，您将 x[2][3]（第三行的第四个元素，因为 C 是从零开始的）寻址为：

x [ (2*length) + 3 ];

此外，在将指向子矩阵的指针传递给您的函数时，您必须制定一个约定来告知您的函数它在哪个子矩阵上运行。例如：

void my_matrix_op(int *X, int quadrant, int length)

我没有检查你所有的代码来确定这种方法是否真的适用于你的算法，但我希望这对你有帮助。