【CUDA并行编程之八】Cuda实现Kmeans算法

「爱情、让人受尽委屈。」 2022-08-01 12:17 238阅读 0赞

本文主要介绍如何使用CUDA并行计算框架编程实现机器学习中的Kmeans算法，Kmeans算法的详细介绍在[这里][Link 1]，本文重点在并行实现的过程。

当然还是简单的回顾一下kmeans算法的串行过程：

**伪代码：**

**\[cpp\]**  [view plain][] [copy][view plain] [![在CODE上查看代码片][CODE]][CODE_CODE] [![派生到我的代码片][ico_fork.svg]][ico_fork.svg 1]

1.  创建k个点作为起始质心(经常是随机选择)  
2.  当任意一个点的簇分配结果发生改变时  
3.      对数据集中的每个数据点  
4.          对每个质心  
5.              计算质心与数据点之间的距离  
6.          将数据点分配到距其最近的簇  
7.      对每一个簇，计算簇中所有点的均值并将均值作为质心

我们可以观察到有两个部分可以并行优化：

①line03-04：将每个数据点到多个质心的距离计算进行并行化

②line05：将数据点到某个执行的距离计算进行并行化

**KMEANS类：**

**\[cpp\]**  [view plain][] [copy][view plain] [![在CODE上查看代码片][CODE]][CODE_CODE] [![派生到我的代码片][ico_fork.svg]][ico_fork.svg 1]

1.  **class** KMEANS  
2.  \{  
3.  **private**:  
4.      **int** numClusters;  
5.      **int** numCoords;  
6.      **int** numObjs;  
7.      **int** \*membership;//\[numObjs\]  
8.      **char** \*filename;   
9.      **float** \*\*objects;//\[numObjs\]\[numCoords\] data objects  
10.     **float** \*\*clusters;//\[numClusters\]\[unmCoords\] cluster center  
11.     **float** threshold;  
12.     **int** loop\_iterations;  
13.   
14. **public**:  
15.     KMEANS(**int** k);  
16.     **void** file\_read(**char** \*fn);  
17.     **void** file\_write();  
18.     **void** cuda\_kmeans();  
19.     **inline** **int** nextPowerOfTwo(**int** n);  
20.     **void** free\_memory();  
21.     **virtual** ~KMEANS();  
22. \};//KMEANS

**成员变量：**

numClusters：中心点的个数

numCoords：每个数据点的维度

numObjs：数据点的个数

membership：每个数据点所属类别的数组，维度为numObjs

filename：读入的文件名

objects：所有数据点，维度为\[numObjs\]\[numCoords\]

clusters：中心点数据，维度为\[numObjs\]\[numCoords\]

threshold：控制循环次数的一个域值

loop\_iterations：循环的迭代次数

**成员函数：**

**KMEANS(int k)：**含参构造函数。初始化成员变量

**file\_read(char \*fn)：**读入文件数据并初始化object以及membership变量

**file\_write()：**将计算结果写回到结果文件中去

**cuda\_kmeans()：**kmeans计算的入口函数

**nextPowerOfTwo(int n)：**它计算大于等于输入参数n的第一个2的幂次数。

**free\_memory()：**释放内存空间

**~KMEANS()：**析构函数

**并行的代码主要三个函数：**

find\_nearest\_cluster(...)

compute\_delta(...)

euclid\_dist\_2(...)

**首先看一下函数euclid\_dist\_2(...)：**

**\[cpp\]**  [view plain][] [copy][view plain] [![在CODE上查看代码片][CODE]][CODE_CODE] [![派生到我的代码片][ico_fork.svg]][ico_fork.svg 1]

1.  \_\_host\_\_ \_\_device\_\_ **inline** **static**   
2.  **float** euclid\_dist\_2(**int** numCoords,**int** numObjs,**int** numClusters,**float** \*objects,**float** \*clusters,**int** objectId,**int** clusterId)  
3.  \{  
4.      **int** i;  
5.      **float** ans = 0;  
6.      **for**( i=0;i<numCoords;i++ )  
7.      \{  
8.          ans += ( objects\[numObjs \* i + objectId\] - clusters\[numClusters\*i + clusterId\] ) \*  
9.                 ( objects\[numObjs \* i + objectId\] - clusters\[numClusters\*i + clusterId\] ) ;  
10.     \}  
11.     **return** ans;  
12. \}

这段代码实际上就是并行的计算向量objects\[objectId\]和clusters\[clusterId\]之间的距离，即第objectId个数据点到第clusterId个中心点的距离。

**再看一下函数compute\_delta(...)：**

**\[cpp\]**  [view plain][] [copy][view plain] [![在CODE上查看代码片][CODE]][CODE_CODE] [![派生到我的代码片][ico_fork.svg]][ico_fork.svg 1]

1.  /\* 
2.  \* numIntermediates:The actual number of intermediates 
3.  \* numIntermediates2:The next power of two 
4.  \*/  
5.  \_\_global\_\_ **static** **void** compute\_delta(**int** \*deviceIntermediates,**int** numIntermediates, **int** numIntermediates2)  
6.  \{  
7.      **extern** \_\_shared\_\_ unsigned **int** intermediates\[\];  
8.    
9.      intermediates\[threadIdx.x\] = (threadIdx.x < numIntermediates) ? deviceIntermediates\[threadIdx.x\] : 0 ;  
10.     \_\_syncthreads();  
11.   
12.     //numIntermediates2 \*must\* be a power of two!  
13.     **for**(unsigned **int** s = numIntermediates2 /2 ; s > 0 ; s>>=1)  
14.     \{  
15.         **if**(threadIdx.x < s)    
16.         \{  
17.             intermediates\[threadIdx.x\] += intermediates\[threadIdx.x + s\];     
18.         \}  
19.         \_\_syncthreads();  
20.     \}  
21.     **if**(threadIdx.x == 0)  
22.     \{  
23.         deviceIntermediates\[0\] = intermediates\[0\];  
24.     \}  
25. \}

这段代码的意义就是将一个线程块中每个线程的对应的intermediates的数据求和最后放到deviceIntermediates\[0\]中去然后拷贝回主存块中去。这个问题的更好的解释在[这里][Link 2]，实际上就是一个数组求和的问题，应用在这里求得的是有改变的membership中所有数据的和，即改变了簇的点的个数。

**最后再看函数finid\_nearest\_cluster(...)：**

**\[cpp\]**  [view plain][] [copy][view plain] [![在CODE上查看代码片][CODE]][CODE_CODE] [![派生到我的代码片][ico_fork.svg]][ico_fork.svg 1]

1.  /\* 
2.  \* objects:\[numCoords\]\[numObjs\] 
3.  \* deviceClusters:\[numCoords\]\[numClusters\] 
4.  \* membership:\[numObjs\] 
5.  \*/  
6.  \_\_global\_\_ **static** **void** find\_nearest\_cluster(**int** numCoords,**int** numObjs,**int** numClusters,**float** \*objects, **float** \*deviceClusters,**int** \*membership ,**int** \*intermediates)  
7.  \{  
8.      **extern** \_\_shared\_\_ **char** sharedMemory\[\];  
9.      unsigned **char** \*membershipChanged = (unsigned **char** \*)sharedMemory;  
10.     **float** \*clusters = deviceClusters;  
11.   
12.     membershipChanged\[threadIdx.x\] = 0;  
13.   
14.     **int** objectId = blockDim.x \* blockIdx.x + threadIdx.x;  
15.     **if**( objectId < numObjs )  
16.     \{  
17.         **int** index;  
18.         **float** dist,min\_dist;  
19.         /\*find the cluster id that has min distance to object\*/  
20.         index = 0;  
21.         min\_dist = euclid\_dist\_2(numCoords,numObjs,numClusters,objects,clusters,objectId,0);  
22.           
23.         **for**(**int** i=0;i<numClusters;i++)  
24.         \{  
25.             dist = euclid\_dist\_2(numCoords,numObjs,numClusters,objects,clusters,objectId,i) ;  
26.             /\* no need square root \*/  
27.             **if**( dist < min\_dist )  
28.             \{  
29.                 min\_dist = dist;  
30.                 index = i;  
31.             \}  
32.         \}  
33.   
34.         **if**( membership\[objectId\]!=index )  
35.         \{  
36.             membershipChanged\[threadIdx.x\] = 1;   
37.         \}  
38.         //assign the membership to object objectId  
39.         membership\[objectId\] = index;  
40.   
41.         \_\_syncthreads(); //for membershipChanged\[\]  
42.   
43. \#if 1  
44.         //blockDim.x \*must\* be a power of two!  
45.         **for**(unsigned **int** s = blockDim.x / 2; s > 0 ;s>>=1)  
46.         \{  
47.             **if**(threadIdx.x < s)    
48.             \{  
49.                 membershipChanged\[threadIdx.x\] += membershipChanged\[threadIdx.x + s\];//calculate all changed values and save result to membershipChanged\[0\]  
50.             \}  
51.             \_\_syncthreads();  
52.         \}  
53.         **if**(threadIdx.x == 0)  
54.         \{  
55.             intermediates\[blockIdx.x\] = membershipChanged\[0\];  
56.         \}  
57. \#endif  
58.     \}  
59. \}//find\_nearest\_cluster

这个函数计算的就是第objectId个数据点到numClusters个中心点的距离，然后根据情况比较更新membership。

这三个函数将所有能够并行的地方都进行了并行，实现了整体算法的并行化~

在此呈上全部代码：

**kmeans.h:**

**\[cpp\]**  [view plain][] [copy][view plain] [![在CODE上查看代码片][CODE]][CODE_CODE] [![派生到我的代码片][ico_fork.svg]][ico_fork.svg 1]

1.  \#ifndef \_H\_KMEANS  
2.  \#define \_H\_KMEANS  
3.    
4.  \#include <assert.h>  
5.    
6.  \#define malloc2D(name, xDim, yDim, type) do \{               \\  
7.      name = (type \*\*)malloc(xDim \* **sizeof**(type \*));          \\  
8.      assert(name != NULL);                                   \\  
9.      name\[0\] = (type \*)malloc(xDim \* yDim \* **sizeof**(type));   \\  
10.     assert(name\[0\] != NULL);                                \\  
11.     **for** (**size\_t** i = 1; i < xDim; i++)                       \\  
12.         name\[i\] = name\[i-1\] + yDim;                         \\  
13. \} **while** (0)  
14.   
15.   
16. **double**  wtime(**void**);  
17.   
18. \#endif

wtime.cu:

**\[cpp\]**  [view plain][] [copy][view plain] [![在CODE上查看代码片][CODE]][CODE_CODE] [![派生到我的代码片][ico_fork.svg]][ico_fork.svg 1]

1.  \#include <sys/time.h>  
2.  \#include <stdio.h>  
3.  \#include <stdlib.h>  
4.    
5.  **double** wtime(**void**)   
6.  \{  
7.      **double**          now\_time;  
8.      **struct** timeval  etstart;  
9.      **struct** timezone tzp;  
10.   
11.     **if** (gettimeofday(&etstart, &tzp) == -1)  
12.         perror("Error: calling gettimeofday() not successful.\\n");  
13.   
14.     now\_time = ((**double**)etstart.tv\_sec) +              /\* in seconds \*/  
15.                ((**double**)etstart.tv\_usec) / 1000000.0;  /\* in microseconds \*/  
16.     **return** now\_time;  
17. \}

cuda\_kmeans.cu：

**\[cpp\]**  [view plain][] [copy][view plain] [![在CODE上查看代码片][CODE]][CODE_CODE] [![派生到我的代码片][ico_fork.svg]][ico_fork.svg 1]

1.  \#include <stdio.h>  
2.  \#include <stdlib.h>  
3.  \#include <string.h>  
4.  \#include <sys/types.h>  
5.  \#include <sys/stat.h>  
6.  \#include <unistd.h>  
7.  \#include <iostream>  
8.  \#include <cassert>  
9.    
10. \#include "kmeans.h"  
11.   
12. **using** **namespace** std;  
13.   
14. **const** **int** MAX\_CHAR\_PER\_LINE = 1024;  
15.   
16. **class** KMEANS  
17. \{  
18. **private**:  
19.     **int** numClusters;  
20.     **int** numCoords;  
21.     **int** numObjs;  
22.     **int** \*membership;//\[numObjs\]  
23.     **char** \*filename;   
24.     **float** \*\*objects;//\[numObjs\]\[numCoords\] data objects  
25.     **float** \*\*clusters;//\[numClusters\]\[unmCoords\] cluster center  
26.     **float** threshold;  
27.     **int** loop\_iterations;  
28.   
29. **public**:  
30.     KMEANS(**int** k);  
31.     **void** file\_read(**char** \*fn);  
32.     **void** file\_write();  
33.     **void** cuda\_kmeans();  
34.     **inline** **int** nextPowerOfTwo(**int** n);  
35.     **void** free\_memory();  
36.     **virtual** ~KMEANS();  
37. \};  
38.   
39. KMEANS::~KMEANS()  
40. \{  
41.     free(membership);  
42.     free(clusters\[0\]);  
43.     free(clusters);  
44.     free(objects\[0\]);  
45.     free(objects);  
46. \}  
47.   
48. KMEANS::KMEANS(**int** k)  
49. \{  
50.     threshold = 0.001;  
51.     numObjs = 0;  
52.     numCoords = 0;  
53.     numClusters = k;  
54.     filename = NULL;  
55.     loop\_iterations = 0;  
56. \}  
57.   
58. **void** KMEANS::file\_write()  
59. \{  
60.     **FILE** \*fptr;  
61.     **char** outFileName\[1024\];  
62.   
63.     //output:the coordinates of the cluster centres  
64.     sprintf(outFileName,"%s.cluster\_centres",filename);  
65.     printf("Writingcoordinates of K=%d cluster centers to file \\"%s\\"\\n",numClusters,outFileName);  
66.     fptr = fopen(outFileName,"w");  
67.     **for**(**int** i=0;i<numClusters;i++)  
68.     \{  
69.         fprintf(fptr,"%d ",i)   ;  
70.         **for**(**int** j=0;j<numCoords;j++)  
71.             fprintf(fptr,"%f ",clusters\[i\]\[j\]);  
72.         fprintf(fptr,"\\n");  
73.     \}  
74.     fclose(fptr);  
75.   
76.     //output:the closest cluster centre to each of the data points  
77.     sprintf(outFileName,"%s.membership",filename);  
78.     printf("writing membership of N=%d data objects to file \\"%s\\" \\n",numObjs,outFileName);  
79.     fptr = fopen(outFileName,"w");  
80.     **for**(**int** i=0;i<numObjs;i++)  
81.     \{  
82.         fprintf(fptr,"%d %d\\n",i,membership\[i\]) ;  
83.     \}  
84.     fclose(fptr);  
85. \}  
86.   
87. **inline** **int** KMEANS::nextPowerOfTwo(**int** n)  
88. \{  
89.     n--;  
90.     n = n >> 1 | n;  
91.     n = n >> 2 | n;  
92.     n = n >> 4 | n;  
93.     n = n >> 8 | n;  
94.     n = n >> 16 | n;  
95.     //n = n >> 32 | n; // for 64-bit ints  
96.     **return** ++n;  
97. \}  
98.   
99. \_\_host\_\_ \_\_device\_\_ **inline** **static**   
100. **float** euclid\_dist\_2(**int** numCoords,**int** numObjs,**int** numClusters,**float** \*objects,**float** \*clusters,**int** objectId,**int** clusterId)  
101. \{  
102.     **int** i;  
103.     **float** ans = 0;  
104.     **for**( i=0;i<numCoords;i++ )  
105.     \{  
106.         ans += ( objects\[numObjs \* i + objectId\] - clusters\[numClusters\*i + clusterId\] ) \*  
107.                ( objects\[numObjs \* i + objectId\] - clusters\[numClusters\*i + clusterId\] ) ;  
108.     \}  
109.     **return** ans;  
110. \}  
111.   
112. /\* 
113. \* numIntermediates:The actual number of intermediates 
114. \* numIntermediates2:The next power of two 
115. \*/  
116. \_\_global\_\_ **static** **void** compute\_delta(**int** \*deviceIntermediates,**int** numIntermediates, **int** numIntermediates2)  
117. \{  
118.     **extern** \_\_shared\_\_ unsigned **int** intermediates\[\];  
119.   
120.     intermediates\[threadIdx.x\] = (threadIdx.x < numIntermediates) ? deviceIntermediates\[threadIdx.x\] : 0 ;  
121.     \_\_syncthreads();  
122.   
123.     //numIntermediates2 \*must\* be a power of two!  
124.     **for**(unsigned **int** s = numIntermediates2 /2 ; s > 0 ; s>>=1)  
125.     \{  
126.         **if**(threadIdx.x < s)    
127.         \{  
128.             intermediates\[threadIdx.x\] += intermediates\[threadIdx.x + s\];     
129.         \}  
130.         \_\_syncthreads();  
131.     \}  
132.     **if**(threadIdx.x == 0)  
133.     \{  
134.         deviceIntermediates\[0\] = intermediates\[0\];  
135.     \}  
136. \}  
137.   
138. /\* 
139. \* objects:\[numCoords\]\[numObjs\] 
140. \* deviceClusters:\[numCoords\]\[numClusters\] 
141. \* membership:\[numObjs\] 
142. \*/  
143. \_\_global\_\_ **static** **void** find\_nearest\_cluster(**int** numCoords,**int** numObjs,**int** numClusters,**float** \*objects, **float** \*deviceClusters,**int** \*membership ,**int** \*intermediates)  
144. \{  
145.     **extern** \_\_shared\_\_ **char** sharedMemory\[\];  
146.     unsigned **char** \*membershipChanged = (unsigned **char** \*)sharedMemory;  
147.     **float** \*clusters = deviceClusters;  
148.   
149.     membershipChanged\[threadIdx.x\] = 0;  
150.   
151.     **int** objectId = blockDim.x \* blockIdx.x + threadIdx.x;  
152.     **if**( objectId < numObjs )  
153.     \{  
154.         **int** index;  
155.         **float** dist,min\_dist;  
156.         /\*find the cluster id that has min distance to object\*/  
157.         index = 0;  
158.         min\_dist = euclid\_dist\_2(numCoords,numObjs,numClusters,objects,clusters,objectId,0);  
159.           
160.         **for**(**int** i=0;i<numClusters;i++)  
161.         \{  
162.             dist = euclid\_dist\_2(numCoords,numObjs,numClusters,objects,clusters,objectId,i) ;  
163.             /\* no need square root \*/  
164.             **if**( dist < min\_dist )  
165.             \{  
166.                 min\_dist = dist;  
167.                 index = i;  
168.             \}  
169.         \}  
170.   
171.         **if**( membership\[objectId\]!=index )  
172.         \{  
173.             membershipChanged\[threadIdx.x\] = 1;   
174.         \}  
175.         //assign the membership to object objectId  
176.         membership\[objectId\] = index;  
177.   
178.         \_\_syncthreads(); //for membershipChanged\[\]  
179.   
180. \#if 1  
181.         //blockDim.x \*must\* be a power of two!  
182.         **for**(unsigned **int** s = blockDim.x / 2; s > 0 ;s>>=1)  
183.         \{  
184.             **if**(threadIdx.x < s)    
185.             \{  
186.                 membershipChanged\[threadIdx.x\] += membershipChanged\[threadIdx.x + s\];//calculate all changed values and save result to membershipChanged\[0\]  
187.             \}  
188.             \_\_syncthreads();  
189.         \}  
190.         **if**(threadIdx.x == 0)  
191.         \{  
192.             intermediates\[blockIdx.x\] = membershipChanged\[0\];  
193.         \}  
194. \#endif  
195.     \}  
196. \}//find\_nearest\_cluster  
197.   
198. **void** KMEANS::cuda\_kmeans()  
199. \{  
200.     **int** index,loop = 0;  
201.     **int** \*newClusterSize;//\[numClusters\]:no.objects assigned in each new cluster  
202.     **float** delta; //% of objects changes their clusters  
203.     **float** \*\*dimObjects;//\[numCoords\]\[numObjs\]  
204.     **float** \*\*dimClusters;  
205.     **float** \*\*newClusters;//\[numCoords\]\[numClusters\]  
206.   
207.     **float** \*deviceObjects; //\[numCoords\]\[numObjs\]  
208.     **float** \*deviceClusters; //\[numCoords\]\[numclusters\]  
209.     **int** \*deviceMembership;  
210.     **int** \*deviceIntermediates;  
211.   
212.     //Copy objects given in \[numObjs\]\[numCoords\] layout to new \[numCoords\]\[numObjs\] layout  
213.     malloc2D(dimObjects,numCoords,numObjs,**float**);  
214.     **for**(**int** i=0;i<numCoords;i++)  
215.     \{  
216.         **for**(**int** j=0;j<numObjs;j++)  
217.         \{  
218.             dimObjects\[i\]\[j\] = objects\[j\]\[i\];     
219.         \}  
220.     \}  
221.     //pick first numClusters elements of objects\[\] as initial cluster centers  
222.     malloc2D(dimClusters, numCoords, numClusters,**float**);  
223.     **for**(**int** i=0;i<numCoords;i++)  
224.     \{  
225.         **for**(**int** j=0;j<numClusters;j++)  
226.         \{  
227.             dimClusters\[i\]\[j\] = dimObjects\[i\]\[j\];  
228.         \}  
229.     \}  
230.     newClusterSize = **new** **int**\[numClusters\];  
231.     assert(newClusterSize!=NULL);  
232.     malloc2D(newClusters,numCoords,numClusters,**float**);  
233.     memset(newClusters\[0\],0,numCoords \* numClusters \* **sizeof**(**float**) );  
234.       
235.     //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.  
236.     **const** unsigned **int** numThreadsPerClusterBlock = 32;  
237.     **const** unsigned **int** numClusterBlocks = (numObjs + numThreadsPerClusterBlock -1)/numThreadsPerClusterBlock;  
238.     **const** unsigned **int** numReductionThreads = nextPowerOfTwo(numClusterBlocks);  
239.   
240.     **const** unsigned **int** clusterBlockSharedDataSize = numThreadsPerClusterBlock \* **sizeof**(unsigned **char**);  
241.   
242.     **const** unsigned **int** reductionBlockSharedDataSize = numReductionThreads \* **sizeof**(unsigned **int**);  
243.   
244.     cudaMalloc(&deviceObjects,numObjs\*numCoords\***sizeof**(**float**));  
245.     cudaMalloc(&deviceClusters,numClusters\*numCoords\***sizeof**(**float**));  
246.     cudaMalloc(&deviceMembership,numObjs\***sizeof**(**int**));  
247.     cudaMalloc(&deviceIntermediates,numReductionThreads\***sizeof**(unsigned **int**));  
248.   
249.     cudaMemcpy(deviceObjects,dimObjects\[0\],numObjs\*numCoords\***sizeof**(**float**),cudaMemcpyHostToDevice);  
250.     cudaMemcpy(deviceMembership,membership,numObjs\***sizeof**(**int**),cudaMemcpyHostToDevice);  
251.   
252.     **do**  
253.     \{  
254.         cudaMemcpy(deviceClusters,dimClusters\[0\],numClusters\*numCoords\***sizeof**(**float**),cudaMemcpyHostToDevice);  
255.   
256.         find\_nearest\_cluster<<<numClusterBlocks,numThreadsPerClusterBlock,clusterBlockSharedDataSize>>>(numCoords,numObjs,numClusters,deviceObjects,deviceClusters,deviceMembership,deviceIntermediates);  
257.   
258.         cudaDeviceSynchronize();  
259.   
260.         compute\_delta<<<1,numReductionThreads,reductionBlockSharedDataSize>>>(deviceIntermediates,numClusterBlocks,numReductionThreads);  
261.   
262.         cudaDeviceSynchronize();  
263.           
264.         **int** d;  
265.         cudaMemcpy(&d,deviceIntermediates,**sizeof**(**int**),cudaMemcpyDeviceToHost);  
266.         delta = (**float**)d;  
267.   
268.         cudaMemcpy(membership,deviceMembership,numObjs\***sizeof**(**int**),cudaMemcpyDeviceToHost);  
269.           
270.         **for**(**int** i=0;i<numObjs;i++)  
271.         \{  
272.             //find the array index of nestest   
273.             index = membership\[i\];  
274.             //update new cluster centers:sum of objects located within  
275.             newClusterSize\[index\]++;  
276.             **for**(**int** j=0;j<numCoords;j++)  
277.             \{  
278.                 newClusters\[j\]\[index\] += objects\[i\]\[j\];  
279.             \}  
280.         \}  
281.         //average the sum and replace old cluster centers with newClusters   
282.         **for**(**int** i=0;i<numClusters;i++)  
283.         \{  
284.             **for**(**int** j=0;j<numCoords;j++)  
285.             \{  
286.                 **if**(newClusterSize\[i\] > 0)      
287.                     dimClusters\[j\]\[i\] = newClusters\[j\]\[i\]/newClusterSize\[i\];  
288.                 newClusters\[j\]\[i\] = 0.0;//set back to 0  
289.             \}  
290.             newClusterSize\[i\] = 0 ; //set back to 0  
291.         \}  
292.         delta /= numObjs;  
293.     \}**while**( delta > threshold && loop++ < 500 );  
294.   
295.     loop\_iterations = loop + 1;  
296.       
297.     malloc2D(clusters,numClusters,numCoords,**float**);  
298.     **for**(**int** i=0;i<numClusters;i++)  
299.     \{  
300.         **for**(**int** j=0;j<numCoords;j++)  
301.         \{  
302.             clusters\[i\]\[j\] = dimClusters\[j\]\[i\];  
303.         \}  
304.     \}  
305.   
306.     cudaFree(deviceObjects) ;  
307.     cudaFree(deviceClusters);  
308.     cudaFree(deviceMembership);  
309.     cudaFree(deviceMembership);  
310.   
311.     free(dimObjects\[0\]);  
312.     free(dimObjects);  
313.     free(dimClusters\[0\]);  
314.     free(dimClusters);  
315.     free(newClusters\[0\]);  
316.     free(newClusters);  
317.     free(newClusterSize);  
318. \}  
319.   
320. **void** KMEANS::file\_read(**char** \*fn)  
321. \{  
322.   
323.     **FILE** \*infile;  
324.     **char** \*line = **new** **char**\[MAX\_CHAR\_PER\_LINE\];  
325.     **int** lineLen = MAX\_CHAR\_PER\_LINE;  
326.   
327.     filename = fn;  
328.     infile = fopen(filename,"r");  
329.     assert(infile!=NULL);  
330.     /\*find the number of objects\*/    
331.     **while**( fgets(line,lineLen,infile) )  
332.     \{  
333.         numObjs++;    
334.     \}  
335.   
336.     /\*find the dimension of each object\*/  
337.     rewind(infile);  
338.     **while**( fgets(line,lineLen,infile)!=NULL )  
339.     \{  
340.         **if**( strtok(line," \\t\\n")!=0 )     
341.         \{  
342.             **while**( strtok(NULL," \\t\\n") )     
343.                 numCoords++;  
344.             **break**;  
345.         \}  
346.     \}  
347.   
348.     /\*allocate space for object\[\]\[\] and read all objcet\*/  
349.     rewind(infile);  
350.     objects = **new** **float**\*\[numObjs\];  
351.     **for**(**int** i=0;i<numObjs;i++)  
352.     \{  
353.         objects\[i\] = **new** **float**\[numCoords\];  
354.     \}  
355.     **int** i=0;  
356.     /\*read all object\*/  
357.     **while**( fgets(line,lineLen,infile)!=NULL )  
358.     \{  
359.         **if**( strtok(line," \\t\\n") ==NULL ) **continue**;  
360.         **for**(**int** j=0;j<numCoords;j++)  
361.         \{  
362.             objects\[i\]\[j\] = atof( strtok(NULL," ,\\t\\n") )   ;  
363.         \}  
364.         i++;  
365.     \}  
366.       
367.     /\* membership: the cluster id for each data object \*/  
368.     membership = **new** **int**\[numObjs\];  
369.     assert(membership!=NULL);  
370.     **for**(**int** i=0;i<numObjs;i++)  
371.         membership\[i\] = -1;  
372.       
373. \}  
374.   
375. **int** main(**int** argc,**char** \*argv\[\])  
376. \{  
377.     KMEANS kmeans(atoi(argv\[1\]));  
378.     kmeans.file\_read(argv\[2\]);  
379.     kmeans.cuda\_kmeans();  
380.     kmeans.file\_write();  
381.     **return** 0;  
382. \}

makefile：

**\[cpp\]**  [view plain][] [copy][view plain] [![在CODE上查看代码片][CODE]][CODE_CODE] [![派生到我的代码片][ico_fork.svg]][ico_fork.svg 1]

1.  target:  
2.      nvcc cuda\_kmeans.cu  
3.      ./a.out  4 ./Image\_data/color100.txt

所有代码和文件数据在这里：[http://yunpan.cn/cKBZMPAJ8tcAs][http_yunpan.cn_cKBZMPAJ8tcAs]（提取码：9476）

**运行代码：**

![Center][]

kmeans的cuda实现代码相对复杂，在阅读的过程中可能会有困难，有问题请留言~

Author：忆之独秀

Email：leaguenew@qq.com

注明出处：[http://blog.csdn.net/lavorange/article/details/41942323][http_blog.csdn.net_lavorange_article_details_41942323]

[Link 1]: http://blog.csdn.net/lavorange/article/details/28854929
[view plain]: http://blog.csdn.net/lavorange/article/details/41942323#
[CODE]: https://code.csdn.net/assets/CODE_ico.png
[CODE_CODE]: https://code.csdn.net/snippets/595143
[ico_fork.svg]: https://code.csdn.net/assets/ico_fork.svg
[ico_fork.svg 1]: https://code.csdn.net/snippets/595143/fork
[Link 2]: http://blog.csdn.net/lavorange/article/details/43031419
[http_yunpan.cn_cKBZMPAJ8tcAs]: http://yunpan.cn/cKBZMPAJ8tcAs
[Center]: /images/20220731/933d1ffe08dd4c4eb99b80428fdb5f13.png
[http_blog.csdn.net_lavorange_article_details_41942323]: http://blog.csdn.net/lavorange/article/details/41942323