When using the CUDA runtime, a single-threaded application can drive
multiple GPUs by calling cudaSetDevice() to specify which GPU will be
operated by the calling CPU thread. This idiom is used in Listing 9-1 to
switch between the source and destination GPUs during the peer-to-peer
memcpy, as well as the single-threaded, multi-GPU implementation of
N-body described in Section 9.5.2.
In the driver API, CUDA maintains a stack of current contexts, so that subroutines can easily change and restore the caller’s current context.
Driver API applications can manage the current context with the
current-context stack: cuCtxPushCurrent() makes a new context
current, pushing it onto the top of the stack;
cuCtxPopCurrent() pops the current context and restores the
previous current context. Listing 9-2 gives a driver API version of
chMemcpyPeerToPeer(), which uses cuCtxPopCurrent() and
cuCtxPushCurrent() to perform a peer-to-peer memcpy between two
contexts.
The current context stack was introduced to CUDA in v2.2, and at the time, the CUDA runtime and driver API could not be used in the same application. That restriction has been relaxed in subsequent versions.
CUresultchMemcpyPeerToPeer( void *_dst, CUcontext dstContext, int dstDevice, const void *_src, CUcontext srcContext, int srcDevice, size_t N ) { CUresult status_cuda; CUdeviceptr dst = (CUdeviceptr) (intptr_t) _dst; CUdeviceptr src = (CUdeviceptr) (intptr_t) _src; int stagingIndex = 0; while ( N ) { size_t thisCopySize = min( N, STAGING_BUFFER_SIZE ); cu(CtxPushCurrent( srcContext ) ); cu(StreamWaitEvent( NULL, g_events[dstDevice][stagingIndex], 0 ) ); cu(MemcpyDtoHAsync( g_hostBuffers[stagingIndex], src, thisCopySize, NULL ) ); cu(EventRecord( g_events[srcDevice][stagingIndex], 0 ) ); cu(CtxPopCurrent( &srcContext ) ); cu(CtxPushCurrent( dstContext ) ); cu(StreamWaitEvent( NULL, g_events[srcDevice][stagingIndex], 0 ) ); cu(MemcpyHtoDAsync( dst, g_hostBuffers[stagingIndex], thisCopySize, NULL ) ); cu(EventRecord( g_events[dstDevice][stagingIndex], 0 ) ); cu(CtxPopCurrent( &dstContext ) ); dst += thisCopySize; src += thisCopySize; N -= thisCopySize; stagingIndex = 1 - stagingIndex; } // Wait until both devices are done cu(CtxPushCurrent( srcContext ) ); cu(CtxSynchronize() ); cu(CtxPopCurrent( &srcContext ) ); cu(CtxPushCurrent( dstContext ) ); cu(CtxSynchronize() ); cu(CtxPopCurrent( &dstContext ) ); Error_cuda: return status_cuda;}chMemcpyPeerToPeer() (driver API version) (source on GitHub)The N-body computation (described in detail in Chapter 14) computes N forces in O(N2) time, and the outputs may be computed independently. On a system with k GPUs, our multi-GPU implementation splits the computation into k parts that each compute \(\frac{N}{k}\) outputs.
Our implementation makes the common assumption that the GPUs are identical, so it divides the computation evenly; applications targeting GPUs of unequal performance or whose workloads have less-predictable runtimes can divide the computation more finely and have the host code submit work items to the GPUs from a queue.
Listing 9-3 gives a modified version of Listing 14-3 that takes two
additional parameters (a base index base and size n of the subarray of
forces), to compute a subset of the output array for an N-body
computation. This __device__ function is invoked by wrapper kernels that
are declared as __global__; it is structured this way to reuse the code
without incurring link errors. When this book was first written,
declaring the function __global__ would have generated a linker error
about duplicate symbols2.
inline __device__ voidComputeNBodyGravitation_Shared_multiGPU( float *force, float *posMass, float softeningSquared, size_t base, size_t n, size_t N ){ extern __shared__ float4 shPosMass[]; for ( int m = blockIdx.x*blockDim.x + threadIdx.x; m < n; m += blockDim.x*gridDim.x ) { size_t i = base+m; float acc[3] = {0}; float4 myPosMass = ((float4 *) posMass)[i];#pragma unroll 4 for ( int j = 0; j < N; j += blockDim.x ) { shPosMass[threadIdx.x] = ((float4 *) posMass)[j+threadIdx.x]; __syncthreads(); for ( size_t k = 0; k < blockDim.x; k++ ) { float fx, fy, fz; float4 bodyPosMass = shPosMass[k]; bodyBodyInteraction( &fx, &fy, &fz, myPosMass.x, myPosMass.y, myPosMass.z, bodyPosMass.x, bodyPosMass.y, bodyPosMass.z, bodyPosMass.w, softeningSquared ); acc[0] += fx; acc[1] += fy; acc[2] += fz; } __syncthreads(); } force[3*m+0] = acc[0]; force[3*m+1] = acc[1]; force[3*m+2] = acc[2]; }}The host code for a single-threaded, multi-GPU version of N-body is
shown in Listing 9-43. The arrays dptrPosMass
and dptrForce track the device pointers for the input and
output arrays for each GPU (the maximum number of GPUs is declared as a
constant in nbody.h - default is 32). Similar to dispatching work into
CUDA streams, the function uses separate loops for different stages of
the computation: the first loop allocates and populates the input array
for each GPU; the second loop launches the kernel and an asynchronous
copy of the output data; and the third loop calls
cudaDeviceSynchronize() on each GPU in turn. Structuring the function
this way maximizes CPU/GPU overlap: during the first loop, asynchronous
host→device memcpy’s to GPUs 0..i-1 can proceed while the CPU
is busy allocating memory for GPU i. If the kernel launch and
asynchronous device→host memcpy were in the first loop, the synchronous
cudaMalloc() calls would decrease performance because they are
synchronous with respect to the current GPU.
floatComputeGravitation_multiGPU_singlethread( float *force, float *posMass, float softeningSquared, size_t N){ cudaError_t status_cudart; float ret = 0.0f; float *dptrPosMass[g_maxGPUs]; float *dptrForce[g_maxGPUs]; int oldDevice; std::chrono::steady_clock::time_point start, end; start = std::chrono::steady_clock::now(); memset( dptrPosMass, 0, sizeof(dptrPosMass) ); memset( dptrForce, 0, sizeof(dptrForce) ); size_t bodiesPerGPU = N / g_numGPUs; if ( (0 != N % g_numGPUs) || (g_numGPUs > g_maxGPUs) ) { return 0.0f; } cuda(GetDevice( &oldDevice ) ); // kick off the asynchronous memcpy's - overlap GPUs pulling // host memory with the CPU time needed to do the memory // allocations. for ( int i = 0; i < g_numGPUs; i++ ) { cuda(SetDevice( i ) ); cuda(Malloc( &dptrPosMass[i], 4*N*sizeof(float) ) ); // we only need 3*N floatsw for the cross-check. otherwise we // would need 3*bodiesPerGPU cuda(Malloc( &dptrForce[i], 3*N*sizeof(float) ) ); cuda(MemcpyAsync( dptrPosMass[i], g_hostAOS_PosMass, 4*N*sizeof(float), cudaMemcpyHostToDevice ) ); } for ( int i = 0; i < g_numGPUs; i++ ) { cuda(SetDevice( i ) ); if ( g_bGPUCrossCheck ) { ComputeNBodyGravitation_multiGPU_onethread<<<300,256,256*sizeof(float4)>>>( dptrForce[i], dptrPosMass[i], softeningSquared, 0, N, N ); cuda(MemcpyAsync( g_hostAOS_gpuCrossCheckForce[i], dptrForce[i], 3*N*sizeof(float), cudaMemcpyDeviceToHost ) ); cuda(MemcpyAsync( g_hostAOS_Force+3*bodiesPerGPU*i, dptrForce[i]+3*bodiesPerGPU*i, 3*bodiesPerGPU*sizeof(float), cudaMemcpyDeviceToHost ) ); } else { ComputeNBodyGravitation_multiGPU_onethread<<<300,256,256*sizeof(float4)>>>( dptrForce[i], dptrPosMass[i], softeningSquared, i*bodiesPerGPU, bodiesPerGPU, N ); cuda(MemcpyAsync( g_hostAOS_Force+3*bodiesPerGPU*i, dptrForce[i], 3*bodiesPerGPU*sizeof(float), cudaMemcpyDeviceToHost ) ); } } // Synchronize with each GPU in turn. for ( int i = 0; i < g_numGPUs; i++ ) { cuda(SetDevice( i ) ); cuda(DeviceSynchronize() ); } end = std::chrono::steady_clock::now(); ret = std::chrono::duration<double>(end - start).count() * 1000.0f; if ( g_fGPUCrosscheckOutput ) { if ( 1 != fwrite( g_hostAOS_Force, 3*N*sizeof(float), 1, g_fGPUCrosscheckOutput ) ) goto Error_cudart; } if ( g_fGPUCrosscheckInput ) { if ( 1 != fread( g_hostAOS_Force_Golden, 3*N*sizeof(float), 1, g_fGPUCrosscheckInput ) ) goto Error_cudart; } Error_cudart: for ( int i = 0; i < g_numGPUs; i++ ) { cudaFree( dptrPosMass[i] ); cudaFree( dptrForce[i] ); } cudaSetDevice( oldDevice ); return ret;}