     for (i=1; i<numClusters; i++)
         newClusters[i] = newClusters[i-1] + numCoords;
 
-    checkCuda(cudaMalloc(&deviceDelta, sizeof(int)));
+    //  To support reduction, numThreadsPerClusterBlock *must* be a power of
+    //  two, and it *must* be no larger than the number of bits that will
+    //  fit into an unsigned char, the type used to keep track of membership
+    //  changes in the kernel.
+    const unsigned int numThreadsPerClusterBlock = 128;
+    const unsigned int numClusterBlocks =
+        (numObjs + numThreadsPerClusterBlock - 1) / numThreadsPerClusterBlock;
+    const unsigned int clusterBlockSharedDataSize =
+        numThreadsPerClusterBlock * sizeof(unsigned char);
+
+    const unsigned int numReductionThreads =
+        nextPowerOfTwo(numClusterBlocks);
+    const unsigned int reductionBlockSharedDataSize =
+        numReductionThreads * sizeof(unsigned int);
+
     checkCuda(cudaMalloc(&deviceObjects, numObjs*numCoords*sizeof(float)));
     checkCuda(cudaMalloc(&deviceClusters, numClusters*numCoords*sizeof(float)));
     checkCuda(cudaMalloc(&deviceMembership, numObjs*sizeof(int)));
+    checkCuda(cudaMalloc(&deviceIntermediates, numReductionThreads*sizeof(unsigned int)));
 //  checkCuda(cudaMalloc(&deviceNewClusterSize, numClusters*sizeof(int)));
 
     checkCuda(cudaMemcpy(deviceObjects, objects[0],

                          float *objects,   //  [numObjs][numCoords]
                          float *clusters,  //  [numObjs][numCoords]
                          int *membership,  //  [numObjs]
                          int *delta)
{
    //  The type chosen here must be large enough to support reductions!
    extern __shared__ unsigned short membershipChanged[];

    int objectId = blockDim.x * blockIdx.x + threadIdx.x;
    float *object = objects + numCoords * objectId;

        //  blockDim.x *must* be a power of two!
        for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
            if (threadIdx.x < s) {
                membershipChanged[threadIdx.x] +=
                    membershipChanged[threadIdx.x + s];
            }
            __syncthreads();
        }

        if (threadIdx.x == 0) {
            intermediates[blockIdx.x] = membershipChanged[0];
        }
    }
}

__global__ static
void compute_delta(int *deviceIntermediates,
                   int numIntermediates,    //  The actual number of intermediates
                   int numIntermediates2)   //  The next power of two
{
    //  The number of elements in this array should be equal to
    //  numIntermediates2, the number of threads launched. It *must* be a power
    //  of two!
    extern __shared__ unsigned int intermediates[];

    //  Copy global intermediate values into shared memory.
    intermediates[threadIdx.x] =
        (threadIdx.x < numIntermediates) ? deviceIntermediates[threadIdx.x] : 0;

    __syncthreads();

    //  numIntermediates2 *must* be a power of two!
    for (unsigned int s = numIntermediates2 / 2; s > 0; s >>= 1) {
        if (threadIdx.x < s) {
            intermediates[threadIdx.x] += intermediates[threadIdx.x + s];
        }
        __syncthreads();
    }

    if (threadIdx.x == 0) {
        deviceIntermediates[0] = intermediates[0];
    }
}

