（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 data 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 memoryh_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 memorycudaMalloc((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 Ifor(int i=0;i<Nx;i++){I[i][i]=1;}//construct Lfor (int i = 0; i < J; i++) {h_L[i] = i;}//construct the Matrix  Afor (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  Bfor(int j=0;j<J;j++){for (int i=0;i<Nx;i++){B[j][i]=i+j;}}// transfer data from host to devicefor(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 equationscudaEvent_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 sidecudaMemcpy(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 memorycudaFree(d_A);cudaFree(d_B);cudaFree(d_L);cudaFree(d_I);cudaFree(d_X);// free host memoryfree(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 matrixfor (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;
}