cuda编程入门示例23

#include <stdio.h>#include <stdlib.h>#include <time.h>#include <cuda.h>#include <cuda_runtime.h>#define BLOCK_SIZE 16static void HandleError(cudaError_t err, const char *file, int line){if (err != cudaSuccess){printf("%s in %s at line %d\n", cudaGetErrorString(err), file, line);exit(EXIT_FAILURE);}}#define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ ))#define HANDLE_NULL( a ) {if ((a) == NULL) { \printf("Host memory failed in %s at line %d\n", \__FILE__, __LINE__); \exit(EXIT_FAILURE); }}static bool InitCUDA(){int count;cudaGetDeviceCount(&count);if (count == 0){fprintf(stderr, "There is no device.\n");return false;}int i;cudaDeviceProp prop = {0};for (i = 0; i < count; i++){if (cudaGetDeviceProperties(&prop, i) == cudaSuccess){if (prop.major >= 1){printf("%s\n", prop.name);break;}}}if (i >= count){fprintf(stderr, "There is no device supporting CUDA 1.x.\n");return false;}cudaSetDevice(i);return true;}static 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);}}}static  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[j + k * ldb];}c[i * ldc + j] = t;}}}static void compare_mat(const float* a, int lda, const float* b, int ldb, int n){float max_err = 0;float average_err = 0;float max_absolute_err = 0;int i, j;for (i = 0; i < n; i++){for (j = 0; j < n; j++){if (b[i * ldb + j] != 0){float tmp = fabs(a[i * lda + j] - b[i * ldb + j]);if (max_absolute_err < tmp){max_absolute_err = tmp;}float err = fabs(tmp / b[i * ldb + j]);if (max_err < err){max_err = err;}average_err += err;}}}printf("Max absolute error: %g\nMax error: %g\nAverage error: %g\n", \max_absolute_err, max_err, average_err / (n * n));}//利用 Kahan's Summation Formula 来提高精确度，并且利用块内共享内存存储矩阵dev_a的一行，一共n块。__global__ static void matMultCUDA(const float* dev_a, size_t lda, const float* dev_b, size_t ldb, float* dev_c, size_t ldc, int n){extern __shared__ float data[];const int tid = threadIdx.x;const int row = blockIdx.x;int i, j;for (i = tid; i < n; i += blockDim.x){data[i] = dev_a[row * lda + i];}__syncthreads();//row行，j列for (j = tid; j < n; j += blockDim.x){float t = 0;float y = 0;for (i = 0; i < n; i++){float r;y -= data[i] * dev_b[i * ldb + j];r = t - y;y = (r - t) + y;t = r;}dev_c[row * ldc + j] = t;}}//显存地址自动对齐，可以提高访问显存的效率static clock_t matmultCUDA(const float* a, int lda, const float* b, int ldb, float* c, int ldc, int n){const int thread_num = 256;float *dev_a, *dev_b, *dev_c;clock_t time;size_t pitch_a, pitch_b, pitch_c;cudaEvent_t     start, stop;float           elapsedTime;time = clock();HANDLE_ERROR(cudaEventCreate(&start));HANDLE_ERROR(cudaEventCreate(&stop));HANDLE_ERROR(cudaMallocPitch((void**)&dev_a, &pitch_a, sizeof(float)* n, n));HANDLE_ERROR(cudaMallocPitch((void**)&dev_b, &pitch_b, sizeof(float)* n, n));HANDLE_ERROR(cudaMallocPitch((void**)&dev_c, &pitch_c, sizeof(float)* n, n));HANDLE_ERROR(cudaMemcpy2D(dev_a, pitch_a, a, sizeof(float)* lda, sizeof(float)* n, n, cudaMemcpyHostToDevice));HANDLE_ERROR(cudaMemcpy2D(dev_b, pitch_b, b, sizeof(float)* ldb, sizeof(float)* n, n, cudaMemcpyHostToDevice));//int blocks = (n * n + thread_num - 1) / thread_num;int blocks = n;HANDLE_ERROR(cudaEventRecord(start, 0));matMultCUDA << <blocks, thread_num, n * sizeof(float) >> >(dev_a, pitch_a / sizeof(float), dev_b, pitch_b / sizeof(float), dev_c, pitch_c / sizeof(float), n);HANDLE_ERROR(cudaEventRecord(stop, 0));HANDLE_ERROR(cudaEventSynchronize(stop));HANDLE_ERROR(cudaEventElapsedTime(&elapsedTime, start, stop));printf("核函数执行时间：%lfms\n", elapsedTime);HANDLE_ERROR(cudaMemcpy2D(c, sizeof(float)* ldc, dev_c, pitch_c, sizeof(float)* n, n, cudaMemcpyDeviceToHost));cudaFree(dev_a);cudaFree(dev_b);cudaFree(dev_c);return clock() - time;}int main(int argc, char *argv[]){const int n = 64 * 4 * 2;float *a, *b, *c, *d;/*HANDLE_NULL(a = (float*)malloc(sizeof(float)* n * n));HANDLE_NULL(b = (float*)malloc(sizeof(float)* n * n));HANDLE_NULL(c = (float*)malloc(sizeof(float)* n * n));HANDLE_NULL(d = (float*)malloc(sizeof(float)* n * n));*/HANDLE_ERROR(cudaHostAlloc((void**)&a, sizeof(float)* n * n, cudaHostAllocDefault));HANDLE_ERROR(cudaHostAlloc((void**)&b, sizeof(float)* n * n, cudaHostAllocDefault));HANDLE_ERROR(cudaHostAlloc((void**)&c, sizeof(float)* n * n, cudaHostAllocDefault));HANDLE_ERROR(cudaHostAlloc((void**)&d, sizeof(float)* n * n, cudaHostAllocDefault));srand(0);matgen(a, n, n);matgen(b, n, n);printf("利用 Kahan's Summation Formula 来提高精确度\n");printf("显存地址自动对齐，可以提高访问显存的效率\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: %.6fs (%.6lf GFLOPS)\n", sec, 2.0 * n * n * n / (sec * 1E9));/*free(a);free(b);free(c);free(d);*/cudaFreeHost(a);cudaFreeHost(b);cudaFreeHost(c);cudaFreeHost(d);//remember to release the devicecudaDeviceReset();return 0;}