Question

A more challenging example is to use CUDA to sum a vector. This is a reducction and requires communicaiton across threads. In the CUDA model, only threads within a block can share state efficiently by using shared memoery as writing to global memory would be disastrously slow. Therefore, we have to do this in stages - if the shared memory size is \(k\) numbers, we will need \(n\) stages to sum \(k^n\) numbers on the GPU.

Coding strategy

• For simplicity, we set up a reduction that only requires 2 stages
• We have an array \(a\) of length \(n\) that we wish to sum
• We consider \(a\) as being made up of a number blocks of size \(n/k\)
• The CPU will launch a kernel to find the \(k\) partial sums of \(a\)
  • Each sum will be of consecutive \(n/k\) elements in the original array
  • The summation of pairs of numbers is performed by a device-only sub-kernel launched by the GPU
  • The summation is done in 3 steps:
    • Each thread in a block writes its values to shared memory in the location corresponding to the thread index
    • Synchronize threads to make sure that all threads have completed writing before proceeding
    • The first thread in the block sums up the values in shared memory (the rest are idle) and stores in the location corresponding to the block index
• Finally, the CPU launches the kernel again to sum the partial sums
• For efficiency, we overwrite partial sums in the original vector

Note that other reductions (e.g. min, max) etc follow the same strategy - just swap the device kernel with another one.

The two strateiges of mapping each operation to a thread and reduction to combine results from several threads are the basic buiding blocks of many CUDA algorithms. It is surprising how many alogrithms can be formulated as combinaitons of mapping and redcution steps - and we will revisit this pattern with Hadoop/SPARK.

@cuda.jit('int32(int32, int32)', device=True)
def dev_sum(a, b):
    return a + b

@cuda.jit('void(int32[:], int32[:])')
def cu_sum(a, b):
    "Simple implementation of reduction kernel"
    # Allocate static shared memory of 512 (max number of threads per block for CC < 3.0)
    # This limits the maximum block size to 512.
    sa = cuda.shared.array(shape=(612,), dtype=int32)
    tx = cuda.threadIdx.x
    bx = cuda.blockIdx.x
    bw = cuda.blockDim.x
    i = tx + bx * bw
    if i < a.shape[0]:
        sa[tx] = a[i]
        if tx == 0:
            # Uses the first thread of each block to perform the actual
            # reduction
            s = sa[tx]
            cuda.syncthreads()
            for j in range(1, bw):
                s = dev_sum(s, sa[j])
            b[bx] = s

k = 4 # numbers to be added in the partial sum (must be less than or equal to 512)
n = 6*4 # total length of vector to be summed

a = np.random.randint(0, n, n).astype(np.int32)

print 'a =', a
print 'a.sum() = ', a.sum()

d_a = cuda.to_device(a)
d_b = cuda.to_device(b, copy=False)

griddim = (k, 1)
blockdim = (a.size//k, 1)

cu_sum[griddim, blockdim](d_a, d_a)

d_a.to_host()

print 'a =', a

cu_sum[1, griddim](d_a[:k], d_a)

d_a.to_host()

print 'a =', a

a = [14 20 16 11 12  4  2  3  3 16  8 12 12 19  6 15 20 20 20 21  8  5 11 17]
a.sum() =  295
a = [77 44 92 82 12  4  2  3  3 16  8 12 12 19  6 15 20 20 20 21  8  5 11 17]
a = [295  44  92  82  12   4   2   3   3  16   8  12  12  19   6  15  20  20
  20  21   8   5  11  17]