The warp shuffle instructions (described in Section 8.6.1) enable
threads to exchange data without writing the data to shared memory. The
__shfl() intrinsic can be used to broadcast one thread’s register value
to all other threads in the warp; as shown in Listing 14-7, instead of
using tiles sized to the threadblock and using shared memory, we can use
tiles of size 32 (corresponding to the warp size) and broadcast the body
description read by each thread to the other threads within the
warp.
Interestingly, this strategy has 25% lower performance than the shared memory implementation (34 billion as opposed to 45.2 billion interactions per second); the gap persists on modern hardware, where a GeForce RTX 3060 runs the shuffle kernel at about 199 billion versus 291 billion for the shared-memory kernel. The warp shuffle instruction takes about as long as a read from shared memory, and the computation is tiled at the warp size (32 threads) rather than a thread block size. So it seems the benefits of warp shuffle are best realized when replacing both a write and a read to shared memory, not just a read; and warp shuffle should only be used if the kernel needs shared memory for other purposes.
__global__ voidComputeNBodyGravitation_Shuffle( float *force, float *posMass, float softeningSquared, size_t N ){ const int laneid = threadIdx.x & 31; for ( int i = blockIdx.x*blockDim.x + threadIdx.x; i < N; i += blockDim.x*gridDim.x ) { float acc[3] = {0}; float4 myPosMass = ((float4 *) posMass)[i]; for ( int j = 0; j < N; j += 32 ) { float4 shufSrcPosMass = ((float4 *) posMass)[j+laneid];#pragma unroll 32 for ( int k = 0; k < 32; k++ ) { float fx, fy, fz; float4 shufDstPosMass; shufDstPosMass.x = __shfl( shufSrcPosMass.x, k ); shufDstPosMass.y = __shfl( shufSrcPosMass.y, k ); shufDstPosMass.z = __shfl( shufSrcPosMass.z, k ); shufDstPosMass.w = __shfl( shufSrcPosMass.w, k ); bodyBodyInteraction( &fx, &fy, &fz, myPosMass.x, myPosMass.y, myPosMass.z, shufDstPosMass.x, shufDstPosMass.y, shufDstPosMass.z, shufDstPosMass.w, softeningSquared); acc[0] += fx; acc[1] += fy; acc[2] += fz; } } force[3*i+0] = acc[0]; force[3*i+1] = acc[1]; force[3*i+2] = acc[2]; }}ComputeNBodyGravitation_Shuffle (source on GitHub)