1

我正在尝试将我的 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？

在网上搜索后，我发现了两种可能的方法：

  1. 包括来自犰狳的图书馆
  2. 包括来自 Octave 的库

但关键是我不知道该怎么做，我的意思是我不知道实现的细节。

概括：

  1. 我想问是否有人知道如何更改 C++ 代码，以便我可以有类似或相同的性能bsxfun？
  2. 有人可以提供一些提示或步骤或示例，以便我可以学习如何包含 Armadillo 或 Octave 来完成此任务。

编辑：

感谢@Peter，我使用选项进行编译，-O3然后问题“解决”了，我的意思是速度与 MATLAB 相同。

4

2 回答 2

6

1-您以错误的顺序运行循环。在 C 和 C++ 中，二维数组以行为主，即在内存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 点）、动态内存分配更少以及缓存利用率更高。

于 2014-01-17T14:33:55.037 回答
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++ 代码版本。我正在使用简单的类来创建向量/矩阵接口。请注意，基础数据基本上存储为具有类似于 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 工作区。然而，这里的问题主要是受内存限制的，因此不太可能看到任何改进（由于数据传输的开销，可能会更慢）。

于 2014-01-18T04:20:37.720 回答