（CUDA）快速GPU开发------Device函数应用

Device函数表示的是仅仅在设备(Device)端能够使用的函数。

Device函数可以是任何的函数，这样可以通过每一个线程来运行一个Device函数来达到并行的目的。

在本文中聚焦软件开发速度，故而不讨论计算性能问题。

---

目录

适用项目

开发方法

内存的分配

Q&A

代码实例

---

适用项目

适用于开发一个明确的单任务被多次运行的工程。举例，一次性运算出n个矩阵乘法，一次性对多个序列进行排序。

开发方法

• 确定需求，并找到合适的device函数
• 确定内存空间
• 确定线程数量
• 直接让每一个线程去运行一个device函数

内存的分配

一定要从一开始就malloc出足够大小的空间出来。需要注意的地方是，每一个device函数都需要malloc出一个内存空间，在核函数中调用了n次device函数，那么就需要在GPU中开辟出n个单device所需要内存大小的显存空间。

Q&A

Q：device函数在定义了之后，还能否被主机(Host)端调用？

A：当然可以，__device__关键字之后再加一段__host__皆可，这样的描述表示该函数既是一个能被kernel调用的函数，也是一个main调用的函数。

__device__  __host__ void function(viod) {
	//implementation
}

Q：device函数需不需要private？

A：如果是单个模块内部使用的话，那么不需要，尤其是PGI编译器，会直接告诉你这是一个仅能被设备调用的函数。

代码实例

下面是一个简单的利用device函数开发的小例子：

主要功能是利用CUDA对多个线性方程组求解。实现方法很简单，device函数是一个矩阵求逆和一个矩阵乘，这样就实现了一个粗粒度的多个线性方程组求解并行方案。

/**************************************************************************************************
*FILE: solveMatrix.cpp
*DESCRIPTION: GPU求解多线性方程组
*MODULE NAME: solveMatrix
*
*CREATE DATE: 2018-11-14
*COPYRIGHT (C): nan
*AUTHOR: CHEN
*
*ENVIRONMENT: VS2015 +CUDA8.0
*
*CURRENT VERSION: V1.0, Quiltman
*
*Note:V1.0 Most of the intermediate variables use registers. If you modify the data size, the dat
           a length needs to be modified.
**************************************************************************************************/
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include <math.h>

__device__ void mvp(float *Q, float *b, int Nx, float *x);
__device__ int inverse(float *A, float *B, int length);
__global__ void MatInvKernel(const float *A, const float *B, float *I, float *L, float *X, int Nx);

int main()
{   
	/*
	**Nx方阵大小 **J 线性方程组数量
	* A：三对角矩阵
	* I：单位矩阵
	* B：b1,b2,b3,b4,b5,b.. 多个向量组成的矩阵
	* L：L1,L2,L3,L4,L5,L.. 多个标量组成的向量
	* X：x1,x2,x3,x4,x5,x.. 多个向量组成的矩阵
	*		（L1*I+A）|x1|=  | b1 |
	*		（L2*I+A）|x2|=  | b2 |
	*		（L3*I+A）|x3|=  | b3 |
	*		（L4*I+A）|x4|=  | b4 |
	*		（L5*I+A）|x..|=  | b.. |
	*/
        const int J = 4;
	const int Nx = 3;
	float dx=1;
        float A[Nx][Nx] = {0};
	float B[J][Nx] = {0};
	float I[Nx][Nx]={0};
	float X[Nx][J]={0};
	
	float *d_A,*d_B,*d_L,*d_I,*d_X;
	float *h_A,*h_B,*h_L,*h_I,*h_X;
	
	// calloc host memory
	h_A=(float *)calloc(Nx*Nx,sizeof(float));
	h_B=(float *)calloc(Nx*J,sizeof(float)); 
	h_L=(float *)calloc(J,sizeof(float));
	h_I=(float *)calloc(Nx*Nx,sizeof(float));
	h_X=(float *)calloc(Nx*J,sizeof(float));

	// malloc device global memory
	cudaMalloc((void**)&d_A,sizeof(float)*Nx*Nx);
	cudaMalloc((void**)&d_B,sizeof(float)*Nx*J);
	cudaMalloc((void**)&d_L,sizeof(float)*J);
	cudaMalloc((void**)&d_I,sizeof(float)*Nx*Nx);
	cudaMalloc((void**)&d_X,sizeof(float)*Nx*J);

	//construct I
	for(int i=0;i<Nx;i++){
		I[i][i]=1;
	}

	//construct L
	for (int i = 0; i < J; i++) {
		h_L[i] = i;
	}

	//construct the Matrix  A
	for (int i=0;i<Nx;i++){
		for(int j=0;j<Nx;j++){
			if(i==j){
				A[i][j]=2/(dx*dx);
			}
			if(1==abs(i-j)){
				A[i][j]=-1/(dx*dx);
			}
		}
	}
	
	//construct the Matrix  B
	for(int j=0;j<J;j++){
		for (int i=0;i<Nx;i++){
			 B[j][i]=i+j;
		}
	}

	// transfer data from host to device
	for(int i=0;i<Nx;i++)
	{
		memcpy(h_A+i*Nx,A[i],Nx*sizeof(float));
		memcpy(h_I+i*Nx,I[i],Nx*sizeof(float));
	}
	for(int j=0;j<J;j++)
	{
		memcpy(h_B+j*Nx,B[j],Nx*sizeof(float));
	}

	cudaMemcpy(d_A,h_A,sizeof(float)*Nx*Nx,cudaMemcpyHostToDevice);
	cudaMemcpy(d_B,h_B,sizeof(float)*Nx*J,cudaMemcpyHostToDevice);
	cudaMemcpy(d_L,h_L,sizeof(float)*J,cudaMemcpyHostToDevice);
	cudaMemcpy(d_I,h_I,sizeof(float)*Nx*Nx,cudaMemcpyHostToDevice);


        // calculate the result of the Linear equations
	cudaEvent_t start, stop;
	float elapsedTime = 0.0;

	cudaEventCreate(&start);
	cudaEventCreate(&stop);
	cudaEventRecord(start, 0);

        MatInvKernel<<<2,2>>>(d_A,d_B,d_I,d_L,d_X,Nx);

	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);

	cudaEventElapsedTime(&elapsedTime, start, stop);
	printf("time:%f\n",elapsedTime);

	 // copy kernel result back to host side
	cudaMemcpy(h_X,d_X,sizeof(float)*Nx*J,cudaMemcpyDeviceToHost);
	for(int i=0;i<Nx;i++)
	{
	      for(int j=0;j<J;j++)
	      {
              X[i][j]=h_X[j*Nx+i];		
		  }
	}
	
	printf("check the result:\n");

	for(int i=0;i<Nx;i++){
		for(int j=0;j<J;j++){
			printf("  %f  ",X[i][j]);
		}
		printf("\n");
	}

	
	// free device global memory
	cudaFree(d_A);
	cudaFree(d_B);
	cudaFree(d_L);
	cudaFree(d_I);
	cudaFree(d_X);
	
	// free host memory
	free(h_A);
	free(h_B);
	free(h_L);
	free(h_I);
	free(h_X);

	system("pause");
    return 0;
}
/************************************
* Function: 核函数
* Parameter: 	
*************************************/
__global__ void MatInvKernel(const float *A, const float *B, float *I, float *L, float *X, int Nx)
{
	
	/*
	* b[Nx]
	* x[Nx]
	* Q[Nx*Nx]
	* invQ[Nx*Nx]
	*/
	float b[3] = { 0 };
	float x[3] = { 0 };
	float Q[9] = { 0 };
	float invQ[9] = { 0 };
	
	int tid = blockIdx.x*blockDim.x + threadIdx.x;

	if (tid<4) {
		for (int i = 0; i<Nx*Nx; i++) {
			Q[i] = L[tid] * I[i] + A[i];
		}

		for (int i = 0; i<Nx; i++) {
			b[i] = B[tid*Nx + i];
		}
		
		inverse(Q, invQ, Nx);
		mvp(invQ, b, Nx, x);

		for (int i = 0; i<Nx; i++) {
			X[tid*Nx + i] = x[i];
		}
	}
}
/************************************
* Function: 矩阵乘向量
* Parameter:Q:（L1*I+A）组成的矩阵
*************************************/
__device__ void mvp(float *Q, float *b, int Nx, float *x) {
	for (int i = 0; i<Nx; i++) {
		float c = 0;
 
		for (int j = 0; j<Nx; j++) {
			c = c + Q[i*Nx + j] * b[j];
		}
		x[i] = c;
	}
}

/************************************
* Function: 矩阵求逆
* Parameter: A：输入矩阵
			 B: 输出矩阵
			 length:方阵长度
*************************************/
__device__ int inverse(float *A, float *B, int length)
{
	int i, j, k;
	float d_max, d_temp;
	float mat_t[3][3]; //temp matrix
					   //save A matrix to t[n][n]中
	for (i = 0; i < length; i++) {
		for (j = 0; j < length; j++)
		{
			mat_t[i][j] = A[i*length + j];
		}
	}
	//Initialization B matrix as identity matrix
	for (i = 0; i < length; i++)
	{
		for (j = 0; j < length; j++)
		{
			B[i*length + j] = (i == j) ? (float)1 : 0;
		}
	}
	for (i = 0; i < length; i++)
	{
		d_max = mat_t[i][i];
		k = i;
		for (j = i + 1; j < length; j++)
		{
			if (fabs(mat_t[j][i]) > fabs(d_max))
			{
				d_max = mat_t[j][i];
				k = j;
			}
		}
		if (k != i)
		{
			for (j = 0; j < length; j++)
			{
				d_temp = mat_t[i][j];
				mat_t[i][j] = mat_t[k][j];
				mat_t[k][j] = d_temp;
				d_temp = B[i*length + j];
				B[i*length + j] = B[k*length + j];
				B[k*length + j] = d_temp;
			}
		}
		if (mat_t[i][i] == 0)
		{
			//printf( "There is no inverse matrix!");
			return 0;
		}
		d_temp = mat_t[i][i];
		for (j = 0; j < length; j++)
		{
			mat_t[i][j] = mat_t[i][j] / d_temp;    
			B[i*length + j] = B[i*length + j] / d_temp; 
		}
		for (j = 0; j < length; j++)
		{
			if (j != i)     
			{
				d_temp = mat_t[j][i];
				for (k = 0; k < length; k++)  
				{
					mat_t[j][k] = mat_t[j][k] - mat_t[i][k] * d_temp;
					B[j*length + k] = B[j*length + k] - B[i*length + k] * d_temp;
				}
			}
		}
	}
	return 1;
}