CUDA计算矩阵幂

**

一.实验内容

本次实验内容为基于CUDA的GPU实现矩阵的幂。要求分别用暴力算法和高效算法实现矩阵的幂。

对于一个 的方阵 ，计算的次幂。首先，生成一个的方阵，保证每行每列元素之和满足(0,1])。

·暴力算法

N个矩阵相乘

·高效算法

利用矩阵乘法的结合律

暴力算法：

#include "cuda_runtime.h"#include "device_launch_parameters.h"#include <stdio.h>#include <stdlib.h> #include <time.h> #include <math.h>const int m = 16; //m*m矩阵const int N = 64; //幂次double v = 0.05;//保证每行的和，每列的和都大于0，小于等于1const int Blocks_Per_grid = 4; const int Threads_Per_block = 4;double **matrix; //原矩阵double **CPU_matrix;//存储CPU串行计算得到的幂结果void init_matrix(){matrix = (double **)malloc(m*sizeof(double*));CPU_matrix = (double **)malloc(m*sizeof(double*));for(int i=0;i<m;i++){matrix[i] = (double*)malloc(m*sizeof(double));CPU_matrix[i] = (double*)malloc(m*sizeof(double));}for(int i=0;i<m;i++){for(int j=0;j<m;j++){matrix[i][j] = v;CPU_matrix[i][j] = 0.0;}}}void CPU_power(){int n = N-1;double **temp; //存储临时结果temp = (double **)malloc(m*sizeof(double*));for(int i=0;i<m;i++){temp[i] = (double*)malloc(m*sizeof(double));}for(int i=0;i<m;i++){for(int j=0;j<m;j++){temp[i][j] = v;}}while(n--){for(int i=0;i<m;i++){for(int j=0;j<m;j++){for(int k=0;k<m;k++){CPU_matrix[i][j] += temp[i][k]*matrix[k][j];}}}for(int i=0;i<m;i++){for(int j=0;j<m;j++){temp[i][j] = CPU_matrix[i][j];}}}free(temp);}__global__ void compute(double *d_matrix1,double *d_a,double *d_b){//每个线程求一行int tid = threadIdx.x + blockDim.x* blockIdx.x; //线程号算对应行for(int i=0;i<N-1;i++){for(int j=0;j<m;j++){double temp =0.0;for(int k=0;k<m;k++){temp += d_a[tid*m+k] * d_matrix1[k*m+j];}d_b[tid*m+j] = temp;}for(int i=0;i<m;i++){d_a[tid*m+i] = d_b[tid*m+i];}}}int main(){double *matrix1;matrix1 = (double*)malloc((m*m)*sizeof(double));double *b; // a保存上次乘完的结果，b为本次乘完的结果double *d_matrix1,*d_b,*d_a; b = (double*)malloc(sizeof(double)*(m*m));for(int i=0;i<(m*m);i++){matrix1[i] = v;}clock_t begin,end;float period;// cudaEvent_t start,stop;//float time;begin = clock();init_matrix(); //初始化矩阵CPU_power(); //CPU串行求矩阵幂end = clock();period = (float)1000*(end - begin)/CLOCKS_PER_SEC;printf("CPU cost time: %.6f ms\n", period);cudaMalloc((void**)&d_matrix1,sizeof(double)*(m*m)); cudaMalloc((void**)&d_b,sizeof(double)*(m*m)); cudaMalloc((void**)&d_a,sizeof(double)*(m*m));cudaMemcpy(d_matrix1,matrix1,sizeof(double)*(m*m),cudaMemcpyHostToDevice); //copy数组，不带地址符号！！！！！//DEBUG： printf！！！ cudaMemcpy(d_a,matrix1,sizeof(double)*(m*m),cudaMemcpyHostToDevice);cudaMemcpy(d_b,matrix1,sizeof(double)*(m*m),cudaMemcpyHostToDevice);//cudaEventCreate(&start);//cudaEventCreate(&stop);//cudaEventRecord(start,0);begin=clock();compute<<<Blocks_Per_grid,Threads_Per_block>>>(d_matrix1,d_a,d_b);end=clock();period = (float)1000*(end - begin)/CLOCKS_PER_SEC;printf("GPU cost time: %.6f ms\n", period);//cudaEventRecord(stop,0);//cudaEventSynchronize(stop);//cudaEventElapsedTime(&time,start,stop);//printf("CPU+GPU cost time: %.6f ms\n", time);cudaMemcpy(b,d_b,sizeof(double)*(m*m),cudaMemcpyDeviceToHost); //copy耗时很久cudaFree(d_b);free(b);free(matrix);free(CPU_matrix);return 0;}在这里插入代码片

高效算法：简单利用结合律减小计算量，即将相邻矩阵合并为平方再算。也可以采用矩阵快速幂。

流程图：

采用__managed统一内存寻址，CUP和GPU之间数据的交互不用人管，机器已做好。

加速比曲线：