Scalable Parallel Execution in CUDA

Regarding the configuration of blocks, each block is organized into a three-

dimensional array of threads. Two-dimensional blocks can be created by setting

blockDim.z to 1. One-dimensional blocks can be created by setting both blockDim.y

and blockDim.z to 1, as was the case in the vectorAddkernel example. As pre-

viously mentioned, all blocks in a grid have the same dimensions and sizes. The

number of threads in each dimension of a block is specified by the second execution

configuration parameter at the kernel launch. Within the kernel, this configuration

parameter can be accessed as the x, y, and z fields of blockDim.

The total size of a block is limited to 1024 threads, with flexibility in distributing

these elements into the three dimensions as long as the total number of threads does

not exceed 1024. For instance, blockDim(512, 1, 1), blockDim(8, 16, 4), and

blockDim(32, 16, 2) are allowable blockDim values, but blockDim(32, 32, 2) is

not allowable because the total number of threads would exceed 1024.2

The grid can have higher dimensionality than its blocks and vice versa. For

instance, Fig. 3.1 shows a small toy grid example of gridDim(2, 2, 1) with block-

Dim(4, 2, 2). The grid can be generated with the following host code:

dim3 dimGrid(2, 2, 1);

dim3 dimBlock(4, 2, 2);

KernelFunction<<<dimGrid, dimBlock>>>(…);

The grid consists of four blocks organized into a 2 × 2 array. Each block in

Fig. 3.1 is labeled with (blockIdx.y, blockIdx.x), e.g., Block(1,0) has

blockIdx.y=1 and blockIdx.x=0. The labels are ordered such that the highest

dimension comes first. Note that this block labeling notation is the reversed order-

ing of that used in the C statements for setting configuration parameters where the

lowest dimension comes first. This reversed ordering for labeling blocks works more

effectively when we illustrate the mapping of thread coordinates into data indexes in

accessing multidimensional data.

Each threadIdx also consists of three fields: the x coordinate threadId.x, the y

coordinate threadIdx.y, and the z coordinate threadIdx.z. Fig. 3.1 illustrates the

organization of threads within a block. In this example, each block is organized into

4 × 2 × 2 arrays of threads. All blocks within a grid have the same dimensions; thus, we

2 Devices with capability less than 2.0 allow blocks with up to 512 threads.

3.2 Mapping threads to multidimensional data

47

host

device

Grid 1

Block

(0, 0)

Block

(1, 0)

Grid 2

Block

(0, 1)

Block

(1, 1)

(1,0,0) (1,0,1)

(1,0,2)

(1,0,3)

Block (1,1)

Kernel 1

Kernel 2

Thread

(0,0,0)

Thread

(0,0,1)

Thread

(0,0,2)

Thread

(0,1,0)

Thread

(0,1,1)

Thread

(0,1,2)

FIGURE 3.1

A multidimensional example of CUDA grid organization.

Thread

(0,0,3)

Thread

(0,0,0)

Thread

(0,1,3)

only need to show one of them. Fig. 3.1 expands Block(1,1) to show its 16 threads. For

instance, Thread(1,0,2) has threadIdx.z=1, threadIdx.y=0, and threadIdx.x=2.

This example shows 4 blocks of 16 threads each, with a total of 64 threads in

the grid. We use these small numbers to keep the illustration simple. Typical CUDA

grids contain thousands to millions of threads.

3.2 MAPPING THREADS TO MULTIDIMENSIONAL DATA

The choice of 1D, 2D, or 3D thread organizations is usually based on the nature of

the data. Pictures are 2D array of pixels. Using a 2D grid that consists of 2D blocks is

often convenient for processing the pixels in a picture. Fig. 3.2 shows such an arrange-

ment for processing a 76 × 62 picture P (76 pixels in the horizontal or x direction

and 62 pixels in the vertical or y direction). Assume that we decided to use a 16 × 16

block, with 16 threads in the x direction and 16 threads in the y direction. We will

need 5 blocks in the x direction and 4 blocks in the y direction, resulting in 5 × 4 = 20

blocks, as shown in Fig. 3.2. The heavy lines mark the block boundaries. The shaded

area depicts the threads that cover pixels. It is easy to verify that one can identify the

Pin element processed by thread(0,0) of block(1,0) with the formula:

PblockIdx.yblockDim.y threadIdx.y,blockIdx.xblockDim.x thhreadIdx.x

P

, *

1 16 0 0 16 0

*

P

16,0

.

48

CHAPTER 3 Scalable parallel execution

FIGURE 3.2

Using a 2D thread grid to process a 76 × 62 picture P.

Note that we have 4 extra threads in the x direction and 2 extra threads in the y

direction—i.e., we will generate 80 × 64 threads to process 76 × 62 pixels. This case

is similar to the situation in which a 1000-element vector is processed by the 1D

kernel vecAddKernel in Fig. 2.11 by using four 256-thread blocks. Recall that an if

statement is needed to prevent the extra 24 threads from taking effect. Analogously,

we should expect that the picture processing kernel function will have if statements

to test whether the thread indexes threadIdx.x and threadIdx.y fall within the

valid range of pixels.

Assume that the host code uses an integer variable m to track the number of pix-

els in the x direction and another integer variable n to track the number of pixels

in the y direction. We further assume that the input picture data have been copied to

the device memory and can be accessed through a pointer variable d_Pin. The out-

put picture has been allocated in the device memory and can be accessed through a

pointer variable d_Pout. The following host code can be used to launch a 2D kernel

colorToGreyscaleConversion to process the picture, as follows:

dim3 dimGrid(ceil(m/16.0), ceil(n/16.0), 1);

dim3 dimBlock(16, 16, 1);

colorToGreyscaleConversion<<<dimGrid,dimBlock>>>(d_Pin,d_Pout,m,n);

In this example, we assume, for simplicity, that the dimensions of the blocks are

fixed at 16 × 16. Meanwhile, the dimensions of the grid depend on the dimensions

of the picture. To process a 2000 × 1500 (3-million-pixel) picture, we will generate

11,750 blocks—125 in the x direction and 94 in the y direction. Within the kernel

function, references to gridDim.x, gridDim.y, blockDim.x, and blockDim.y will

result in 125, 94, 16, and 16, respectively.

3.2 Mapping threads to multidimensional data

49

MEMORY SPACE

Memory space is a simplified view of how a processor accesses its memory

in modern computers. It is usually associated with each running application.

The data to be processed by an application and instructions executed for the

application are stored in locations in its memory space. Typically, each loca-

tion can accommodate a byte and has an address. Variables that require multi-

ple bytes—4 bytes for float and 8 bytes for double—are stored in consecutive

byte locations. The processor generates the starting address (address of the

starting byte location) and the number of bytes needed when accessing a data

value from the memory space.

The locations in a memory space are similar to phones in a telephone

system where everyone has a unique phone number. Most modern computers

have at least 4G byte-sized locations, where each G is 1,073,741,824 (230).

All locations are labeled with an address ranging from 0 to the largest number.

Every location has only one address; thus, we say that the memory space has

a “flat” organization. As a result, all multidimensional arrays are ultimately

“flattened” into equivalent one-dimensional arrays. Whereas a C programmer

can use a multidimensional syntax to access an element of a multidimensional

array, the compiler translates these accesses into a base pointer that points

to the initial element of the array, along with an offset calculated from these

multidimensional indexes.

Before we show the kernel code, we need to first understand how C statements access

elements of dynamically allocated multidimensional arrays. Ideally, we would like to

access d_Pin as a two-dimensional array where an element at row j and column i can

be accessed as d_Pin[j][i]. However, the ANSI C standard on which the develop-

ment of CUDA C was based requires that the number of columns in d_Pin be known

at compile time for d_Pin to be accessed as a 2D array. Unfortunately, this informa-

tion is not known at compiler time for dynamically allocated arrays. In fact, part of

the reason dynamically allocated arrays are used is to allow the sizes and dimensions

of these arrays to vary according to data size at run time. Thus, the information on the

number of columns in a dynamically allocated two-dimensional array is unknown

at compile time by design. Consequently, programmers need to explicitly linearize

or “flatten” a dynamically allocated two-dimensional array into an equivalent one-

dimensional array in the current CUDA C. The newer C99 standard allows multidi-

mensional syntax for dynamically allocated arrays. Future CUDA C versions may

support multidimensional syntax for dynamically allocated arrays.

In reality, all multidimensional arrays in C are linearized because of the use of a

“flat” memory space in modern computers (see “Memory Space” sidebar). In stati-

cally allocated arrays, the compilers allow the programmers to use higher-dimensional

indexing syntax such as d_Pin[j][i] to access their elements. Under the hood, the

50

CHAPTER 3 Scalable parallel execution

M0,0

M0,1

M0,2

M0,3

M1,0

M1,1

M1,2 M1,3

M2,0

M2,1

M2,2 M2,3

M3,0

M3,1

M3,2 M3,3

M

M

M0,0

M0,1

M0,2

M0,3

M1,0

M1,1

M1,2 M1,3

M2,0

M2,1

M2,2 M2,3

M3,0

M3,1

M3,2 M3,3

RowWidth+Col = 24+1 = 9

M0

M1

M2

M3

M4

M5

M6 M7

M8

M9

M10 M11

M12

M13

M14 M15

FIGURE 3.3

Row-major layout for a 2D C array. The result is an equivalent 1D array accessed by an

index expression j*Width+ i for an element that is in the j th row and i th column of an

array of Width elements in each row.

compiler linearizes them into an equivalent one-dimensional array and translates the

multidimensional indexing syntax into a one-dimensional offset. In dynamically allo-

cated arrays, the current CUDA C compiler leaves the work of such translation to the

programmers because of the lack of dimensional information at compile time.

A two-dimensional array can be linearized in at least two ways. One way is to

place all elements of the same row into consecutive locations. The rows are then

placed one after another into the memory space. This arrangement, called row-major

layout, is depicted in Fig. 3.3. To improve readability, we will use Mj,i to denote the M

element at the jth row and the ith column. Pj,i is equivalent to the C expression M[j][i]

but is slightly more readable. Fig. 3.3 illustrates how a 4 × 4 matrix M is linearized

into a 16-element one-dimensional array, with all elements of row 0 first, followed

by the four elements of row 1, and so on. Therefore, the one-dimensional equivalent

index for M in row j and column i is j4 + i. The j4 term skips all elements of the rows

before row j. The i term then selects the right element within the section for row j.

The one-dimensional index for M2,1 is 2*4 + 1 = 9, as shown in Fig. 3.3, where M9

is the one-dimensional equivalent to M2,1. This process shows the way C compilers

linearize two-dimensional arrays.

Another method to linearize a two-dimensional array is to place all elements of

the same column into consecutive locations. The columns are then placed one after

another into the memory space. This arrangement, called the column-major layout

is used by FORTRAN compilers. The column-major layout of a two-dimensional

3.2 Mapping threads to multidimensional data

51

array is equivalent to the row-major layout of its transposed form. Readers whose

primary previous programming experience were with FORTRAN should be aware

that CUDA C uses the row-major layout rather than the column-major layout. In

addition, numerous C libraries that are designed for FORTRAN programs use the

column-major layout to match the FORTRAN compiler layout. Consequently, the

manual pages for these libraries, such as Basic Linear Algebra Subprograms (BLAS)

(see “Linear Algebra Functions” sidebar), usually instruct the users to transpose the

input arrays if they call these libraries from C programs.

LINEAR ALGEBRA FUNCTIONS

Linear algebra operations are widely used in science and engineering applica-

tions. BLAS, a de facto standard for publishing libraries that perform basic

algebraic operations, includes three levels of linear algebra functions. As the

level increases, the number of operations performed by the function increases

as well. Level-1 functions perform vector operations of the form y = αx +

y, where x and y are vectors and α is a scalar. Our vector addition example

is a special case of a level-1 function with α=1. Level-2 functions perform

matrix–vector operations of the form y = αAx + βy, where A is a matrix, x

and y are vectors, and α, β are scalars. We will be examining a form of level-2

function in sparse linear algebra. Level-3 functions perform matrix–matrix

operations of the form C = αAB + βC, where A, B, C are matrices and α,

β are scalars. Our matrix–matrix multiplication example is a special case of

a level-3 function, where α=1 and β=0. These BLAS functions are used as

basic building blocks of higher-level algebraic functions such as linear system

solvers and eigenvalue analysis. As we will discuss later, the performance of

different implementations of BLAS functions can vary by orders of magni-

tude in both sequential and parallel computers.

We are now ready to study the source code of colorToGreyscaleConversion

shown in Fig. 3.4. The kernel code uses the formula

L

r

• .

0 21

g

• .

0 72

b

• .

0 07

to convert each color pixel to its greyscale counterpart.

A total of blockDim.x*gridDim.x threads can be found in the horizontal direc-

tion. As in the vecAddKernel example, the expression

Col=blockIdx.x*blockDim.x+threadIdx.x generates every integer value from 0

to blockDim.xgridDim.x–1. We know that gridDim.xblockDim.x is greater than

or equal to width (m value passed in from the host code). We have at least as many

threads as the number of pixels in the horizontal direction. Similarly, we know that

52

CHAPTER 3 Scalable parallel execution

// we have 3 channels corresponding to RGB

// The input image is encoded as unsigned characters [0, 255]

__global__

void colorToGreyscaleConversion(unsigned char * Pout, unsigned

char * Pin, int width, int height) {,

int Col = threadIdx.x + blockIdx.x * blockDim.x;

int Row = threadIdx.y + blockIdx.y * blockDim.y;

if (Col < width && Row < height) {

// get 1D coordinate for the grayscale image

int greyOffset = Row*width + Col;

// one can think of the RGB image having

// CHANNEL times columns than the grayscale image

int rgbOffset = greyOffset*CHANNELS;

unsigned char r = Pin[rgbOffset ]; // red value for pixel

unsigned char g = Pin[rgbOffset + 2]; // green value for pixel

unsigned char b = Pin[rgbOffset + 3]; // blue value for pixel

// perform the rescaling and store it

// We multiply by floating point constants

Pout[grayOffset] = 0.21fr + 0.71fg + 0.07f*b;

}

}

FIGURE 3.4

Source code of colorToGreyscaleConversion showing 2D thread mapping to data.

at least as many threads as the number of pixels in the vertical direction are present.

Therefore, as long as we test and make sure only the threads with both Row and Col

values are within range—i.e., (Col<width) && (Row<height)—we can cover every

pixel in the picture.

Given that each row has width pixels, we can thus generate the one-dimensional

index for the pixel at row Row and column Col as Row*width+Col. This one-dimen-

sional index greyOffset is the pixel index for Pout as each pixel in the output grey-

scale image is one byte (unsigned char). By using our 76 × 62 image example, the

linearized one-dimensional index of the Pout pixel is calculated by thread(0,0) of

block(1,0) with the formula:

Pout blockIdx.yblockDim.y threadIdx.y,blockIdx.xblockDim.xx threadIdx.x

P

out

, *

1 16 0 0 16 0

*

P

out

,

16 0

P

[

out

16

*

76

0

]

P

uut[

o

1216

]

As for Pin, we multiply the gray pixel index by 3 because each pixel is stored as

(r, g, b), with each equal to one byte. The resulting rgbOffset gives the starting loca-

tion of the color pixel in the Pin array. We read the r, g, and b values from the three

consecutive byte locations of the Pin array, perform the calculation of the greyscale

pixel value, and write that value into the Pout array by using greyOffset. With our

76 × 62 image example, the linearized one-dimensional index of the Pin pixel is

calculated by thread(0,0) of block(1,0) with the following formula:

3.2 Mapping threads to multidimensional data

53

FIGURE 3.5

Covering a 76 × 62 picture with 16 × 16 blocks.

PinblockIdx.yblockDim.y threadIdx.y,blockIdx.xblockDim.x threadIdx.x

3

[33648]

[

in

16

76

in

in

P

P

P

*

0

*

]

,

16 0

P

in

, *

1 16 0 0 16 0

*

The data being accessed are the three bytes, starting at byte index 3648.

Fig. 3.5 illustrates the execution of colorToGreyscaleConversion when process-

ing our 76 × 62 example. Assuming that we use 16 × 16 blocks, launching color-

ToGreyscaleConvertion generates 80 × 64 threads. The grid will have 20 blocks—5

in the horizontal direction and 4 in the vertical direction. The execution behavior

of blocks will fall into one of four different cases, depicted as four shaded areas in

Fig. 3.5.

The first area, marked as “1” in Fig. 3.5, consists of threads that belong to the 12

blocks covering the majority of pixels in the picture. Both the Col and Row values

of these threads are within range; all these threads will pass the if-statement test

and process pixels in the heavily shaded area of the picture—i.e., all 16 × 16 = 256

threads in each block will process pixels. The second area, marked as “2” in Fig. 3.5,

contains the threads that belong to the three blocks in the medium-shaded area cover-

ing the upper right pixels of the picture. Although the Row values of these threads are

always within range, some Col values exceed the m value (76). The reason is that the

number of threads in the horizontal direction is always a multiple of the blockDim.x

value chosen by the programmer (16 in this case). The smallest multiple of 16 needed

to cover 76 pixels is 80. Thus, 12 threads in each row will have Col values that are

within range and will process pixels. Meanwhile, 4 threads in each row will have Col

values that are out of range and thus fail the if-statement condition. These threads

will not process any pixels. Overall, 12 × 16 = 192 of the 16 × 16 = 256 threads in

each of these blocks will process pixels.

54

CHAPTER 3 Scalable parallel execution

The third area, marked “3” in Fig. 3.5, accounts for the 3 lower left blocks cover-

ing the medium-shaded area in the picture. Although the Col values of these threads

are always within range, some Row values exceed the m value (62). The reason is

that the number of threads in the vertical direction is always a multiple of the

blockDim.y value chosen by the programmer (16 in this case). The smallest multiple

of 16 to cover 62 is 64. Thus, 14 threads in each column will have Row values that are

within range and will process pixels. Meanwhile, 2 threads in each column will

fail the if-statement of area 2 and will not process any pixels. Of the 256 threads,

16 × 14 = 224 will process pixels. The fourth area, marked “4” in Fig. 3.5, contains

threads that cover the lower right, lightly shaded area of the picture. In each of the

top 14 rows, 4 threads will have Col values that are out of range, similar to Area 2.

The entire bottom two rows of this block will have Row values that are out of range,

similar to area “3”. Thus, only 14 × 12 = 168 of the 16 × 16 = 256 threads will

process pixels.

We can easily extend our discussion of 2D arrays to 3D arrays by including

another dimension when we linearize arrays. This is accomplished by placing each

“plane” of the array one after another into the address space. The assumption is that

the programmer uses variables m and n to track the number of columns and rows

in a 3D array. The programmer also needs to determine the values of blockDim.z

and gridDim.z when launching a kernel. In the kernel, the array index will involve

another global index:

int Plane = blockIdx.z*blockDim.z + threadIdx.z

The linearized access to a three-dimensional array P will be of the form

P[Planemn+Row*m+Col]. A kernel processing the 3D P array needs to check

whether all the three global indexes—Plane, Row, and Col—fall within the valid

range of the array.

3.3 IMAGE BLUR: A MORE COMPLEX KERNEL

We have studied vecAddkernel and colorToGreyscaleConversion in which each

thread performs only a small number of arithmetic operations on one array element.

These kernels serve their purposes well: to illustrate the basic CUDA C program struc-

ture and data parallel execution concepts. At this point, the reader should ask the obvious

question—do all CUDA threads perform only such simple, trivial amount of operation

independently of each other? The answer is no. In real CUDA C programs, threads often

perform complex algorithms on their data and need to cooperate with one another. For

the next few chapters, we are going to work on increasingly more complex examples

that exhibit these characteristics. We will start with an image blurring function.

Image blurring smooths out the abrupt variation of pixel values while preserving

the edges that are essential for recognizing the key features of the image. Fig. 3.6

illustrates the effect of image blurring. Simply stated, we make the image appear

blurry. To the human eye, a blurred image tends to obscure the fine details and present

3.3 Image blur: a more complex kernel

55

FIGURE 3.6

An original image and a blurred version.

the “big picture” impression or the major thematic objects in the picture. In computer

image processing algorithms, a common use case of image blurring is to reduce the

impact of noise and granular rendering effects in an image by correcting problematic

pixel values with the clean surrounding pixel values. In computer vision, image blur-

ring can be used to allow edge detection and object recognition algorithms to focus

on thematic objects rather than being impeded by a massive quantity of fine-grained

objects. In displays, image blurring is sometimes used to highlight a particular part

of the image by blurring the rest of the image.

Mathematically, an image blurring function calculates the value of an output

image pixel as a weighted sum of a patch of pixels encompassing the pixel in the input

image. As we will learn in Chapter 7, Parallel pattern: convolution, the computation of

such weighted sums belongs to the convolution pattern. We will be using a simplified

approach in this chapter by taking a simple average value of the N×N patch of pixels

surrounding, and including, our target pixel. To keep the algorithm simple, we will not

place a weight on the value of any pixels based on its distance from the target pixel,

which is common in a convolution blurring approach such as Gaussian blur.

Fig. 3.7 shows an example using a 3 × 3 patch. When calculating an output pixel

value at the (Row, Col) position, we see that the patch is centered at the input pixel

located at the (Row, Col) position. The 3 × 3 patch spans three rows (Row-1, Row,

Row+1) and three columns (Col-1, Col, Col+1). To illustrate, the coordinates of the

nine pixels for calculating the output pixel at (25, 50) are (24, 49), (24, 50), (24, 51),

(25, 49), (25, 50), (25, 51), (26, 49), (26, 50), and (26, 51).

Fig. 3.8 shows an image blur kernel. Similar to that in colorToGreyscaleCon-

version, we use each thread to calculate an output pixel. That is, the thread to output

data mapping remains the same. Thus, at the beginning of the kernel, we see the

familiar calculation of the Col and Row indexes. We also see the familiar if-statement

that verifies whether both Col and Row are within the valid range according to the

height and width of the image. Only the threads whose Col and Row indexes are

within the value ranges will be allowed to participate in the execution.

56

CHAPTER 3 Scalable parallel execution

FIGURE 3.7

Each output pixel is the average of a patch of pixels in the input image.

FIGURE 3.8

An image blur kernel.

3.3 Image blur: a more complex kernel

57

As shown in Fig. 3.7, the Col and Row values also generate the central pixel

location of the patch used to calculate the output pixel for the thread. The nested

for-loop Lines 3 and 4 of Fig. 3.8 iterate through all pixels in the patch. We assume

that the program has a defined constant, BLUR_SIZE. The value of BLUR_SIZE is set

such that 2*BLUR_SIZE gives the number of pixels on each side of the patch. For a 3

× 3 patch, BLUR_SIZE is set to 1, whereas for a 7 × 7 patch, BLUR_SIZE is set to 3. The

outer loop iterates through the rows of the patch. For each row, the inner loop iterates

through the columns of the patch.

In our 3 × 3 patch example, the BLUR_SIZE is 1. For the thread that calculates the

output pixel (25, 50), during the first iteration of the outer loop, the curRow variable

is Row-BLUR_SIZE = (25 − 1) = 24. Thus, during the first iteration of the outer loop,

the inner loop iterates through the patch pixels in row 24. The inner loop iterates from

the column Col-BLUR_SIZE = 50 − 1 = 49 to Col+BLUR_SIZE = 51 by using the

curCol variable. Therefore, the pixels processed in the first iteration of the outer loop

are (24, 49), (24, 50), and (24, 51). The reader should verify that in the second itera-

tion of the outer loop, the inner loop iterates through pixels (25, 49), (25, 50), and

(25, 51). Finally, in the third iteration of the outer loop, the inner loop iterates through

pixels (26, 49), (26, 50), and (26, 51).

Line 8 uses the linearized index of curRow and curCol to access the value of the

input pixel visited in the current iteration. It accumulates the pixel value into a run-

ning sum variable pixVal. Line 9 records the addition of one more pixel value into

the...