您的代码和/或您的理解存在多个严重问题。让我试着解释一下。
矩阵乘法受到处理器加载和存储值到内存的速率的限制。大多数当前架构使用 cache 来帮助解决这个问题。数据以块的形式从内存移动到缓存,从缓存移动到内存。为了最大限度地利用缓存,您需要确保使用该块中的所有数据。为此,请确保您按顺序访问内存中的数据。
在 C 中,多维数组在 row-major order 中指定。表示最右边的索引在内存中是连续的;即a[i][k] 和a[i][k+1] 在内存中是连续的。
根据架构,处理器等待数据从 RAM 移动到缓存(反之亦然)所花费的时间可能包含在 CPU 时间中,也可能不包含在 CPU 时间中(即例如clock() 测量,尽管分辨率很差)。对于这种测量(“microbenchmark”),最好同时测量和报告 CPU 和实际(或挂钟)使用的时间;特别是如果微基准测试在不同的机器上运行,以便更好地了解更改的实际影响。
会有很多变化,因此通常情况下,您会测量数百次重复所花费的时间(每次重复可能会进行多次操作;足以轻松测量),存储每次的持续时间,并报告它们的中位数。为什么是中位数,而不是最小值、最大值、平均值?因为总会偶尔出现故障(由于外部事件或其他原因导致的不合理测量),通常会产生比正常情况高得多的值;这使得最大值变得无趣,并且除非被移除,否则会扭曲平均值(平均值)。最小值通常是一种过度乐观的情况,一切都恰巧完美;这在实践中很少发生,所以只是一种好奇心,没有实际意义。另一方面,中值时间为您提供了一个实际的测量结果:您可以预期 50% 的测试用例运行所花费的时间不会超过测量的中值时间。
在 POSIXy 系统(Linux、Mac、BSD)上,您应该使用clock_gettime() 来测量时间。 struct timespec 格式具有纳秒精度(1 秒 = 1,000,000,000 纳秒),但分辨率可能更小(即,时钟变化超过 1 纳秒,无论何时发生变化)。我个人使用
#define _POSIX_C_SOURCE 200809L
#include <time.h>
static struct timespec cpu_start, wall_start;
double cpu_seconds, wall_seconds;
void timing_start(void)
{
clock_gettime(CLOCK_REALTIME, &wall_start);
clock_gettime(CLOCK_THREAD_CPUTIME_ID, &cpu_start);
}
void timing_stop(void)
{
struct timespec cpu_end, wall_end;
clock_gettime(CLOCK_REALTIME, &wall_end);
clock_gettime(CLOCK_THREAD_CPUTIME_ID, &cpu_end);
wall_seconds = (double)(wall_end.tv_sec - wall_start.tv_sec)
+ (double)(wall_end.tv_nsec - wall_start.tv_nsec) / 1000000000.0;
cpu_seconds = (double)(cpu_end.tv_sec - cpu_start.tv_sec)
+ (double)(cpu_end.tv_nsec - cpu_start.tv_nsec) / 1000000000.0;
}
你在操作之前调用timing_start(),在操作之后调用timing_stop();然后,cpu_seconds 包含占用的 CPU 时间量和 wall_seconds 占用的实际挂钟时间(均以秒为单位,例如使用 %.9f 打印所有有意义的小数)。
以上内容不适用于 Windows,因为 Microsoft 不希望您的 C 代码可移植到其他系统。它更喜欢开发自己的“标准”。 (与例如 POSIX getline() 或除 Windows 之外的所有系统上的宽字符支持状态相比,那些 C11 “安全”_s() I/O 函数变体是愚蠢的骗局。)
矩阵乘法是
c[r][c] = a[r][0] * b[0][c]
+ a[r][1] * b[1][c]
: :
+ a[r][L] * b[L][c]
其中a 有L+1 列,b 有L+1 行。
为了使求和循环使用连续元素,我们需要转置b。如果B[c][r] = b[r][c],那么
c[r][c] = a[r][0] * B[c][0]
+ a[r][1] * B[c][1]
: :
+ a[r][L] * B[c][L]
请注意,a 和 B 在内存中是连续的,但分开(可能彼此“远离”)就足够了,这样处理器在这种情况下可以有效地利用缓存。
OP 使用类似于以下伪代码的简单循环来转置b:
For r in rows:
For c in columns:
temporary = b[r][c]
b[r][c] = b[c][r]
b[c][r] = temporary
End For
End For
上面的问题是每个元素都参与了两次交换。例如,如果 b 有 10 行和 10 列,r = 3, c = 5 交换 b[3][5] 和 b[5][3],但后来,r = 5, c = 3 再次交换 b[5][3] 和 b[3][5]!本质上,双循环最终将矩阵恢复到原始顺序;它不做转置。
考虑以下条目和实际转置:
b[0][0] b[0][1] b[0][2] b[0][0] b[1][0] b[2][0]
b[1][0] b[1][1] b[1][2] ⇔ b[0][1] b[1][1] b[2][1]
b[2][0] b[2][1] b[2][2] b[0][2] b[1][2] b[2][2]
不交换对角线条目。您只需在上三角部分(c > r)或下三角部分(r > c)进行交换,即可交换所有条目,因为每次交换都会将一个条目从上三角交换到下三角,反之亦然。
所以,回顾一下:
是不是做错了什么?
是的。你的转置什么都不做。您还没有理解为什么要转置第二个矩阵的原因。您的时间测量依赖于低精度 CPU 时间,这可能无法反映在 RAM 和 CPU 缓存之间移动数据所花费的时间。在第二个测试用例中,使用m2“转置”(除了它不是,因为您将每个元素对交换两次,将它们返回到原来的状态),您的最内层循环位于最左边的数组索引之上,这意味着它计算错误的结果。 (此外,由于最内层循环的连续迭代会访问内存中彼此相距较远的项目,因此它是反优化:它使用在速度方面最差的模式.)
以上所有内容可能听起来很刺耳,但实际上并非如此,。我不认识你,也不想评价你;我只是在您当前的理解中指出此特定答案中的错误,并且仅希望它可以帮助您以及在类似情况下遇到此问题的任何其他人学习。