Why MATLAB is faster than C++ in creating random numbers?

Question

I have been using MATLAB for a while for my projects and I have almost never had an experience in C++.

I needed speed and I heard that C++ can be more efficient and faster than MATLAB. So I tried this:

I created a matrix of random numbers using rand(5000,5000) on MATLAB.

And in C++, I have initialized a 2D vector created 2 for loops each of them looping for 5000 times and each time. MATLAB was 4-5x faster, so I thought it is because matlab executes vectorized codes in parallel, then I written the C++ code using parallel_for. Here is the code:

#include "stdafx.h"
#include <iostream>
#include <vector>
#include <fstream>
#include <ppl.h>
using namespace std;
using namespace concurrency;
int main();
{
    int a = 5000, b = 5000, j, k;
    vector< vector<int> > vec(a, vector<imt>(b));
    parallel_for(int(0), a, [&](int i) {
        for (j = 0; j <b; j++)
        {
            vec[i][j] = rand();
        }
    });
}

So the code above is about 25% faster than MATLAB's rand(5000,5000) Yet C++ is using 100% of the CPU while MATLAB is using 30% of CPU.

So I forced MATLAB to use all of the CPU by running 3 instances of MATLAB using rand(5000,5000) and divided the time it takes by 3. It made MATLAB twice as fast as C++.

I wonder what am I missing? I know this is a tiny example but I need an answer to be sure to port my code to C++.

Current status:

When I write C++ code without parallel_for I get half of the MATLAB's speed with the same CPU usage. Yet people who gave answers say that they are almost the same. I do not understand what I am missing

here is a snapshot of the optimization menu

Answer 1

This is maybe no answer, but a litle hint. The comparison might be a bit unfair due to the usage of vectors.

Here is a comparison I've written. Both take up roughly 100% of one of the four threads available. In both cases I create 5000x5000 random numbers and do this 100 times for timing

Matlab

function stackoverflow

tic
for i=1:100
    A =rand(5000);
end
toc

Runtime: ~27.9 sec

C++

#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <ctime>

using namespace std;


int main(){

    int N = 5000;
    double ** A = new double*[N];
    for (int i=0;i<N;i++)
        A[i] = new double[N];


    srand(time(NULL));

    clock_t start = clock();
    for (int k=0;k<100;k++){
        for (int i=0;i<N;i++){
            for (int j=0;j<N;j++){
                A[i][j] = rand();
            }
        }
    }

    cout << "T="<< (clock()-start)/(double)(CLOCKS_PER_SEC/1000)<< "ms " << endl;

}

Runtime: ~28.7 sec

So both examples run almost equally fast.

Answer 2

When you call rand(5000,5000) in Matlab, Matlab executes the command by calling Intel MKL library, which is a highly optimized library written in C/C++ with lots of hand-coded assembly.

MKL should be faster than any straightforward C++ implementation, but there is an overhead for Matlab to call external library. The net result is that, for random number generation in smaller sizes (less than 1K for instance), plain C/C++ implementation will be faster, but for larger sizes, Matlab will benefit from super optimized MKL.

Answer 3

After looking at @sonystarmap's answer, I added a few types of containers: double*, vector<double> and vector<vector<double> >. I also added tests where the "pointer-containers" are memset, since vector initialises all memory.

The C++ code was compiled with these optimization flag: -O3 -march=native

The results:

Matlab: Elapsed time is 28.457788 seconds.

C++:

T=23844.2ms

T=25161.5ms

T=25154ms

T=24197.3ms

T=24235.2ms

T=24166.1ms

I can essentially not find the large gain you mention.

#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <ctime>
#include <vector>
#include <cstring>

using namespace std;


int main(){

    const int N = 5000;

    {
        vector<double> A(N*N);

        srand(0);

        clock_t start = clock();
        for (int k=0;k<100;k++){
            for (int i=0;i<N;i++){
                for (int j=0;j<N;j++){
                    A[i*N+j] = rand();
                }
            }
        }

        cout << "T="<< (clock()-start)/(double)(CLOCKS_PER_SEC/1000)<< "ms " << endl;
    }

    {
        vector<vector<double> > A(N);
        for (int i=0;i<N;i++)
            A[i] = vector<double>(N);

        srand(0);

        clock_t start = clock();
        for (int k=0;k<100;k++){
            for (int i=0;i<N;i++){
                for (int j=0;j<N;j++){
                    A[i][j] = rand();
                }
            }
        }

        cout << "T="<< (clock()-start)/(double)(CLOCKS_PER_SEC/1000)<< "ms " << endl;
    }

    {
        double ** A = new double*[N];
        for (int i=0;i<N;i++)
            A[i] = new double[N];

        srand(0);

        clock_t start = clock();
        for (int k=0;k<100;k++){
            for (int i=0;i<N;i++){
                for (int j=0;j<N;j++){
                    A[i][j] = rand();
                }
            }
        }

        cout << "T="<< (clock()-start)/(double)(CLOCKS_PER_SEC/1000)<< "ms " << endl;
    }

    {
        double ** A = new double*[N];
        for (int i=0;i<N;i++) {
            A[i] = new double[N];
            memset(A[i], 0, sizeof(double) * N);
        }

        srand(0);

        clock_t start = clock();
        for (int k=0;k<100;k++){
            for (int i=0;i<N;i++){
                for (int j=0;j<N;j++){
                    A[i][j] = rand();
                }
            }
        }

        cout << "T="<< (clock()-start)/(double)(CLOCKS_PER_SEC/1000)<< "ms " << endl;
    }

    {
        double * A = new double[N * N];

        srand(0);

        clock_t start = clock();
        for (int k=0;k<100;k++){
            for (int i=0;i<N;i++){
                for (int j=0;j<N;j++){
                    A[i*N + j] = rand();
                }
            }
        }

        cout << "T="<< (clock()-start)/(double)(CLOCKS_PER_SEC/1000)<< "ms " << endl;
    }

    {
        double * A = new double[N * N];
        memset(A, 0, sizeof(double) * N * N);

        srand(0);

        clock_t start = clock();
        for (int k=0;k<100;k++){
            for (int i=0;i<N;i++){
                for (int j=0;j<N;j++){
                    A[i*N + j] = rand();
                }
            }
        }

        cout << "T="<< (clock()-start)/(double)(CLOCKS_PER_SEC/1000)<< "ms " << endl;
    }
}

Answer 4

#include <vector>
#include <iostream>
#include <cstdlib>
#include <ctime>
#include <cstring>

int main() {
  const int N = 5000;
  std::vector<int> A(N*N);
  srand(0);
  clock_t start = clock();
  for(int k = 0; k < 100; ++k){
    for(int i = 0; i < N * N; ++i) {
        A[i] = rand();
    }
  }
  std::cout << (clock()-start)/(double)(CLOCKS_PER_SEC/1000) << "ms" << "\n";
  return 0;
}

Went from 25-27 seconds on my workstation without any optimization flag on the compiler to 21 seconds with

-O3 -g -Wall -ftree-vectorizer-verbose=5 -msse -msse2 -msse3 -march=native -mtune=native -ffast-math