标准矩阵乘法算法扩容（2^10+阶）时性能超预期下降

实现代码 #

// matrix.h

typedef struct Matrix {
    int rows;
    int columns;
    double *elements;
} Matrix;

double *addressAt(Matrix *matrix, int row, int column);
double *rowAt(Matrix *matrix, int row);
double dotProductColumnArray(Matrix *matrix, const int column, double *array);
void mmProduct(Matrix *first, Matrix *second, Matrix *result);

// matrix.c

// Helper function to find address of elements at row and column.
inline double *addressAt(Matrix *matrix, int row, int column) {
    return &(matrix->elements[row * matrix->columns + column]);
}

// Helper function to find row array
inline double *rowAt(Matrix *matrix, int row) {
    return &(matrix->elements[row * matrix->columns]);
}

// Finds the dot product of a column of a matrix and an array.
double dotProductColumnArray(Matrix *matrix, const int column, double *array) {
    double sum = 0.0;
    const int rows = matrix->rows, columns = matrix->columns;
    int j = column;
    const double *elements = matrix->elements;
    for (int i = 0; i < rows; i++) {
        sum += array[i] * elements[j];
        j += columns;
    }
    return sum;
}

void mmProduct(Matrix *first, Matrix *second, Matrix *result) {
    const int rows = result->rows;
    const int columns = result->columns;
    for (int i = 0; i < rows; i++) {
        double *row = rowAt(first, i);
        for (int j = 0; j < columns; j++) {
            *addressAt(result, i, j) = dotProductColumnArray(second, j, row);
        }
    }
}

// Snippet of main.c

// Fills the matrix with random numbers between -1 and 1
void randomFill(Matrix *matrix);

int main() {
    struct timeval timestamp;
    long start, now;
    Matrix first, second, result;
    
    // 2^10
    first = createMatrix((int)pow(2, 10), (int)pow(2, 10));
    second = createMatrix((int)pow(2, 10), (int)pow(2, 10));
    randomFill(&first);
    randomFill(&second);
    result = createMatrix((int)pow(2, 10), (int)pow(2, 10));

    gettimeofday(&timestamp, NULL);
    start = timestamp.tv_sec * 1000 + timestamp.tv_usec / 1000;
    mmProduct(&first, &second, &result);
    gettimeofday(&timestamp, NULL);
    now = timestamp.tv_sec * 1000 + timestamp.tv_usec / 1000;
    printf("Time taken: %ldms\n", now - start);
    
    deleteMatrix(&first);
    deleteMatrix(&second);
    deleteMatrix(&result);
    
    // Same code as above but for 2^11, 2^12, replacing pow() when necessary.
    // ...
    // ...
}

测试结果 #

• 2^10×210矩阵乘法耗时约2秒；
• 2^11×211矩阵预期耗时约16秒，实际耗时36秒；
• 2^12×212矩阵预期耗时约128/288秒，实际耗时391秒。

运行环境 #

• 编译环境：Apple Clang 14.0.0（测试了-O1/-Ofast/无优化，均存在类似缩放问题）
• 硬件：Apple M1芯片

问题背景 #

矩阵乘法运算量为O(n³)，理论上矩阵规模扩容2倍应耗时8倍，但实际性能下降幅度远超预期。同时实现的Strassen算法（mmStrassen）符合预期的7倍缩放规律。

已排查方向 #

1. addressAt/rowAt辅助函数：手动内联后无性能提升，且调用量缩放远低于运算量；
2. 非正规数（subnormal numbers）：调整随机填充值后问题依旧；
3. 怀疑dotProductColumnArray函数存在CPU预取失效问题，但不确定原因。

寻求问题原因分析及优化建议。

内容的提问来源于stack exchange，提问作者kmeh