     } while (delta > threshold && loop++ < 500);
 
     *loop_iterations = loop + 1;
+
+    /* allocate a 2D space for returning variable clusters[] (coordinates
+       of cluster centers) */
+    malloc2D(clusters, numClusters, numCoords, float);
+    for (i = 0; i < numClusters; i++) {
+        for (j = 0; j < numCoords; j++) {
+            clusters[i][j] = dimClusters[j][i];
+        }
+    }
 
     checkCuda(cudaFree(deviceObjects));
     checkCuda(cudaFree(deviceClusters));

    float    delta;          /* % of objects change their clusters */
    float  **dimObjects;
    float  **clusters;       /* out: [numClusters][numCoords] */
    float  **newClusters;    /* [numClusters][numCoords] */

    float *deviceObjects;
    float *deviceClusters;
    int *deviceMembership;
    int *deviceIntermediates;

    //  Copy objects given in [numObjs][numCoords] layout
    //  to new [numCoords][numObjs] layout
    malloc2D(dimObjects, numCoords, numObjs, float);
    for (i = 0; i < numCoords; i++) {
        for (j = 0; j < numObjs; j++) {

              numObjs*sizeof(int), cudaMemcpyHostToDevice));

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

        find_nearest_cluster
