June 24, 2013 Jason Su

1. June 24, 2013 Jason Su

2. Technologies for C/C++/Fortran • Single machine, multi-core • P(OSIX) threads: bare metal multi-threading • OpenMP: compiler directives that implement various constructs like parallel-for • Single machine, GPU • CUDA/OpenCL: bare metal GPU coding • Thrust: algorithms library for CUDA inspired by C++ STL • Like MATLAB, where you use a set of fast core functions and data structures to implement your program • Multi-machine • Message Passing Interface (MPI): a language-independent communication protocol to coordinate and program a cluster of machines

3. Challenges of Parallel Programming • Correctness • Race conditions • Synchronization/deadlock • Floating point arithmetic is not associative or distributive • Debugging • How do you fix problems that are not reproducible? • How do you assess the interaction of many threads running simultaneously? • Performance • Management and algorithm overhead • Amdahl’s Law

4. Challenges: Race Conditions As a simple example, let us assume that two threads each want to increment the value of a global integer variable by one. What we expected: What we got:

5. Challenges: Deadlock A real world example would be an illogical statute passed by the Kansas legislature in the early 20th century, which stated: “When two trains approach each other at a crossing, both shall come to a full stop and neither shall start up again until the other has gone.”

6. Challenges: Deadlock • Consider a program that manages bank accounts: BankAccount: string owner float balance withdraw(float amount) deposit(float amount) transfer(self, Account to, float amount): lock(self) lock(to) from.withdraw(amount) to.deposit(amount) release(to) release(from)

7. Challenges: Race Conditions As a simple example, let us assume that two threads each want to increment the value of a global integer variable by one. What we expected: What we got:

8. Challenges: Floating Point Arithmetic What we expected: What we got: Say we had a floating point format that has 1 base-10 exponent and 1 decimal place: □e□

9. Challenges: Floating Point Arithmetic • Floating point integers have limited precision • Any operation implicitly rounds, which means that order matters • Float arithmetic is not associative or distributive • Do operations on values of similar precision first • Actually an issue in all code, but becomes apparent when serial and parallel answers differ “Many programmers steeped in sequential programming for so many years make the assumption that there is only one right answer for their algorithm. After all, their code has always delivered the same answer every time it was run…however, all the answers are equally correct”

10. Challenges: Floating Point Arithmetic • The answer changes significantly with the number of threads • Because the order of operations has changed • Indicates an unstable numerical algorithm • Parallel programming revealed a flaw in the algorithm which is otherwise unapparent with serial code

11. OpenMP • Uses a fork-join programming paradigm

13. OpenMP • Declare private and shared variables for a thread • Atomic operations • Ensure that code is uninterrupted by other threads, important for shared variables • Parallel for-loop: see MATLAB’s parfor • #pragma omp parallel for num_threads(nthreads)for (inti=0; i < N; i++) { ... } • Reduction: ℝn_threads→ℝ • #pragma omp parallel for reduction(+:val)for (uint i=0; i < v.size(); i++) val = v[i];

14. Moore’s Law • Serial scaling performance has reached its peak. • Processors are not getting faster, but wider

15. GPUs

16. GPUs

17. CPU vs GPU • GPU devotes more transistors to data processing

18. CPU vs GPU • GPU devotes more transistors to data processing

19. CPU vs GPU • CPU minimizes time to complete a given task: latency of a thread • GPU maximizes number of tasks in a fixed time: throughput of all threads • Which is better? It depends on the problem “If you were plowing a field, which would you rather use: two strong oxen or 1024 chicken?” Seymour Cray

20. NVIDIA CUDA

21. NVIDIA CUDA

22. NVIDIA GPU Architecture: Processors • GPU contains many streaming multiprocessors (MP or SM) • These do the work • Groups threads into “warps” and executes these in lock-step • Different warps are swapped in and out constantly to hide latency = “parallelism”

23. What is a warp? • A warp is a group of 32 threads that are executed in lock-step (SIMT) • Threads in a warp demand obedience, can either: be idle or do the same instruction as its siblings • Divergent code (if/else, loops, return, etc.) can wastefully leave some threads idle

24. NVIDIA GPU Architecture: Memory

25. NVIDIA GPU Architecture • Imagine a pool of threads representing all the work there is to do • For example, a thread for each pixel in the image • These threads need to be broken up into blocks that fit neatly onto each MP, this is chosen by you • There are several limitations that affect how well these blocks fill up an MP, called its occupancy • Total threads < 2048 • Total blocks < 16 • Total shared memory < 48KB • Total registers < 64K

26. CUDA Program Structure

27. CUDA Program Structure • Memory and processors on the GPU exist in big discrete chunks • Designing for GPU is not only developing a parallel algorithm • But also shaping it so that it optimally maps to GPU hardware structures

28. CUDA Program Structure

29. Finally, we get to the problem.

30. Connected Component Labeling • Say you had a segmented mask of lesions: • How would you count the total number of (connected) lesions in the image? • How would you track an individual lesion over time? 1 2

32. Connected labels share a similar property or state • In this case color, where we’ve simplified to black and white

33. Stars labeled 1 through 5467

34. Initialize with Raster-Scan Numbering Input Image Initial Labels • Goal is to label each connected region with a unique number, min(a, b) is used to combine labels. • That means 1, 7 , 9, and 54 should remain.

35. Kernel A – Label Propagation • Iterate until nothing changes: • Assign a pixel for each thread • For each thread: • Look at my neighboring pixels • If neighbor is smaller than me, re-label myself to the lowest adjacent label • K. Hawick, A. Leist, D. Playne, Parallel graph component labelling with GPUs • and CUDA, Parallel Computing 36 (12) (2010) 655–678.

36. Kernel A – Label Propagation • K. Hawick, A. Leist, D. Playne, Parallel graph component labelling with GPUs • and CUDA, Parallel Computing 36 (12) (2010) 655–678.

37. Kernel A – Label Propagation • K. Hawick, A. Leist, D. Playne, Parallel graph component labelling with GPUs • and CUDA, Parallel Computing 36 (12) (2010) 655–678.

38. Kernel A – Label Propagation • A label can only propagate itself by a maximum of one cell per iteration • Many iterations are required • Very slow for large clusters in the image • We’re getting killed by many accesses to global memory, each taking O(100) cycles • Even having many parallel threads is not enough to hide it, need 32 cores/MP*400 cycles = 12,800 threads active/MP • K. Hawick, A. Leist, D. Playne, Parallel graph component labelling with GPUs • and CUDA, Parallel Computing 36 (12) (2010) 655–678.

39. Kernel B – Local Label Propagation • Take advantage of L1-speed shared memory • Iterate until global convergence: • Assign a pixel for each thread • For each thread: • Load myself and neighbor labels into shared memory • Iterate until local convergence: • Look at my neighboring pixels • If neighbor is smaller than me, re-label myself to the lowest adjacent label • K. Hawick, A. Leist, D. Playne, Parallel graph component labelling with GPUs • and CUDA, Parallel Computing 36 (12) (2010) 655–678.

40. Kernel B – Local Label Propagation Initial Labels

41. Kernel B – Local Label Propagation • K. Hawick, A. Leist, D. Playne, Parallel graph component labelling with GPUs • and CUDA, Parallel Computing 36 (12) (2010) 655–678.

42. Kernel B – Local Label Propagation • K. Hawick, A. Leist, D. Playne, Parallel graph component labelling with GPUs • and CUDA, Parallel Computing 36 (12) (2010) 655–678.

43. Kernel B – Local Label Propagation • A label can only propagate itself by a maximum of one block per iteration • Iterations are reduced but still can be a lot • Tiling the data into shared memory blocks is a very common technique • Appears in many applications, including FDTD

44. Kernel D – Label Equivalence • Completely change the algorithm • How do we make it so labels can propagate extensively in 1 iteration? • Track which labels are equivalent and resolve equivalency chains Iterate until convergence:

45. Kernel D – Label Equivalence • K. Hawick, A. Leist, D. Playne, Parallel graph component labelling with GPUs • and CUDA, Parallel Computing 36 (12) (2010) 655–678.

46. CPU – MATLAB’s bwlabel

47. Performance

48. Results: CPU vs GPU

49. Figure Credits • Christopher Cooper, Boston University • Eric Darve, Stanford University • Gordon Erlebacher, Florida State University • CUDA C Computing Guide