 
     msg("Loop iterations: %d\n", loop)
 
-    checkCuda(cudaFree(deviceDelta));
     checkCuda(cudaFree(deviceObjects));
     checkCuda(cudaFree(deviceClusters));
     checkCuda(cudaFree(deviceMembership));
+    checkCuda(cudaFree(deviceIntermediates));
 //  checkCuda(cudaFree(deviceNewClusterSize));
 
     free(newClusters[0]);

        checkCuda(cudaMemcpy(deviceClusters, clusters[0],
                  numClusters*numCoords*sizeof(float), cudaMemcpyHostToDevice));

        const unsigned int numThreadsPerBlock = 128;
//      const unsigned int numThreadsPerBlock = nextPowerOfTwo(numObjs);
        const unsigned int numBlocks = (numObjs + numThreadsPerBlock - 1) / numThreadsPerBlock;
//      const unsigned int numBlocks = 1;
        const unsigned int sharedDataSize = numThreadsPerBlock * sizeof(unsigned short);

        msg("--- CUDA time! ---\n")
        find_nearest_cluster <<< numBlocks, numThreadsPerBlock, sharedDataSize >>>
            (numCoords, numObjs, numClusters,
             deviceObjects, deviceClusters, deviceMembership, deviceDelta);

        cudaThreadSynchronize();
        checkLastCudaError();

        int d;
        checkCuda(cudaMemcpy(&d, deviceDelta,
                  sizeof(int), cudaMemcpyDeviceToHost));
        delta = (float)d;
        msg("Delta = %d\n", d)

        checkCuda(cudaMemcpy(clusters[0], deviceClusters,
              numClusters*numCoords*sizeof(float), cudaMemcpyDeviceToHost));
        checkCuda(cudaMemcpy(membership, deviceMembership,
                  numObjs*sizeof(int), cudaMemcpyDeviceToHost));
//      checkCuda(cudaMemcpy(newClusterSize, deviceNewClusterSize,

        //  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];
    }
}

