Transcript PPT - Dr. Bo Yuan
Open Multiprocessing
Dr. Bo Yuan E-mail: [email protected]
OpenMP
• An API for shared memory multiprocessing (parallel) programming in C, C++ and Fortran.
– Supports multiple platforms (processor architectures and operating systems).
– Higher level implementation (a block of code that should be executed in parallel).
• A method of parallelizing whereby a master thread forks a number of slave threads and a task is divided among them.
• Based on preprocessor directives ( Pragma ) – Requires compiler support.
– omp.h
• References – http://openmp.org/ – https://computing.llnl.gov/tutorials/openMP/ – http://supercomputingblog.com/openmp/ 2
Hello, World!
#include
Definitions
# pragma omp parallel [clauses] { code_block } implicit barrier text to modify the directive Error Checking #ifdef _OPENMP # include
The Trapezoidal Rule
/* Input: a, b, n */ h=(b-a)/n; approx=(f(a)+f(b))/2.0; for (i=1; i<=n-1; i++) { x_i=a+i*h; approx+=f(x_i); } approx=h*approx; Thread 0 Shared Memory Shared Variables Race Condition Thread 2 # pragma omp critical global_result+=my_result; 5
The critical Directive
# pragma omp critical y=f(x); ...
double f(double x) { # pragma omp critical z=g(x); ... } Cannot be executed simultaneously!
# pragma omp critical(one) y=f(x); ...
double f(double x) { # pragma omp critical(two) z=g(x); ... } 6
The atomic Directive
# pragma omp atomic x
• Only the load and store of x critical is protected.
directive.
# pragma omp atomic x+=f(y); x++ ++x x- --x # pragma omp critical x=g(x); Can be executed simultaneously!
7
Locks
/* Executed by one thread */ Initialize the lock data structure; ...
/* Executed by multiple threads */ Attempt to lock or set the lock data structure; Critical section; Unlock or unset the lock data structure; ...
/* Executed by one thread */ Destroy the lock data structure; void omp_init_lock(omp_lock_t* lock_p); void omp_set_lock(omp_lock_t* lock_p); void omp_unset_lock(omp_lock_t* lock_p); void omp_destroy_lock(omp_lock_t* lock_p); 8
Trapezoidal Rule in OpenMP
#include
Trapezoidal Rule in OpenMP
void Trap(double a, double b, int n, double* global_result_p) { double h, x, my_result; double local_a, local_b; int int int i, local_n; my_rank=omp_get_thread_num(); thread_count=omp_get_num_threads(); h=(b-a)/n; local_n=n/thread_count; local_a=a+my_rank*local_n*h; local_b=local_a+local_n*h; my_result=(f(local_a)+f(local_b))/2.0; for(i=1; i<=local_n-1; i++) { x=local_a+i*h; my_result+=f(x); } my_result=my_result*h; # pragma omp critical *global_result_p+=my_result; } 10
Scope of Variables
In serial programming: • Function-wide scope • File-wide scope Shared Scope • Accessible by all threads in a team • Declared before a parallel directive Private Scope • Only accessible by a single thread • Declared in the code block • a, b, n • global_result • thread_count • my_rank • my_result • global_result_p • *global_result_p 11
Another Trap Function
double Local_trap(double a, double b, int n); global_result=0.0; # pragma omp parallel num_threads(thread_count) { # pragma omp critical global_result+=Local_trap(a, b, n); } global_result=0.0; # pragma omp parallel num_threads(thread_count) { double my_result=0.0; # pragma omp critical /* Private */ my_result=Local_trap(a, b, n); global_result+=my_result; } 12
The Reduction Clause
• Reduction
: A computation (binary operation) that repeatedly applies the same reduction operator (e.g., addition or multiplication) to a sequence of operands in order to get a single result.
reduction(
– A private variable is created for each thread in the team.
– The private variables are initialized to 0 for addition operator.
global_result=0.0; # pragma omp parallel num_threads(thread_count)\ reduction(+: global_result) global_result=Local_trap(a, b, n); 13
The parallel for Directive
h=(b-a)/n; approx=(f(a)+f(b))/2.0; # pragma omp parallel for num_threads(thread_count)\ reduction(+: approx) for (i=1; i<=n-1; i++) { approx+=f(a+i*h); } approx=h*approx; h=(b-a)/n; approx=(f(a)+f(b))/2.0; for (i=1; i<=n-1; i++) { approx+=f(a+i*h); } approx=h*approx; • The code block must be a for loop.
• Iterations of the for loop are divided among threads.
• approx is a reduction variable.
• i is a private variable.
14
The parallel for Directive
• • • Sounds like a truly wonderful approach to parallelizing serial programs.
Does not work with while or do-while loops.
– How about converting them into for loops?
The number of iterations must be determined in advance.
for (; ;) { ...
} for (i=0; i } int x, y; # pragma omp parallel for num_threads(thread_count) for(x=0; x < width; x++) { for(y=0; y < height; y++) { finalImage[x][y] = f(x, y); } } private(y) 15 4 1 1 3 1 5 1 7 4 k 0 2 ( 1 ) k K 1 double factor=1.0; double sum=0.0; for(k=0; k Loop-carried dependence double factor=1.0; double sum=0.0; # pragma omp parallel for\ num_threads(thread_count)\ reduction(+: sum) for(k=0; k if(k%2 == 0) factor=1.0; else factor=-1.0; sum+=factor/(2*k+1); factor=(k%2 == 0)?1.0: -1.0; sum+=factor/(2*k+1); double factor=1.0; double sum=0.0; # pragma omp parallel for num_threads(thread_count)\ reduction(+: sum) private(factor) for(k=0; k double factor=1.0; double sum=0.0; # pragma omp parallel for num_threads(thread_count)\ default(none) reduction(+: sum) private(k, factor) shared(n) for(k=0; k } pi_approx=4.0*sum; • With the default (none) clause, we need to specify the scope of each variable that we use in the block that has been declared outside the block. • The value of a variable with private scope is unspecified at the beginning (and after completion) of a parallel or parallel for block. 18 for (len=n; len>=2; len--) for (i=0; i • Can we parallelize the outer loop? • Can we parallelize the inner loop? 19 Phase 0 1 2 3 0 9 7 7 7 7 6 6 6 Subscript in Array 1 7 9 9 6 6 7 7 7 2 8 6 6 9 9 8 8 8 Any opportunities for parallelism? 3 6 8 8 8 8 9 9 9 20 void Odd_even_sort (int a[], int n) { int phase, i, temp; for (phase=0; phase for (phase=0; phase # pragma omp parallel num_thread(thread_count) \ default(none) shared(a, n) private(i, tmp, phase) for (phase=0; phase Iterations 0 1 2 3 4 5 6 7 8 Block Threads 0 1 2 Iterations 0 1 2 3 4 5 6 7 8 Cyclic Threads 0 1 2 24 sum=0.0; for (i=0; i<=n; i++) sum+=f(i); double f(int i) { int j, start=i*(i+1)/2, finish=start+i; double return_val=0.0; } for (j=start; j<=finish; j++) { return_val+=sin(j); } return return_val; 25 sum=0.0; # pragma omp parallel for num_threads(thread_count) \ reduction(+:sum) schedule(static, 1) for (i=0; i n=12, t=3 schedule(static, 1) schedule(static, 2) schedule(static, 4) Thread 0 : 0, 3, 6, 9 Thread 1 : 1, 4, 7, 10 Thread 2 : 2, 5, 8, 11 Thread 0 : 0, 1, 6, 7 Thread 1 : 2, 3, 8, 9 Thread 2 : 4, 5, 10, 11 Thread 0 : 0, 1, 2, 3 Thread 1 : 4, 5, 6, 7 Thread 2 : 8, 9, 10, 11 schedule(static, total_iterations/thread_count) 26 • In a schedule: – Iterations are broken into chunks of chunksize consecutive iterations. – Default chunksize value: 1 – Each thread executes a chunk. – When a thread finishes a chunk, it requests another one. • In a schedule: – Each thread executes a chunk. – When a thread finishes a chunk, it requests another one. – As chunks are completed, the size of the new chunks decreases. – Approximately equals to the number of iterations remaining divided by the number of threads. – The size of chunks decreases down to chunksize or 1 (default). 27 • The optimal schedule depends on: – The type of problem – The number of iterations – The number of threads • Overhead – guided > dynamic > static – If you are getting satisfactory results (e.g., close to the theoretically maximum speedup) without a schedule clause, go no further. • The Cost of Iterations – If it is roughly the same, use the default schedule. – If it decreases or increases linearly as the loop executes, a static schedule with small chunksize values will be good. – If it cannot be determined in advance, try to explore different options. 28 A x y # pragma omp parallel for num_threads(thread_count) \ default(none) private(i,j) shared(A, x, y, m, n) for(i=1; i Number of Threads 1 2 3 8,000,000 x 8 Time Efficiency 0.322 1.000 0.219 0.141 0.735 0.571 Matrix Dimension 8,000 x 8,000 Time Efficiency 0.264 0.189 0.119 1.000 0.698 0.555 8 x 8,000,000 Time Efficiency 0.333 1.000 0.300 0.303 0.555 0.275 30 • 8,000,000-by-8 – y has 8,000,000 elements Potentially large number of write misses • 8-by-8,000,000 – x has 8,000,000 elements Potentially large number of read misses • 8-by-8,000,000 – y has 8 elements (8 doubles) Could be stored in the same cache line (64 bytes). – Potentially serious false sharing effect for multiple processors • 8000-by-8000 – y has 8,000 elements (8,000 doubles). – Thread 2: 4000 to 5999 Thread 3: 6000 to 7999 – {y[5996], y[5997], y[5998], y[5999], y[6000], y[6001], y[6002], y[6003] } – The effect of false sharing is highly unlikely. 31 • How to generate random numbers in C? – First, call srand() – Second, call rand() with an integer seed. to create a sequence of random numbers. • • Pseudorandom Number Generator (PRNG) X n 1 aX n c mod m Is it thread safe? – Can it be simultaneously executed by multiple threads without causing problems? 32 • Partitioning – Divide the computation and the data into small tasks. – Identify tasks that can be executed in parallel. • Communication – Determine what communication needs to be carried out. – Local Communication vs. Global Communication • Agglomeration – Group tasks into larger tasks. – Reduce communication. – Task Dependence • Mapping – Assign the composite tasks to processes/threads. 33 34 • To predict the motion of a group of objects that interact with each other gravitationally over a period of time. – Inputs: Mass, Position and Velocity • Astrophysicist – The positions and velocities of a collection of stars • Chemist – The positions and velocities of a collection of molecules 35 f qk s q Gm q m k s k 3 s q ( t ) s k ( t ) F q Gm q k k n 1 0 q s q ( t ) m k s k ( t ) 3 s q ( t ) s k ( t ) a q G k k n 1 0 q s q ( t ) m k s k ( t ) 3 s q ( t ) s k ( t ) 36 Get input data; for each timestep { if (timestep output) Print positions and velocities of particles; for each particle q Compute total force on q; for each particle q Compute position and velocity of q; } for each particle q { forces[q][0]=forces[q][1]=0; for each particle k!=q { x_diff=pos[q][0]-pos[k][0]; y_diff=pos[q][1]-pos[k][1]; dist=sqrt(x_diff*x_diff+y_diff*y_diff); dist_cubed=dist*dist*dist; forces[q][0]-=G*masses[q]*masses[k]/dist_cubed*x_diff; forces[q][1]-=G*masses[q]*masses[k]/dist_cubed*y_diff; } } 37 for each particle q forces[q][0]=forces[q][1]=0; for each particle q { for each particle k>q { x_diff=pos[q][0]-pos[k][0]; y_diff=pos[q][1]-pos[k][1]; dist=sqrt(x_diff*x_diff+y_diff*y_diff); dist_cubed=dist*dist*dist; force_qk[0]=-G*masses[q]*masses[k]/dist_cubed*x_diff; force_qk[1]=-G*masses[q]*masses[k]/dist_cubed*y_diff; forces[q][0]+=force_qk[0]; forces[q][1]+=force_qk[1]; forces[k][0]-=force_qk[0]; forces[k][1]-=force_qk[1]; } } 38 y ( t 0 t ) y ( t 0 ) y ( t 0 )( t 0 t t 0 ) y ( t 0 ) y ( t 0 ) t 39 s q ( t ) s q ( 0 ) ts q ' ( 0 ) s q ( 0 ) tv q ( 0 ) v q ( t ) v q ( 0 ) tv q ' ( 0 ) v q ( 0 ) ta q ( 0 ) v q ( 0 ) t F q ( 0 ) m q s q ( 2 t ) s q ( t ) ts q ' ( t ) s q ( t ) tv q ( t ) v q ( 2 t ) v q ( t ) tv q ' ( t ) v q ( t ) ta q ( t ) v q ( t ) t F q ( t ) m q for each particle q { pos[q][0]+=delta_t*vel[q][0]; pos[q][1]+=delta_t*vel[q][1]; vel[q][0]+=delta_t*forces[q][0]/masses[q]; vel[q][1]+=delta_t*forces[q][1]/masses[q]; } 40 s q (t) v q (t) F q (t) s r (t) s q (t + △ t) v q (t + △ t) v r (t) s r (t + △ t) v r (t + △ t) F r (t) F q (t+ △ t) F r (t+ △ t) 41 t s q s q , v q , F q s r s r , v r , F r t + △ t s q , v q , F q s q s r s r , v r , F r 42 # pragma omp parallel for each timestep { if (timestep output){ # pragma omp single nowait Print positions and velocities of particles; } # pragma omp for for each particle q Compute total force on q; # pragma omp for for each particle q Compute position and velocity of q; } 43 # pragma omp for for each particle q forces[q][0]=forces[q][1]=0; # pragma omp for for each particle q { for each particle k>q { x_diff=pos[q][0]-pos[k][0]; y_diff=pos[q][1]-pos[k][1]; dist=sqrt(x_diff*x_diff+y_diff*y_diff); dist_cubed=dist*dist*dist; force_qk[0]=-G*masses[q]*masses[k]/dist_cubed*x_diff; force_qk[1]=-G*masses[q]*masses[k]/dist_cubed*y_diff; forces[q][0]+=force_qk[0]; forces[q][1]+=force_qk[1]; forces[k][0]-=force_qk[0]; forces[k][1]-=force_qk[1]; } 44 } • Consider 2 threads and 4 particles. • Thread 1 is assigned particle 0 and particle 1. • Thread 2 is assigned particle 2 and particle 3. • F 3 =-f 03 -f 13 -f 23 • Who will calculate f 03 and f 13 ? • Who will calculate f 23 ? • Any race conditions? 45 Thread 0 1 2 Particle 0 1 2 3 4 5 0 f 01 +f 02 +f 03 +f 04 +f 05 -f 01 +f 12 +f 13 +f 14 +f 15 -f 02 -f 12 -f 03 -f 13 -f 04 -f 14 -f 05 -f 15 Thread 1 0 0 f 23 +f 24 +f 25 -f 23 +f 34 +f 35 -f 24 –f 34 -f 25 –f 35 2 0 0 0 0 f 45 -f 45 3 Threads, 6 Particles, Block Partition 46 # pragma omp for for each particle q { for each particle k>q { x_diff=pos[q][0]-pos[k][0]; y_diff=pos[q][1]-pos[k][1]; dist=sqrt(x_diff*x_diff+y_diff*y_diff); dist_cubed=dist*dist*dist; force_qk[0]=-G*masses[q]*masses[k]/dist_cubed*x_diff; force_qk[1]=-G*masses[q]*masses[k]/dist_cubed*y_diff; loc_forces[my_rank][q][0]+=force_qk[0]; loc_forces[my_rank][q][1]+=force_qk[1]; loc_forces[my_rank][k][0]-=force_qk[0]; loc_forces[my_rank][k][1]-=force_qk[1]; } } 47 # pragma omp for for (q=0; q • In the second phase, the thread that has been assigned particle the contributions that have been computed by different threads. q will add 48 • In the reduced code: – Loop 1: Initialization of the loc_forces array – Loop 2: The first phase of the computation of forces – Loop 3: The second phase of the computation of forces – Loop 4: The updating of positions and velocities • Which schedule should be used? Threads 1 2 4 8 Basic 7.71 3.87 1.95 0.99 Reduced Default 3.90 2.94 1.73 0.95 Reduced Forces Cyclic 3.90 1.98 1.01 0.54 Reduced All Cyclic 3.90 2.01 1.08 0.61 49 • What are the major differences between MPI and OpenMP? • What is the scope of a variable? • What is a reduction variable? • How to ensure mutual exclusion in a critical section? • What are the common loop scheduling options? • What factors may potentially affect the performance of an OpenMP program? • What is a thread safe function? 50Estimating π
Estimating π
Scope Matters
Bubble Sort
Odd-Even Sort
Odd-Even Sort
Odd-Even Sort in OpenMP
Odd-Even Sort in OpenMP
Data Partitioning
Scheduling Loops
The schedule clause
The dynamic and guided Types
dynamic
guided
Which schedule?
Performance Issue
Performance Issue
Performance Issue
Thread Safety
Foster’s Methodology
Foster’s Methodology
The n-body Problem
Newton’s Law
The Basic Algorithm
The Reduced Algorithm
Euler Method
Position and Velocity
Communications
Agglomeration
Parallelizing the Basic Solver
Parallelizing the Reduced Solver
Does it work properly?
Thread Contributions
First Phase
Second Phase
Evaluating the OpenMP Codes
Review