A2. Parallel Programming in C

Embarrassingly parallel computational problems are the easiest to parallelize and you can achieve impressive speedups if you have a computer with many cores. Even if you have just two cores, you can get close to a two-times speedup. An example of an embarrassingly parallel problem is when you need to run a preprocessing pipeline on datasets collected for 15 subjects. Each subject's data can be processed independently of the others. In other words, the computations involved in processing one subject's data do not in any way depend on the results of the computations for processing some other subject's data.

As an example, a grad student in my lab (Heather) figured out how to distribute her FSL preprocessing pipeline for 24 fMRI subjects across multiple cores on her Mac Pro desktop (it has 8) and as a result what used to take about 48 hours to run, now takes "just" over 6 hours.

In contrast to embarrassingly parallel problems, there is a class of problems that cannot be split into independent sub-problems, we can call them inherently sequential or serial problems. For these types of problems, the computation at one stage does depend on the results of a computation at an earlier stage, and so it is not so easy to parallelize across independent processing units. In these kinds of problems, there is a need for some communication or coordination between sub-tasks.

An example of a serial problem is a simulation of an arm movement. We run simulations of arm movements like reaching, that use detailed mathematical models of muscle mechanics, activation dynamics, musculoskeletal dynamics and spinal reflexes. Differential equations govern the relationship between muscle stimulation (the input) and the resulting arm movement (the output). These equations are "unwrapped" in time by a differential equation integrator, that takes small steps (like 1 millisecond at a time) to generate a simulation of a whole movement (e.g. 1 second of simulated time). On each step the current state of the system depends on both the current input (muscle command) and on the previous state of the system. With a 1 ms step, it takes (at least) 1000 computations to simulate a 1 sec arm movement … but we cannot simply split up those 1000 computations and distribute them to a set of independent processing units. This is an inherently serial problem where the current computation cannot be carried out without the results of the previous computation.

As an aside, the way to take advantage of parallelism even with a serial problem, is to parallelize meta-computations. So for example, we typically want to run "sensitivity analyses" where we vary some parameter(s) of the arm model and re-run the simulation. If we have 100 such simulations to run, even though each one of them is a serial problem on its own, they are independent of each other, and so we can parallelize the sensitivity analysis by distributing those 100 simulations to multiple processing units.

A good example of a problem that has both embarrassingly parallel properties as well as serial dependency properties, is the computations involved in training and running an artificial neural network (ANN). An ANN is made up of several layers of neuron-like processing units, each layer having many (even hundreds or thousands) of these units. If the ANN is a pure feedforward architecture, then computations within each layer are embarrassingly parallel, while computations between layers are serial.

By using the private() and shared() directives, you can specify variables within the parallel region as being shared, i.e. visible and accessible by all threads simultaneously, or private, i.e. private to each thread, meaning each thread will have its own local copy. In the code example below for parallelizing a for loop, you can see that we specify the thread_id and nloops variables as private.

Parallelizing for loops is really simple (see code below). By default, loop iteration counters in OpenMP loop constructs (in this case the i variable) in the for loop are set to private variables.

// gcc -fopenmp -o go go.c
// ./go

#include <stdio.h>
#include <omp.h>

int main(int argc, char **argv)
{
  int i, thread_id, nloops;

#pragma omp parallel private(thread_id, nloops)
  {
    nloops = 0;

#pragma omp for
    for (i=0; i<1000; ++i)
      {
        ++nloops;
      }

    thread_id = omp_get_thread_num();

    printf("Thread %d performed %d iterations of the loop.\n",
           thread_id, nloops );
  }

  return 0;
}

plg@wildebeest:~/Desktop$ gcc -fopenmp -o go go.c
plg@wildebeest:~/Desktop$ ./go
Thread 4 performed 125 iterations of the loop.
Thread 7 performed 125 iterations of the loop.
Thread 2 performed 125 iterations of the loop.
Thread 6 performed 125 iterations of the loop.
Thread 5 performed 125 iterations of the loop.
Thread 0 performed 125 iterations of the loop.
Thread 3 performed 125 iterations of the loop.
Thread 1 performed 125 iterations of the loop.

Using OpenMP you can specify something called a "critical" section of code. This is code that is performed by all threads, but is only performed one thread at a time (i.e. in serial). This provides a convenient way of letting you do things like updating a global variable with local results from each thread, and you don't have to worry about things like other threads writing to that global variable at the same time (a collision).

#include <stdio.h>
#include <omp.h>

int main(int argc, char *argv[])
{
    int i, thread_id;
    int glob_nloops, priv_nloops;
    glob_nloops = 0;

    // parallelize this chunk of code
    #pragma omp parallel private(priv_nloops, thread_id)
    {
        priv_nloops = 0;
        thread_id = omp_get_thread_num();

        // parallelize this for loop
        #pragma omp for
        for (i=0; i<100000; ++i)
        {
            ++priv_nloops;
        }

        // make this a "critical" code section
        #pragma omp critical
        {
            printf("Thread %d is adding its iterations (%d) to sum (%d), ",
                   thread_id, priv_nloops, glob_nloops);
            glob_nloops += priv_nloops;
            printf(" total nloops is now %d.\n", glob_nloops);
        }
    }
    printf("Total # loop iterations is %d\n",
           glob_nloops);
    return 0;
}

plg@wildebeest:~/Desktop$ gcc -fopenmp -o go go.c
plg@wildebeest:~/Desktop$ ./go
Thread 1 is adding its iterations (12500) to sum (0),  total nloops is now 12500.
Thread 4 is adding its iterations (12500) to sum (12500),  total nloops is now 25000.
Thread 0 is adding its iterations (12500) to sum (25000),  total nloops is now 37500.
Thread 5 is adding its iterations (12500) to sum (37500),  total nloops is now 50000.
Thread 3 is adding its iterations (12500) to sum (50000),  total nloops is now 62500.
Thread 6 is adding its iterations (12500) to sum (62500),  total nloops is now 75000.
Thread 2 is adding its iterations (12500) to sum (75000),  total nloops is now 87500.
Thread 7 is adding its iterations (12500) to sum (87500),  total nloops is now 100000.
Total # loop iterations is 100000

Reduction refers to the process of combining the results of several sub-calculations into a final result. This is a very common paradigm (and indeed the so-called "map-reduce" framework used by Google and others is very popular). Indeed we used this paradigm in the code example above, where we used the "critical code" directive to accomplish this. The map-reduce paradigm is so common that OpenMP has a specific directive that allows you to more easily implement this.

#include <stdio.h>
#include <omp.h>

int main(int argc, char *argv[])
{
    int i, thread_id;
    int glob_nloops, priv_nloops;
    glob_nloops = 0;

    // parallelize this chunk of code
    #pragma omp parallel private(priv_nloops, thread_id) reduction(+:glob_nloops)
    {
        priv_nloops = 0;
        thread_id = omp_get_thread_num();

        // parallelize this for loop
        #pragma omp for
        for (i=0; i<100000; ++i)
        {
            ++priv_nloops;
        }
        glob_nloops += priv_nloops;
    }
    printf("Total # loop iterations is %d\n",
           glob_nloops);
    return 0;
}

plg@wildebeest:~/Desktop$ gcc -fopenmp -o go go.c
plg@wildebeest:~/Desktop$ ./go
Total # loop iterations is 100000

There are a host of other directives you can issue using OpenMP, see here for a list (wikipedia). Some other clauses of interest are:

• barrier: each thread will wait until all threads have reached this point in the code, before proceeding
• nowait: threads will not wait until everybody is finished
• schedule(type, chunk) allows you to specify how tasks are spawned out to threads in a for loop. There are three types of scheduling you can specify
• if: allows you to parallelize only if a certain condition is met
• … and a host of others

We saw with OpenMP that we can use a reduce directive to sum values across all threads. A similar function exists in MPI called MPI_Reduce().

// Estimating pi using the dartboard algorithm
// All processes contribute to the calculation, with the
// master process averaging the values for pi.
// We then use mpc_reduce to collect the results
//
// mpicc -o go mpi_pi_reduce.c
// mpirun -n 8 go

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

#define MASTER 0         // task ID of master task
#define NDARTS 1000      // # dart throws per round
#define NROUNDS 10     // # of rounds of dart throwing

// our function for throwing darts and estimating pi
double dartboard(int ndarts)
{
  double x, y, r, pi; 
  int n, hits;
  hits = 0;

  // throw darts
  for (n = 1; n <= ndarts; n++)  {
    // (x,y) are random between -1 and 1
    r = (double)random()/RAND_MAX;
    x = (2.0 * r) - 1.0;
    r = (double)random()/RAND_MAX;
    y = (2.0 * r) - 1.0;
    // if our random dart landed inside the unit circle, increment the score
    if (((x*x) + (y*y)) <= 1.0) {
      hits++;
    }
  }

  // estimate pi
  pi = 4.0 * (double)hits / (double)ndarts;
  return(pi);
} 

// the main program
int main (int argc, char *argv[])
{
  double my_pi, pi_sum, pi_est, mean_pi, err;
  int task_id, n_tasks, rc, i;
  MPI_Status status;

  // Obtain number of tasks and task ID
  MPI_Init(&argc,&argv);
  MPI_Comm_size(MPI_COMM_WORLD,&n_tasks);
  MPI_Comm_rank(MPI_COMM_WORLD,&task_id);
  //  printf ("task %d of %d reporting for duty...\n", task_id, n_tasks);

  // different seed for random number generator for each task
  srandom (task_id);

  mean_pi = 0.0;
  for (i=0; i<NROUNDS; i++) {
    // all tasks will execute dartboard() to calculate their own estimate of pi
    my_pi = dartboard(NDARTS);

    // now we use MPI_Reduce() to sum values of my_pi across all tasks
    // the master process (id=MASTER) will store the accumulated value
    // in pi_sum. We tell MPI_Reduce() to sum by passing it
    // the MPI_SUM value (define in mpi.h)
    rc = MPI_Reduce(&my_pi, &pi_sum, 1, MPI_DOUBLE, MPI_SUM,
                    MASTER, MPI_COMM_WORLD);

    // now, IF WE ARE THE MASTER process, we will compute the mean
    if (task_id == MASTER) {
      pi_est = pi_sum / n_tasks;
      mean_pi = ( (mean_pi * i) + pi_est ) / (i + 1); // running average
      err = mean_pi - 3.14159265358979323846;
      printf("%d throws: mean_pi %.12f: error %.12f\n",
             (NDARTS * (i + 1)), mean_pi, err);
    }
  }
  if (task_id == MASTER) 
    printf ("PS, the real value of pi is about 3.14159265358979323846\n");

  MPI_Finalize();
  return 0;
}

Here we run it with just one parallel process:

plg@wildebeest:~/Desktop/mpi$ time mpirun -n 1 go
1000 throws: mean_pi 3.088000000000: error -0.053592653590
2000 throws: mean_pi 3.104000000000: error -0.037592653590
3000 throws: mean_pi 3.101333333333: error -0.040259320256
4000 throws: mean_pi 3.120000000000: error -0.021592653590
5000 throws: mean_pi 3.124800000000: error -0.016792653590
6000 throws: mean_pi 3.127333333333: error -0.014259320256
7000 throws: mean_pi 3.134285714286: error -0.007306939304
8000 throws: mean_pi 3.128500000000: error -0.013092653590
9000 throws: mean_pi 3.132444444444: error -0.009148209145
10000 throws: mean_pi 3.119600000000: error -0.021992653590
PS, the real value of pi is about 3.14159265358979323846

real	0m0.032s
user	0m0.020s
sys	0m0.012s

Now let's run it with 4:

plg@wildebeest:~/Desktop/mpi$ time mpirun -n 4 go
1000 throws: mean_pi 3.105000000000: error -0.036592653590
2000 throws: mean_pi 3.122500000000: error -0.019092653590
3000 throws: mean_pi 3.122000000000: error -0.019592653590
4000 throws: mean_pi 3.137750000000: error -0.003842653590
5000 throws: mean_pi 3.143600000000: error 0.002007346410
6000 throws: mean_pi 3.140166666667: error -0.001425986923
7000 throws: mean_pi 3.142000000000: error 0.000407346410
8000 throws: mean_pi 3.140250000000: error -0.001342653590
9000 throws: mean_pi 3.136666666667: error -0.004925986923
10000 throws: mean_pi 3.135000000000: error -0.006592653590
PS, the real value of pi is about 3.14159265358979323846

real	0m0.034s
user	0m0.044s
sys	0m0.024s

We see the final error is much reduced. Each of the 4 processes (which are parallelized across the cores of my CPU) contributes an estimate of pi, which are then averaged by the master process to come up with the final estimate of pi.