高质量C++进阶[2]：如何让线性代数加速1000倍？

一直没搞明白传参、传值、传指针哪个快。同时想验证一下一些可行的加速，于是做了一个试验。

先说结论

线性代数相关的加速，有优先级地采用下面处理，可以使计算速度加快1000倍。

1. 调用Intel MKL 线性代数库加速 （据说针对Intel的处理器，可另外触发硬件加速）。
2. 尽可能采用Eigen原生函数计算。
3. 适当采用传引用的方法来控制输入输出。

问题构建

自定义一个struct，内含一个矩阵变量Eigen::Matrix。 定义一个作用于结构体的加法，使结构体成员矩阵的每一行加到第一列上，并返回该列向量，计时后cout该列向量的和以验证。

1. include :

#include <iostream>
#include <time.h>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;

2. struct定义

struct mm
    MatrixXd A;
    mm(int D);
mm::mm(int D)
    A.setZero(D, D);

3. 时长获取（包含<time.h>就可以用，返回的是毫秒）

size_t begin,end;
begin = clock();
// do something
end = clock();
cout << "time = " << end - begin << endl;

输入值：传值、传引用与传指针的对比

1. 输入传值【No.1】

VectorXd sum1(mm M)
	int n = M.A.rows();	
	VectorXd sum(n);
	sum.setZero();
	double s = 0.0;
	for (int ii = 0; ii < n; ii++)
		s = 0.0;
		for (int jj = 0; jj < n; jj++)
			s +=  M.A(ii, jj);
		sum(ii) = s;
	return sum;

2. 输入传引用【No.2】

VectorXd sum2(mm& M)
	int n = M.A.rows();	
	VectorXd sum(n);
	sum.setZero();
	double s = 0.0;
	for (int ii = 0; ii < n; ii++)
		s = 0.0;
		for (int jj = 0; jj < n; jj++)
			s += M.A(ii, jj);
		sum(ii) = s;
	return sum;

3. 传指针代码【N0.4】

VectorXd sum4(mm* M)
	int n = M->A.rows();
	VectorXd re(n);
	re.setZero();
	double s = 0.0;
	for (int ii = 0; ii < n; ii++)
		s = 0.0;
		for (int jj = 0; jj < n; jj++)
			s += M->A(ii, jj);
		re(ii) = s;
	return re;

4. main函数：

void main()
	int N = 3000;     // 矩阵大小：N×N        
	mm MA(N);         
	MA.A.setRandom();  // 矩阵填满随机值  
	size_t bg, ed;
	VectorXd s1, s2, s3, s4, s5, s6, s7;
	bg = clock();
	s1 = sum1(MA);
	ed = clock();
	cout << "sum1: " << s1.sum() << endl;
	cout << "time1: " << ed - bg << endl;
	cout << endl;
	bg = clock();
	s2 = sum2(MA);
	ed = clock();
	cout << "sum2: " << s2.sum() << endl;
	cout << "time2: " << ed - bg << endl;
	cout << endl;
	bg = clock();
	s4 = sum4(&MA);
	ed = clock();
	cout << "sum4: " << s4.sum() << endl;
	cout << "time4: " << ed - bg << endl;
	cout << endl;

5. 对比结果：

sum1: -207.832
time1: 4627
sum2: -207.832
time2: 4565
sum4: -207.832
time4: 4593

传值、传指针和传引用在我这个问题里影响不大。如果一定要比的话： 传引用最快 ，传值最慢，传指针居中。

输出值：return与传引用的对比

上面部分的输出都是通过return作用的，这里就用同样的问题对传引用试一试。

1. 输出值传引用+输入值传引用【No.3】

void sum3(VectorXd& re, mm& M)
	int n = M.A.rows();
	re.setZero(n);
	double s = 0.0;
	for (int ii = 0; ii < n; ii++)
		s = 0.0;
		for (int jj = 0; jj < n; jj++)
			s += M.A(ii, jj);
		re(ii) = s;

2. 输出值传引用+输入值传指针【No.5】

void sum5(VectorXd& re, mm* M)
	int n = M->A.rows();
	re.setZero(n);
	double s = 0.0;
	for (int ii = 0; ii < n; ii++)
		double s = 0.0;
		for (int jj = 0; jj < n; jj++)
			s += M->A(ii, jj);
		re(ii) = s;

3. main函数与上面是一样的，调用这两个函数，并评估时间

void main()
	int N = 3000;
	mm MA(N);
	MA.A.setRandom();
	size_t bg, ed;
	VectorXd s1, s2, s3, s4, s5, s6, s7;
	bg = clock();
	sum3(s3, MA);
	ed = clock();
	cout << "sum3: " << s3.sum() << endl;
	cout << "time3: " << ed - bg << endl;
	cout << endl;
	bg = clock();
	sum5(s5, &MA);
	ed = clock();
	cout << "sum5: " << s5.sum() << endl;
	cout << "time5: " << ed - bg << endl;
	cout << endl;

4. 对比结果

sum3: -207.832
time3: 4583
sum5: -207.832
time5: 4571

结论：将输出做成传引用，对我这个问题并没有带来什么改善。

提速大利器之一：Eigen原生函数

Eigen下定义的矩阵，自带很多函数，比如前面用到的Eigen::Matrix::sum(void)函数。按传值输入输出和传引用输入输出分别写一遍：

1. 利用Eigen::Matrix::sum(void)传引用输入输出【No.6】

void sum6(VectorXd& re, mm& M)
	int n = M.A.rows();
	re.setZero(n);
	for (int i = 0; i < n; i++)
		re(i) = M.A.col(i).sum();

2. 利用Eigen::Matrix::sum(void)传值输入输出【No.7】

VectorXd sum7(mm M)
	int n = M.A.rows();
	VectorXd re(n);
	for (int i = 0; i < n; i++)
		re(i) = M.A.col(i).sum();
	return re;

3. 与前面同样的问题，结果对比

sum6: -207.832
time6: 547
sum7: -207.832
time7: 590

结论是： 利用原生函数能极大地加快计算（快了10倍） 。但是传引用对我的问题虽然有改善，但幅度不大。

提速大利器之二：MKL加速

如果说Eigen原生函数的加速令人惊讶，那MKL的加速就要令人目瞪口呆了。

1. 使用方法，极其简单，只需要在#include<>代码之后的预处理区加上，其他什么也不用处理：

#ifndef EIGEN_USE_MKL_ALL
#define EIGEN_USE_MKL_ALL
#endif
#ifndef EIGEN_VECTORIZE_SSE4_2
#define EIGEN_VECTORIZE_SSE4_2
#endif

【1】直接#define EIGEN_USE_MKL_ALL就行，我加了一个宏判断是防止大型程序中被重复定义了。
【2】#define EIGEN_VECTORIZE_SSE4_2 起辅助作用，具体如何运行的目前还没搞清楚。

使用完全相同的代码跑出来的结果：

采用MKL加速和输入输出都传引用之后，一个简单的矩阵计算问题，运行速度能相差1000倍！！！太令人惊叹了，

试一试比起Matlab如何？

close all
clear all
N=3000;
M = rand(N,N);
V = zeros(N,1);
t1 = clock;
for i=1:N
    V(i) = sum(M(i,:));
t2 = clock;
sum8 = sum(V)
time8 = 1000*etime(t2,t1)   // 化为毫秒单位

运行结果：

sum8 = 4.4974e+06
time8 = 75.0000

可以看到，MKL优化后的普通代码（双for循环）就跟Matlab接近了，而开启全部优化之后，速度吊打Matlab! Cpp + Eigen + MKL 真香！！！！！

矩阵太小了，开到10000*10000对飚一下看看。
【1】Cpp + Eigen +MKL + 输入输出传引用

sum6: 10690.3
time6: 58
sum7: 10690.3
time7: 334

【2】Matlab

sum8 = 4.9995e+07
time8 = 1236.0

与MatLab对比的结论

1. 矩阵大了之后，传引用输入输出的对比非常明显（对比 time6 和 time7 ），再结合前面3000*3000的结果来看，能稳定地差个6倍左右。
2. 相同的程序设计思路下，Matlab的速度被秒杀，对比两次计算结果（ time8 与 time6、time7 ）能差到15-20倍。

写在后面

想记录下这个，是因为看到结果真的太Amazing了。另一方面还是我太垃圾了。要学的还很多。希望看到的有缘人多指教。

更新补充【1】

（1）如下代码比之前最快的【No.6】方法稍微快一些，但是相差不大，虽然【No.6】方法使用了for循环，但内部似乎将该for循环优化得很好【No.9】。

void sum9(VectorXd& sum, mm& M)
	sum.setZero(M.A.rows());
	sum = M.A.colwise().sum().transpose();

（2）下面代码比更新（1）中的代码慢很多，估计大概慢5倍。colwise更快是因为计算机的存储位置连续？【No.10】

void sum9(VectorXd& sum, mm& M)