Конвертирайте matlab в c++, bsxfun

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- Обичайно е да съхранявате 2D масиви в C++ като 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.

Не е лошо подобрение (почти x4 по-бързо) само чрез добавяне на един ред към кода (макросът pargma omp).

Това последното бие 37 ms, които получавам в MATLAB (R2013b). Кодът е компилиран с помощта на GCC 4.8.1 (MinGW-w64, работещ на Windows 8, лаптоп Core i7).

Ако наистина искате да надхвърлите ограниченията тук за C++ кода, ще трябва да добавите векторизация (SSE/AVX intrinsics) в допълнение към многонишковостта, постигната с OpenMP.

Може също да обмислите използването на програмиране на 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. Проблемът тук обаче до голяма степен е обвързан с паметта, така че вероятно няма да видите подобрение (вероятно дори по-бавно поради претоварването на прехвърлянето на данни).