What is OpenMP

What is OpenMP?

• Used for multi-threaded parallel processing
• Used on shared-memory multi-processor (core) computers
• Part of program is a single thread and part is multi-threaded
• Has 3 components
• directives that can be put into C/C++ or Fortran programs
• runtime library for setting and querying parallel parameters (ex. # of threads)
• environment variables that define runtime parallel parameters (ex. # of threads)

c Fortran example
 call omp_set_num_threads(nthread) !requests "nthread" threads
!$OMP PARALLEL DO
 DO i=1,N
 DO j=1,M
 .
 .
 .
 END DO
 END DO
!$OMP END PARALLEL DO

omp_set_num_threads(nthread); /* requests nthread threads */
#pragma omp parallel for
{
 for (i=0; i<n; i++) {
 for (j=0; j<m; j++) {
 .
 .
 .
 }
 }
}

Memory Architectures and Parallel Programming

Pros and Cons Of OpenMP

Approaches to Parallelism Using OpenMP

Example OpenMP Hello World

#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
int main (int argc, char *argv[]) {

int nthreads, tid;

/* Fork a team of threads giving them their own copies of variables */
#pragma omp parallel private(nthreads, tid)
 {

 /* Obtain thread number */
 tid = omp_get_thread_num();
 printf("Hello World from thread = %d\n", tid);

 /* Only master thread does this */
 if (tid == 0) 
 {
  nthreads = omp_get_num_threads();
 printf("Number of threads = %d\n", nthreads);
 }

 } /* All threads join master thread and disband */

}

 PROGRAM HELLO
 
 INTEGER NTHREADS, TID, OMP_GET_NUM_THREADS,
 + OMP_GET_THREAD_NUM
 
C Fork a team of threads giving them their own copies of variables
!$OMP PARALLEL PRIVATE(NTHREADS, TID)


C Obtain thread number
 TID = OMP_GET_THREAD_NUM()
 PRINT *, 'Hello World from thread = ', TID

C Only master thread does this
 IF (TID .EQ. 0) THEN
  NTHREADS = OMP_GET_NUM_THREADS()
 PRINT *, 'Number of threads = ', NTHREADS
 END IF

C All threads join master thread and disband
!$OMP END PARALLEL

 END

Loop level Parallelizaion

• no dependencies between loop indicies
• an element of an array is assigned to by at most one iteration
• no loop iteration reads array elements modified by any other dependency
• due to overhead of parallelization - use only on  loops where individual iterations take a long time 

!     Fortran example

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

/* C/C++ Example */

  #pragma omp parallel for
   for(i=1; i<=n; i++)
      a[i] = b[i] + c[i]  

Example of Code With Data Dependencies

! Fortran example

do i = 2, 5
   a(i) = a(i) + a(i-1)
enddo

/* C/C++ example */

for(i=2; i<=5; i++)
   a[i] = a[i] + a[i-1];

Shared vs. Private Variables

• By default all variables in a loop share the same address space 
• All threads can modify and access all variables (except the loop index)
• Can result in incorrect results
• Can use shared and private clauses with parallel for or parallel do

!Fortran example 

!$omp parallel do private(temp) shared(n,a,b,c)
do i = 1, n
   temp = 2.0*a(i)
   a(i) = temp
   b(i) = c(i)/temp
enddo

/* C/C++ Example */

#pragma omp parallel for private(temp) shared(n,a,b,c)
{
   for(i=1; i<=n; i++){
      temp = 2.0*a[i];
      a[i] = temp;
      b[i] = c[i]/temp;
   }
}

Example of Parallelizing A Loop

#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#define N       100

int main (int argc, char *argv[]) {

int nthreads, tid, i;
float a[N], b[N], c[N];

/* Some initializations */
for (i=0; i < N; i++)
  a[i] = b[i] = i;


#pragma omp parallel shared(a,b,c,nthreads) private(i,tid)
  {
  tid = omp_get_thread_num();
  if (tid == 0)
    {
    nthreads = omp_get_num_threads();
    printf("Number of threads = %d\n", nthreads);
    }

  printf("Thread %d starting...\n",tid);

  #pragma omp for 
  for (i=0; i<N; i++)
    {
    c[i] = a[i] + b[i];
    printf("Thread %d: c[%d]= %f\n",tid,i,c[i]);
    }

  }  /* end of parallel section */

}

 
      PROGRAM WORKSHARE1

      INTEGER NTHREADS, TID, OMP_GET_NUM_THREADS,
     +  OMP_GET_THREAD_NUM, N, I
      PARAMETER (N=100) 
      REAL A(N), B(N), C(N)

!     Some initializations
      DO I = 1, N
        A(I) = I 
        B(I) = A(I)
      ENDDO

!$OMP PARALLEL SHARED(A,B,C) PRIVATE(I,TID)

      TID = OMP_GET_THREAD_NUM()
      IF (TID .EQ. 0) THEN
        NTHREADS = OMP_GET_NUM_THREADS()
        PRINT *, 'Number of threads =', NTHREADS
      END IF
      PRINT *, 'Thread',TID,' starting...'

!$OMP DO 
      DO I = 1, N
        C(I) = A(I) + B(I)
        WRITE(*,100) TID,I,C(I)
 100    FORMAT(' Thread',I2,': C(',I3,')=',F8.2)
      ENDDO
!$OMP END DO NOWAIT

      PRINT *, 'Thread',TID,' done.'

!$OMP END PARALLEL

      END

Basic OpenMP Functions

#include <omp.h>
int omp_get_thread_num();
# pragma omp parallel
{
  printf("Thread rank: %d\n", omp_get_thread_num());
}

INTEGER FUNCTION OMP_GET_THREAD_NUM()
!$OMP PARALLEL
      write(*,*)'Thread rank: ', OMP_GET_THREAD_NUM()
!$OMP END PARALLEL

Thread rank:  2
Thread rank:  0
Thread rank:  3
Thread rank:  1

How to Compile and Run an OpenMP Program

$ gcc -o omp_helloc -fopenmp omp_hello.c
$ export OMP_NUM_THREADS=2
$ ./omp_helloc
Hello World from thread = 0
Hello World from thread = 1
Number of threads = 2
$
$ gfortran -o omp_hellof -fopenmp omp_hello.f
$  export OMP_NUM_THREADS=4
$ ./omp_hellof
 Hello World from thread =            2
 Hello World from thread =            1
 Hello World from thread =            3
 Hello World from thread =            0
 Number of threads =            4

$ icc -o omp_helloc -openmp omp_hello.c
omp_hello.c(22): (col. 1) remark: OpenMP DEFINED REGION WAS PARALLELIZED.
$ export OMP_NUM_THREADS=3
$ ./omp_helloc
Hello World from thread = 0
Hello World from thread = 2
Hello World from thread = 1
Number of threads = 3
$
$ ifort -o omp_hellof -openmp omp_hello.f
omp_hello.f(20): (col. 7) remark: OpenMP DEFINED REGION WAS PARALLELIZED.
$ export OMP_NUM_THREADS=2
$ ./omp_hellof
 Hello World from thread =            0
 Number of threads =            2
 Hello World from thread =            1

[sas@discovery intro_openmp]$ pgcc -o omp_helloc -mp omp_hello.c
[sas@discovery intro_openmp]$ export OMP_NUM_THREADS=2
[sas@discovery intro_openmp]$ ./omp_helloc
Hello World from thread = 0
Number of threads = 2
Hello World from thread = 1
$
$ pgf90 -o omp_hellof -mp omp_hello.f
$ export OMP_NUM_THREADS=2
$ ./omp_hellof
 Hello World from thread =             0
 Number of threads =             2
 Hello World from thread =             1

Other OpenMP Clauses

• firstprivate
• lastprivate
• ordered

• initialize a variable from the serial part of the code
• private clause doesn't initialize the variable

!     Fortran example

j = jstart
!$omp parallel do firstprivate(j)
do i = 1, n
   if(i == 1  .or. i == n) then
      j = j + 1
   endif
   a(i) = a(i) + j
enddo

/* C/C++ Example */

 j = jstart;
#pragma omp parallel for firstprivate(j)
{
   for(i=1; i<=n; i++){
      if(i == 1 || i == n)
         j = j + 1;
      a[i] = a[i] + j;
   }
}

• thread that executes the ending loop index copies its value to the master (serial) thread
• this gives the same result as serial execution

Example of  lastprivate clause

! Fortran example

!$omp parallel do lastprivate(x)
do i = 1, n
   x = sin(pi*dx*real(i))
   a(i) = exp(x)
enddo
lastx = x

/* C/C++ example */

#pragma omp parallel for lastprivate(x)
{
   for(i=1; i<=n; i++){
      x = sin( pi * dx * (float)i );
      a[i] = exp(x);
   }
}
lastx = x;

• used when part of the loop must execute in serial order
• ordered clause plus an ordered directive

! Fortran example

!$omp parallel do private(myval) ordered
do i = 1, n
   myval = do_lots_of_work(i)
   !$omp ordered
   print*, i, myval
   !$omp end ordered
enddo
lastx = x

/* C/C++ example */

#pragma omp parallel for private(myval) ordered
{
   for(i=1; i<=n; i++){
      myval = do_lots_of_work(i);
      #pragma omp ordered
      {
         printf("%d %d\n", i, myval);
      }
   }
}

Reduction Operations

!     Fortran example

do i = 1, n
   sum = sum + a(i)
enddo

/* C/C++ Example */

for(i=1; i<=n; i++){
   sum = sum + a[i];
}

• sum is the reduction variable
• cannot be declared shared
• threads would overwrite the value of sum
• cannot be declared private
• private variables don't persist outside of parallel region
• specified reduction operation performed on individual values from each thread

Example of  reduction clause

! Fortran example

!$omp parallel do reduction(+:sum)
do i = 1, n
   sum = sum + a(i)
enddo

/* C/C++ example */

#pragma omp parallel for reduction(+:sum)
{
   for(i=1; i<=n; i++){
      sum = sum + a[i];
   }
}

Fortran Reduction Operands

Thread Control

!     Fortran example

!$omp parallel private(myid,istart,iend)
call myrange(myid, nthreads, istart, iend)
do i = istart, iend
   a(i) = a(i) - b(i)
enddo
!$omp barrier
call dowork(a)
!$omp end parallel

/* C/C++ Example */

#pragma omp parallel private(myid, istart, iend)
{
  myrange(myid, nthreads,  &istart, &iend);
  for(i=istart; i<=iend; i++){
     a[i] = a[i] - b[i];
  }
  #pragma omp barrier
  dowork(a);
}

! Fortran example

!$omp parallel private(myid, istart, iend)
call myrange(myid, nthreads, global_start, global_end, istart, iend)
do i = istart, iend
   a(i) = b(i)
enddo
!$omp barrier
!$omp master
write(21) a
!$omp end master
call do_work(istart, iend)
!$omp end parallel

/* C/C++ example */

#pragma omp parallel private(myid, istart, iend)
{
  myrange(myid, nthreads, global_start, global_end, &istart, &iend);
  for(i=istart; i<=iend; i++){
     a[i] = b[i];
  }
  #pragma omp barrier
  #pragma omp master
  {
  n = global_end - global_start + 1;
  write_size = fwrite(a, 1, n, file_pointer);
  }
  do_work(istart, iend);
}

Thread Control (continued) 

• Only one thread executes a specified section of the code at a time
• Threads can execute in any order
• Similar to ORDERED directive except ordered specifies that threads go in numerical order    

!     Fortran example

the_max = 0.0
!$omp parallel private(myid, istart, iend)
  call myrange(myid, nthreads, global_start, global_end, istart, iend)
  call compute_a(a(istart:iend))
    !$omp critical
    the_max = max( maxval(a(istart:iend), the_max )
    !$omp end critical
  call more_work_on_a(a)
!$omp end parallel

/* C/C++ Example */

the_max = 0.0;
#pragma omp parallel private(myid, istart, iend)
{
  myrange(myid, nthreads, global_start, global_end, &istart, &iend);
  nvals = iend-istart+1;
  compute_a(a[istart],nvals);
    #pragma omp critical
    the_max = max( maxval(a[istart],nvals), the_max );
    #pragma omp end critical
  call more_work_on_a(a)
}

• A section of code that is run by only one thread
• Sections are performed in parallel

! Fortran example

!$omp parallel
!$omp sections

  !$omp section
  call init_field(field)

  !$omp section
  call check_grid(grid)

!$omp end sections
!$omp end parallel

/* C/C++ example */

#pragma omp parallel
{
  #pragma omp sections
  {
    #pragma omp section
    init_field(field);

    #pragma omp section
    check_grid(grid);
  }
}

Profiling Your Code

• Compile with "-g -gp" (no optimization) options
• Run the program in a serial mode and produce file gmon.out
• Use gprof  to display the profile data  and the lines where the program spends the most time.

[sas@discovery]$ icc -g -pg -openmp-stubs -o matmul_par matmul_par.c
matmul_par.c(60): warning #161: unrecognized #pragma
 #pragma omp parallel for private(tmp, i, j, k) 
 ^
[sas@discovery]$ ./matmul_par
 Order 1000 multiplication in 13.498900 seconds 
 1 threads
 Hey, it worked
 all done 

[sas@discovery] ls -ls gmon.out
4 -rw-rw-r-- 1 sas sas 1317 Feb 17 11:27 gmon.out

[sas@discovery] gprof -l matmul_par >gprof.out

[sas@discovery] more gprof.out

Flat profile:

Each sample counts as 0.01 seconds.
 % cumulative self self total 
 time seconds seconds calls Ts/call Ts/call name 
 94.39 18.06 18.06 main (matmul_par.c:68 @ 400aea)
 5.72 19.15 1.09 main (matmul_par.c:66 @ 400b22)
 0.05 19.16 0.01 main (matmul_par.c:50 @ 4009f8)
 0.05 19.17 0.01 main (matmul_par.c:54 @ 400a57)
 0.05 19.18 0.01 main (matmul_par.c:70 @ 400b30)
 0.05 19.19 0.01 main (matmul_par.c:85 @ 400c39)
.
.

60 #pragma omp parallel for private(tmp, i, j, k) 
61 for (i=0; i<Ndim; i++){
62 for (j=0; j<Mdim; j++){
63 
64 tmp = 0.0;
65
66 for(k=0;k<Pdim;k++){
67 /* C(i,j) = sum(over k) A(i,k) * B(k,j) */
68 tmp += *(A+(i*Ndim+k)) * *(B+(k*Pdim+j));
69 }
70 *(C+(i*Ndim+j)) = tmp;
71 }
72 }

OpenMP v3.0

• OpenMP v2.0 for C: 2002
• OpenMP v2.0 for Fortran: 2000
• OpenMP v2.5: 2005
• OpenMP v3.0: 2008 -  tasks most important addition

• allow you to parallelize irregular problems
• unbounded loops (ex. while loops)
• recursive algorithms
• unstructured parallelism
• dynamically generated units of work
• task can be executed by any thread in the team , in parallel with others
• execution can be immediate or deferred until later
• execution might be suspended and continued later by the same or different thread
• parallel tasks can be nested within parallel loops or sections

c  C example  
   
 #pragma omp parallel
{
    #pragma omp single
    {
       p=listhead;
       while (p) {
            /* create a task for each element of the list */
             #pragma omp task
             process(p);  /* porcess list element p */
             p=next(P);
       }
    }
}