Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
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
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <cuda.h>
#define THREADS_PER_DIM 16
#define BLOCKS_PER_DIM 16
#define THREADS_PER_BLOCK THREADS_PER_DIM*THREADS_PER_DIM
#include "kmeans_cuda_kernel.cu"
//#define BLOCK_DELTA_REDUCE
//#define BLOCK_CENTER_REDUCE
#define CPU_DELTA_REDUCE
#define CPU_CENTER_REDUCE
extern "C"
int setup(int argc, char** argv); /* function prototype */
// GLOBAL!!!!!
unsigned int num_threads_perdim = THREADS_PER_DIM; /* sqrt(256) -- see references for this choice */
unsigned int num_blocks_perdim = BLOCKS_PER_DIM; /* temporary */
unsigned int num_threads = num_threads_perdim*num_threads_perdim; /* number of threads */
unsigned int num_blocks = num_blocks_perdim*num_blocks_perdim; /* number of blocks */
/* _d denotes it resides on the device */
int *membership_new; /* newly assignment membership */
float *feature_d; /* inverted data array */
float *feature_flipped_d; /* original (not inverted) data array */
int *membership_d; /* membership on the device */
float *block_new_centers; /* sum of points in a cluster (per block) */
float *clusters_d; /* cluster centers on the device */
float *block_clusters_d; /* per block calculation of cluster centers */
int *block_deltas_d; /* per block calculation of deltas */
/* -------------- allocateMemory() ------------------- */
/* allocate device memory, calculate number of blocks and threads, and invert the data array */
extern "C"
void allocateMemory(int npoints, int nfeatures, int nclusters, float **features)
{
num_blocks = npoints / num_threads;
if (npoints % num_threads > 0) /* defeat truncation */
num_blocks++;
num_blocks_perdim = sqrt((double) num_blocks);
while (num_blocks_perdim * num_blocks_perdim < num_blocks) // defeat truncation (should run once)
num_blocks_perdim++;
num_blocks = num_blocks_perdim*num_blocks_perdim;
/* allocate memory for memory_new[] and initialize to -1 (host) */
membership_new = (int*) malloc(npoints * sizeof(int));
for(int i=0;i<npoints;i++) {
membership_new[i] = -1;
}
/* allocate memory for block_new_centers[] (host) */
block_new_centers = (float *) malloc(nclusters*nfeatures*sizeof(float));
/* allocate memory for feature_flipped_d[][], feature_d[][] (device) */
cudaMalloc((void**) &feature_flipped_d, npoints*nfeatures*sizeof(float));
cudaMemcpy(feature_flipped_d, features[0], npoints*nfeatures*sizeof(float), cudaMemcpyHostToDevice);
cudaMalloc((void**) &feature_d, npoints*nfeatures*sizeof(float));
/* invert the data array (kernel execution) */
invert_mapping<<<num_blocks,num_threads>>>(feature_flipped_d,feature_d,npoints,nfeatures);
/* allocate memory for membership_d[] and clusters_d[][] (device) */
cudaMalloc((void**) &membership_d, npoints*sizeof(int));
cudaMalloc((void**) &clusters_d, nclusters*nfeatures*sizeof(float));
#ifdef BLOCK_DELTA_REDUCE
// allocate array to hold the per block deltas on the gpu side
cudaMalloc((void**) &block_deltas_d, num_blocks_perdim * num_blocks_perdim * sizeof(int));
//cudaMemcpy(block_delta_d, &delta_h, sizeof(int), cudaMemcpyHostToDevice);
#endif
#ifdef BLOCK_CENTER_REDUCE
// allocate memory and copy to card cluster array in which to accumulate center points for the next iteration
cudaMalloc((void**) &block_clusters_d,
num_blocks_perdim * num_blocks_perdim *
nclusters * nfeatures * sizeof(float));
//cudaMemcpy(new_clusters_d, new_centers[0], nclusters*nfeatures*sizeof(float), cudaMemcpyHostToDevice);
#endif
}
/* -------------- allocateMemory() end ------------------- */
/* -------------- deallocateMemory() ------------------- */
/* free host and device memory */
extern "C"
void deallocateMemory()
{
free(membership_new);
free(block_new_centers);
cudaFree(feature_d);
cudaFree(feature_flipped_d);
cudaFree(membership_d);
cudaFree(clusters_d);
#ifdef BLOCK_CENTER_REDUCE
cudaFree(block_clusters_d);
#endif
#ifdef BLOCK_DELTA_REDUCE
cudaFree(block_deltas_d);
#endif
}
/* -------------- deallocateMemory() end ------------------- */
////////////////////////////////////////////////////////////////////////////////
// Program main //
int
main( int argc, char** argv)
{
// make sure we're running on the big card
cudaSetDevice(1);
// as done in the CUDA start/help document provided
setup(argc, argv);
}
// //
////////////////////////////////////////////////////////////////////////////////
/* ------------------- kmeansCuda() ------------------------ */
extern "C"
int // delta -- had problems when return value was of float type
kmeansCuda(float **feature, /* in: [npoints][nfeatures] */
int nfeatures, /* number of attributes for each point */
int npoints, /* number of data points */
int nclusters, /* number of clusters */
int *membership, /* which cluster the point belongs to */
float **clusters, /* coordinates of cluster centers */
int *new_centers_len, /* number of elements in each cluster */
float **new_centers /* sum of elements in each cluster */
int delta = 0; /* if point has moved */
int i,j; /* counters */
/* copy membership (host to device) */
cudaMemcpy(membership_d, membership_new, npoints*sizeof(int), cudaMemcpyHostToDevice);
/* copy clusters (host to device) */
cudaMemcpy(clusters_d, clusters[0], nclusters*nfeatures*sizeof(float), cudaMemcpyHostToDevice);
cudaChannelFormatDesc chDesc0 = cudaCreateChannelDesc<float>();
t_features.filterMode = cudaFilterModePoint;
t_features.normalized = false;
t_features.channelDesc = chDesc0;
if(cudaBindTexture(NULL, &t_features, feature_d, &chDesc0, npoints*nfeatures*sizeof(float)) != (cudaError_t) CUDA_SUCCESS)
printf("Couldn't bind features array to texture!\n");
cudaChannelFormatDesc chDesc1 = cudaCreateChannelDesc<float>();
t_features_flipped.filterMode = cudaFilterModePoint;
t_features_flipped.normalized = false;
t_features_flipped.channelDesc = chDesc1;
if(cudaBindTexture(NULL, &t_features_flipped, feature_flipped_d, &chDesc1, npoints*nfeatures*sizeof(float)) != (cudaError_t) CUDA_SUCCESS)
printf("Couldn't bind features_flipped array to texture!\n");
cudaChannelFormatDesc chDesc2 = cudaCreateChannelDesc<float>();
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
t_clusters.filterMode = cudaFilterModePoint;
t_clusters.normalized = false;
t_clusters.channelDesc = chDesc2;
if(cudaBindTexture(NULL, &t_clusters, clusters_d, &chDesc2, nclusters*nfeatures*sizeof(float)) != (cudaError_t) CUDA_SUCCESS)
printf("Couldn't bind clusters array to texture!\n");
/* copy clusters to constant memory */
cudaMemcpyToSymbol("c_clusters",clusters[0],nclusters*nfeatures*sizeof(float),0,cudaMemcpyHostToDevice);
/* setup execution parameters.
changed to 2d (source code on NVIDIA CUDA Programming Guide) */
dim3 grid( num_blocks_perdim, num_blocks_perdim );
dim3 threads( num_threads_perdim*num_threads_perdim );
/* execute the kernel */
kmeansPoint<<< grid, threads >>>( feature_d,
nfeatures,
npoints,
nclusters,
membership_d,
clusters_d,
block_clusters_d,
block_deltas_d);
cudaThreadSynchronize();
/* copy back membership (device to host) */
cudaMemcpy(membership_new, membership_d, npoints*sizeof(int), cudaMemcpyDeviceToHost);
#ifdef BLOCK_CENTER_REDUCE
/*** Copy back arrays of per block sums ***/
float * block_clusters_h = (float *) malloc(
num_blocks_perdim * num_blocks_perdim *
nclusters * nfeatures * sizeof(float));
cudaMemcpy(block_clusters_h, block_clusters_d,
num_blocks_perdim * num_blocks_perdim *
nclusters * nfeatures * sizeof(float),
cudaMemcpyDeviceToHost);
#endif
#ifdef BLOCK_DELTA_REDUCE
int * block_deltas_h = (int *) malloc(
num_blocks_perdim * num_blocks_perdim * sizeof(int));
cudaMemcpy(block_deltas_h, block_deltas_d,
num_blocks_perdim * num_blocks_perdim * sizeof(int),
cudaMemcpyDeviceToHost);
#endif
/* for each point, sum data points in each cluster
and see if membership has changed:
if so, increase delta and change old membership, and update new_centers;
otherwise, update new_centers */
delta = 0;
for (i = 0; i < npoints; i++)
{
int cluster_id = membership_new[i];
new_centers_len[cluster_id]++;
if (membership_new[i] != membership[i])
{
#ifdef CPU_DELTA_REDUCE
delta++;
#endif
membership[i] = membership_new[i];
}
#ifdef CPU_CENTER_REDUCE
for (j = 0; j < nfeatures; j++)
{
new_centers[cluster_id][j] += feature[i][j];
}
#endif
}
#ifdef BLOCK_DELTA_REDUCE
/*** calculate global sums from per block sums for delta and the new centers ***/
//debug
//printf("\t \t reducing %d block sums to global sum \n",num_blocks_perdim * num_blocks_perdim);
for(i = 0; i < num_blocks_perdim * num_blocks_perdim; i++) {
//printf("block %d delta is %d \n",i,block_deltas_h[i]);
delta += block_deltas_h[i];
}
#endif
#ifdef BLOCK_CENTER_REDUCE
for(int j = 0; j < nclusters;j++) {
for(int k = 0; k < nfeatures;k++) {
block_new_centers[j*nfeatures + k] = 0.f;
}
}
for(i = 0; i < num_blocks_perdim * num_blocks_perdim; i++) {
for(int j = 0; j < nclusters;j++) {
for(int k = 0; k < nfeatures;k++) {
block_new_centers[j*nfeatures + k] += block_clusters_h[i * nclusters*nfeatures + j * nfeatures + k];
}
}
}
#ifdef CPU_CENTER_REDUCE
//debug
/*for(int j = 0; j < nclusters;j++) {
for(int k = 0; k < nfeatures;k++) {
if(new_centers[j][k] > 1.001 * block_new_centers[j*nfeatures + k] || new_centers[j][k] < 0.999 * block_new_centers[j*nfeatures + k]) {
printf("\t \t for %d:%d, normal value is %e and gpu reduced value id %e \n",j,k,new_centers[j][k],block_new_centers[j*nfeatures + k]);
}
}
}*/
#endif
#ifdef BLOCK_CENTER_REDUCE
for(int j = 0; j < nclusters;j++) {
for(int k = 0; k < nfeatures;k++)
new_centers[j][k]= block_new_centers[j*nfeatures + k];
}
#endif
#endif
return delta;
}
/* ------------------- kmeansCuda() end ------------------------ */