Multicore Programming Models: Threads vs OpenMP

The Problems With Threads

#include <pthread.h>

volatile int delay=0;
void doWork(const char *caller) {
	std::cout<<caller<<"\n";
	for (int i=0;i<1000*1000*1;i++) delay++;
}
void *fnA(void *arg) {
	for (int i=0;i<10;i++) doWork("A");
	return 0;
}
void *fnB(void *arg) {
	for (int i=0;i<10;i++) doWork("B");
	return 0;
}
int foo(void) {
	pthread_t B; 
	pthread_create(&B,0,fnB,0); 
	fnA(0); 
	void *Bret;
	pthread_join(B,&Bret); 
	return delay;
}

(Try this in NetRun now!)

• Moving the per-thread work into a separate function is a pain.
• You need to remember to both create *and* join threads.
• The A and B printouts come in some unpredictable order, depending on the speed of the machine, OS, etc.
  
• "delay" ends up with a weird random answer (read the full story here).
  

A New Hope: OpenMP

    for (int i=0;i<n;i++) do_fn(i);

#pragma omp parallel for
    for (int i=0;i<n;i++) do_fn(i);

• The compiler decides how many threads to create (typically the number of CPUs)
• The compiler moves the guts of your loop off into a separate function.
• The compiler creates and deallocates the threads.
• The compiler can help you find and fix memory race conditions.
  

OpenMP Gotchas

volatile int sum=0;
int foo(void) {
	int i,j;
	sum=0;
	for (i=0;i<1000;i++) 
	for (j=0;j<1000;j++) {
		sum++;
	}
	return sum;
}

(Try this in NetRun now!)

volatile int sum=0;
int foo(void) {
	int i,j;
	sum=0;
#pragma omp parallel for
	for (i=0;i<1000;i++) 
	for (j=0;j<1000;j++) {
		sum++;
	}
	return sum;
}

(Try this in NetRun now!)

• "j" is shared between the threads.
• "sum" is shared between the threads.
  

volatile int sum=0;
int foo(void) {
	int i,j;
	sum=0;
#pragma omp parallel for private(j)
	for (i=0;i<1000;i++) 
	for (j=0;j<1000;j++) {
		sum++;
	}
	return sum;
}

(Try this in NetRun now!)

volatile int sum=0;
int foo(void) {
	int i,j;
	sum=0;
#pragma omp parallel for private(j) reduction(+:sum)
	for (i=0;i<1000;i++) 
	for (j=0;j<1000;j++) {
		sum++;
	}
	return sum;
}

(Try this in NetRun now!)

This is now 0.44ns/iteration, a solid 4x speedup over the original code, and it gets the right answer!

Even in OpenMP, you can still get multithreaded race conditions, especially when calling library code.  The compiler can help you fix these, but often you need to explicitly ask it for help.

False Sharing

#include <omp.h>

const int space=16;
volatile int delay[4*space]={0};
int foo(void) {
	int i;
	for (int core=0;core<4;core++) delay[core*space]=0;
#pragma omp parallel for //reduction(+:delay) 
	for (i=0;i<1000;i++) 
	{
		int d=(i/250)*space;
		for (int j=0;j<1000;j++) 
		{
			delay[d]++;
		}
	}
	return delay[0*space];
}

(Try this in NetRun now!)

OpenMP Features

• "if(n>100)" only creates parallel threads if the number of iterations is big enough; otherwise it's a normal sequential loop.  You can use any condition there.
  
• "private(n)" would give each CPU a separate read-write copy of the variable "n" (or any other variable).  By default, the new copy is uninitialized.  WARNING: there's only one shared copy of all variables declared outside the parallel loop, unless you explicitly declare them private!
• "default(none)" is for error-checking your loops: it verifies that every variable used in the loop is either local, or explicitly private or shared.
• "shared(n)" declares that we only want a single shared copy of n.  This is the default for variables declared outside the loop.
  
• "copyin(n)" initializes n from the master's version.
  
• "firstprivate(n)" is like private, but initializes the new copy from the master.  It's basically private and copyin together.
  
• "lastprivate(n)" is also like private, but initializes each copy from the last-used version.
• "reduction(+:k)" would total up each CPU's (private) copy of the variable "k" (or any other variable) using the "+" operation (or any other operation).  The compiler generates atomic ("lock" prefix) instructions to do this.
  
• "schedule(static,10)" tells each thread to do blocks of ten iterations of the loop, with blocks distributed beforehand in round-robin fashion.  The default scheduling is basically "schedule(static,n/num_threads)".
  
• "schedule(dynamic,10)" tells each thread to do ten iterations of the loop at a time, but each thread takes the next available set of iterations.
• "schedule(guided,10)" takes at least ten iterations at a time, dynamically, but usually more.
  
• "schedule(runtime)" says to guess the scheduling at runtime.
  

• Just "pragma omp parallel" around a random block creates separate threads; call omp_get_thread_num() to figure out who you are.
  
• "barrier" tells each thread to not to continue until all threads have reached this point. 
  
• "critical" only allows one thread at a time to execute the following block (wraps a lock around the block).
• "atomic" makes a simple ++ or += expression atomic (using atomic instructions).
  
• "master" makes a block of code only run on the master.
• "single" makes a block of code only run on a single (random) thread.
• "ordered" makes a part of the loop execute in sequential order.  This might be more cache efficient than just running a second sequential loop from the master thread.
  
• "parallel sections" lets you statically create several parallel pieces of work.

OpenMP Problems

OpenMP Sorting

#include <omp.h>
#include <cassert>
#include <algorithm> /* for std::sort */

enum{n=8000};
float arr[n+1]; /* silly: farray_fill uses n+1 elements */
float darr[n+1];
int foo(void) 
{
  farray_fill(arr,n,0.1); /* fill with random values */

/* Parallel stage 1: sort sub-pieces */
  int nth=1; /* runtime thread count */
  enum{MAXTHREADS=1024}; /* STUPID maximum thread count */
  float *subpieces[MAXTHREADS]; /* subarray for each thread */
  int subn[MAXTHREADS]; /* number of elements in each subarray */
#pragma omp parallel  /*  each thread sorts a subrange of the array */
  {
     int tid=omp_get_thread_num();
     nth=omp_get_num_threads();
     assert(nth<=MAXTHREADS);
     float *marr=&arr[tid*n/nth];
     subpieces[tid]=marr;
     int mn=n/nth; /* FIXME: what if nth doesn't divide n evenly?! */
     subn[tid]=mn;
     printf("Thread %d: %d to %d more\n",tid,tid*n/nth,mn);
     std::sort(marr,marr+mn);
  }

/* Stage 2: merge sub-pieces */
  int di=0;
  int subi[MAXTHREADS]; /* index into this thread's subarray */
  for (int i=0;i<nth;i++) subi[i]=0;
  for (di=0;di<n;di++) {
     int ts=-1; /* thread with smallest unconsumed number so far */
     float smallest=1.0e30;
     for (int t=0;t<nth;t++) {
       if (subi[t]<subn[t]) /* this thread still has values to offer */
         if (subpieces[t][subi[t]]<smallest) 
         { /* this thread has the new smallest value */
             ts=t;
             smallest=subpieces[t][subi[t]];
         }
     }
     subi[ts]++; /* consume that value */
     darr[di]=smallest; /* copy to output array */
  }

  farray_print(arr,n);
  farray_print(darr,n);
  for (int i=0;i<n-1;i++) if (darr[i]>darr[i+1]) {printf("Not sorted at place %d!\n",i); return i;}
  printf("Array sorted properly\n");
  return 0;
}

(Try this in NetRun now!)