I have the following code that I am trying to optimize. I am running it on a Linux machine with 24 cores. I thought I could use multithreading to make it faster, but it's somehow making it way slower than the single thread. I believe it is because of cache misses but I am not sure how to structure this code to avoid them.

#include <stdio.h>              
#include <stdlib.h>             

#include "mex.h"
#include "string.h"
#include "pthread.h"

/* NOT sure why multi-threading is worse than single thread */
#define NTHREADS 1

/* An argument list structure */
struct argList{
    double *X;
    double *p;
    int cycles;
    int N1;
    int N2;
typedef struct argList aList;

long int bino_rand(double p, double n)
  long int    x = 0;             
  int i;
  for (i=0; i<n; i++)
    if (((double)rand())/RAND_MAX < p) x++;


void *bino_process(void* arglist){
    double* X   = ((aList*)arglist)->X;
    double* p   = ((aList*)arglist)->p;
    int cycles  = ((aList*)arglist)->cycles;
    int N1      = ((aList*)arglist)->N1;
    int N2      = ((aList*)arglist)->N2;
    int i,j;

/* mexPrintf("N1=%d, N2=%d, cycles=%d\n", N1, N2, cycles); */

for(i = N1;i<N2;i++){
    double tmp = X[i];
    double tmp_p = p[i];
        tmp = tmp + bino_rand(tmp_p,tmp);
    X[i] = tmp;

void mexFunction( int nlhs, mxArray *plhs[],
              int nrhs, const mxArray *prhs[])
double* n;              /* input scalar */
double* p;
int cycles;
mwSize N;               /* size of matrix */
double *X;              /* output matrix */
int i;

pthread_t threadA[NTHREADS];
int iterA[NTHREADS];
aList arglistA[NTHREADS];


/* get the value of the scalar input  */
n = mxGetPr(prhs[0]);
p = mxGetPr(prhs[1]);
cycles = (int)mxGetScalar(prhs[2]);

N = mxGetM(prhs[0]);

/* create the output matrix */
plhs[0] = mxCreateDoubleMatrix(N,1,mxREAL);

/* get a pointer to the real data in the output matrix */
X = mxGetPr(plhs[0]);


    N1[i] = i*(N/NTHREADS);
    N2[i] = (i+1)*(N/NTHREADS);
    arglistA[i].X = X;
    arglistA[i].p = p;
    arglistA[i].cycles = cycles;

    arglistA[i].X = (double*)mxCalloc(N2[i]-N1[i],sizeof(double));
    memcpy(arglistA[i].X, &X[N1[i]], (N2[i]-N1[i])*sizeof(double));

    arglistA[i].N2 = N2[i]-N1[i];        
    arglistA[i].N1 = 0;


/* Spawn the threads */
    iterA[i] = pthread_create(&threadA[i],NULL,bino_process,(void*)&arglistA[i]);       

/* Join back the threads */

/* Copy back the data from the chunks into X */
    memcpy(&X[N1[i]], arglistA[i].X, (N2[i]-N1[i])*sizeof(double));


  • 2
    \$\begingroup\$ As we all want to make our code more efficient or improve it in one way or another, try to write a title that summarizes what your code does, not what you want to get out of a review. \$\endgroup\$
    – Jamal
    Jun 29, 2016 at 16:14
  • 3
    \$\begingroup\$ For context, you might want to consider also including additional information such as which code / section of code you're timing, the size of the arrays you're processing, and approximate timings for single/multi threaded runs. As far as the title goes, I think @Jamal was suggesting something more along the lines of 'Multi-threaded function for calculating XXXX' where XXXX is whatever it is your code is actually doing. At the moment, your title is still a description of what you want from the review, not what the code does. \$\endgroup\$
    – forsvarir
    Jun 29, 2016 at 17:40
  • \$\begingroup\$ I still cannot quite tell what this code is supposed to do. You need to provide more context. \$\endgroup\$
    – Jamal
    Jul 7, 2016 at 3:16

2 Answers 2


How did you time your program?

I ran a modified version of your program (I didn't have "mex.h" so I rewrote the parts that needed it). When I ran your program with 1, 2, and 4 threads, I got the following results:

1 thread  = 9.38 sec
2 threads = 4.93 sec
4 threads = 2.81 sec

So it seems to be working for me. How did you time your program? You may have been using a function such as clock() that measures cpu time across all cores, instead of measuring wall clock time. See this Stackoverflow question for more information on that.

  • 1
    \$\begingroup\$ mex.h is provided by MATLAB for its native-code interface. \$\endgroup\$ Jun 29, 2016 at 18:27
  • \$\begingroup\$ I recompiled the code on my laptop and it worked as expected - time goes down as threads increase. Previously, I was running the code on a 24 core linux cluster. May be there is something about the way that server is set up. I can check with the system's admin. I am not experienced multithreading and thought I may be doing something bad, which caused the poor performance. I use MATLAB's tic-toc for timing. But also, the multithreading is so slow that it is just obvious without precise timing. The results for some input were 95secs without multithreading and 1111secs with 8 threads. \$\endgroup\$ Jun 29, 2016 at 18:48
  • \$\begingroup\$ @200_success Ah thanks, I've never used matlab before. \$\endgroup\$
    – JS1
    Jun 30, 2016 at 3:21

After some more Google searching, I found that the fundamental problem is the call to the rand() function. It is not thread-safe as it says here, and then I found this as well.

So, the solution is to use a thread-safe rand_r() function instead. This issue seems to be only on a Linux machine. There is no issue on my MacBook.


Not the answer you're looking for? Browse other questions tagged or ask your own question.