Work Sharing Constructs (Loop)

Serial version discussion

To begin to understand the work-sharing constructs, we need to learn how to parallelise the for (C/C++) or do (FORTRAN) loop. For this, we will learn simple vector addition examples.

As we can see from the above figure, the two vectors should be added to get a single vector. This is done by iterating over the elements and adding them together. Since there are no data dependencies, the loop indexes do not have any data dependency on the other indexes. Therefore, it is easy to parallelise.

Examples: Loop

Serial (C/C++)

for(int i = 0; i < n; i ++)
  {
    c[i] = a[i] + b[i];
  }

Serial (FORTRAN)

do i = 1, n
  c(i) = a(i) + b(i)
end do

Note

FORTRAN has a column-major order and C/C++ has a row-major order

Fortran array index starts from 1
C/C++ array index starts from 0

Parallel version discussion

Now we will look into how to parallelise the for (C/C++) or do (FORTRAN) loops. We just need to add the OpenMP directives below.

Functionality	Syntax in C/C++	Syntax in FORTRAN
Distribute iterations over threads	#pragma omp for	!$omp do

With the help of the above syntax, the loops can be easily parallelised. The figure below shows an example of how the loops are parallelised. Using 5 threads would divide 10 iterations by two — each thread handles 2 iterations.

Examples: Loop parallelisation

C/C++

#pragma omp parallel for
for(int i = 0; i < n; i ++)
  {
    c[i] = a[i] + b[i];
  }

// or

#pragma omp parallel
#pragma omp for
for(int i = 0; i < n; i ++)
  {
    c[i] = a[i] + b[i];
  }

FORTRAN

!$omp parallel do
do i = 1, n
  c(i) = a(i) + b(i)
end do
!$omp end parallel do

! or

!$omp parallel
!$omp do
do i = 1, n
  c(i) = a(i) + b(i)
end do
!$omp end do
!$omp end parallel

Questions and Solutions

Examples: Vector Addition

Serial (C/C++)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <time.h>

#define N 5120
#define MAX_ERR 1e-6

float * Vector_Add(float *a, float *b, float *c, int n)
{
  for(int i = 0; i < n; i ++)
    { c[i] = a[i] + b[i]; }
  return c;
}

int main()
{
  float *a, *b, *c;
  a = (float*)malloc(sizeof(float) * N);
  b = (float*)malloc(sizeof(float) * N);
  c = (float*)malloc(sizeof(float) * N);

  for(int i = 0; i < N; i++)
    { a[i] = 1.0f; b[i] = 2.0f; }

  clock_t start = clock();
  Vector_Add(a, b, c, N);
  clock_t end = clock();
  printf("Time measured: %.3f seconds.\n", (double)(end - start)/CLOCKS_PER_SEC);

  for(int i = 0; i < N; i++)
    { assert(fabs(c[i] - a[i] - b[i]) < MAX_ERR); }

  printf("c[0] = %f\n", c[0]);
  printf("PASSED\n");

  free(a); free(b); free(c);
  return 0;
}

Template (C/C++)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <time.h>

#define N 5120
#define MAX_ERR 1e-6

float * Vector_Add(float *a, float *b, float *c, int n)
{
// ADD YOUR PARALLEL REGION FOR THE LOOP
  for(int i = 0; i < n; i ++)
    { c[i] = a[i] + b[i]; }
  return c;
}

int main()
{
  float *a, *b, *c;
  a = (float*)malloc(sizeof(float) * N);
  b = (float*)malloc(sizeof(float) * N);
  c = (float*)malloc(sizeof(float) * N);

  for(int i = 0; i < N; i++)
    { a[i] = 1.0f; b[i] = 2.0f; }

  clock_t start = clock();
  // ADD YOUR PARALLEL REGION HERE
  Vector_Add(a, b, c, N);
  clock_t end = clock();
  printf("Time measured: %.3f seconds.\n", (double)(end - start)/CLOCKS_PER_SEC);

  for(int i = 0; i < N; i++)
    { assert(fabs(c[i] - a[i] - b[i]) < MAX_ERR); }

  printf("c[0] = %f\n", c[0]);
  printf("PASSED\n");

  free(a); free(b); free(c);
  return 0;
}

Solution (C/C++)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <omp.h>

#define N 5120
#define MAX_ERR 1e-6

float * Vector_Add(float *a, float *b, float *c, int n)
{
#pragma omp for
  for(int i = 0; i < n; i ++)
    { c[i] = a[i] + b[i]; }
  return c;
}

int main()
{
  float *a, *b, *c;
  a = (float*)malloc(sizeof(float) * N);
  b = (float*)malloc(sizeof(float) * N);
  c = (float*)malloc(sizeof(float) * N);

  for(int i = 0; i < N; i++)
    { a[i] = 1.0f; b[i] = 2.0f; }

  omp_set_num_threads(omp_get_max_threads());
  double start = omp_get_wtime();
  #pragma omp parallel
  Vector_Add(a, b, c, N);
  double end = omp_get_wtime();
  printf("Time measured: %.3f seconds.\n", end - start);

  for(int i = 0; i < N; i++)
    { assert(fabs(c[i] - a[i] - b[i]) < MAX_ERR); }

  printf("c[0] = %f\n", c[0]);
  printf("PASSED\n");

  free(a); free(b); free(c);
  return 0;
}

Compilation and Output

Serial (C/C++)

// compilation
$ gcc Vector-addition-Serial.c -o Vector-addition-Serial

// execution
$ ./Vector-addition-Serial

Solution (C/C++)

// compilation
$ gcc -fopenmp Vector-addition-OpenMP-solution.c -o Vector-addition-Solution

// execution
$ ./Vector-addition-Solution

Solution (FORTRAN)

// compilation
$ gfortran -fopenmp Vector-addition-OpenMP-solution.f90 -o Vector-addition-Solution

// execution
$ ./Vector-addition-Solution

Questions

• Can you measure the performance speedup for parallelising loop? Do you see any speedup?

• Can you create more threads to speed up the computation? If yes or no, why?