1. Kmeans 步驟

常規的 Kmeans 步驟：?
1. 初始化聚類中心?
2. 迭代?
1. 計算每個樣本與聚類中心的歐式距離?
2. 根據樣本與聚類中心的歐式距離更新每個樣本的類標簽?
3. 根據類標簽更新聚類中心

本文中并行化的 Kmeans 的步驟：

初始化每個樣本的類標簽
迭代?
統計每一類的樣本和
統計每一類的樣本個數
計算每一類的聚類中心：樣本和 / 樣本個數
計算每個樣本與每個聚類中心的歐式距離
根據每個樣本與聚類中心的歐式距離更新樣本的類標簽
如下所示：

/* 樣本聚類索引的初始化*/ KmeansCUDA_Init_ObjClusterIdx<<<>>>();for (int i = 0; i < maxKmeansIter; i++) {/* 統計每一類的樣本和 */KmeansCUDA_Update_Cluster<<<>>>();/* 統計每一類的樣本個數 */KmeansCUDA_Count_objNumInCluster<<<>>>();/* 聚類中心平均 = 樣本和 / 樣本個數 */KmeansCUDA_Scale_Cluster<<<>>>();/* 計算每個樣本與每個聚類中心的歐式距離 */KmeansCUDA_distOfObjAndCluster<<<>>>();/* 根據每個樣本與聚類中心的歐式距離更新樣本的類標簽 */KmeansCUDA_Update_ObjClusterIdx<<<>>>(); }

2. 完全分解

不同的數據，可能會有不同的優化方式，所以本文首先給出本文的實驗條件：

數據?
將一個 256 * 256 * 3 的三維圖像進行塊劃分，塊的大小為 5 * 5 * 3，取塊的間隔為 1 個像素，所以總的塊數（樣本數）為（256 - 5 + 1）*（256 - 5 + 1）= 63504 個，樣本的維度為 5 * 5 * 3 = 75。聚類的個數為 80。
環境?
GPU：NVIDIA GTX980，2048 個核?
CUDA版本：6.5
第一部分對應的幾個函數的詳細分解方式如下所示：

KmeansCUDA_Init_ObjClusterIdx函數
由于只執行一次，且計算量非常小，故不作為優化的重點。

線程塊維度：256
線程格維度：(63504 + 256 - 1) / 256 = 249

每個線程負責初始化一個樣本的類標簽。

KmeansCUDA_Update_Cluster函數
如果每個類的樣本都放在連續空間，那么此函數可以使用類似于規約求和的方式就可以實現，且效率很高。但是此處每個類是分散的，所以換一種方式：

線程塊維度：（16，16）
線程格維度：((75 + 16 - 1) / 16, (63504 + 16 - 1) / 16) = (5, 3969)

每個線程負責一個樣本中的一個數據，此處要使用原子操作，因為多個線程可能同時寫一個聚類中心的數據。當然，也可以按下列的方式劃分線程塊和線程格：

線程塊維度：256
線程格維度：(63504 + 256 - 1) / 256 = 249

每個線程負責更新一個樣本，但是此時的效率不高：

第一種方式運行 14 次的時間為： 6.686ms
第二種方式運行 14 次的時間為：14.142ms

KmeansCUDA_Count_objNumInCluster函數
與KmeansCUDA_Init_ObjClusterIdx函數類似，計算量很少，每個負責處理一個樣本的計數，只有一個原子操作的加法，所以依舊采用一維的線程方式：

線程塊維度：256
線程格維度：(63504 + 256 - 1) / 256 = 249

由于會有 63504 個線程寫 80（類數）個位置，所以會存在許多沖突，所以此處對其進行優化，開辟線程數變為：

線程塊維度：1024
線程格維度：(63504 + 1024 - 1) / 1024 = 63

此時，在每個塊內申請 80 個整數大小的共享內存，先在塊內進行統計，最后再寫到全局內存，而不是像第一種方式那樣，直接寫全局內存，這樣就避免了很多的沖突，兩種方式的時間對比為：

第一種方式運行 14 次的時間為：733.3us
第二種方式運行 14 次的時間為：61.44us

從上可以看出，通過減少沖突，使用共享內存，可以加速 12 倍左右。

KmeansCUDA_Scale_Cluster函數?
此函數用于對每個樣本求和之后的取平均操作，計算量極其少，不是優化的重點，開辟二維線程塊

線程塊維度：（16，16）?
線程格維度：((75 + 16 - 1) / 16, (80 + 16 - 1) / 16) = (5, 5)

每個線程負責更新聚類中心中的一個數。

KmeansCUDA_distOfObjAndCluster函數
此函數是優化的重點，因為 Kmeans 的絕大部分計算量都集中在求每個樣本與每個聚類中心的歐式距離，在此也開辟二維線程塊：

線程塊維度：（16，16）
線程格維度：((80 + 16 - 1) / 16, (63504 + 16 - 1) / 16) = (5, 3969)

每個線程負責計算一個樣本與一個聚類中心的歐式距離。考慮到此函數的計算過程中聚類中心是不變的，可以考慮使用常量內存，以此加快對聚類中心的訪問速度（此方案稱為方案二）。另外，考慮到此函數的計算方式與矩陣乘法類似，可以考慮使用共享內存，每個線程塊的任務是計算 16 個樣本與 16 個聚類中心的距離（此方案稱為方案三）將16 個樣本與 16 個聚類中心的數據存到共享內存中（每一個線程塊的共享內存為：32 * 75 * 4 / 1024 = 9.375kb 小于 每塊最大可以使用量 48kb，每個 SM 最多調度 2048 個線程，也就是最多 8 個線程塊，所以每個 SM 使用的共享內存為 75kb，小于 GTX980 每塊 SM 的 96kb，綜上，共享內存是夠用的）。三種方案的執行14次的時間為：

方案一：183.562ms
方案二：243.298ms
方案三：25.213ms

從上可以看出，使用共享內存的計算時間遠遠小于未使用共享內存的計算時間，方案三相較于方案一加速 7 倍多，考慮到是計算量最大的模塊加速 7 倍多，所以加速效果相當理想。但是方案二使用常量內存的計算時間反而比不使用常量內存的時間還長，這是因為聚類中心的數據大小為：80 *75*4/1024 = 23.4375 kb，雖然 GTX980 的最大可用常量內存為 64 kb，但是在每個 SM 共共常量內存使用的高速緩存只有 10 kb（常量內存是借助于高速緩存來實現快速訪問，對于計算能力小于 5.0 的設備，此值為 8 kb），也就是說不夠緩存所有的聚類中心，導致反復緩存，緩存命中的效率極低。

__shared__ float objShared[BLOCKSIZE][OBJLENGTH]; // 存樣本 __shared__ float cluShared[BLOCKSIZE][OBJLENGTH]; // 存聚類中心for (int xidx = threadIdx.x; xidx < OBJLENGTH; xidx += BLOCKSIZE) {int index = myParameter.objLength * threadIdx.y + xidx;objShared[threadIdx.y][xidx] = objects[index];cluShared[threadIdx.y][xidx] = clusters[index]; }

上述代碼中將數據從全局內存拷貝到共享內存中的方式如下圖所示，每次讀同一種顏色的塊，經多次循環之后全部讀入共享內存中。

上述訪問方式會人為的將連續的內存區域分成不連續的訪問，將無法實現合并訪問，進而無法隱藏內存訪問的延遲，導致效率不高，因此對其進行優化。

__shared__ float objShared[OBJLENGTH * BLOCKSIZE]; // 存樣本 __shared__ float cluShared[OBJLENGTH * BLOCKSIZE]; // 存聚類中心for (int index = BLOCKSIZE * threadIdx.y + threadIdx.x; index < OBJLENGTH * BLOCKSIZE; index = BLOCKSIZE * BLOCKSIZE + index) {objShared[index] = objects[index];cluShared[index] = clusters[index]; }

讀取方式如下圖所示：

優化前后的時間對比為：

優化前：25.213ms
優化后：24.173ms

雖然效率沒有明顯提升，但是還是很可觀的。

KmeansCUDA_Update_ObjClusterIdx函數
用于尋找每個樣本與所有聚類中心中距離最小的，將當前樣本劃歸到該類。可以開辟一維線程：

線程塊維度：256
線程格維度：(63504 + 256 - 1) / 256 = 249

每個線程用于查找當前樣本對應的最短距離（共 80 個）。因為求最小與規約類似，而上述方式的每個樣本卻是完全串行的方式，所以對其進行優化，開辟二維線程：

線程塊維度：（16，16）
線程格維度：((1, (63504 + 16 - 1) / 16) = (1, 3969)

每個線程塊用于計算 16 個樣本的最小距離，也就是說用 16 個線程來計算原來 1 個線程的工作（查找 80 個數的最小值）。

最后，把剩余兩個元素的較小值作為最短距離，優化前后的時間對比為：

優化前：12.725ms
優化后：1.887ms

優化后加速 7 倍左右，效果相當可觀。

3. 結果

如下圖所示，是優化的最終版本的 visual profile 結果，可見，依舊是計算樣本與聚類中心的距離最耗時，也必然還有優化的空間。

另外，本文對 Matlab 和 GPU 上的運行時間進行對比，以衡量加速效果（采用相同數據，相同的迭代次數）：

GPU 運行時間為：32.1842 ms
Matlab 運行時間為：11102.919 ms
最終的加速比為：345 倍

下圖所示的為聚類結果，按每一類中的樣本數從大到小排列，由于初始聚類中心的不同，有部份計算誤差。

代碼：

kernel.cu

#include "KmeansCUDA.h"#include <cuda_runtime.h> #include <device_launch_parameters.h>#include "ClassParameter.h" #include "ReadSaveImage.h"#include <iostream>using std::cout; using std::endl;int main() {sParameter myParameter{63504, 75, 80, 40, 150, 14};float *objData = (float*)malloc(myParameter.objNum * myParameter.objLength * sizeof(float));float *centerData = (float*)malloc(myParameter.clusterNum * myParameter.objLength * sizeof(float));int *objClassIdx = (int*)malloc(myParameter.objNum * sizeof(int));ReadData(objData, myParameter);KmeansCUDA(objData, objClassIdx, centerData, myParameter);SaveData(objClassIdx, myParameter);cudaDeviceReset();return 0; }

ReadSaveImage.h

#ifndef READSAVEIMAGE_H #define READSAVEIMAGE_H#include "ClassParameter.h"void ReadData(float *fileData, sParameter myParameter); void SaveData(int *fileData, sParameter myParameter);#endif // !READSAVEIMAGE_H

ReadSaveImage.cpp

#include "ReadSaveImage.h" #include <iostream> #include <fstream>/** * 功能：從txt文件中讀取數據 * 輸出：fileData 輸出數據的頭指針 * 輸入：fileName 讀取的文本文件的文件名 * 輸入：dataNum 讀取的數據個數 */ void ReadFile(float *fileData, std::string fileName, int dataNum) {std::fstream file;file.open(fileName, std::ios::in);if (!file.is_open()){std::cout << "不能打開文件" << std::endl;return;}// 讀入數據到內存中for (int i = 0; i < dataNum; i++) file >> fileData[i];file.close(); }void ReadData(float *fileData, sParameter myParameter) {std::string str = "D:\\Document\\vidpic\\CUDA\\Kmeans\\objData.txt"; //文件路徑ReadFile(fileData, str, myParameter.objNum*myParameter.objLength); }/** * 功能：把數據保存到文本文件中 * 輸入：fileData 輸入數據的頭指針 * 輸入：fileName 保存的文本文件的文件名 * 輸入：dataNum 保存的數據個數 */ void SaveFile(int *fileData, std::string fileName, int dataNum) {std::fstream file;file.open(fileName, std::ios::out);if (!file.is_open()){std::cout << "不能打開文件" << std::endl;return;}// 讀入數據到內存中for (int i = 0; i < dataNum; i++) file << fileData[i] << std::endl;file.close(); }void SaveData(int *fileData, sParameter myParameter) {std::string str = "D:\\Document\\vidpic\\CUDA\\Kmeans\\objClusterIdxPre.txt"; //文件路徑SaveFile(fileData, str, myParameter.objNum); }

KmeansCUDA.h

#ifndef KMEANSCUDA_H #define KMEANSCUDA_H#include "ClassParameter.h"/** * 功能：并行 Kmeans 聚類 * 輸入：objData 樣本數據 * 輸出：objClusterIdx 每個樣本的類別索引 * 輸出：clusterData 聚類中心 * 輸入：myPatameter 輸入參數 */ void KmeansCUDA(float *objData, int *objClusterIdx, float*clusterData, sParameter myParameter);#endif // KMEANSCUDA_H

KmeansCUDA.cu

#include "KmeansCUDA.h"#include <cuda_runtime.h> #include <device_launch_parameters.h> #include <curand_kernel.h>#include <iostream>#include <stdlib.h>#define BLOCKSIZE_16 16 #define BLOCKSIZE_32 32 #define OBJLENGTH 75/** * 功能：初始化每個樣本的類索引 * 輸出：objClusterIdx_Dev 每個樣本的類別索引 * 輸入：objNum 樣本個數 * 輸入：maxIdx 索引的最大值 */ __global__ void KmeansCUDA_Init_ObjClusterIdx(int *objClusterIdx_Dev, int objNum, int maxIdx) {int index = blockDim.x * blockIdx.x + threadIdx.x; curandState s;curand_init(index, 0, 0, &s);if (index < objNum) objClusterIdx_Dev[index] = (int(curand_uniform(&s) * maxIdx)); }/** * 功能：更新 Kmeans 的聚類中心 * 輸入：objData_Dev 樣本數據 * 輸入：objClusterIdx_Dev 每個樣本的類別索引 * 輸出：clusterData_Dev 聚類中心 * 輸入：myPatameter 輸入參數 */ __global__ void KmeansCUDA_Update_Cluster(float *objData_Dev, int *objClusterIdx_Dev, float *clusterData_Dev, sParameter myParameter) {int x_id = blockDim.x * blockIdx.x + threadIdx.x; // 列坐標int y_id = blockDim.y * blockIdx.y + threadIdx.y; // 行坐標if (x_id < myParameter.objLength && y_id < myParameter.objNum){int index = y_id * myParameter.objLength + x_id;int clusterIdx = objClusterIdx_Dev[y_id];atomicAdd(&clusterData_Dev[clusterIdx * myParameter.objLength + x_id], objData_Dev[index]);} }/** *功能：更新 Kmeans 的聚類中心 * 輸入：objClusterIdx_Dev 每個樣本的類別索引 * 輸出：objNumInCluster 每個聚類中的樣本數 * 輸入：myPatameter 輸入參數 */ __global__ void KmeansCUDA_Count_objNumInCluster(int *objClusterIdx_Dev, int *objNumInCluster, sParameter myParameter) {int index = blockDim.x * blockIdx.x + threadIdx.x;if (index < myParameter.objNum){int clusterIdx = objClusterIdx_Dev[index];atomicAdd((int*)&objNumInCluster[clusterIdx], 1); // 計數} }/** *功能：更新 Kmeans 的聚類中心 * 輸入：objClusterIdx_Dev 每個樣本的類別索引 * 輸出：objNumInCluster 每個聚類中的樣本數 * 輸入：myPatameter 輸入參數 */ __global__ void KmeansCUDA_Count_objNumInCluster1(int *objClusterIdx_Dev, int *objNumInCluster, sParameter myParameter) {int index = blockDim.x * blockIdx.x + threadIdx.x;__shared__ int sData[80];if (threadIdx.x < myParameter.clusterNum)sData[threadIdx.x] = 0;__syncthreads();if (index < myParameter.objNum){int clusterIdx = objClusterIdx_Dev[index];atomicAdd((int*)&sData[clusterIdx], 1);}__syncthreads();if (threadIdx.x < myParameter.clusterNum)atomicAdd((int*)&objNumInCluster[threadIdx.x], sData[threadIdx.x]); // 計數 }/** *功能：平均 Kmeans 的聚類中心 * 輸出：clusterData_Dev 聚類中心 * 輸出：objNumInCluster 每個聚類中的樣本數 * 輸入：myPatameter 輸入參數 */ __global__ void KmeansCUDA_Scale_Cluster(float *clusterData_Dev, int *objNumInCluster, sParameter myParameter) {int x_id = blockDim.x * blockIdx.x + threadIdx.x; // 列坐標int y_id = blockDim.y * blockIdx.y + threadIdx.y; // 行坐標if (x_id < myParameter.objLength && y_id < myParameter.clusterNum){int index = y_id * myParameter.objLength + x_id;clusterData_Dev[index] /= float(objNumInCluster[y_id]);} }/** * 功能：計算兩個向量的歐拉距離 * 輸入：objects 樣本數據 * 輸出：clusters 聚類中心數據 * 輸入：objLength 樣本長度 */ __device__ inline static float EuclidDistance(float *objects, float *clusters, int objLength) {float dist = 0.0f;for (int i = 0; i < objLength; i++){float onePoint = objects[i] - clusters[i];dist = onePoint * onePoint + dist;}return(dist); }/** * 功能：計算所有樣本與聚類中心的歐式距離 * 輸入：objData_Dev 樣本數據 * 輸入：objClusterIdx_Dev 每個樣本的類別索引 * 輸入：clusterData_Dev 聚類中心 * 輸出：distOfObjAndCluster_Dev 每個樣本與聚類中心的歐式距離 * 輸入：objNumInCluster_Dev 每個聚類中的樣本數 * 輸入：iter 迭代次數 * 輸入：myPatameter 輸入參數 */ __global__ void KmeansCUDA_distOfObjAndCluster(float *objData_Dev, int *objClusterIdx_Dev, float *clusterData_Dev, float *distOfObjAndCluster_Dev, int *objNumInCluster_Dev, int iter, sParameter myParameter) {int x_id = blockDim.x * blockIdx.x + threadIdx.x; // 列坐標int y_id = blockDim.y * blockIdx.y + threadIdx.y; // 行坐標const int oneBlockData = OBJLENGTH * BLOCKSIZE_16;__shared__ float objShared[oneBlockData]; // 存樣本__shared__ float cluShared[oneBlockData]; // 存聚類中心/* 數據讀入共享內存 */if (y_id < myParameter.objNum){float *objects = &objData_Dev[myParameter.objLength * blockDim.y * blockIdx.y]; // 當前塊需要樣本對應的首地址float *clusters = &clusterData_Dev[myParameter.objLength * blockDim.x * blockIdx.x]; // 當前塊需要聚類中心對應的首地址for (int index = BLOCKSIZE_16 * threadIdx.y + threadIdx.x; index < oneBlockData; index = BLOCKSIZE_16 * BLOCKSIZE_16 + index){objShared[index] = objects[index];cluShared[index] = clusters[index];}__syncthreads();}if (x_id < myParameter.clusterNum && y_id < myParameter.objNum){//if (objNumInCluster_Dev[x_id] < myParameter.minObjInClusterNum && iter >= myParameter.maxKmeansIter - 2)// distOfObjAndCluster_Dev[y_id * myParameter.clusterNum + x_id] = 3e30;//elsedistOfObjAndCluster_Dev[y_id * myParameter.clusterNum + x_id] = EuclidDistance(&objShared[myParameter.objLength * threadIdx.y], &cluShared[myParameter.objLength * threadIdx.x], myParameter.objLength);} }/** * 功能：計算所有樣本與聚類中心的歐式距離 * 輸入：objData_Dev 樣本數據 * 輸入：objClusterIdx_Dev 每個樣本的類別索引 * 輸入：clusterData_Dev 聚類中心 * 輸出：distOfObjAndCluster_Dev 每個樣本與聚類中心的歐式距離 * 輸入：objNumInCluster_Dev 每個聚類中的樣本數 * 輸入：iter 迭代次數 * 輸入：myPatameter 輸入參數 */ __global__ void KmeansCUDA_distOfObjAndCluster1(float *objData_Dev, int *objClusterIdx_Dev, float *clusterData_Dev, float *distOfObjAndCluster_Dev, int *objNumInCluster_Dev, int iter, sParameter myParameter) {int x_id = blockDim.x * blockIdx.x + threadIdx.x; // 列坐標int y_id = blockDim.y * blockIdx.y + threadIdx.y; // 行坐標__shared__ float objShared[BLOCKSIZE_16][OBJLENGTH]; // 存樣本__shared__ float cluShared[BLOCKSIZE_16][OBJLENGTH]; // 存聚類中心float *objects = &objData_Dev[myParameter.objLength * blockDim.y * blockIdx.y]; // 當前塊需要樣本對應的首地址float *clusters = &clusterData_Dev[myParameter.objLength * blockDim.x * blockIdx.x]; // 當前塊需要聚類中心對應的首地址/* 數據讀入共享內存 */if (y_id < myParameter.objNum){for (int xidx = threadIdx.x; xidx < OBJLENGTH; xidx += BLOCKSIZE_16){int index = myParameter.objLength * threadIdx.y + xidx;objShared[threadIdx.y][xidx] = objects[index];cluShared[threadIdx.y][xidx] = clusters[index];}__syncthreads();}if (x_id < myParameter.clusterNum && y_id < myParameter.objNum){if (objNumInCluster_Dev[x_id] < myParameter.minObjInClusterNum && iter >= myParameter.maxKmeansIter - 2)distOfObjAndCluster_Dev[y_id * myParameter.clusterNum + x_id] = 3e30;elsedistOfObjAndCluster_Dev[y_id * myParameter.clusterNum + x_id] = EuclidDistance(objShared[threadIdx.y], cluShared[threadIdx.x], myParameter.objLength);} }/** * 功能：計算所有樣本與聚類中心的歐式距離 * 輸出：objClusterIdx_Dev 每個樣本的類別索引 * 輸入：distOfObjAndCluster_Dev 每個樣本與聚類中心的歐式距離 * 輸入：myPatameter 輸入參數 */ __global__ void KmeansCUDA_Update_ObjClusterIdx1(int *objClusterIdx_Dev, float *distOfObjAndCluster_Dev, sParameter myParameter) {int index = blockDim.x * blockIdx.x + threadIdx.x;if (index < myParameter.objNum){float *objIndex = &distOfObjAndCluster_Dev[index * myParameter.clusterNum];int idx = 0;float dist = objIndex[0];for (int i = 1; i < myParameter.clusterNum; i++){if (dist > objIndex[i]){dist = objIndex[i];idx = i;}}objClusterIdx_Dev[index] = idx;} }/** * 功能：計算所有樣本與聚類中心的歐式距離（優化后的） * 輸出：objClusterIdx_Dev 每個樣本的類別索引 * 輸入：distOfObjAndCluster_Dev 每個樣本與聚類中心的歐式距離 * 輸入：myPatameter 輸入參數 */ __global__ void KmeansCUDA_Update_ObjClusterIdx(int *objClusterIdx_Dev, float *distOfObjAndCluster_Dev, sParameter myParameter) {int y_id = blockDim.y * blockIdx.y + threadIdx.y; // 行坐標__shared__ float sData[BLOCKSIZE_16][BLOCKSIZE_16]; // 樣本與聚類中心距離__shared__ int sIndx[BLOCKSIZE_16][BLOCKSIZE_16]; // 距離對應的類標號sData[threadIdx.y][threadIdx.x] = 2e30;sIndx[threadIdx.y][threadIdx.x] = 0;__syncthreads();if (y_id < myParameter.objNum){float *objIndex = &distOfObjAndCluster_Dev[y_id * myParameter.clusterNum];sData[threadIdx.y][threadIdx.x] = objIndex[threadIdx.x];sIndx[threadIdx.y][threadIdx.x] = threadIdx.x;__syncthreads();/* 每 BLOCKSIZE_16 個進行比較 */for (int index = threadIdx.x + BLOCKSIZE_16; index < myParameter.clusterNum; index += BLOCKSIZE_16){float nextData = objIndex[index];if (sData[threadIdx.y][threadIdx.x] > nextData){sData[threadIdx.y][threadIdx.x] = nextData;sIndx[threadIdx.y][threadIdx.x] = index;}}/* BLOCKSIZE_16 個的內部規約，到只剩 2 個 */for (int step = BLOCKSIZE_16 / 2; step > 1; step = step >> 1){int idxStep = threadIdx.x + step;if (threadIdx.x < step && sData[threadIdx.y][threadIdx.x] > sData[threadIdx.y][idxStep]){sData[threadIdx.y][threadIdx.x] = sData[threadIdx.y][idxStep];sIndx[threadIdx.y][threadIdx.x] = sIndx[threadIdx.y][idxStep];}//__syncthreads();}if (threadIdx.x == 0){objClusterIdx_Dev[y_id] = sData[threadIdx.y][0] < sData[threadIdx.y][1] ? sIndx[threadIdx.y][0] : sIndx[threadIdx.y][1];}} }/** * 功能：并行 Kmeans 聚類 * 輸入：objData_Host 樣本數據 * 輸出：objClassIdx_Host 每個樣本的類別索引 * 輸出：centerData_Host 聚類中心 * 輸入：myPatameter 輸入參數 */ void KmeansCUDA(float *objData_Host, int *objClassIdx_Host, float*centerData_Host, sParameter myParameter) {/* 開辟設備端內存 */float *objData_Dev, *centerData_Dev;cudaMalloc((void**)&objData_Dev, myParameter.objNum * myParameter.objLength * sizeof(float));cudaMalloc((void**)&centerData_Dev, myParameter.clusterNum * myParameter.objLength * sizeof(float));cudaMemcpy(objData_Dev, objData_Host, myParameter.objNum * myParameter.objLength * sizeof(float), cudaMemcpyHostToDevice);int *objClassIdx_Dev;cudaMalloc((void**)&objClassIdx_Dev, myParameter.objNum * sizeof(int));float *distOfObjAndCluster_Dev; // 每個樣本與聚類中心的歐式距離cudaMalloc((void**)&distOfObjAndCluster_Dev, myParameter.objNum * myParameter.clusterNum * sizeof(float));int *objNumInCluster_Dev; // 每個聚類中的樣本數cudaMalloc((void**)&objNumInCluster_Dev, myParameter.clusterNum * sizeof(int));/* 線程塊和線程格 */dim3 dimBlock1D_16(BLOCKSIZE_16 * BLOCKSIZE_16);dim3 dimBlock1D_32(BLOCKSIZE_32 * BLOCKSIZE_32);dim3 dimGrid1D_16((myParameter.objNum + BLOCKSIZE_16 * BLOCKSIZE_16 - 1) / dimBlock1D_16.x);dim3 dimGrid1D_32((myParameter.objNum + BLOCKSIZE_32 * BLOCKSIZE_32 - 1) / dimBlock1D_32.x);dim3 dimBlock2D(BLOCKSIZE_16, BLOCKSIZE_16);dim3 dimGrid2D_Cluster((myParameter.objLength + BLOCKSIZE_16 - 1) / dimBlock2D.x, (myParameter.clusterNum + BLOCKSIZE_16 - 1) / dimBlock2D.y);dim3 dimGrid2D_ObjNum_Objlen((myParameter.objLength + BLOCKSIZE_16 - 1) / dimBlock2D.x, (myParameter.objNum + BLOCKSIZE_16 - 1) / dimBlock2D.y);dim3 dimGrid2D_ObjCluster((myParameter.clusterNum + BLOCKSIZE_16 - 1) / dimBlock2D.x, (myParameter.objNum + BLOCKSIZE_16 - 1) / dimBlock2D.y);dim3 dimGrid2D_ObjNum_BLOCKSIZE_16(1, (myParameter.objNum + BLOCKSIZE_16 - 1) / dimBlock2D.y);// 記錄時間cudaEvent_t start_GPU, end_GPU;float elaspsedTime;cudaEventCreate(&start_GPU);cudaEventCreate(&end_GPU);cudaEventRecord(start_GPU, 0);/* 樣本聚類索引的初始化*/KmeansCUDA_Init_ObjClusterIdx<<<dimGrid1D_16, dimBlock1D_16>>>(objClassIdx_Dev, myParameter.objNum, myParameter.clusterNum);for (int i = 0; i < myParameter.maxKmeansIter; i++){cudaMemset(centerData_Dev, 0, myParameter.clusterNum * myParameter.objLength * sizeof(float));cudaMemset(objNumInCluster_Dev, 0, myParameter.clusterNum * sizeof(int));/* 統計每一類的樣本和 */KmeansCUDA_Update_Cluster<<<dimGrid2D_ObjNum_Objlen, dimBlock2D>>>(objData_Dev, objClassIdx_Dev, centerData_Dev, myParameter);/* 統計每一類的樣本個數 *///KmeansCUDA_Count_objNumInCluster1<<<dimGrid1D_16, dimBlock1D_16>>>(objClassIdx_Dev, objNumInCluster_Dev, myParameter);KmeansCUDA_Count_objNumInCluster<<<dimGrid1D_32, dimBlock1D_32>>>(objClassIdx_Dev, objNumInCluster_Dev, myParameter);/* 聚類中心平均 = 樣本和 / 樣本個數 */KmeansCUDA_Scale_Cluster<<<dimGrid2D_Cluster, dimBlock2D>>>(centerData_Dev, objNumInCluster_Dev, myParameter);/* 計算每個樣本與每個聚類中心的歐式距離 */KmeansCUDA_distOfObjAndCluster<<<dimGrid2D_ObjCluster, dimBlock2D>>>(objData_Dev, objClassIdx_Dev, centerData_Dev, distOfObjAndCluster_Dev, objNumInCluster_Dev, i, myParameter);/* 根據每個樣本與聚類中心的歐式距離更新樣本的類標簽 *///KmeansCUDA_Update_ObjClusterIdx1<<<dimGrid1D_16, dimBlock1D_16>>>(objClassIdx_Dev, distOfObjAndCluster_Dev, myParameter);KmeansCUDA_Update_ObjClusterIdx<<<dimGrid2D_ObjNum_BLOCKSIZE_16, dimBlock2D>>>(objClassIdx_Dev, distOfObjAndCluster_Dev, myParameter);}// 計時結束cudaEventRecord(end_GPU, 0);cudaEventSynchronize(end_GPU);cudaEventElapsedTime(&elaspsedTime, start_GPU, end_GPU);std::cout << "Kmeans 的運行時間為：" << elaspsedTime << "ms." << std::endl;/* 輸出從設備端拷貝到內存 */cudaMemcpy(objClassIdx_Host, objClassIdx_Dev, myParameter.objNum * sizeof(int), cudaMemcpyDeviceToHost);cudaMemcpy(centerData_Host, centerData_Dev, myParameter.objNum * myParameter.objLength * sizeof(float), cudaMemcpyDeviceToHost);/* 釋放設備端內存 */cudaFree(objData_Dev);cudaFree(objClassIdx_Dev);cudaFree(centerData_Dev);cudaFree(distOfObjAndCluster_Dev);cudaFree(objNumInCluster_Dev); }

ClassParameter.h

#ifndef CLASSPARAMETER_H #define CLASSPARAMETER_H// 參數 class sParameter { public:int objNum; // 樣本數int objLength; // 樣本維度int clusterNum; // 聚類數int minClusterNum; // 最少的聚類數int minObjInClusterNum; // 每個聚類中的最少樣本數int maxKmeansIter; // 最大迭代次數 };#endif // !CLASSPARAMETER_H

轉載至：https://github.com/muzichao/Kmeans_CUDA

總結

以上是生活随笔為你收集整理的Kmeans CUDA的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。