取两个3x3矩阵的乘积A*B=C。天真的,使用标准算法需要27次乘法。如果一个人很聪明,您可以只使用23个乘法来完成此操作,这是1973年Laderman发现的结果。该技术涉及保存中间步骤,并以正确的方式组合它们。
A*B=C
现在,让我们修复一种语言和一种类型,例如使用元素的C ++ double。如果Laderman算法是硬编码的,而不是简单的双循环的,那么我们是否可以期望现代编译器的性能超越算法的差异?
double
关于此问题的注释: 这是一个 编程 站点,在时间紧迫的内部循环的最佳实践中提出该问题;过早的优化不是这样。有关实施的提示非常受好评。
我自己进行了时序测试,结果令我惊讶(因此,为什么我首先问这个问题)。简而言之,在标准编译下,它的laderman速度快〜225%,但是在带有-03优化标志的情况下,它的 速度慢了50% !每次在-O3标志期间我都必须向矩阵中添加一个随机元素,否则编译器会 完全优化 简单乘法,并在时钟精度范围内将时间设为零。由于该laderman算法很难检查/仔细检查,因此我将在下面张贴完整的代码以供后代参考。
laderman
-03
-O3
规格:Ubuntu 12.04,Dell Prevision T1600,gcc。时间差异百分比:
g++ [2.22, 2.23, 2.27]
g++ -O3 [-0.48, -0.49, -0.48]
g++ -funroll-loops -O3 [-0.48, -0.48, -0.47]
基准测试代码以及Laderman实现:
#include <iostream> #include <ctime> #include <cstdio> #include <cstdlib> using namespace std; void simple_mul(const double a[3][3], const double b[3][3], double c[3][3]) { int i,j,m,n; for(i=0;i<3;i++) { for(j=0;j<3;j++) { c[i][j] = 0; for(m=0;m<3;m++) c[i][j] += a[i][m]*b[m][j]; } } } void laderman_mul(const double a[3][3], const double b[3][3], double c[3][3]) { double m[24]; // not off by one, just wanted to match the index from the paper m[1 ]= (a[0][0]+a[0][1]+a[0][2]-a[1][0]-a[1][1]-a[2][1]-a[2][2])*b[1][1]; m[2 ]= (a[0][0]-a[1][0])*(-b[0][1]+b[1][1]); m[3 ]= a[1][1]*(-b[0][0]+b[0][1]+b[1][0]-b[1][1]-b[1][2]-b[2][0]+b[2][2]); m[4 ]= (-a[0][0]+a[1][0]+a[1][1])*(b[0][0]-b[0][1]+b[1][1]); m[5 ]= (a[1][0]+a[1][1])*(-b[0][0]+b[0][1]); m[6 ]= a[0][0]*b[0][0]; m[7 ]= (-a[0][0]+a[2][0]+a[2][1])*(b[0][0]-b[0][2]+b[1][2]); m[8 ]= (-a[0][0]+a[2][0])*(b[0][2]-b[1][2]); m[9 ]= (a[2][0]+a[2][1])*(-b[0][0]+b[0][2]); m[10]= (a[0][0]+a[0][1]+a[0][2]-a[1][1]-a[1][2]-a[2][0]-a[2][1])*b[1][2]; m[11]= a[2][1]*(-b[0][0]+b[0][2]+b[1][0]-b[1][1]-b[1][2]-b[2][0]+b[2][1]); m[12]= (-a[0][2]+a[2][1]+a[2][2])*(b[1][1]+b[2][0]-b[2][1]); m[13]= (a[0][2]-a[2][2])*(b[1][1]-b[2][1]); m[14]= a[0][2]*b[2][0]; m[15]= (a[2][1]+a[2][2])*(-b[2][0]+b[2][1]); m[16]= (-a[0][2]+a[1][1]+a[1][2])*(b[1][2]+b[2][0]-b[2][2]); m[17]= (a[0][2]-a[1][2])*(b[1][2]-b[2][2]); m[18]= (a[1][1]+a[1][2])*(-b[2][0]+b[2][2]); m[19]= a[0][1]*b[1][0]; m[20]= a[1][2]*b[2][1]; m[21]= a[1][0]*b[0][2]; m[22]= a[2][0]*b[0][1]; m[23]= a[2][2]*b[2][2]; c[0][0] = m[6]+m[14]+m[19]; c[0][1] = m[1]+m[4]+m[5]+m[6]+m[12]+m[14]+m[15]; c[0][2] = m[6]+m[7]+m[9]+m[10]+m[14]+m[16]+m[18]; c[1][0] = m[2]+m[3]+m[4]+m[6]+m[14]+m[16]+m[17]; c[1][1] = m[2]+m[4]+m[5]+m[6]+m[20]; c[1][2] = m[14]+m[16]+m[17]+m[18]+m[21]; c[2][0] = m[6]+m[7]+m[8]+m[11]+m[12]+m[13]+m[14]; c[2][1] = m[12]+m[13]+m[14]+m[15]+m[22]; c[2][2] = m[6]+m[7]+m[8]+m[9]+m[23]; } int main() { int N = 1000000000; double A[3][3], C[3][3]; std::clock_t t0,t1; timespec tm0, tm1; A[0][0] = 3/5.; A[0][1] = 1/5.; A[0][2] = 2/5.; A[1][0] = 3/7.; A[1][1] = 1/7.; A[1][2] = 3/7.; A[2][0] = 1/3.; A[2][1] = 1/3.; A[2][2] = 1/3.; t0 = std::clock(); for(int i=0;i<N;i++) { // A[0][0] = double(rand())/RAND_MAX; // Keep this in for -O3 simple_mul(A,A,C); } t1 = std::clock(); double tdiff_simple = (t1-t0)/1000.; cout << C[0][0] << ' ' << C[0][1] << ' ' << C[0][2] << endl; cout << C[1][0] << ' ' << C[1][1] << ' ' << C[1][2] << endl; cout << C[2][0] << ' ' << C[2][1] << ' ' << C[2][2] << endl; cout << tdiff_simple << endl; cout << endl; t0 = std::clock(); for(int i=0;i<N;i++) { // A[0][0] = double(rand())/RAND_MAX; // Keep this in for -O3 laderman_mul(A,A,C); } t1 = std::clock(); double tdiff_laderman = (t1-t0)/1000.; cout << C[0][0] << ' ' << C[0][1] << ' ' << C[0][2] << endl; cout << C[1][0] << ' ' << C[1][1] << ' ' << C[1][2] << endl; cout << C[2][0] << ' ' << C[2][1] << ' ' << C[2][2] << endl; cout << tdiff_laderman << endl; cout << endl; double speedup = (tdiff_simple-tdiff_laderman)/tdiff_laderman; cout << "Approximate speedup: " << speedup << endl; return 0; }