c - 在 C 中超过 1024x1024 矩阵时，总线错误/核心转储

问题描述

在 N > 1024 上运行此代码时，我收到总线错误/核心转储错误。我正在使用远程 HPC 和 gcc/8.1。这是一个矩阵乘法 NxN。我不明白错误来自哪里。特别是为什么较小的 N 没有任何问题。我之前有其他代码运行到 2^20。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

#define N 2048

float *A[N], *B[N];
int i, j, k, count = 0;

float** matrix_create(int n){
  float** M = malloc(n * sizeof(float*));
  for (i = 0; i < n; i++)
    M[i] = (float*)malloc(n * sizeof(float));
  return M;
}

float** add(float* M1[], float* M2[], int n){
  float** M3 = matrix_create(n);
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      M3[i][j] = M1[i][j] + M2[i][j];
  return M3;
}

float** sub(float* M1[], float* M2[], int n){
  float** M3 = matrix_create(n);
  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      M3[i][j] = M1[i][j] - M2[i][j];
  return M3;
}

void print(float* M[], int n){
  for (i = 0; i <  n; i++){
    for (j = 0; j < n; j++)
      printf("%f\t", M[i][j] );
    printf("\n");
    }
}

float** strassen_multiply(float* A[], float* B[], int n){
  if(n == 1 ){
    float** C = matrix_create(n);
    C[0][0] = A[0][0] * B[0][0];
    return C;
  }
  count++;
  float** C = matrix_create(n);
  int k = n/2;
  /** Creating sub matrecies**/
  float** A11 = matrix_create(k);
  float** A12 = matrix_create(k);
  float** A21 = matrix_create(k);
  float** A22 = matrix_create(k);
  float** B11 = matrix_create(k);
  float** B12 = matrix_create(k);
  float** B21 = matrix_create(k);
  float** B22 = matrix_create(k);

  /**Dividing the Data Matrecies A & B**/
  for(i = 0; i < k; i++)
      for(j = 0; j < k; j++){
          A11[i][j] = A[i][j];
          A12[i][j] = A[i][k+j];
          A21[i][j] = A[k+i][j];
          A22[i][j] = A[k+i][k+j];
          B11[i][j] = B[i][j];
          B12[i][j] = B[i][k+j];
          B21[i][j] = B[k+i][j];
          B22[i][j] = B[k+i][k+j];
      }

  float** P1 = strassen_multiply(A11, sub(B12, B22, k), k);
  float** P2 = strassen_multiply(add(A11, A12, k), B22, k);
  float** P3 = strassen_multiply(add(A21, A22, k), B11, k);
  float** P4 = strassen_multiply(A22, sub(B21, B11, k), k);
  float** P5 = strassen_multiply(add(A11, A22, k), add(B11, B22, k), k);
  float** P6 = strassen_multiply(sub(A12, A22, k), add(B21, B22, k), k);
  float** P7 = strassen_multiply(sub(A11, A21, k), add(B11, B12, k), k);

  float** C11 = sub(add(add(P5, P4, k), P6, k), P2, k);
  float** C12 = add(P1, P2, k);
  float** C21 = add(P3, P4, k);
  float** C22 = sub(sub(add(P5, P1, k), P3, k), P7, k);

  for(i = 0; i < k; i++)
      for(j = 0; j < k; j++){
          C[i][j] = C11[i][j];
          C[i][j+k] = C12[i][j];
          C[k+i][j] = C21[i][j];
          C[k+i][k+j] = C22[i][j];
      }

      for(i = 0; i < k; i++){
          free( A11[i]);
          free( A12[i]);
          free( A21[i]);
          free( A22[i]);
          free( B11[i]);
          free( B12[i]);
          free( B21[i]);
          free( B22[i]);
          free( P1[i]);
          free( P2[i]);
          free( P3[i]);
          free( P4[i]);
          free( P5[i]);
          free( P6[i]);
          free( P7[i]);
          free( C11[i]);
          free( C12[i]);
          free( C21[i]);
          free( C22[i]);
      }

      free( A11);
      free( A12);
      free( A21);
      free( A22);
      free( B11);
      free( B12);
      free( B21);
      free( B22);
      free( P1);
      free( P2);
      free( P3);
      free( P4);
      free( P5);
      free( P6);
      free( P7);
      free( C11);
      free( C12);
      free( C21);
      free( C22);

      return C;

}


int main(){

  int i,j, k;
  struct timeval begin, end;



  for (i = 0; i < N; i++)
    A[i] = (float*)malloc(N * sizeof(float));
  for (i = 0; i < N; i++)
    B[i] = (float*)malloc(N * sizeof(float));


   for (i = 0; i <  N; i++)
     for (j = 0; j < N; j++){
        A[i][j] = -1+2*((float)rand())/RAND_MAX;
        B[i][j] = -1+2*((float)rand())/RAND_MAX;
        }

    float** C = matrix_create(N);
    gettimeofday(&begin, 0);
    C = strassen_multiply(A, B, N);
    gettimeofday(&end, 0);
    long seconds = end.tv_sec - begin.tv_sec;
    long microseconds = end.tv_usec - begin.tv_usec;
    double elapsed = seconds + microseconds*1e-6;
    printf("number of recursion: %d\n\n", count);
    printf("Total wall time: %f\n", elapsed);
}

标签： c

解决方案

^{从聊天中转移一些评论。}

诊断

您没有检查您的内存是否已成功分配。你不知道是否一切正常。
您从两个 2048x2048float矩阵开始。然后，您的strassen_multiply()函数 (1) 创建 8 个矩阵，每个矩阵的大小只有一半（就行数和列数而言），加载它们，然后连续递归 7 次。这些递归中的每一个都会创建大量矩阵——我还没有坐下来计算所需的总空间，但它会是相当可观的。您确实需要检查您的内存分配是否有效。可能是您的 64 位机器有足够的空间，这不是问题（两个初始矩阵需要 32 MiB 的数据，这可能没问题）。
你有这样的电话
```
float** P1 = strassen_multiply(A11, sub(B12, B22, k), k);
float** P2 = strassen_multiply(add(A11, A12, k), B22, k);
```
sub()您无法释放对and的嵌套调用返回的矩阵add()。你不能不释放那段记忆。因此，您正在泄漏大量内存。你需要一个函数来释放你的矩阵——并且可以说是一个记录矩阵大小的矩阵结构类型，因为你需要函数中的大小来释放矩阵。
您可以通过检查返回的空指针来检查内存是否已分配malloc()。在大多数系统上，这是可靠的。在 Linux 上，它具有 OOM（内存不足）管理器，并且往往会返回一个非空指针，然后当您尝试使用它告诉您可用但实际上不可用的内存时会崩溃。我认为这是非常不受欢迎的行为，但是……如果您未能分配其中一个行，请不要忘记释放该矩阵中任何先前分配的行。
你不能使用全局矩阵；你必须从函数返回矩阵，并且你有递归函数，所以全局矩阵不起作用。您需要将矩阵（都是方阵）转换为如下结构：
```
struct Matrix
{
    int     size;
    float **data;
};
```
你现有的两个全局指针数组float应该被替换——否则，你需要特殊的代码来释放分配给它们的内存。

在main()你有：

float** C = matrix_create(N);
…
C = strassen_multiply(A, B, N);

所以你泄漏了一个全尺寸的矩阵。

返回矩阵的函数将返回一个矩阵结构，而接受两个矩阵参数的函数将采用两个指向（常量）矩阵结构的指针作为参数。概述的矩阵结构非常小，返回指向矩阵结构的指针没有任何好处。
在您当前的代码中main()，您应该有：
```
float **A = matrix_create(N);
float **B = matrix_create(N);
```
您在 main() 中的矩阵C应使用以下命令创建：
```
float **C = strassen_multiply(A, B, N);
```
矩阵C从来都不是全局的。
matrix_create()像现在一样使用。请记住在调用add()or的函数中释放返回值sub()，这也意味着您需要将这些中间结果保存在局部变量中，以便您可以释放它们。
您正在为数组索引使用全局变量i, j。k所有的地狱都会崩溃。数组索引必须是局部变量，特别是如果您使用递归。
这意味着您必须在每个函数中声明循环变量。你应该写
```
for (int i = 0; i < n; i++)
```
或每个循环的等效项。这将比使用全局变量更有效；它还使您的代码更有可能正确。就目前而言，您没有丝毫机会使代码正确。

处方

将这些点放在一起会产生如下代码：

#include <assert.h>
#include <errno.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>

#ifndef N
#define N 128
#endif

typedef struct Matrix
{
    int    size;
    float **data;
} Matrix;

static int count = 0;
static size_t cnt_create = 0;
static size_t cnt_destroy = 0;
static size_t cnt_add = 0;
static size_t cnt_sub = 0;

static void err_nomemory(const char *file, const char *func, int line, size_t size)
{
    fprintf(stderr, "%s:%s():%d: out of memory attempting to allocate %zu bytes "
            "(%d: %s)\n", file, func, line, size, errno, strerror(errno));
    exit(EXIT_FAILURE);
}

static void matrix_destroy(Matrix *M)
{
    cnt_destroy++;
    for (int i = 0; i < M->size; i++)
        free(M->data[i]);
    free(M->data);
}

static Matrix matrix_create(int n)
{
    cnt_create++;
    Matrix M = { .size = n, .data = malloc(n * sizeof(float *)) };
    if (M.data == NULL)
        err_nomemory(__FILE__, __func__, __LINE__, n * sizeof(float *));
    for (int i = 0; i < n; i++)
    {
        if ((M.data[i] = (float *)malloc(n * sizeof(float))) == NULL)
            err_nomemory(__FILE__, __func__, __LINE__, n * sizeof(float));
    }
    return M;
}

static Matrix add(const Matrix *M1, const Matrix *M2)
{
    cnt_add++;
    assert(M1->size == M2->size);
    Matrix M3 = matrix_create(M1->size);
    for (int i = 0; i < M1->size; i++)
    {
        for (int j = 0; j < M1->size; j++)
            M3.data[i][j] = M1->data[i][j] + M2->data[i][j];
    }

    return M3;
}

static Matrix sub(const Matrix *M1, const Matrix *M2)
{
    cnt_sub++;
    assert(M1->size == M2->size);
    Matrix M3 = matrix_create(M1->size);
    for (int i = 0; i < M1->size; i++)
    {
        for (int j = 0; j < M1->size; j++)
            M3.data[i][j] = M1->data[i][j] - M2->data[i][j];
    }

    return M3;
}

static void matrix_print(const char *tag, const Matrix *M)
{
    printf("%s (%dx%d):\n", tag, M->size, M->size);
    if (M->size > 128)
    {
        printf("Printing suppressed - matrix too large\n");
        return;
    }
    char buffer[32];
    int len = snprintf(buffer, sizeof(buffer), "%d", M->size);
    for (int i = 0; i < M->size; i++)
    {
        printf("[%*d]: ", len, i);
        const char *pad = "";
        for (int j = 0; j < M->size; j++)
        {
            printf("%s%f", pad, M->data[i][j]);
            pad = "\t";
        }
        printf("\n");
    }
}

static Matrix strassen_multiply(const Matrix *A, const Matrix *B)
{
    assert(A->size == B->size);
    if (A->size == 1)
    {
        Matrix C = matrix_create(A->size);
        C.data[0][0] = A->data[0][0] * B->data[0][0];
        return C;
    }
    count++;
    Matrix C = matrix_create(A->size);
    int k = A->size / 2;
    /** Creating sub matrices**/
    Matrix A11 = matrix_create(k);
    Matrix A12 = matrix_create(k);
    Matrix A21 = matrix_create(k);
    Matrix A22 = matrix_create(k);
    Matrix B11 = matrix_create(k);
    Matrix B12 = matrix_create(k);
    Matrix B21 = matrix_create(k);
    Matrix B22 = matrix_create(k);

    /** Dividing the Data Matrices A & B **/
    for (int i = 0; i < k; i++)
    {
        for (int j = 0; j < k; j++)
        {
            A11.data[i][j] = A->data[i + 0][j + 0];
            A12.data[i][j] = A->data[i + 0][k + j];
            A21.data[i][j] = A->data[k + i][j + 0];
            A22.data[i][j] = A->data[k + i][k + j];
            B11.data[i][j] = B->data[i + 0][j + 0];
            B12.data[i][j] = B->data[i + 0][k + j];
            B21.data[i][j] = B->data[k + i][j + 0];
            B22.data[i][j] = B->data[k + i][k + j];
        }
    }

    Matrix T1 = sub(&B12, &B22);
    Matrix P1 = strassen_multiply(&A11, &T1);
    matrix_destroy(&T1);

    Matrix T2 = add(&A11, &A12);
    Matrix P2 = strassen_multiply(&T2, &B22);
    matrix_destroy(&T2);

    Matrix T3 = add(&A21, &A22);
    Matrix P3 = strassen_multiply(&T3, &B11);
    matrix_destroy(&T3);

    Matrix T4 = sub(&B21, &B11);
    Matrix P4 = strassen_multiply(&A22, &T4);
    matrix_destroy(&T4);

    Matrix T5A = add(&A11, &A22);
    Matrix T5B = add(&B11, &B22);
    Matrix P5 = strassen_multiply(&T5A, &T5B);
    matrix_destroy(&T5A);
    matrix_destroy(&T5B);

    Matrix T6A = sub(&A12, &A22);
    Matrix T6B = add(&B21, &B22);
    Matrix P6 = strassen_multiply(&T6A, &T6B);
    matrix_destroy(&T6A);
    matrix_destroy(&T6B);

    Matrix T7A = sub(&A11, &A21);
    Matrix T7B = add(&B11, &B12);
    Matrix P7 = strassen_multiply(&T7A, &T7B);
    matrix_destroy(&T7A);
    matrix_destroy(&T7B);

    matrix_destroy(&A11);
    matrix_destroy(&A12);
    matrix_destroy(&A21);
    matrix_destroy(&A22);
    matrix_destroy(&B11);
    matrix_destroy(&B12);
    matrix_destroy(&B21);
    matrix_destroy(&B22);

    Matrix C1A = add(&P5, &P4);
    Matrix C1B = add(&C1A, &P6);
    Matrix C11 = sub(&C1B, &P2);
    Matrix C12 = add(&P1, &P2);
    Matrix C21 = add(&P3, &P4);
    Matrix C2A = add(&P5, &P1);
    Matrix C2B = sub(&C2A, &P3);
    Matrix C22 = sub(&C2B, &P7);
    matrix_destroy(&C1A);
    matrix_destroy(&C1B);
    matrix_destroy(&C2A);
    matrix_destroy(&C2B);

    matrix_destroy(&P1);
    matrix_destroy(&P2);
    matrix_destroy(&P3);
    matrix_destroy(&P4);
    matrix_destroy(&P5);
    matrix_destroy(&P6);
    matrix_destroy(&P7);

    for (int i = 0; i < k; i++)
    {
        for (int j = 0; j < k; j++)
        {
            C.data[i + 0][j + 0] = C11.data[i][j];
            C.data[i + 0][j + k] = C12.data[i][j];
            C.data[k + i][j + 0] = C21.data[i][j];
            C.data[k + i][k + j] = C22.data[i][j];
        }
    }

    matrix_destroy(&C11);
    matrix_destroy(&C12);
    matrix_destroy(&C21);
    matrix_destroy(&C22);

    return C;
}

int main(void)
{
    struct timeval begin, end;
    Matrix A = matrix_create(N);
    Matrix B = matrix_create(N);

    for (int i = 0; i <  N; i++)
    {
        for (int j = 0; j < N; j++)
        {
            A.data[i][j] = -1.0 + 2.0 * ((double)rand()) / RAND_MAX;
            B.data[i][j] = -1.0 + 2.0 * ((double)rand()) / RAND_MAX;
        }
    }

    gettimeofday(&begin, 0);
    Matrix C = strassen_multiply(&A, &B);
    gettimeofday(&end, 0);

    matrix_print("A", &A);
    matrix_print("B", &B);
    matrix_print("C", &C);

    matrix_destroy(&A);
    matrix_destroy(&B);
    matrix_destroy(&C);

    long seconds = end.tv_sec - begin.tv_sec;
    long microseconds = end.tv_usec - begin.tv_usec;
    double elapsed = seconds + microseconds * 1e-6;

    printf("Number of non-minimal recursive calls: %d\n", count);
    printf("Number of matrices created:    %zu\n", cnt_create);
    printf("Number of matrices destroyed:  %zu\n", cnt_destroy);
    printf("Number of matrix additions:    %zu\n", cnt_add);
    printf("Number of matrix subtractions: %zu\n", cnt_sub);
    printf("Total wall time: %f\n", elapsed);
    return 0;
}

这通过调用一个简单退出的函数来检测内存分配错误，而不是释放任何成功分配的内存并返回给调用者。

代码可以用2 的-DN=256任何其他幂进行编译。目前尚不清楚如果大小不是 2 的幂会发生什么。

表现

各种大小的一些采样时间和其他统计信息：

N=8

Number of non-minimal recursive calls: 57
Number of matrices created:    1884
Number of matrices destroyed:  1884
Number of matrix additions:    627
Number of matrix subtractions: 399
Total wall time: 0.000480

N=16

Number of non-minimal recursive calls: 400
Number of matrices created:    13203
Number of matrices destroyed:  13203
Number of matrix additions:    4400
Number of matrix subtractions: 2800
Total wall time: 0.003723

N=32

Number of non-minimal recursive calls: 2801
Number of matrices created:    92436
Number of matrices destroyed:  92436
Number of matrix additions:    30811
Number of matrix subtractions: 19607
Total wall time: 0.025097

N=64

Number of non-minimal recursive calls: 19608
Number of matrices created:    647067
Number of matrices destroyed:  647067
Number of matrix additions:    215688
Number of matrix subtractions: 137256
Total wall time: 0.161971

N=128

Number of non-minimal recursive calls: 137257
Number of matrices created:    4529484
Number of matrices destroyed:  4529484
Number of matrix additions:    1509827
Number of matrix subtractions: 960799
Total wall time: 1.164555

N=256

Number of non-minimal recursive calls: 960800
Number of matrices created:    31706403
Number of matrices destroyed:  31706403
Number of matrix additions:    10568800
Number of matrix subtractions: 6725600
Total wall time: 7.632881

N=512

Number of non-minimal recursive calls: 6725601
Number of matrices created:    221944836
Number of matrices destroyed:  221944836
Number of matrix additions:    73981611
Number of matrix subtractions: 47079207
Total wall time: 53.730002

N=1024

Number of non-minimal recursive calls: 47079208
Number of matrices created:    1553613867
Number of matrices destroyed:  1553613867
Number of matrix additions:    517871288
Number of matrix subtractions: 329554456
Total wall time: 373.596480

N=2048

Number of non-minimal recursive calls: 329554457
Number of matrices created:    10875297084
Number of matrices destroyed:  10875297084
Number of matrix additions:    3625099027
Number of matrix subtractions: 2306881199
Total wall time: 2737.750096

请注意，创建的矩阵数量与销毁的数量相同；这让人放心。还要注意，有大量的矩阵被创建和销毁。

然而，将要相乘的矩阵的大小加倍并不是三次时间；它比 O(N³) 好，而朴素算法是 O(N³)。

提高性能

提高代码速度的一种方法是特殊情况下的 2x2 矩阵乘法。实施后，结果如下：

N=16

Number of large multiplications: 57
Number of 1x1 multiplications:   0
Number of 2x2 multiplications:   343
Number of matrices created:      1884
Number of matrices destroyed:    1884
Number of matrix additions:      627
Number of matrix subtractions:   399
Total wall time: 0.001045

N=32

Number of large multiplications: 400
Number of 1x1 multiplications:   0
Number of 2x2 multiplications:   2401
Number of matrices created:      13203
Number of matrices destroyed:    13203
Number of matrix additions:      4400
Number of matrix subtractions:   2800
Total wall time: 0.006532

N=64

Number of large multiplications: 2801
Number of 1x1 multiplications:   0
Number of 2x2 multiplications:   16807
Number of matrices created:      92436
Number of matrices destroyed:    92436
Number of matrix additions:      30811
Number of matrix subtractions:   19607
Total wall time: 0.038640

N=128

Number of large multiplications: 19608
Number of 1x1 multiplications:   0
Number of 2x2 multiplications:   117649
Number of matrices created:      647067
Number of matrices destroyed:    647067
Number of matrix additions:      215688
Number of matrix subtractions:   137256
Total wall time: 0.263008

N=256

Number of large multiplications: 137257
Number of 1x1 multiplications:   0
Number of 2x2 multiplications:   823543
Number of matrices created:      4529484
Number of matrices destroyed:    4529484
Number of matrix additions:      1509827
Number of matrix subtractions:   960799
Total wall time: 1.796228

N=512

Number of large multiplications: 960800
Number of 1x1 multiplications:   0
Number of 2x2 multiplications:   5764801
Number of matrices created:      31706403
Number of matrices destroyed:    31706403
Number of matrix additions:      10568800
Number of matrix subtractions:   6725600
Total wall time: 12.383302

作为比较，使用 1x1 和 2x2 特殊情况创建和销毁的矩阵数量为：

N              1x1           2x2
 16         13,203         1,884
 32         92,436        13,203
 64        647,067        92,436
128      4,528,484       647,067
256     31,706,403     4,529,484
512    221,944,836    31,706,403

观察为乘以 NxN 矩阵的 1x1 最小情况创建的矩阵数量与使用 2Nx2N 矩阵的 2x2 最小情况相同。它还提供了相当显着的加速（比较 N=512 和 1x1 为 53.73 秒，而 N=512 和 2x2 为 12.38 秒）。很多原始成本是创建 1x1 矩阵以相乘。

其他建议

诽谤者莫妮卡建议：

只有在大量使用子矩阵时才应复制子矩阵——以提高缓存局部性。否则，“子矩阵”就不是一个变量，而是一个概念。这意味着当您执行任何矩阵运算时，您应该传递一些描述索引范围和索引步幅的对象，例如矩阵乘法。这样你就不会创建那些子矩阵。一般来说，要使这段代码合理，还有很多工作要做。

这将使矩阵结构更加复杂，但也会从根本上提高性能。您最终可能会matrix_create()返回 a Matrix *，并且该结构将包含额外的元素：int tl_r; int tl_c; int br_r; int br_c;（左上的行和列，左下的行和列）。您将有另一个函数将矩阵拆分为 4 个四分之一矩阵，这些矩阵都将引用未拆分矩阵的数据，但子矩阵的左上角和右下角坐标具有不同的值。如果您继续使用指向浮点组织数组的当前指针数组，则不需要记录“步幅”（原始数组中每行的宽度，这也是高度，因为这仅处理方阵） . 您必须小心内存管理。结果数组将被重新创建。您不会发布季度矩阵中的数据——只会发布那些重新创建的数据。

带翅膀的小行星评论：

为什么将指针数组用于方数组？这是无缘无故的大量开销。只需创建一个 N*N 浮点数组！然后你就可以开始简化所有这些疯狂的内存管理了。

尽管需要谨慎，但其中有一些正义。我仍然认为您将使用结构，但数据元素将float *代替float **，并且您将计算数组索引 ( row * width + col) 而不是使用两个下标。如果您放弃结构，则可以改用“可变长度数组”（VLA）表示法。需要小心。数组仍将是动态分配的。

进一步的实验和建议

我已经尝试过 4x4 和 8x8 的特殊情况，由于减少了内存管理开销（更少的矩阵分配和破坏），两者都提供了相当大的好处。将具有不同最小尺寸的 1024x1024 矩阵相乘得到：

Size    Time        #Matrices
1x1   6m 32s    1,553,613,867
2x2   1m 31s      221,944,836
4x4      23s       31,706,403
8x8       7s        4,529,484

我还编写了一个执行直接原始矩阵乘法（O(N³) 算法 - 使用与我用于 NxN 的 8x8 乘法相同的代码）的版本，它比 Strassen 算法快很多，主要是因为有几乎不需要内存管理。

Size             Time
128x128           3ms
256x256          25ms
512x512         280ms
1024x1024     1,802ms
2048x2048    84,686ms
4096x4096   850,860ms

请注意，1024 处的 1.80s 和 2048 处的 84.7s 之间的时间乘积大于三次关系（或多或少适用于其他情况的 8 倍）——我没有调查原因。

我认为从这里加快速度的关键不是复制矩阵——使用 Unslander Monica 的建议。我注意到您可能不需要我之前建议的 4 个坐标；2个坐标和大小就足够了（因为矩阵都是正方形的）。这减少了复制和内存分配——这将对性能产生重大影响。

我不认为施特拉森被证明是失败的。我认为这表明您（我们）正在使用的幼稚内存管理不是要走的路。但我也认为，除非您的矩阵大于 1024x1024，否则天真的乘法算法可能足够快。如果您必须处理大数组，那么 Strassen 算法可能是有益的。需要更多的编码和测试！

c - 在 C 中超过 1024x1024 矩阵时，总线错误/核心转储

问题描述

解决方案

诊断

处方

表现

N=8

N=16

N=32

N=64

N=128

N=256

N=512

N=1024

N=2048

提高性能

N=16

N=32

N=64

N=128

N=256

N=512

其他建议

进一步的实验和建议

推荐阅读