正如我在评论中提到的,您的 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++ 代码版本。我正在使用简单的类来创建向量/矩阵接口。请注意,基础数据基本上存储为具有类似于 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
只需在代码(宏)中添加一行,就可以进行不错的改进(几乎快 4 倍)。
最后一个超过了我在 MATLAB (R2013b) 中获得的 37 毫秒。该代码是使用 GCC 4.8.1 编译的(MinGW-w64 在 Windows 8、Core i7 笔记本电脑上运行)。
如果您真的想在此处突破 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 工作区。然而,这里的问题主要是受内存限制的,因此不太可能看到任何改进(由于数据传输的开销,可能会更慢)。