博主CUDA学习系列汇总传送门(持续更新):编程语言|CUDA入门
本文章为 《 GPU编程与优化 大众高性能计算》的读书笔记,例子也都取自书中教程。
矩阵乘法的运算量和数据量的关系不再是线性关系,对应的运算量为
O
(
n
3
)
O(n^3)
O(n3)。
矩阵乘法的数学表达式如下:
本质上看,矩阵陈发是向量内积的集合,
C
C
C矩阵的每一个元素都是
A
A
A矩阵一行和
B
B
B矩阵一列的向量内积。
在C语言中,数据是按行存储的,所以
B
B
B矩阵的访存不连续,导致catch命中失效。
一、CPU下一般矩阵乘法
其中
A
A
A矩阵大小为
m
m
m行
l
l
l列,
B
B
B矩阵的大小为
l
l
l行
n
n
n列,相乘后得到
C
C
C矩阵大小为
m
m
m行
n
n
n列。
这种方法的思路是一次完成
C
C
C矩阵中一个元素的计算,然后通过循环迭代的方法完成
C
C
C矩阵所有元素的运算。
平均耗时为 4458ms
#include <stdio.h>
#include <iostream>
#include <time.h>
#include <sys/time.h>
long get_time()
{
struct timeval t_start;
gettimeofday(&t_start, NULL);
long end = ((long)t_start.tv_sec) * 1000 + (long)t_start.tv_usec / 1000;;
return end;
}
/*
A * B = C
A 矩阵大小为 m * l, 对应[0, 1, 2, ...1023],[0, 1, 2, ...1023]
B 矩阵大小为 l * n 对应[1, 1, 1, 1, 1, 1, 1]
C 矩阵大小为 m * n
*/
const int M = 2048;
const int L = 1024;
const int N = 512;
/* 一、cpu 上矩阵相乘 一般矩阵乘法*/
template <typename T>
void mat_dot_cpu_1(T *a, T *b, T *c, int m, int n, int l)
{
double dTemp = 0.0;
for (int i = 0; i < m; ++i)
{
for (int j = 0; j < n; ++j)
{
dTemp = 0.0;
for (int k = 0; k < l; ++k)
{
dTemp += a[i*l+k] * b[n*k+j];
}
c[ i * n + j] = dTemp;
}
}
}
int main()
{
float a[M * L], b[L * N];
float c[M * N];
for (int i=0; i<M; ++i) // 为数组赋值
{
for (int j=0; j< L; ++j)
{
a[i*L+j] = j;
}
}
for (int i=0; i<L; ++i)
{
for(int j=0; j<N; ++j)
{
b[i*N + j] = 1;
}
}
mat_dot_cpu_1(a, b, c, M, N, L);
long start = get_time();
for (int i=0; i<10; ++i)
{
mat_dot_cpu_1(a, b, c, M, N, L);
}
long end = get_time();
printf("cost time is %f\n", 1.0*(end - start) / 10.0);
printf("c[0]~c[10]");
for (int i=0; i<10;++i)
{
printf(" = %.1f, ", c[i]);
}
printf("\nhello world\n");
return 0;
}
二、CPU下循环交换矩阵乘法
本例子平均耗时 2479ms;
注意,因为这种方法没有采用中间变量,所以当矩阵数字很大时,超出float能表示的范围会导致结果不正确。
#include <stdio.h>
#include <iostream>
#include <time.h>
#include <sys/time.h>
long get_time()
{
struct timeval t_start;
gettimeofday(&t_start, NULL);
long end = ((long)t_start.tv_sec) * 1000 + (long)t_start.tv_usec / 1000;;
return end;
}
/*
A * B = C
A 矩阵大小为 m * l, 对应[0, 1, 2, ...1023],[0, 1, 2, ...1023]
B 矩阵大小为 l * n 对应[1, 1, 1, 1, 1, 1, 1]
C 矩阵大小为 m * n
*/
const int M = 2048;
const int L = 1024;
const int N = 512;
/* 二、cpu 上矩阵相乘 循环交换矩阵乘法
A数据的内存是连续访问的
每次读取A的一个元素,与B矩阵对应的一行做乘积运算,累加乘积结果到C矩阵的对应行*/
template <typename T>
void mat_dot_cpu_2(T *a, T *b, T *c, int m, int n, int l)
{
double temp = 0.0;
for (int i = 0; i < m; ++i)
{
for (int k = 0; k < l; ++k)
{
temp = a[i * l + k];
for (int j = 0; j < n; ++j)
{
c[i * n + j] += temp * b[k * n + j];
}
}
}
}
int main()
{
float a[M * L], b[L * N];
float c[M * N];
for (int i=0; i<M; ++i) // 为数组赋值
{
for (int j=0; j< L; ++j)
{
a[i*L+j] = j;
}
}
for (int i=0; i<L; ++i)
{
for(int j=0; j<N; ++j)
{
b[i*N + j] = 1;
}
}
for (int i = 0; i<M*N; ++i)
{
c[i] = 0.0;
}
mat_dot_cpu_2(a, b, c, M, N, L);
long start = get_time();
for (int i=0; i<10; ++i)
{
for (int i = 0; i<M*N; ++i) // 得初始化一下
{
c[i] = 0.0;
}
mat_dot_cpu_2(a, b, c, M, N, L);
}
long end = get_time();
printf("cost time is %f\n", 1.0*(end - start) / 10.0);
printf("c[0]~c[10]");
for (int i=0; i<10;++i)
{
printf(" = %.1f, ", c[i]);
}
printf("\nhello world\n");
return 0;
}
三、CPU下转置矩阵乘法
为了解决访存不连续,无法向量化的问题,除了内部两层循环交换的方法,还可以采用转置矩阵的方法。
即对
B
B
B矩阵转置,转置后,就编程
A
A
A矩阵的一行和
B
B
B矩阵的一行做向量内积运算。因为
B
B
B矩阵是按行访问的,因此具有cache命中率、数据连续访存、可向量化等有点。
#include <stdio.h>
#include <iostream>
#include <time.h>
#include <sys/time.h>
/* 三、cpu 上矩阵相乘 转置矩阵乘法
要解决矩阵乘法中访存不连续、无法向量化的问题,除了内部两次循环交换的方法,还有矩阵转置的方法
对B矩阵转置,转置后,原来A矩阵一行和B矩阵一列的向量内积运算就变成了A矩阵的一行和B矩阵的一行的向量内积运算*/
template <typename T>
void mat_dot_cpu_3(T *a, T *b, T *c, int m, int n, int l)
{
T *b_t = NULL;
b_t = (T *)malloc(L * N * sizeof(T));
// 将b矩阵转置成 b_t矩阵
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < l; ++j)
{
b_t[i * l + j] = b[j * n + i];
}
}
for (int i = 0; i < m; ++i)
{
for (int j = 0; j < n; ++j)
{
double temp = 0.0;
for (int k = 0; k < l; ++k)
{
temp += a[i * l + k] * b_t[j * l + k];
}
c[i * n + j] = temp;
}
}
if(NULL != b_t)
{
free(b_t);
b_t = NULL;
}
}
未完待续。。。