Blame view

stim/cuda/kmeans.cuh 13 KB
6911de0d   Davar   added acknowledge...
1
2
3
4
5
6
7
8
9
10
11
12
13
  //This software is dervied from Professor Wei-keng Liao's parallel k-means
  //clustering code obtained on November 21, 2010 from
  // http://users.eecs.northwestern.edu/~wkliao/Kmeans/index.html
  //(http://users.eecs.northwestern.edu/~wkliao/Kmeans/simple_kmeans.tar.gz).
  //
  //With his permission, Serban Giuroiu is publishing his CUDA implementation based on his code
  //under the open-source MIT license. See the LICENSE file for more details.
  
  // The original code can be found on Github ( https://github.com/serban/kmeans )
  // Here I have just made a few changes to get it to work
  
  
  
7f02091a   Davar   added cuda kmeans...
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
  
  #define malloc2D(name, xDim, yDim, type) do {               \
      name = (type **)malloc(xDim * sizeof(type *));          \
      assert(name != NULL);                                   \
      name[0] = (type *)malloc(xDim * yDim * sizeof(type));   \
      assert(name[0] != NULL);                                \
      for (size_t i = 1; i < xDim; i++)                       \
          name[i] = name[i-1] + yDim;                         \
  } while (0)
  
  
  
  static void handleError(cudaError_t error, const char* file, int line){
  
  	if(error != cudaSuccess){
  		cout << cudaGetErrorString(error) << " in " << file <<  " at line " << line << endl;
  		exit(1);
  	}
  }
  
  #define  handle_error(error)  handleError(error, __FILE__ , __LINE__)
  
  
  
  static inline int nextPowerOfTwo(int n) {
      n--;
  
      n = n >>  1 | n;
      n = n >>  2 | n;
      n = n >>  4 | n;
      n = n >>  8 | n;
      n = n >> 16 | n;
     // n = n >> 32 | n;    //  For 64-bit ints
  
      return ++n;
  }
  
  /*----< euclid_dist_2() >----------------------------------------------------*/
  /* square of Euclid distance between two multi-dimensional points            */
  __host__ __device__ inline static
  float euclid_dist_2(int    numCoords,
                      int    numObjs,
                      int    numClusters,
                      float *objects,     // [numCoords][numObjs]
                      float *clusters,    // [numCoords][numClusters]
                      int    objectId,
                      int    clusterId)
  {
      int i;
      float ans=0.0;
  
      for (i = 0; i < numCoords; i++) {
          ans += (objects[numObjs * i + objectId] - clusters[numClusters * i + clusterId]) *
                 (objects[numObjs * i + objectId] - clusters[numClusters * i + clusterId]);
      }
  
      return(ans);
  }
  
  /*----< find_nearest_cluster() >---------------------------------------------*/
  __global__ static
  void find_nearest_cluster(int numCoords,
                            int numObjs,
                            int numClusters,
                            float *objects,           //  [numCoords][numObjs]
                            float *deviceClusters,    //  [numCoords][numClusters]
                            int *membership,          //  [numObjs]
                            int *intermediates)
  {
      extern __shared__ char sharedMemory[];
  
      //  The type chosen for membershipChanged must be large enough to support
      //  reductions! There are blockDim.x elements, one for each thread in the
      //  block. See numThreadsPerClusterBlock in cuda_kmeans().
      unsigned char *membershipChanged = (unsigned char *)sharedMemory;
  #ifdef BLOCK_SHARED_MEM_OPTIMIZATION
      float *clusters = (float *)(sharedMemory + blockDim.x);
  #else
      float *clusters = deviceClusters;
  #endif
  
      membershipChanged[threadIdx.x] = 0;
  
  #ifdef BLOCK_SHARED_MEM_OPTIMIZATION
      //  BEWARE: We can overrun our shared memory here if there are too many
      //  clusters or too many coordinates! For reference, a Tesla C1060 has 16
      //  KiB of shared memory per block, and a GeForce GTX 480 has 48 KiB of
      //  shared memory per block.
      for (int i = threadIdx.x; i < numClusters; i += blockDim.x) {
          for (int j = 0; j < numCoords; j++) {
              clusters[numClusters * j + i] = deviceClusters[numClusters * j + i];
          }
      }
      __syncthreads();
  #endif
  
      int objectId = blockDim.x * blockIdx.x + threadIdx.x;
  
      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, numObjs, numClusters,
                                   objects, clusters, objectId, 0);
  
          for (i=1; i<numClusters; i++) {
              dist = euclid_dist_2(numCoords, numObjs, numClusters,
                                   objects, clusters, objectId, 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) {
              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];
      }
  }
  
  /*----< cuda_kmeans() >-------------------------------------------------------*/
  //
  //  ----------------------------------------
  //  DATA LAYOUT
  //
  //  objects         [numObjs][numCoords]
  //  clusters        [numClusters][numCoords]
  //  dimObjects      [numCoords][numObjs]
  //  dimClusters     [numCoords][numClusters]
  //  newClusters     [numCoords][numClusters]
  //  deviceObjects   [numCoords][numObjs]
  //  deviceClusters  [numCoords][numClusters]
  //  ----------------------------------------
  //
  /* return an array of cluster centers of size [numClusters][numCoords]       */
  float** cuda_kmeans(float **objects,      /* in: [numObjs][numCoords] */
            unsigned int     numCoords,    /* no. features */
            unsigned int     numObjs,      /* no. objects */
            unsigned int     numClusters,  /* no. clusters */
                     float   threshold,    /* % objects change membership */
                     int    *membership,   /* out: [numObjs] */
                     int	   loops)
  {
      int      i, j, index, loop=0;
      int     *newClusterSize; /* [numClusters]: no. objects assigned in each
                                  new cluster */
      float    delta;          /* % of objects change their clusters */
      float  **dimObjects;
      float  **clusters;       /* out: [numClusters][numCoords] */
      float  **dimClusters;
      float  **newClusters;    /* [numCoords][numClusters] */
  
      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++) {
              dimObjects[i][j] = objects[j][i];
          }
      }
  
      /* pick first numClusters elements of objects[] as initial cluster centers*/
      malloc2D(dimClusters, numCoords, numClusters, float);
      for (i = 0; i < numCoords; i++) {
          for (j = 0; j < numClusters; j++) {
              dimClusters[i][j] = dimObjects[i][j];
  		}
      }
  
      /* initialize membership[] */
      for (i=0; i<numObjs; i++) membership[i] = -1;
  
      /* need to initialize newClusterSize and newClusters[0] to all 0 */
      newClusterSize = (int*) calloc(numClusters, sizeof(int));
      assert(newClusterSize != NULL);
  
      malloc2D(newClusters, numCoords, numClusters, float);
      memset(newClusters[0], 0, numCoords * numClusters * sizeof(float));
  
      //  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.
  	cudaDeviceProp props;
  	handle_error(cudaGetDeviceProperties(&props, 0));
      const unsigned int numThreadsPerClusterBlock = props.maxThreadsPerBlock;
      const unsigned int numClusterBlocks =
          ceil(numObjs / (double)numThreadsPerClusterBlock);
  
  #ifdef BLOCK_SHARED_MEM_OPTIMIZATION
      const unsigned int clusterBlockSharedDataSize =
          numThreadsPerClusterBlock * sizeof(unsigned char) +
          numClusters * numCoords * sizeof(float);
  
      cudaDeviceProp deviceProp;
      int deviceNum;
      cudaGetDevice(&deviceNum);
      cudaGetDeviceProperties(&deviceProp, deviceNum);
  
      if (clusterBlockSharedDataSize > deviceProp.sharedMemPerBlock) {
  		std::cout << "ERROR: insufficient shared memory. Please don't use the definition 'BLOCK_SHARED_MEM_OPTIMIZATION'" << endl;
  		exit(1);
      }
  #else
      const unsigned int clusterBlockSharedDataSize =
          numThreadsPerClusterBlock * sizeof(unsigned char);
  #endif
  
      const unsigned int numReductionThreads =
          nextPowerOfTwo(numClusterBlocks);
      const unsigned int reductionBlockSharedDataSize =
          numReductionThreads * sizeof(unsigned int);
  
      handle_error(cudaMalloc((void**)&deviceObjects, numObjs*numCoords*sizeof(float)));
      handle_error(cudaMalloc((void**)&deviceClusters, numClusters*numCoords*sizeof(float)));
      handle_error(cudaMalloc((void**)&deviceMembership, numObjs*sizeof(int)));
      handle_error(cudaMalloc((void**)&deviceIntermediates, numReductionThreads*sizeof(unsigned int)));
  
      handle_error(cudaMemcpy(deviceObjects, dimObjects[0],
                numObjs*numCoords*sizeof(float), cudaMemcpyHostToDevice));
      handle_error(cudaMemcpy(deviceMembership, membership,
                numObjs*sizeof(int), cudaMemcpyHostToDevice));
  
      do {
          handle_error(cudaMemcpy(deviceClusters, dimClusters[0],
                    numClusters*numCoords*sizeof(float), cudaMemcpyHostToDevice));
  
          find_nearest_cluster
              <<< numClusterBlocks, numThreadsPerClusterBlock, clusterBlockSharedDataSize >>>
              (numCoords, numObjs, numClusters,
               deviceObjects, deviceClusters, deviceMembership, deviceIntermediates);
  
          cudaDeviceSynchronize(); 
  
          compute_delta <<< 1, numReductionThreads, reductionBlockSharedDataSize >>>
              (deviceIntermediates, numClusterBlocks, numReductionThreads);
  
          cudaDeviceSynchronize(); 
  
          int d;
          handle_error(cudaMemcpy(&d, deviceIntermediates,
                    sizeof(int), cudaMemcpyDeviceToHost));
          delta = (float)d;
  
          handle_error(cudaMemcpy(membership, deviceMembership,
                    numObjs*sizeof(int), cudaMemcpyDeviceToHost));
  
          for (i=0; i<numObjs; i++) {
              /* find the array index of nestest cluster center */
              index = membership[i];
  
              /* update new cluster centers : sum of objects located within */
              newClusterSize[index]++;
              for (j=0; j<numCoords; j++)
                  newClusters[j][index] += objects[i][j];
          }
  
          //  TODO: Flip the nesting order
          //  TODO: Change layout of newClusters to [numClusters][numCoords]
          /* average the sum and replace old cluster centers with newClusters */
          for (i=0; i<numClusters; i++) {
              for (j=0; j<numCoords; j++) {
                  if (newClusterSize[i] > 0)
                      dimClusters[j][i] = newClusters[j][i] / newClusterSize[i];
                  newClusters[j][i] = 0.0;   /* set back to 0 */
              }
              newClusterSize[i] = 0;   /* set back to 0 */
          }
  
          delta /= numObjs;
      } while (delta > threshold && loop++ < loops);
  
      
  
      /* 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];
          }
      }
  
      handle_error(cudaFree(deviceObjects));
      handle_error(cudaFree(deviceClusters));
      handle_error(cudaFree(deviceMembership));
      handle_error(cudaFree(deviceIntermediates));
  
      free(dimObjects[0]);
      free(dimObjects);
      free(dimClusters[0]);
      free(dimClusters);
      free(newClusters[0]);
      free(newClusters);
      free(newClusterSize);
  
      return clusters;
  }