我有一个标准的矩阵乘法算法.
// 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(×tamp, NULL);
start = timestamp.tv_sec * 1000 + timestamp.tv_usec / 1000;
mmProduct(&first, &second, &result);
gettimeofday(×tamp, 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.
// ...
// ...
}
以下是您自己运行它的完整代码: https://gist.github.com/kevinMEH/c37bda728bf5ee53c32a405ef611e5a7个
当我用两个210x 210矩阵运行该算法时，乘法大约需要2秒.
当我使用两个211x211矩阵运行该算法时，我预计该算法会在8*2秒~16秒内将它们相乘.然而，它需要36秒才能进行倍增.
当我使用两个212x212矩阵运行该算法时，我预计该算法将它们相乘所需的时间为64*2000毫秒~128秒.最坏的情况是，我预计8*36=288秒会成倍增加.然而，它需要391秒才能倍增.
(使用Apple Clang 14.0.0和-O1编译，但使用-Ofast并且没有优化，它的伸缩性类似，仍然不是8的倍数，但更高.运行在Apple M1芯片上.)
我很困惑为什么会这样.执行的加法和乘法运算的次数正好是(2n)3，当我将矩阵的宽度增加2倍时，我预计后续的乘法运算每次只需要2=8倍的时间.
嫌疑人:
addressAt和rowAt帮助器函数是导致速度减慢的原因:这没有意义；加法和乘法的数量按n3zoom ，但调用addressAt的时候仅按n2zoom ，而rowAt的调用仅按n2zoom ，因此我们实际上应该看到速度的相对提高.(8*(n3+n2+n)>；8n3+4n2+2n)
我还手动内联了函数，发现没有明显的速度提高，所以应该不是因为函数调用开销.
非正态数:我的随机数生成器得到一个随机的U64，并将其除以2*0x8..0u，得到一个介于0和1之间的数字，因此应该不存在非正态数.尽管如此，在随机填充函数中，我将所有内容乘以100000，结果仍然是相同的.
我怀疑这就是罪魁祸首，但我找不到原因.
我最好的猜测是，速度变慢是因为CPU预取了错误的列元素.当我们按列数递增j时，CPU对此感到困惑，不会预取matrix->elements[j + n * columns]处的数组元素，因为它不希望j按列递增(所以它预取j、j + 1或其他任何值，但不是j + columns)，但这对我来说没有意义，因为columns是常量，所以CPU应该知道得更清楚.
Side note:
我还编写了实现Strassen matrix multiplication algorithm的mmStrassen()，当我测试它时，它每次都出乎意料地长7倍.它是在常规的mmProduct函数之后一起测试的，所以它不是运行时间问题.
另外，欢迎任何优化建议，谢谢.
对于这个长长的问题，我要提前表示感谢和歉意.