cuda程式設計與gpu平行計算（六）：圖稀疏矩陣轉為CSR結構並傳入gpu

雖然sepgraph有這部分程式碼，還是自己先試著實現一下，這樣讀起來也方便

行壓縮格式Compressed Sparse Row (CSR)
CSR需要三種資料來表達：數值、列號、行轉移。CSR不是三元組，而是整體的編碼方式。

CSR

編碼：
行優先遍歷矩陣Matrix
values陣列中儲存矩陣中非零元素。
column indices陣列儲存values陣列中對應位置非零元素的列索引。
row offsets陣列的下標表示每一行第一個非零元素的行索引，元素值為values陣列的下標，最後一個元素值為非零元素的個數。

解碼:
遍歷values陣列，對其中的元素值x（下標記為index_X）去column indices陣列中對應位置取出x在原矩陣中的列索引（記為col_index），然後在row offsets陣列查詢index_X，如果index_X在row offsets中，其在row offsets中的下標即為x在原矩陣中的行索引（記為row_index），如果index_X不在row offsets中，則使用上一個row_index值（因為編碼時是逐行編碼，且row offsets僅儲存每行的第一個非零元素，所以當前x與上一x在同一行）。

此為前言，清楚邏輯，程式碼就好辦了

#include<stdlib.h>
#include<time.h>
#include<iostream>
#include<stdio.h>
#include<string.h>
using namespace std;


__global__ void csr(int* d_in_1,int* d_in_2,float* d_in_3){
  
        int i = blockIdx.x * blockDim.x + threadIdx.x; 

	int start = d_in_1[i];
	int end = d_in_1[i + 1];
	
	for(int j = start;j < end ;j++){
		printf("(%d,%d) ",i,d_in_2[j]);
 	}

}


//隨機生成一個圖的鄰接矩陣
void fill_random(float*data, int m, int n) {
	srand((unsigned)(time(NULL))); //每次生成的隨機數不一樣
	for (int i = 0; i < m*n; i++) {
		data[i] = rand() % 100;  //生成100以內的隨機數
		if (data[i] < 80) 
			data[i] = 0;
	}
}

void print_matrix(float*data, int m, int n) {
	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) {
			int l = i + j * m;
			cout << data[l] << " ";
		}
		cout << endl;
	}
	cout << endl;
}

void print_matrix(int*data, int m, int n) {
	for (int i = 0; i < m; i++) {
		for (int j = 0; j < n; j++) {
			int l = i + j * m;
			cout << data[l] << " ";
		}
		cout << endl;
	}

	cout << endl;
}

void dense2csr(float*data, int*&rowPtr, int*&colInd, float*&val, int m, int n) {
	rowPtr = (int*)malloc(sizeof(int)*(m + 1));

	int* tcolInd = (int*)malloc(sizeof(int)*(m *n));
	float* tval = (float*)malloc(sizeof(float)*(m *n));
	int towtal = m * n;
	int nnv = 0;

	for (int i = 0; i < m; i++) {
		rowPtr[i] = nnv;//記錄行偏移，其實也是前面一共已經有多少邊
		for (int j = 0; j < n; j++) {
			int l = i + j * m;
			if (data[l] != 0) {
				tcolInd[nnv] = j;//記錄列索引
				tval[nnv] = data[l];//記錄邊權
				nnv++;//找到一個邊，行偏移加一
			}
		}
	}
	rowPtr[m] = nnv;

	colInd = (int*)malloc(sizeof(int)*(nnv));
	val = (float*)malloc(sizeof(float)*(nnv));

	memcpy(colInd, tcolInd, sizeof(float)*nnv);
	memcpy(val, tval, sizeof(float)*nnv);

	free(tcolInd);
	free(tval);
}



int main() {
	int m = 5;
	int n = 5;
	float*A = (float*)malloc(sizeof(float)*m*n);

	fill_random(A, m, n);
	print_matrix(A, m, n);

	int*csrRowPtr;
	int*csrColInd;
	float*csrVal;
	dense2csr(A, csrRowPtr, csrColInd, csrVal, m, n);

        print_matrix(csrRowPtr, 1, m + 1);
	print_matrix(csrColInd, 1, csrRowPtr[m]);
	print_matrix(csrVal, 1, csrRowPtr[m]);


	const int ARRAY_SIZE = m + 1;
  	const int ARRAY_BYTES_1 = ARRAY_SIZE * sizeof(int);
	const int ARRAY_BYTES_2 = csrRowPtr[ARRAY_SIZE - 1] * sizeof(int);
	const int ARRAY_BYTES_3 = csrRowPtr[ARRAY_SIZE - 1] * sizeof(float);

	// 生成cpu上的csr結構
	int h_in_1[m + 1];
	int h_in_2[csrRowPtr[m]];
	float h_in_3[csrRowPtr[m]];
	for(int i = 0;i < m + 1 ; i++){
    	        h_in_1[i] = csrRowPtr[i];
  	}
	for(int i = 0;i < csrRowPtr[m]; i++){
		h_in_2[i] = csrColInd[i];
		h_in_3[i] = csrVal[i];
	}

  	// 宣告GPU上的CSR結構
  	int* d_in_1;
  	int* d_in_2;
	float* d_in_3;

  	// 分配空間
  	cudaMalloc((void**) &d_in_1,ARRAY_BYTES_1);
  	cudaMalloc((void**) &d_in_2,ARRAY_BYTES_2);
	cudaMalloc((void**) &d_in_3,ARRAY_BYTES_3);

  	// 轉移全部
  	cudaMemcpy(d_in_1,h_in_1,ARRAY_BYTES_1,cudaMemcpyHostToDevice);
	cudaMemcpy(d_in_2,h_in_2,ARRAY_BYTES_2,cudaMemcpyHostToDevice);
	cudaMemcpy(d_in_3,h_in_3,ARRAY_BYTES_3,cudaMemcpyHostToDevice);

  	// 操作全部
  	csr<<<1,ARRAY_SIZE>>>(d_in_1,d_in_2,d_in_3);


	// 釋放記憶體
  	cudaFree(d_in_1);
	cudaFree(d_in_2);
	cudaFree(d_in_3);

	return 0;
}

執行的結果如下

先輸出一個隨機生成的帶權重有向稀疏圖，在cpu中將其壓縮為csr格式並輸出，再傳遞到gpu上，並輸出每一條邊的起點和終點。