CUDA 编程学习 <2> - 7 种还原算法优化方法 - 代码部分
最编程
2024-04-25 22:12:32
...
3.1、优化1：交错寻址（Interleaved addressing divergent branches）
3.1.1、核函数（优化1、优化2、优化3）：
__global__ void reduce0(unsigned int *g_idata, long long int *g_odata, int N, clock_t* time,int GuiYueWay) //
{//包含归约1、优化2、优化3的程序。
//__shared__ int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
//__shared__ unsigned int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
extern __shared__ long long int sdata[];//动态申请共享内存，空间大小在三个尖括号的第三个参数进行传入<<<1，2，3>>>开辟2的指数个空间，2、
clock_t start = clock();
// each thread loads one element from global to shared mem
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
//printf("tid:%d,i:%d\n",tid, i);
//方法1，要求一次性把线程数量建完，这样不确定有没有所有最大的线程限定，特点是速度快
sdata[tid] = (i < N) ? g_idata[i] : 0;
//printf("CUDA1 i:%d\n", i);
方法2，可以一个线程多次处理，但是在数据较大的时候结果容易出错
//unsigned int temp = 0;
//while (i < N)//N
//{//这种方式可以一个线程多次处理
// temp += g_idata[i];
// i += blockDim.x * gridDim.x;
// //printf("tid:%d\n", tid);
//}
设置cache中相应位置上的值
//sdata[tid] = temp;
__syncthreads();//同步共享变量
switch (GuiYueWay)
{
case 1:
GuiYue_Opt1(sdata, tid);//初始版本的归约函数
break;
case 2:
GuiYue_Opt2(sdata, tid);//优化的第2版本的归约函数
break;
case 3:
GuiYue_Opt3(sdata, tid);//优化的第3版本的归约函数
break;
default:
break;
}
// //do reduction in shared mem
//for (unsigned int s = 1; s < blockDim.x; s *= 2) {
// if (tid % (2 * s) == 0) {
// sdata[tid] += sdata[tid + s];
// }
// __syncthreads();
//}
// write result for this block to global mem
//printf("tid2:%d,i2:%d\n", tid, i);
if (tid == 0)
{
g_odata[blockIdx.x] = sdata[0];
//printf("sdata[0]:%u\n", sdata[0]);
//printf("g_odata[%d]:%u\n", blockIdx.x, sdata[0]);
}
*time = clock() - start;//计时
}
核函数说明
需要说明的是里面用到了共享内存，共享内存的开辟有两种方式，静态定义和动态开辟，可以参考如下链接：
CUDA之静态、动态共享内存分配详解
简要介绍共享内存
// __shared__ 为共享内存标识
//动态定义如下，在核函数内部定义，
extern __shared__ unsigned int sdata[];//动态申请共享内存，空间大小在核函数三个尖括号的第三个参数进行传入<<<1，2，3>>>，如需开辟20个空间，则第三个参数应该是20*sizeof(unsigned int)
//静态定义
__shared__ unsigned int sdata[20];//静态开辟20个静态空间
3.1.2、交错寻址函数代码：
__device__ void GuiYue_Opt1(long long int* sdata, unsigned int tid)
{//归约的方式1，交错寻址
for (unsigned int s = 1; s < blockDim.x; s *= 2) {
if (tid % (2 * s) == 0) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
}
3.2、优化2：跨步寻址(Interleaved addressing bank conflicts)
3.2.1、核函数
核函数使用3.1交错寻址的核函数的switch语句调用
3.2.2、跨步寻址函数代码：
__device__ void GuiYue_Opt2(long long int* sdata, unsigned int tid)
{//归约的方式2，跨步寻址
for (unsigned int s = 1; s < blockDim.x; s *= 2) {
int index = 2 * s*tid;
if (index<blockDim.x) {
sdata[index] += sdata[index + s];
}
__syncthreads();
}
}
3.3、优化3：连续寻址（Sequential Addressing）
3.1、核函数
核函数使用3.1交错寻址的核函数的switch语句调用
3.2、连续寻址函数代码：
__device__ void GuiYue_Opt3(long long int* sdata, unsigned int tid)
{//归约的方式3，连续寻址
for (unsigned int s = blockDim.x / 2;s>0; s >>= 1) //注意：>>右移1位，假设s=5,那么二进制为0101 s>>1 右移1位，即把最右边的1位删掉，变为010 此时结果为2;如果s=5,s>>2,即右移2位，变成01，结果为1
{
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
}
3.4、优化4：把一次归约放在加载时进行（first add during load）
3.4.1、核函数：
主要还是用优化3的代码算法，更改核函数的赋值共享内存里，核函数代码如下
__global__ void reduce4( unsigned int *g_idata, long long int *g_odata, int N, clock_t* time, int GuiYueWay) //
{//在优化3的基础上进行优化4和优化5的程序。
//__shared__ int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
//__shared__ unsigned int sdata[threadsPerBlock];//全局变量定义共享内存空间，开辟2的指数个空间，2、4、8、32 ...
extern __shared__ long long int sdata[];//动态申请共享内存，空间大小在三个尖括号的第三个参数进行传入<<<1，2，3>>>开辟2的指数个空间。当数据过大容易越界，需定义大点的空间
clock_t start = clock();
// each thread loads one element from global to shared mem
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*(blockDim.x*2) + threadIdx.x;//优化位置1
//printf("tid:%d,i:%d\n",tid, i);
//sdata[tid] = (i < N) ? g_idata[i] : 0;
sdata[tid] = (i < N) ? g_idata[i]+ g_idata[i + blockDim.x] : 0;//优化位置2，把第一次归约放在加载数据的时候进行，减少了一次CUDA内部操作
//printf("CUDA4 i:%d,2i:%d\n", i, i + blockDim.x);
__syncthreads();//同步共享变量
switch (GuiYueWay)
{
case 4:
GuiYue_Opt3(sdata, tid);//优化的第3版本的归约函数
break;
case 5:
GuiYue_Opt5(sdata, tid);//优化的第5版本的归约函数
break;
default:
break;
}
//GuiYue_Opt3(sdata, tid);//优化的第3版本的归约函数
if (tid == 0)
{
g_odata[blockIdx.x] = sdata[0];
//printf("g_odata[%d]:%u\n", blockIdx.x,sdata[0]);
}
*time = clock() - start;//计时
}
3.4.2 算法优化部分
unsigned int i = blockIdx.x*(blockDim.x*2) + threadIdx.x;//优化位置1
sdata[tid] = (i < N) ? g_idata[i]+ g_idata[i + blockDim.x] : 0;//优化位置2，把第一次归约放在加载数据的时候进行，减少了一次CUDA内部操作
如PPT
3.5、优化5：最后线程束展开（Unroll the last warp）
5.1 核函数
核函数使用3.4优化4的核函数的switch语句调用
5.2 最后线束展开函数代码：
__device__ void warpReduce(volatile long long int* sdata,int tid)
{//该函数注意volatile，如果没有这个关键字，编译器可能会将sdata变量优化掉，导致结果出错的
sdata[tid] += sdata[tid + 32];
sdata[tid] += sdata[tid + 16];
sdata[tid] += sdata[tid + 8];
sdata[tid] += sdata[tid + 4];
sdata[tid] += sdata[tid + 2];
sdata[tid] += sdata[tid + 1];
}
__device__ void GuiYue_Opt5(long long int* sdata, unsigned int tid)
{//归约的方式5，该函数在归约方式3的基础上进行优化，省去warp过程中不必要的操作
for (unsigned int s = blockDim.x / 2; s > 32; s >>= 1) //注意：>>右移1位，假设s=5,那么二进制为0101 s>>1 右移1位，即把最右边的1位删掉，变为010 此时结果为2;如果s=5,s>>2,即右移2位，变成01，结果为1
{
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
if (tid < 32)
{
warpReduce(sdata, tid);
}
}
3.6、优化6：完全展开（completely Unrolled）
3.6.1、核函数：
把核函数优化成了模板方法增加通用性，另外增加了warpReduceT模板函数
//下面函数优化程序6使用
template<unsigned int blockSize>
__device__ void warpReduceT(volatile long long int* sdata, unsigned int tid)
{//该函数注意volatile，如果没有这个关键字，编译器可能会将sdata变量优化掉，导致结果出错的
if (blockSize >= 64)sdata[tid] += sdata[tid + 32];
if (blockSize >= 32)sdata[tid] += sdata[tid + 16];
if (blockSize >= 16)sdata[tid] += sdata[tid + 8];
if (blockSize >= 8)sdata[tid] += sdata[tid + 4];
if (blockSize >= 4)sdata[tid] += sdata[tid + 2];
if (blockSize >= 2)sdata[tid] += sdata[tid + 1];
}
template <unsigned int blockSize>
__global__ void reduce6( unsigned int *g_idata, long long int *g_odata, int N, clock_t* time) //
{//优化6程序
//__shared__ int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
//__shared__ unsigned int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
extern __shared__ long long int sdata[];//动态申请共享内存，空间大小在三个尖括号的第三个参数进行传入<<<1，2，3>>>开辟2的指数个空间，2、4、8、32 ...
clock_t start = clock();
// each thread loads one element from global to shared mem
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*(blockDim.x * 2) + threadIdx.x;
//printf("tid:%d,i:%d\n",tid, i);
//sdata[tid] = (i < N) ? g_idata[i] : 0;
sdata[tid] = (i < N) ? g_idata[i] + g_idata[i + blockDim.x] : 0;//把第一次归约放在加载数据的时候进行，减少了一次CUDA内部操作
__syncthreads();//同步共享变量
//优化的第6版本的归约函数
if (blockSize >= 512)
{
if (tid < 256) { sdata[tid] += sdata[tid + 256]; }__syncthreads();
}
if (blockSize >= 256)
{
if (tid < 128) { sdata[tid] += sdata[tid + 128]; }__syncthreads();
}
if (blockSize >= 128)
{
if (tid < 64) { sdata[tid] += sdata[tid + 64]; }__syncthreads();
}
if (tid < 32)warpReduceT<blockSize>(sdata, tid);
if (tid == 0)
{
g_odata[blockIdx.x] = sdata[0];
//printf("sdata[0]:%d\n", sdata[0]);
}
*time = clock() - start;//计时
}
7、优化7:每个线程有多个元素（Multiple elements per thread）
与第6个优化版本相比，在核函数的时候就进行了合并，另外使用了优化6中的warpReduceT模板函数
template <unsigned int blockSize>
__global__ void reduce7(unsigned int *g_idata,long long int *g_odata, int N, clock_t* time) //
{//优化6程序
//__shared__ int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
//__shared__ unsigned int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
extern __shared__ long long int sdata[];//动态申请共享内存，空间大小在三个尖括号的第三个参数进行传入<<<1，2，3>>>开辟2的指数个空间，2、4、8、32 ...
clock_t start = clock();
// each thread loads one element from global to shared mem
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*(blockSize * 2) + threadIdx.x;//更改位置1
unsigned int gridSize = blockSize * 2 * gridDim.x;
sdata[tid] = 0;
//printf("tid:%d,i:%d\n",tid, i);
//sdata[tid] = (i < N) ? g_idata[i] : 0;
//sdata[tid] = (i < N) ? g_idata[i] + g_idata[i + blockDim.x] : 0;//把第一次归约放在加载数据的时候进行，减少了一次CUDA内部操作
while (i<N)
{
sdata[tid] += g_idata[i] + g_idata[i + blockSize];
i += gridSize;
}
__syncthreads();//同步共享变量
//优化的第6版本的归约函数
if (blockSize >= 512)
{
if (tid < 256) { sdata[tid] += sdata[tid + 256]; }__syncthreads();
}
if (blockSize >= 256)
{
if (tid < 128) { sdata[tid] += sdata[tid + 128]; }__syncthreads();
}
if (blockSize >= 128)
{
if (tid < 64) { sdata[tid] += sdata[tid + 64]; }__syncthreads();
}
if (tid < 32)warpReduceT<blockSize>(sdata, tid);//
if (tid == 0)
{
g_odata[blockIdx.x] = sdata[0];
//printf("sdata[0]:%d\n", sdata[0]);
}
*time = clock() - start;//计时
}
英文pdf上有几处错误和需要注意的地方：
优化7最后的这个代码，因为是调用模板方法，所以需要增加模板值
8、汇总以上所有测试代码，包括主函数
/*
归约样例，增加时间统计，关注宽带占用情况
参考链接：
中文：https://wenku.baidu.com/link?url=LOvDSVU9fyqQko7h4rIWEuxwmKCgCOzrFSGoE1saaNWpGlmfl2wR9Qu--IcdwM9D0H1cIT_0mUFLrUazYPFuUrc57O_RbD8_IbS27MjFmtS
英文原版：https://developer.download.nvidia.cn/assets/cuda/files/reduction.pdf
优化方式：
归约的方式1，交错寻址
归约的方式2，跨步寻址
归约的方式3，连续寻址
归约的方式4，在连续寻址基础上进行改进。由于第一次循环迭代中，一半线程处于空闲状态，太浪费了，所以优化。比前三种优化，每个共享内存中的数值大约是
归约的方式5，在最后进行了swarp展开，实测速度不如第四次快
归约的方式6，完全展开，把核函数写成模板方法进行调用，在优化部分进行了模板方法的完全展开。模板方法会在编译的时候进行评估计算
归约的方式7，把第4种方法和第6种方法结合。
*/
///*
#include "device_launch_parameters.h"
#include "cuda_runtime.h"
#include "device_functions.h"
#include <stdio.h>
#include <iostream>
#include <time.h>
#define imin(a,b) (a<b?a:b)
//#define N 131072*2//1048576 //数组长度,总任务 1024//
//#define threadsPerBlock 128//后面执行归约并行计算，该值必须是2的指数值
//#define blocksPerGrid 32//规约后的长度
__device__ void GuiYue_Opt1(long long int* sdata, unsigned int tid)
{//归约的方式1，交错寻址
for (unsigned int s = 1; s < blockDim.x; s *= 2) {
if (tid % (2 * s) == 0) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
}
__device__ void GuiYue_Opt2(long long int* sdata, unsigned int tid)
{//归约的方式2，跨步寻址
for (unsigned int s = 1; s < blockDim.x; s *= 2) {
int index = 2 * s*tid;
if (index<blockDim.x) {
sdata[index] += sdata[index + s];
}
__syncthreads();
}
}
__device__ void GuiYue_Opt3(long long int* sdata, unsigned int tid)
{//归约的方式3，连续寻址
for (unsigned int s = blockDim.x / 2;s>0; s >>= 1) //注意：>>右移1位，假设s=5,那么二进制为0101 s>>1 右移1位，即把最右边的1位删掉，变为010 此时结果为2;如果s=5,s>>2,即右移2位，变成01，结果为1
{
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
}
__device__ void warpReduce(volatile long long int* sdata,int tid)
{//该函数注意volatile，如果没有这个关键字，编译器可能会将sdata变量优化掉，导致结果出错的
sdata[tid] += sdata[tid + 32];
sdata[tid] += sdata[tid + 16];
sdata[tid] += sdata[tid + 8];
sdata[tid] += sdata[tid + 4];
sdata[tid] += sdata[tid + 2];
sdata[tid] += sdata[tid + 1];
}
__device__ void GuiYue_Opt5(long long int* sdata, unsigned int tid)
{//归约的方式5，该函数在归约方式3的基础上进行优化，省去warp过程中不必要的操作
for (unsigned int s = blockDim.x / 2; s > 32; s >>= 1) //注意：>>右移1位，假设s=5,那么二进制为0101 s>>1 右移1位，即把最右边的1位删掉，变为010 此时结果为2;如果s=5,s>>2,即右移2位，变成01，结果为1
{
if (tid < s) {
sdata[tid] += sdata[tid + s];
}
__syncthreads();
}
if (tid < 32)
{
warpReduce(sdata, tid);
}
}
//下面函数优化程序6使用
template<unsigned int blockSize>
__device__ void warpReduceT(volatile long long int* sdata, unsigned int tid)
{//该函数注意volatile，如果没有这个关键字，编译器可能会将sdata变量优化掉，导致结果出错的
if (blockSize >= 64)sdata[tid] += sdata[tid + 32];
if (blockSize >= 32)sdata[tid] += sdata[tid + 16];
if (blockSize >= 16)sdata[tid] += sdata[tid + 8];
if (blockSize >= 8)sdata[tid] += sdata[tid + 4];
if (blockSize >= 4)sdata[tid] += sdata[tid + 2];
if (blockSize >= 2)sdata[tid] += sdata[tid + 1];
}
__global__ void reduce0(unsigned int *g_idata, long long int *g_odata, int N, clock_t* time,int GuiYueWay) //
{//包含归约1、优化2、优化3的程序。
//__shared__ int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
//__shared__ unsigned int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
extern __shared__ long long int sdata[];//动态申请共享内存，空间大小在三个尖括号的第三个参数进行传入<<<1，2，3>>>开辟2的指数个空间，2、
clock_t start = clock();
// each thread loads one element from global to shared mem
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
//printf("tid:%d,i:%d\n",tid, i);
//方法1，要求一次性把线程数量建完，这样不确定有没有所有最大的线程限定，特点是速度快
sdata[tid] = (i < N) ? g_idata[i] : 0;
//printf("CUDA1 i:%d\n", i);
方法2，可以一个线程多次处理，但是在数据较大的时候结果容易出错
//unsigned int temp = 0;
//while (i < N)//N
//{//这种方式可以一个线程多次处理
// temp += g_idata[i];
// i += blockDim.x * gridDim.x;
// //printf("tid:%d\n", tid);
//}
设置cache中相应位置上的值
//sdata[tid] = temp;
__syncthreads();//同步共享变量
switch (GuiYueWay)
{
case 1:
GuiYue_Opt1(sdata, tid);//初始版本的归约函数
break;
case 2:
GuiYue_Opt2(sdata, tid);//优化的第2版本的归约函数
break;
case 3:
GuiYue_Opt3(sdata, tid);//优化的第3版本的归约函数
break;
default:
break;
}
// //do reduction in shared mem
//for (unsigned int s = 1; s < blockDim.x; s *= 2) {
// if (tid % (2 * s) == 0) {
// sdata[tid] += sdata[tid + s];
// }
// __syncthreads();
//}
// write result for this block to global mem
//printf("tid2:%d,i2:%d\n", tid, i);
if (tid == 0)
{
g_odata[blockIdx.x] = sdata[0];
//printf("sdata[0]:%u\n", sdata[0]);
//printf("g_odata[%d]:%u\n", blockIdx.x, sdata[0]);
}
*time = clock() - start;//计时
}
__global__ void reduce4( unsigned int *g_idata, long long int *g_odata, int N, clock_t* time, int GuiYueWay) //
{//在优化3的基础上进行优化4和优化5的程序。
//__shared__ int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
//__shared__ unsigned int sdata[threadsPerBlock];//全局变量定义共享内存空间，开辟2的指数个空间，2、4、8、32 ...
extern __shared__ long long int sdata[];//动态申请共享内存，空间大小在三个尖括号的第三个参数进行传入<<<1，2，3>>>开辟2的指数个空间。当数据过大容易越界，需定义大点的空间
clock_t start = clock();
// each thread loads one element from global to shared mem
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*(blockDim.x*2) + threadIdx.x;//优化位置1
//printf("tid:%d,i:%d\n",tid, i);
//sdata[tid] = (i < N) ? g_idata[i] : 0;
sdata[tid] = (i < N) ? g_idata[i]+ g_idata[i + blockDim.x] : 0;//优化位置2，把第一次归约放在加载数据的时候进行，减少了一次CUDA内部操作
//printf("CUDA4 i:%d,2i:%d\n", i, i + blockDim.x);
__syncthreads();//同步共享变量
switch (GuiYueWay)
{
case 4:
GuiYue_Opt3(sdata, tid);//优化的第3版本的归约函数
break;
case 5:
GuiYue_Opt5(sdata, tid);//优化的第5版本的归约函数
break;
default:
break;
}
//GuiYue_Opt3(sdata, tid);//优化的第3版本的归约函数
if (tid == 0)
{
g_odata[blockIdx.x] = sdata[0];
//printf("g_odata[%d]:%u\n", blockIdx.x,sdata[0]);
}
*time = clock() - start;//计时
}
template <unsigned int blockSize>
__global__ void reduce6( unsigned int *g_idata, long long int *g_odata, int N, clock_t* time) //
{//优化6程序
//__shared__ int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
//__shared__ unsigned int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
extern __shared__ long long int sdata[];//动态申请共享内存，空间大小在三个尖括号的第三个参数进行传入<<<1，2，3>>>开辟2的指数个空间，2、4、8、32 ...
clock_t start = clock();
// each thread loads one element from global to shared mem
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*(blockDim.x * 2) + threadIdx.x;
//printf("tid:%d,i:%d\n",tid, i);
//sdata[tid] = (i < N) ? g_idata[i] : 0;
sdata[tid] = (i < N) ? g_idata[i] + g_idata[i + blockDim.x] : 0;//把第一次归约放在加载数据的时候进行，减少了一次CUDA内部操作
__syncthreads();//同步共享变量
//优化的第6版本的归约函数
if (blockSize >= 512)
{
if (tid < 256) { sdata[tid] += sdata[tid + 256]; }__syncthreads();
}
if (blockSize >= 256)
{
if (tid < 128) { sdata[tid] += sdata[tid + 128]; }__syncthreads();
}
if (blockSize >= 128)
{
if (tid < 64) { sdata[tid] += sdata[tid + 64]; }__syncthreads();
}
if (tid < 32)warpReduceT<blockSize>(sdata, tid);
if (tid == 0)
{
g_odata[blockIdx.x] = sdata[0];
//printf("sdata[0]:%d\n", sdata[0]);
}
*time = clock() - start;//计时
}
template <unsigned int blockSize>
__global__ void reduce7(unsigned int *g_idata,long long int *g_odata, int N, clock_t* time) //
{//优化6程序
//__shared__ int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
//__shared__ unsigned int sdata[threadsPerBlock];//共享内存，开辟2的指数个空间，2、4、8、32 ...
extern __shared__ long long int sdata[];//动态申请共享内存，空间大小在三个尖括号的第三个参数进行传入<<<1，2，3>>>开辟2的指数个空间，2、4、8、32 ...
clock_t start = clock();
// each thread loads one element from global to shared mem
unsigned int tid = threadIdx.x;
unsigned int i = blockIdx.x*(blockSize * 2) + threadIdx.x;//更改位置1
unsigned int gridSize = blockSize * 2 * gridDim.x;
sdata[tid] = 0;
//printf("tid:%d,i:%d\n",tid, i);
//sdata[tid] = (i < N) ? g_idata[i] : 0;
//sdata[tid] = (i < N) ? g_idata[i] + g_idata[i + blockDim.x] : 0;//把第一次归约放在加载数据的时候进行，减少了一次CUDA内部操作
while (i<N)
{
sdata[tid] += g_idata[i] + g_idata[i + blockSize];
i += gridSize;
}
__syncthreads();//同步共享变量
//优化的第6版本的归约函数
if (blockSize >= 512)
{
if (tid < 256) { sdata[tid] += sdata[tid + 256]; }__syncthreads();
}
if (blockSize >= 256)
{
if (tid < 128) { sdata[tid] += sdata[tid + 128]; }__syncthreads();
}
if (blockSize >= 128)
{
if (tid < 64) { sdata[tid] += sdata[tid + 64]; }__syncthreads();
}
if (tid < 32)warpReduceT<blockSize>(sdata, tid);//
if (tid == 0)
{
g_odata[blockIdx.x] = sdata[0];
//printf("sdata[0]:%d\n", sdata[0]);
}
*time = clock() - start;//计时
}
int main(int argc, char*argv[])
{
printf("CUDA测试程序:\n");
int N = 131072;//总任务数
int GuiYueWayCPU = 1;//归约的方法，共有7种优化方式
int threadsPerBlock = 128;//后面执行归约并行计算，该值必须是2的指数值
printf("argc:%d\n", argc);
if (argc == 2)
{
N = atoi(argv[1]);
}
else if (argc == 3)
{
N = atoi(argv[1]);
GuiYueWayCPU = atoi(argv[2]);
}
else if (argc >= 4)
{
N = atoi(argv[1]);
GuiYueWayCPU = atoi(argv[2]);
threadsPerBlock = atoi(argv[3]);
}
printf("N:%d\tGuiYueWay:%d\n", N, GuiYueWayCPU);
int blocksPerGrid = 32;//规约后的长度，重新计算
blocksPerGrid = (N+ threadsPerBlock-1) / threadsPerBlock;//这样可以并发一次搞定
printf("blocksPerGrid:%d,threadsPerBlock:%d\n", blocksPerGrid, threadsPerBlock);
unsigned int * inData = new unsigned int[N] {0};
//int * OutData = new int[N] {0};
long long int * OutData = new long long int[blocksPerGrid] {0};//输出开辟空间为blocksPerGrid
clock_t time_used;//CPU 上的
//long int * inData = new long int[N] {0};
int * OutData = new int[N] {0};
//long int * OutData = new long int[blocksPerGrid] {0};//输出开辟空间为blocksPerGrid
for (int i = 1; i <= N; i++)
{//赋初始值
inData[i - 1] = i;
}
unsigned int *inData_dev;
long long int *OutData_dev;//进GPU上的
clock_t* time;//进GPU上的
cudaMalloc((void**)&inData_dev, N * sizeof(unsigned int));
//cudaMalloc((void**)&OutData_dev, N * sizeof(int));
cudaMalloc((void**)&OutData_dev, blocksPerGrid * sizeof(long long int));
cudaMalloc((void**)&time, sizeof(clock_t));//开辟空间
cudaMemcpy(inData_dev, inData, N * sizeof(unsigned int), cudaMemcpyHostToDevice);
//reduce0 << <3, 8 >> > (inData_dev, OutData_dev);
int smemSize = threadsPerBlock * sizeof(long long int);//动态分配的共享内存大小
if (GuiYueWayCPU <= 3)
{//尖括号中第三个参数为动态分配共享内存大小
reduce0 << <blocksPerGrid, threadsPerBlock, smemSize >> > (inData_dev, OutData_dev, N, time, GuiYueWayCPU);
}
else if (GuiYueWayCPU == 4|| GuiYueWayCPU == 5)
{
reduce4 << <blocksPerGrid, threadsPerBlock, smemSize >> > (inData_dev, OutData_dev, N, time, GuiYueWayCPU);
}
else if (GuiYueWayCPU == 6)
{
switch (threadsPerBlock)
{
case 512:
reduce6 <512> << <blocksPerGrid, threadsPerBlock, smemSize >> > (inData_dev, OutData_dev, N, time);
break;
case 256:
reduce6 <256> << <blocksPerGrid, threadsPerBlock, smemSize >> > (inData_dev, OutData_dev, N, time);
break;
case 128:
reduce6 <128> << <blocksPerGrid, threadsPerBlock, smemSize >> > (inData_dev, OutData_dev, N, time);
break;
case 64:
reduce6 <64> << <blocksPerGrid, threadsPerBlock, smemSize >> > (inData_dev, OutData_dev, N, time);
break;
case 32:
reduce6 <32> << <blocksPerGrid, threadsPerBlock, smemSize >> > (inData_dev, OutData_dev, N, time);
break;
case 16:
reduce6 <16> << <blocksPerGrid, threadsPerBlock, smemSize >> > (inData_dev, OutData_dev, N, time);
break;
case 8:
reduce6 <8> << <blocksPerGrid, threadsPerBlock, smemSize >> > ( inData_dev, OutData_dev, N, time);
break;
case 4:
reduce6 <4> << <blocksPerGrid, threadsPerBlock, smemSize >> > ( inData_dev, OutData_dev, N, time);
break;
case 2:
reduce6 <2> << <blocksPerGrid, threadsPerBlock, smemSize >> > (inData_dev, OutData_dev, N, time);
break;
case 1:
reduce6 <1> << <blocksPerGrid, threadsPerBlock, smemSize >> > (inData_dev, OutData_dev, N, time);
break;
}
}
else if (GuiYueWayCPU == 7)
{
reduce7 <128> << <blocksPerGrid, threadsPerBlock, smemSize >> > (inData_dev, OutData_dev, N, time);
}
else
{
printf("GuiYueWayCPU:%d define error!\n", GuiYueWayCPU);
return -1;
}
cudaMemcpy(OutData, OutData_dev, blocksPerGrid * sizeof(long long int), cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t), cudaMemcpyDeviceToHost);
long long int sum = 0;
for (int i = 0; i < blocksPerGrid; i++)
{
//printf("%d=>%u\n", i, OutData[i]);
sum += OutData[i];
}
printf("N:%d,result sum:%lld,GuiYueWay:%d,time: %d\n", N,sum, GuiYueWayCPU,time_used);
//释放GPU上的内存
cudaFree(inData_dev);
cudaFree(OutData_dev);
cudaFree(time);
//释放CPU上的内存
delete[]inData;
delete[]OutData;
return 0;
}
//*/
推荐阅读
-
CUDA 编程学习 <2> - 7 种还原算法优化方法 - 代码部分
-
F#探险之旅（二）：函数式编程（上）-函数式编程范式简介 F#主要支持三种编程范式：函数式编程（Functional Programming，FP）、命令式编程（Imperative Programming）和面向对象（Object-Oriented，OO）的编程。回顾它们的历史，FP是最早的一种范式，第一种FP语言是IPL，产生于1955年，大约在Fortran一年之前。第二种FP语言是Lisp，产生于1958，早于Cobol一年。Fortan和Cobol都是命令式编程语言，它们在科学和商业领域的迅速成功使得命令式编程在30多年的时间里独领风骚。而产生于1970年代的面向对象编程则不断成熟，至今已是最流行的编程范式。有道是“*代有语言出，各领风骚数十年”。 尽管强大的FP语言（SML，Ocaml，Haskell及Clean等）和类FP语言（APL和Lisp是现实世界中最成功的两个）在1950年代就不断发展，FP仍停留在学院派的“象牙塔”里；而命令式编程和面向对象编程则分别凭着在商业领域和企业级应用的需要占据领先。今天，FP的潜力终被认识——它是用来解决更复杂的问题的（当然更简单的问题也不在话下）。 纯粹的FP将程序看作是接受参数并返回值的函数的集合，它不允许有副作用（side effect，即改变了状态），使用递归而不是循环进行迭代。FP中的函数很像数学中的函数，它们都不改变程序的状态。举个简单的例子，一旦将一个值赋给一个标识符，它就不会改变了，函数不改变参数的值，返回值是全新的值。 FP的数学基础使得它很是优雅，FP的程序看起来往往简洁、漂亮。但它无状态和递归的天性使得它在处理很多通用的编程任务时没有其它的编程范式来得方便。但对F#来说这不是问题，它的优势之一就是融合了多种编程范式，允许开发人员按照需要采用最好的范式。 关于FP的更多内容建议阅读一下这篇文章：Why Functional Programming Matters（中文版）。F#中的函数式编程 从现在开始，我将对F#中FP相关的主要语言结构逐一进行介绍。标识符（Identifier） 在F#中，我们通过标识符给值（value）取名字，这样就可以在后面的程序中引用它。通过关键字let定义标识符，如： let x = 42 这看起来像命令式编程语言中的赋值语句，两者有着关键的不同。在纯粹的FP中，一旦值赋给了标识符就不能改变了，这也是把它称为标识符而非变量（variable）的原因。另外，在某些条件下，我们可以重定义标识符；在F#的命令式编程范式下，在某些条件下标识符的值是可以修改的。 标识符也可用于引用函数，在F#中函数本质上也是值。也就是说，F#中没有真正的函数名和参数名的概念，它们都是标识符。定义函数的方式与定义值是类似的，只是会有额外的标识符表示参数： let add x y = x + y 这里共有三个标识符，add表示函数名，x和y表示它的参数。关键字和保留字关键字是指语言中一些标记，它们被编译器保留作特殊之用。在F#中，不能用作标识符或类型的名称（后面会讨论“定义类型”）。它们是： abstract and as asr assert begin class default delegate do donedowncast downto elif else end exception extern false finally forfun function if in inherit inline interface internal land lazy letlor lsr lxor match member mod module mutable namespace new nullof open or override private public rec return sig static structthen to true try type upcast use val void when while with yield 保留字是指当前还不是关键字，但被F#保留做将来之用。可以用它们来定义标识符或类型名称，但编译器会报告一个警告。如果你在意程序与未来版本编译器的兼容性，最好不要使用。它们是： atomic break checked component const constraint constructor continue eager event external fixed functor global include method mixinobject parallel process protected pure sealed trait virtual volatile 文字值（Literals） 文字值表示常数值，在构建计算代码块时很有用，F#提供了丰富的文字值集。与C#类似，这些文字值包括了常见的字符串、字符、布尔值、整型数、浮点数等，在此不再赘述，详细信息请查看F#手册。 与C#一样，F#中的字符串常量表示也有两种方式。一是常规字符串（regular string），其中可包含转义字符；二是逐字字符串（verbatim string），其中的（"）被看作是常规的字符，而两个双引号作为双引号的转义表示。下面这个简单的例子演示了常见的文字常量表示： let message = "Hello World"r"n!" // 常规字符串let dir = @"C:"FS"FP" // 逐字字符串let bytes = "bytes"B // byte 数组let xA = 0xFFy // sbyte, 16进制表示let xB = 0o777un // unsigned native-sized integer，8进制表示let print x = printfn "%A" xlet main = print message; print dir; print bytes; print xA; print xB; main Printf函数通过F#的反射机制和.NET的ToString方法来解析“%A”模式，适用于任何类型的值，也可以通过F#中的print_any和print_to_string函数来完成类似的功能。值和函数（Values and Functions） 在F#中函数也是值，F#处理它们的语法也是类似的。 let n = 10let add a b = a + blet addFour = add 4let result = addFour n printfn "result = %i" result 可以看到定义值n和函数add的语法很类似，只不过add还有两个参数。对于add来说a + b的值自动作为其返回值，也就是说在F#中我们不需要显式地为函数定义返回值。对于函数addFour来说，它定义在add的基础上，它只向add传递了一个参数，这样对于不同的参数addFour将返回不同的值。考虑数学中的函数概念，F(x, y) = x + y，G(y) = F(4, y)，实际上G(y) = 4 + y，G也是一个函数，它接收一个参数，这个地方是不是很类似？这种只向函数传递部分参数的特性称为函数的柯里化（curried function）。 当然对某些函数来说，传递部分参数是无意义的，此时需要强制提供所有参数，可是将参数括起来，将它们转换为元组（tuple）。下面的例子将不能编译通过： let sub(a, b) = a - blet subFour = sub 4 必须为sub提供两个参数，如sub(4, 5)，这样就很像C#中的方法调用了。 对于这两种方式来说，前者具有更高的灵活性，一般可优先考虑。 如果函数的计算过程中需要定义一些中间值，我们应当将这些行进行缩进： let halfWay a b = let dif = b - a let mid = dif / 2 mid + a 需要注意的是，缩进时要用空格而不是Tab，如果你不想每次都按几次空格键，可以在VS中设置，将Tab字符自动转换为空格；虽然缩进的字符数没有限制，但一般建议用4个空格。而且此时一定要用在文件开头添加#light指令。作用域（Scope）作用域是编程语言中的一个重要的概念，它表示在何处可以访问（使用）一个标识符或类型。所有标识符，不管是函数还是值，其作用域都从其声明处开始，结束自其所处的代码块。对于一个处于最顶层的标识符而言，一旦为其赋值，它的值就不能修改或重定义了。标识符在定义之后才能使用，这意味着在定义过程中不能使用自身的值。 let defineMessage = let message = "Help me" print_endline message // error 对于在函数内部定义的标识符，一般而言，它们的作用域会到函数的结束处。 但可使用let关键字重定义它们，有时这会很有用，对于某些函数来说，计算过程涉及多个中间值，因为值是不可修改的，所以我们就需要定义多个标识符，这就要求我们去维护这些标识符的名称，其实是没必要的，这时可以使用重定义标识符。但这并不同于可以修改标识符的值。你甚至可以修改标识符的类型，但F#仍能确保类型安全。所谓类型安全，其基本意义是F#会避免对值的错误操作，比如我们不能像对待字符串那样对待整数。这个跟C#也是类似的。 let changeType = let x = 1 let x = "change me" let x = x + 1 print_string x 在本例的函数中，第一行和第二行都没问题，第三行就有问题了，在重定义x的时候，赋给它的值是x + 1，而x是字符串，与1相加在F#中是非法的。 另外，如果在嵌套函数中重定义标识符就更有趣了。 let printMessages = let message = "fun value" printfn "%s" message; let innerFun = let message = "inner fun value" printfn "%s" message innerFun printfn "%s" message printMessages 打印结果： fun value inner fun valuefun value 最后一次不是inner fun value，因为在innerFun仅仅将值重新绑定而不是赋值，其有效范围仅仅在innerFun内部。递归（Recursion）递归是编程中的一个极为重要的概念，它表示函数通过自身进行定义，亦即在定义处调用自身。在FP中常用于表达命令式编程的循环。很多人认为使用递归表示的算法要比循环更易理解。 使用rec关键字进行递归函数的定义。看下面的计算阶乘的函数： let rec factorial x = match x with | x when x < 0 -> failwith "value must be greater than or equal to 0" | 0 -> 1 | x -> x * factorial(x - 1) 这里使用了模式匹配（F#的一个很棒的特性），其C#版本为： public static long Factorial(int n) { if (n < 0) { throw new ArgumentOutOfRangeException("value must be greater than or equal to 0"); } if (n == 0) { return 1; } return n * Factorial (n - 1); } 递归在解决阶乘、Fibonacci数列这样的问题时尤为适合。但使用的时候要当心，可能会写出不能终止的递归。匿名函数（Anonymous Function） 定义函数的时候F#提供了第二种方式：使用关键字fun。有时我们没必要给函数起名，这种函数就是所谓的匿名函数，有时称为lambda函数，这也是C#3.0的一个新特性。比如有的函数仅仅作为一个参数传给另一个函数，通常就不需要起名。在后面的“列表”一节中你会看到这样的例子。除了fun，我们还可以使用function关键字定义匿名函数，它们的区别在于后者可以使用模式匹配（本文后面将做介绍）特性。看下面的例子： let x = (fun x y -> x + y) 1 2let x1 = (function x -> function y -> x + y) 1 2let x2 = (function (x, y) -> x + y) (1, 2) 我们可优先考虑fun，因为它更为紧凑，在F#类库中你能看到很多这样的例子。 注意：本文中的代码均在F# 1.9.4.17版本下编写，在F# CTP 1.9.6.0版本下可能不能通过编译。 F#系列随笔索引页面