CUDA: Introduction

Christian Trefftz / Greg Wolffe Grand Valley State University Supercomputing 2008 Education Program

Terms
 What is GPGPU?  General-Purpose computing on a Graphics GeneralProcessing Unit  Using graphic hardware for non-graphic noncomputations  What is CUDA?  Compute Unified Device Architecture  Software architecture for managing data-parallel dataprogramming
Supercomputing 2008 Education Program 2

Motivation

Supercomputing 2008 Education Program

CPU vs. GPU



CPU
  

Fast caches Branching adaptability High performance Multiple ALUs Fast onboard memory High throughput on parallel tasks
Executes program on each fragment/vertex

GPU
  

CPUs are great for task parallelism  GPUs are great for data parallelism

Supercomputing 2008 Education Program 4

CPU vs. GPU - Hardware

More transistors devoted to data processing

Supercomputing 2008 Education Program 5

Traditional Graphics Pipeline

Vertex processing Rasterizer Fragment processing Renderer (textures)
Supercomputing 2008 Education Program 6

Pixel / Thread Processing

Supercomputing 2008 Education Program

GPU Architecture

Supercomputing 2008 Education Program

Processing Element

Processing element = thread processor = ALU

Supercomputing 2008 Education Program 9

Memory Architecture
  

Constant Memory Texture Memory Device Memory

Supercomputing 2008 Education Program

10

DataData-parallel Programming
 Think

of the CPU as a massively-threaded massivelycoco-processor  Write kernel functions that execute on the device -- processing multiple data elements in parallel
 Keep

it busy! > massive threading  Keep your data close! > local memory
Supercomputing 2008 Education Program 11

Hardware Requirements
   

CUDACUDA-capable video card Power supply Cooling PCIPCI-Express

Supercomputing 2008 Education Program

12

Supercomputing 2008 Education Program

13

Acknowledgements
 NVidia

Corporation

[Link]/CUDA
 NVidia

Technical Brief Architecture Overview CUDA Programming Guide

 ACM


Queue
Supercomputing 2008 Education Program 14

[Link]

A Gentle Introduction to CUDA Programming

Supercomputing 2008 Education Program

15

Credits
 The

code used in this presentation is based on code available in:

   

the Tutorial on CUDA in Dr. Dobbs Journal Andrew Bellenirs code for matrix multiplication Igor Majdandzics code for Voronoi diagrams NVIDIAs CUDA programming guide

Supercomputing 2008 Education Program

16

Software Requirements/Tools
 CUDA


device driver  CUDA Software Development Kit

Emulator
 CUDA

Toolkit

 Occupancy

calculator  Visual profiler

Supercomputing 2008 Education Program

17

To compute, we need to:

   

Allocate memory that will be used for the computation (variable declaration and allocation) Read the data that we will compute on (input) Specify the computation that will be performed Write to the appropriate device the results (output)

Supercomputing 2008 Education Program

18

A GPU is a specialized computer

 

 

We need to allocate space in the video cards memory for the variables. The video card does not have I/O devices, hence we need to copy the input data from the memory in the host computer into the memory in the video card, using the variable allocated in the previous step. We need to specify code to execute. execute. Copy the results back to the memory in the host computer.
Supercomputing 2008 Education Program 19

Initially:

array Hosts Memory GPU Cards Memory

Supercomputing 2008 Education Program

20

Allocate Memory in the GPU card

array Hosts Memory

array_d GPU Cards Memory

Supercomputing 2008 Education Program

21

Copy content from the hosts memory to the GPU card memory

array Hosts Memory

array_d GPU Cards Memory

Supercomputing 2008 Education Program

22

Execute code on the GPU

GPU MPs

array Hosts Memory

array_d GPU Cards Memory

Supercomputing 2008 Education Program

23

Copy results back to the host memory

array Hosts Memory

array_d GPU Cards Memory

Supercomputing 2008 Education Program

24

The Kernel


It is necessary to write the code that will be executed in the stream processors in the GPU card That code, called the kernel, will be downloaded and executed, simultaneously and in lock-step fashion, in several lock(all?) stream processors in the GPU card How is every instance of the kernel going to know which piece of data it is working on?

Supercomputing 2008 Education Program

25

Grid Size and Block Size

 Programmers


need to specify:

The grid size: The size and shape of the data that the program will be working on The block size: The block size indicates the subsub-area of the original grid that will be assigned to an MP (a set of stream processors that share local memory)

Supercomputing 2008 Education Program

26

Block Size
 Recall

that the stream processors of the GPU are organized as MPs (multi(multiprocessors) and every MP has its own set of resources:
 

Registers Local memory

 The

block size needs to be chosen such that there are enough resources in an MP to execute a block at a time.
Supercomputing 2008 Education Program 27

In the GPU:
Processing Elements

Array Elements
Block 0 Block 1
28

Supercomputing 2008 Education Program

Lets look at a very simple example

 The
 

code has been divided into two files:

simple.c [Link]

 simple.c

is ordinary code in C  It allocates an array of integers, initializes it to values corresponding to the indices in the array and prints the array.  It calls a function that modifies the array  The array is printed again.
Supercomputing 2008 Education Program 29

simple.c


#include <stdio.h> #define SIZEOFARRAY 64 extern void fillArray(int *a,int size); /* The main program */ int main(int argc,char *argv[]) { /* Declare the array that will be modified by the GPU */ int a[SIZEOFARRAY]; int i; /* Initialize the array to 0s */ for(i=0;i < SIZEOFARRAY;i++) { a[i]=i; } /* Print the initial array */ printf("Initial state of the array:\n"); array:\ for(i = 0;i < SIZEOFARRAY;i++) { printf("%d ",a[i]); } printf("\ printf("\n"); /* Call the function that will in turn call the function in the GPU that will fill the array */ fillArray(a,SIZEOFARRAY); /* Now print the array after calling fillArray */ printf("Final state of the array:\n"); array:\ for(i = 0;i < SIZEOFARRAY;i++) { printf("%d ",a[i]); } printf("\ printf("\n"); return 0; }

Supercomputing 2008 Education Program

30

[Link]
 [Link]


contains two functions

fillArray(): A function that will be executed on the host and which takes care of:
Allocating variables in the global GPU memory Copying the array from the host to the GPU memory Setting the grid and block sizes Invoking the kernel that is executed on the GPU Copying the values back to the host memory Freeing the GPU memory

Supercomputing 2008 Education Program

31

fillArray (part 1)
#define BLOCK_SIZE 32 extern "C" void fillArray(int *array,int arraySize){ /* a_d is the GPU counterpart of the array that exists on the host memory */ int *array_d; cudaError_t result; /* allocate memory on device */ /* cudaMalloc allocates space in the memory of the GPU card */ result = cudaMalloc((void**)&array_d,sizeof(int)*arraySize); /* copy the array into the variable array_d in the device */ /* The memory from the host is being copied to the corresponding variable in the GPU global memory */ result = cudaMemcpy(array_d,array,sizeof(int)*arraySize, cudaMemcpyHostToDevice);

Supercomputing 2008 Education Program

32

fillArray (part 2)
/* execution configuration... */ /* Indicate the dimension of the block */ dim3 dimblock(BLOCK_SIZE); /* Indicate the dimension of the grid in blocks */ dim3 dimgrid(arraySize/BLOCK_SIZE); /* actual computation: Call the kernel, the function that is */ /* executed by each and every processing element on the GPU card */ cu_fillArray<<<dimgrid,dimblock>>>(array_d); /* read results back: */ /* Copy the results from the GPU back to the memory on the host */ result = cudaMemcpy(array,array_d,sizeof(int)*arraySize,cudaMemcpyDevice ToHost); /* Release the memory on the GPU card */ cudaFree(array_d); }

Supercomputing 2008 Education Program

33

[Link] (cont.)
 The


other function in [Link] is

cu_fillArray()
This is the kernel that will be executed in every stream processor in the GPU It is identified as a kernel by the use of the keyword: __global__ This function uses the built-in variables built 

blockIdx.x and threadIdx.x

to identify a particular position in the array

Supercomputing 2008 Education Program

34

cu_fillArray
__global__ void cu_fillArray(int *array_d){ int x; /* blockIdx.x is a built-in variable in CUDA builtthat returns the blockId in the x axis of the block that is executing this block of code threadIdx.x is another built-in variable in CUDA builtthat returns the threadId in the x axis of the thread that is being executed by this stream processor in this particular block */ x=blockIdx.x*BLOCK_SIZE+threadIdx.x; array_d[x]+=array_d[x]; }

Supercomputing 2008 Education Program

35

To compile:
 nvcc

simple.c [Link] o simple  The compiler generates the code for both the host and the GPU  Demo on [Link]

Supercomputing 2008 Education Program

36

What are those blockIds and threadIds?

 With

a minor modification to the code, we can print the blockIds and threadIds  We will use two arrays instead of just one.
 

One for the blockIds One for the threadIds

 The

code in the kernel:

x=blockIdx.x*BLOCK_SIZE+threadIdx.x; block_d[x] = blockIdx.x; thread_d[x] = threadIdx.x;

Supercomputing 2008 Education Program

37

In the GPU:
Processing Elements
Thread 0 Thread 1 Thread 2 Thread 3 Thread 0 Thread 1 Thread 2 Thread 3

Array Elements
Block 0 Block 1
38

Supercomputing 2008 Education Program

HandsHands-on Activity


Compile with (one single line) nvcc blockAndThread.c -o blockAndThread Run the program ./blockAndThread

[Link]

Edit the file [Link]  Modify the constant BLOCK_SIZE. The current value is 8, try replacing it with 4.  Recompile as above  Run the program and compare the output with the previous run.

Supercomputing 2008 Education Program 39

This can be extended to 2 dimensions



See files:
 

blockAndThread2D.c [Link]

The gist in the kernel

x = blockIdx.x*BLOCK_SIZE+threadIdx.x; y = blockIdx.y*BLOCK_SIZE+threadIdx.y; pos = x*sizeOfArray+y; block_dX[pos] = blockIdx.x;

Compile and run blockAndThread2D



nvcc blockAndThread2D.c [Link] -o blockAndThread2D ./blockAndThread2D

Supercomputing 2008 Education Program 40

When the kernel is called:

dim3 dimblock(BLOCK_SIZE,BLOCK_SIZE); nBlocks = arraySize/BLOCK_SIZE; dim3 dimgrid(nBlocks,nBlocks); cu_fillArray<<<dimgrid,dimblock>>> ( params);

Supercomputing 2008 Education Program

41

Another Example: saxpy

 SAXPY


(Scalar Alpha X Plus Y) loop iteration thread

A common operation in linear algebra

 CUDA:

Supercomputing 2008 Education Program

42

Traditional Sequential Code

void saxpy_serial(int n, float alpha, float *x, float *y) { for(int i = 0;i < n;i++) y[i] = alpha*x[i] + y[i]; }
Supercomputing 2008 Education Program 43

CUDA Code
__global__ void saxpy_parallel(int n, float alpha, float *x, float *y) { int i = blockIdx.x*blockDim.x+threadIdx.x; if (i<n) y[i] = alpha*x[i] + y[i]; }

Supercomputing 2008 Education Program

44

Keeping Multiprocessors in mind

Each hardware multiprocessor has the ability to actively process multiple blocks at one time.  How many depends on the number of registers per thread and how much shared memory per block is required by a given kernel.  The blocks that are processed by one multiprocessor at one time are referred to as active.  If a block is too large, then it will not fit into the resources of an MP.


Supercomputing 2008 Education Program

45

Warps


 

Each active block is split into SIMD ("Single Instruction Multiple Data") groups of threads called "warps". Each warp contains the same number of threads, called the "warp size", which are executed by the multiprocessor in a SIMD fashion. On if statements, or while statements (control transfer) the threads may diverge. Use: __syncthreads()

Supercomputing 2008 Education Program

46

A Real Application


The Voronoi Diagram: A fundamental data structure in Computational Geometry

Supercomputing 2008 Education Program

47

Definition
 Definition

: Let S be a set of n sites in Euclidean space of dimension d. For each site p of S, the Voronoi cell V(p) of p is the set of points that are closer to p than to other sites of S. The Voronoi diagram V(S) is the space partition induced by Voronoi cells.

Supercomputing 2008 Education Program

48

Algorithms
 The

classical sequential algorithm has complexity O(n log n) where n is the number of sites (seeds).  If one only needs an approximation, on a grid of points (e.g. digital display):
 

Assign a different color to each seed Calculate the distance from every point in the grid to all seeds Color each point with the color of its closest seed
Supercomputing 2008 Education Program 49

Lends itself to implementation on a GPU

 The

calculation for every pixel is a good candidate to be carried out in parallel  Notice that the locations of the seeds are readread-only in the kernel  Thus we can use the texture map area in the GPU card, which is a fast read-only readcache to store the seeds: __device__ __constant__
Supercomputing 2008 Education Program 50

Demo on cuda

Supercomputing 2008 Education Program

51

Tips for improving performance

 Special

thanks to Igor Majdandzic.

Supercomputing 2008 Education Program

52

Memory Alignment
 Memory

access on the GPU works much better if the data items are aligned at 64 byte boundaries.  Hence, allocating 2D arrays so that every row starts at a 64-byte boundary address 64will improve performance.  But that is difficult to do for a programmer

Supercomputing 2008 Education Program

53

Allocating 2D arrays with pitch

 CUDA


offers special versions of:

Memory allocation of 2D arrays so that every row is padded (if necessary). The function determines the best pitch and returns it to the program. The function name is cudaMallocPitch() Memory copy operations that take into account the pitch that was chosen by the memory allocation operation. The function name is cudaMemcpy2D()

Supercomputing 2008 Education Program

54

Pitch
Columns Padding

Rows

Pitch
Supercomputing 2008 Education Program 55

A simple example:
 See

[Link]  A matrix of 30 rows and 10 columns  The work is divided into 3 blocks of 10 rows:
 

Block size is 10 Grid size is 3

Supercomputing 2008 Education Program

56

Key portions of the code (1)

result = cudaMallocPitch( (void **)&devPtr, &pitch, width*sizeof(int), height);

Supercomputing 2008 Education Program

57

Key portions of the code (2)

result = cudaMemcpy2D( devPtr, pitch, mat, width*sizeof(int), width*sizeof(int), height, cudaMemcpyHostToDevice);
Supercomputing 2008 Education Program 58

In the kernel:
__global__ void myKernel(int int int int *devPtr, pitch, width, height)

{ int c; int thisRow; thisRow = blockIdx.x * 10 + threadIdx.x; int *row = (int *)((char *)devPtr + thisRow*pitch); for(c = 0;c < width;c++) row[c] = row[c] + 1; }
\

Supercomputing 2008 Education Program

59

The call to the kernel

myKernel<<<3,10>>>( devPtr, pitch, width, height);

Supercomputing 2008 Education Program

60

pitch
 Notice

Divide work by rows

that when using pitch, we divide the work by rows.  Instead of using the 2D decomposition of 2D blocks, we are dividing the 2D matrix into blocks of rows.

Supercomputing 2008 Education Program

61

Dividing the work by blocks:

Columns Block 0

Rows

Block 1

Block 2

Pitch
Supercomputing 2008 Education Program 62

An application that uses pitch: Mandelbrot



The Mandelbrot set: A set of points in the complex plane, the boundary of which forms a fractal. A complex number, c, is in the Mandelbrot set if, when starting with x0=0 and applying the iteration
xn+1 = xn2 + c

repeatedly, the absolute value of xn never exceeds a certain number (that number depends on c) however large n gets.

Supercomputing 2008 Education Program

63

Demo: Comparison
 We
 

can compare the execution times of:

The sequential version The CUDA version

Supercomputing 2008 Education Program

64

Performance Tip: Block Size

   

 

Critical for performance Recommended value is 192 or 256 Maximum value is 512 Should be a multiple of 32 since this is the warp size for Series 8 GPUs and thus the native execution size for multiprocessors Limited by number of registers on the MP Series 8 GPU MPs have 8192 registers which are shared between all the threads on an MP
Supercomputing 2008 Education Program 65

Performance Tip: Grid Size

    

Critical for scalability Recommended value is at least 100, but 1000 would scale for many generations of hardware Actual value depends on size of the problem data It should be a multiple of the number of MPs for an even distribution of work (not a requirement though) Example: 24 blocks


Grid will work efficiently on Series 8 (12 MPs), but it will waste resources on new GPUs with 32MPs

Supercomputing 2008 Education Program

66

Performance Tip: Code Divergance

    

Control flow instructions diverge (threads take different paths of execution) Example: if, for, while Diverged code prevents SIMD execution it forces serial execution (kills efficiency) One approach is to invoke a simpler kernel multiple times Liberal use of __syncthreads()

Supercomputing 2008 Education Program

67

Performance Tip: Memory Latency

    

4 clock cycles for each memory read/write plus additional 400-600 cycles for latency 400Memory latency can be hidden by keeping a large number of threads busy Keep number of threads per block (block size) and number of blocks per grid (grid size) as large as possible Constant memory can be used for constant data (variables that do not change). Constant memory is cached.

Supercomputing 2008 Education Program

68

Performance Tip: Memory Reads

  

Device is capable of reading a 32, 64 or 128-bit 128number from memory with a single instruction Data has to be aligned in memory (this can be accomplished by using cudaMallocPitch() calls) If formatted properly, multiple threads from a warp can each receive a piece of memory with a single read instruction

Supercomputing 2008 Education Program

69

Watchdog timer


  

Operating system GUI may have a "watchdog" timer that causes programs using the primary graphics adapter to time out if they run longer than the maximum allowed time. Individual GPU program launches are limited to a run time of less than this maximum. Exceeding this time limit usually causes a launch failure. Possible solution: run CUDA on a GPU that is NOT attached to a display.
Supercomputing 2008 Education Program 70

Resources on line
    

[Link] Content&pa=showpage&pid=532 [Link] [Link] [Link] Computation of Voronoi diagrams using a graphics processing unit by Igor Majdandzic et al. available through IEEE Digital Library, DOI: 10.1109/EIT.2008.4554342
Supercomputing 2008 Education Program 71