我正在尝试将我的MATLAB代码转换为C++，我发现在以下情况下存在问题：
MATLAB
A = rand(1000,40000);
b = rand(1000,1);
tic;
ans = bsxfun(@ne,b,A);
toc
c++
std::vector<std::vector<int> > A;
std::vector<int> b;
std::vector<int> ans(10000);
// initial A and b
const clock_t begin_time = clock();
for(int i = 0; i < 40000; ++i){
for(int j = 0; j < 1000; ++j){
if(A[i][j] != b[j])
ans[i]++;
}
}
double run_time = static_cast<double>((clock() - begin_time)) / CLOCKS_PER_SEC;
我发现C++的速度是MATLAB的三倍。我想问一下，是否有人知道如何更改C++代码，以便我可以获得与bsxfun类似或相同的性能？
在我搜索网络后，我找到了两种可能的方法：
- 包含来自Armadillo的库
- 包括Octave中的库
但关键是我不确定如何做到这一点，我的意思是我不知道实现的细节。
摘要：
- 我想问是否有人知道如何更改C++代码，以便我可以获得与
bsxfun类似或相同的性能 - 有人能提供一些提示、步骤或例子吗？这样我就可以学习如何包括Armadillo或Octave来完成这项任务
编辑：
多亏了@Peter，我用选项-O3编译，然后问题就"解决"了，我的意思是速度和MATLAB一样。
1-您的循环运行顺序错误。在C和C++中，2D阵列以行为主存储，这意味着A[j][i]和A[j][i+1]在存储器中相邻。(可以这样想：A[j]是第一个下标运算，返回对另一个向量的引用，然后用[i]再次为其下标)。
将数据保存在缓存中以进行尽可能多的操作是现代处理器性能的关键之一，这意味着您希望在可能的时候访问相邻的元素。所以切换循环的顺序：
for(int j = 0; j < 1000; ++j){
for(int i = 0; i < 40000; ++i){
2-编译器选项非常重要。确保你是在"发布"模式下构建的，或者在上进行优化
3-将C++中的2D数组存储为1D数组是很常见的，使用乘法对行/列进行索引。也就是说，A将是大小为1000*40000的向量，而A[j][i]将改为A[j*row_length + i]。这样做的好处是有更多的连续内存(参见第1点)、更少的动态内存分配和更好的缓存利用率。
正如我在评论中提到的，您的MATLAB代码缺少对sum函数的调用(否则，这两个代码计算的内容不同！)。所以它应该是：
MATLAB
A = rand(1000,40000);
B = rand(1000,1);
tic
count = sum(bsxfun(@ne, A, B));
toc
在我的机器上我得到：
Elapsed time is 0.036931 seconds.
请记住，上面的语句是矢量化的(想想SIMD并行化)。如果大小足够大，MATLAB也可以自动运行这个多线程。
这是我用C++编写的代码版本。我使用简单的类来创建向量/矩阵接口。请注意，底层数据基本上存储为1D阵列，其列主序类似于MATLAB。
C++
#include <iostream>
#include <cstdlib> // rand
#include <ctime> // time
#include <sys/time.h> // gettimeofday
class Timer
{
private:
timeval t1, t2;
public:
Timer() {}
~Timer() {}
void start() { gettimeofday(&t1, NULL); }
void stop() { gettimeofday(&t2, NULL); }
double elapsedTime() { return (t2.tv_sec - t1.tv_sec)*1000.0 + (t2.tv_usec - t1.tv_usec)/1000; }
};
template<typename T>
class Vector
{
private:
T *data;
const size_t num;
public:
Vector(const size_t num) : num(num) { data = new T[num]; }
~Vector() { delete[] data; }
inline T& operator() (const size_t i) { return data[i]; }
inline const T& operator() (const size_t i) const { return data[i]; }
size_t size() const { return num; }
};
template<typename T>
class Matrix
{
private:
T *data;
const size_t nrows, ncols;
public:
Matrix(const size_t nr, const size_t nc) : nrows(nr), ncols(nc) { data = new T[nrows * ncols]; }
~Matrix() { delete[] data; }
inline T& operator() (const size_t r, const size_t c) { return data[c*nrows + r]; }
inline const T& operator() (const size_t r, const size_t c) const { return data[c*nrows + r]; }
size_t size1() const { return nrows; }
size_t size2() const { return ncols; }
};
inline double rand_double(double min=0.0, double max=1.0)
{
return (max - min) * (static_cast<double>(rand()) / RAND_MAX) + min;
}
int main() {
// seed random number generator
srand( static_cast<unsigned int>(time(NULL)) );
// intialize data
const int m = 1000, n = 40000;
Matrix<double> A(m,n);
Vector<double> B(m);
for(size_t i=0; i<A.size1(); i++) {
B(i) = rand_double();
for(size_t j=0; j<A.size2(); j++) {
A(i,j) = rand_double();
}
}
// measure timing
Timer timer;
timer.start();
// in MATLAB: count = sum(bsxfun(@ne, A, B))
Vector<double> count(n);
#pragma omp parallel for
for(int j=0; j<n; ++j) {
count(j) = 0.0;
for(int i=0; i<m; i++) {
count(j) += (A(i,j) != B(i));
}
}
timer.stop();
// elapsed time in milliseconds
std::cout << "Elapsed time is " << timer.elapsedTime() << " milliseconds." << std::endl;
return 0;
}
结果：
$ g++ -Wall -O3 test.cpp -o test
$ ./test
Elapsed time is 63 milliseconds.
如果我在启用OpenMP支持的情况下编译和运行它，我会得到：
$ g++ -Wall -O3 -fopenmp test.cpp -o test_omp
$ ./test_omp
Elapsed time is 16 milliseconds.
只需在代码中添加一行(pargma omp宏)，这是一个不错的改进(几乎快了x4)。
最后一个超过了我在MATLAB(R2013b)中获得的37毫秒。该代码是使用GCC 4.8.1(运行在Windows 8 Core i7笔记本电脑上的MinGW-w64)编译的。
如果你真的想在这里突破C++代码的限制，除了使用OpenMP实现的多线程外，你还必须添加矢量化(SSE/AVX内部)。
您可能还想考虑使用GPGPU编程(CUDA、OpenCL)。在MATLAB中，这很容易做到：
AA = gpuArray(A);
BB = gpuArray(B);
CC = sum(bsxfun(@ne, AA, BB));
C = gather(CC);
gpuArray(.)将矩阵传输到GPU，之后对其进行的所有操作都在GPU设备上执行，而不是在CPU上执行。gather(.)将把数组传输回MATLAB工作空间。然而，这里的问题很大程度上是内存限制的，因此不太可能有任何改进(由于数据传输的开销，可能会更慢)。