利用CUDA的矩阵乘法1 <利用 Kahan's Summation Formula 来提高CUDA 的浮点数运算精确度>

//矩阵乘法

#include<stdio.h>

#include<stdlib.h>

#include<cuda_runtime.h>

#include<math.h>

#include<time.h>

 

#define NUM_THREADS 256

 

bool InitCUDA()

{

   int count;

   cudaGetDeviceCount(&count);

   if(count == 0)

   {

      fprintf(stderr,"There is no device.\n");

      return false;

   }

   int i;

   for(i = 0;i < count;i++)

   {

      cudaDeviceProp prop;

      if(cudaGetDeviceProperties(&prop, i) == cudaSuccess)

      {

        if(prop.major >= 1)

           break;

      }

   }

   if(i == count)

   {

      fprintf(stderr,"There is no device supporting CUDA 1.X.\n");

      return false;

   }

   cudaSetDevice(i);

   return true;

}

 

/*利用随机数生成器把矩阵填满 0 ~ 1 之间的数字。

  特别注意到因为 C 语言中无法声明变动大小的二维矩阵，

  所以我们使用 i * lda + j 的方式*/

//产生矩阵

void matgen(float *a,int lda,int n)

{

   int i,j;

   for(i = 0;i < n;i++)

   {

      for(j = 0;j < n;j++)

        a[i*lda + j] = (float)rand()/RAND_MAX +

           (float)rand()/(RAND_MAX*RAND_MAX);

   }

}

 

//矩阵乘法

void matmult(const float* a,int lda,const float* b,int ldb,

            float* c,int ldc,int n)

{

   int i,j,k;

   for(i = 0;i < n;i++)

   {

      for(j = 0;j < n;j++)

      {

        double t = 0;

        for(k = 0;k < n;k++)

           t += a[i*lda + k]*b[k*ldb + j];

        c[i*ldc + j] = t;

      }

   }

}

 

/*这是以 CPU 进行矩阵乘法、用来进行验证答案正确与否的程序。

 

  特别注意到它用double 来储存暂时的计算结果，以提高精确度*/

 

//验证结果：

//这个函式计算两个矩阵的最大相对误差和平均相对误差，并把结果印出来

void compare_mat(const float* a,int lda,

                const float* b,int ldb,int n)

 

{

 

   float max_err = 0;

   float average_err = 0;

   int i,j;

   for(i = 0;i < n;i++)

   {

     for(j = 0;j < n;j++)

      {

        if(b[i*ldb + j] != 0)

        {

           float err = fabs((a[i*lda + j] -

                      b[i*ldb + j])/b[i*ldb + j]);

           if(max_err < err)

            max_err = err;

           average_err += err;

        }

      }

   }

   printf("Max error: %g Average error: %g\n",

         max_err,average_err/(n*n));

 

}

 

//进行计算的 kernel

/*一开始先从 bid 和 tid 计算出这个 thread 应该计算的 row 和 column,在判断

row 和 column 在范围内之后，就直接进行计算，并把结果写到 c 矩阵中*/

__global__ static void matMultCUDA(const float* a,size_t lda,

  const float* b,size_t ldb,float* c,size_t ldc,int n)

{

   const int tid = threadIdx.x;

   const int bid = blockIdx.x;

   const int idx = bid*blockDim.x + tid;

   const int row = idx/n;

   const int column = idx%n;

   int i;

 

   if(row < n && column < n)

   {

      float t = 0;

      float y = 0;

      for(i = 0;i < n;i++)

      {

        float r;

           y -= a[row*lda + i]*b[i*ldb + column];

        r = t - y;

        y = (r - t) + y;

        t = r;

      }

      //c[row*ldc + column] = t;

   }

}

 

//CUDA 的矩阵乘法的部份

clock_t matmultCUDA(const float* a,int lda,const float*b,

        int ldb,float* c,int ldc,int n)

{

   float *ac,*bc,*cc;

   clock_t start,end;

 

   start = clock();

   cudaMalloc((void**)&ac,sizeof(float)*n*n);

   cudaMalloc((void**)&bc,sizeof(float)*n*n);

   cudaMalloc((void**)&cc,sizeof(float)*n*n);

   cudaMemcpy2D(ac,sizeof(float)*n,a,sizeof(float)*lda,

           sizeof(float)*n,n,cudaMemcpyHostToDevice);

   cudaMemcpy2D(bc,sizeof(float)*n,b,sizeof(float)*ldb,

           sizeof(float)*n,n,cudaMemcpyHostToDevice);

   int blocks = (n + NUM_THREADS - 1)/NUM_THREADS;

   matMultCUDA<<<blocks*n,NUM_THREADS>>>

               (ac,n,bc,n,cc,n,n);

   cudaMemcpy2D(c,sizeof(float)*ldc,cc,sizeof(float)*n,

            sizeof(float)*n,n,cudaMemcpyDeviceToHost);

 

   cudaFree(ac);

   cudaFree(bc);

   cudaFree(cc);

   end = clock();

   return end - start;

}

 

/*就是在显卡内存中配置存放矩阵的内存，然后把主内存中的矩阵数据复

 

  制到显卡内存上*/

int main()

{

   float *a,*b,*c,*d;

   int n = 1000;

   if(!InitCUDA())

      return 0;

   a = (float*)malloc(sizeof(float)*n*n);

   b = (float*)malloc(sizeof(float)*n*n);

   c = (float*)malloc(sizeof(float)*n*n);

   d = (float*)malloc(sizeof(float)*n*n);

   srand(0);

   matgen(a,n,n);

   matgen(b,n,n);

   clock_t time = matmultCUDA(a,n,b,n,c,n,n);

   matmult(a,n,b,n,d,n,n);

   compare_mat(c,n,d,n,n);

 

   double sec = (double)time/CLOCKS_PER_SEC;

   printf("Time used:%.2f(%.21f GFLOPS)\n",sec,

         2.0*n*n*n/(sec*1E9));

   int i;

   scanf("%d",&i);

   return 0;

}