 }
 
 /*----< find_nearest_cluster() >---------------------------------------------*/
-__inline static
-int find_nearest_cluster(int     numClusters, /* no. clusters */
-                         int     numCoords,   /* no. coordinates */
-                         float  *object,      /* [numCoords] */
-                         float **clusters)    /* [numClusters][numCoords] */
+__global__ static
+void find_nearest_cluster(int numCoords,
+                          int numObjs,
+                          int numClusters,
+                          float *objects,   //  [numObjs][numCoords]
+                          float *clusters,  //  [numObjs][numCoords]
+                          int *membership,  //  [numObjs]
+                          int *delta)
 {
-    int   index, i;
-    float dist, min_dist;
-
-    /* find the cluster id that has min distance to object */
-    index    = 0;
-    min_dist = euclid_dist_2(numCoords, object, clusters[0]);
-
-    for (i=1; i<numClusters; i++) {
-        dist = euclid_dist_2(numCoords, object, clusters[i]);
-        /* no need square root */
-        if (dist < min_dist) { /* find the min and its array index */
-            min_dist = dist;
-            index    = i;
+    //  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;
+
+    membershipChanged[threadIdx.x] = 0;
+
+    if (objectId < numObjs) {
+        int   index, i;
+        float dist, min_dist;
+
+        /* find the cluster id that has min distance to object */
+        index    = 0;
+        min_dist = euclid_dist_2(numCoords, object, clusters);
+
+        for (i=1; i<numClusters; i++) {
+            dist = euclid_dist_2(numCoords, object, clusters + numCoords * i);
+            /* no need square root */
+            if (dist < min_dist) { /* find the min and its array index */
+                min_dist = dist;
+                index    = i;
+            }
+        }
+
+        if (membership[objectId] != index) {
+            membershipChanged[threadIdx.x] = 1;
+        }
+
+        /* assign the membership to object objectId */
+        membership[objectId] = index;
+
+        __syncthreads();    //  For membershipChanged[]
+
+        //  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) {
+            *delta = membershipChanged[0];
         }
     }
-    return(index);
 }
 
 /*----< cuda_kmeans() >-------------------------------------------------------*/

        delta /= numObjs;
    } while (delta > threshold && loop++ < 500);

    free(newClusters[0]);
    free(newClusters);
    free(newClusterSize);

                   float   threshold,    /* % objects change membership */
                   int    *membership)   /* out: [numObjs] */
{
    int      i, j, index, loop=1;
    int     *newClusterSize; /* [numClusters]: no. objects assigned in each
                                new cluster */
    float    delta;          /* % of objects change their clusters */
    float  **clusters;       /* out: [numClusters][numCoords] */
    float  **newClusters;    /* [numClusters][numCoords] */

    int *deviceDelta;
    float *deviceObjects;
    float *deviceClusters;
    int *deviceMembership;
//  int *deviceNewClusterSize;

    /* allocate a 2D space for returning variable clusters[] (coordinates
       of cluster centers) */
