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,小于等于1
const 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之间数据的交互不用人管,机器已做好。
加速比曲线: