归约算法1(reductionInner)

归约可以分为两个层次。一个层次是kernel层面的。在下图中需要对一维的数据求和，假设整个数组是在同一个group中，自然可以用barrier来做同步，每次迭代问题规模减半。这是一个层次。过程如下图。

此外不论是opencl还是cuda的GPU编程(硬件)模型都不支持全局线程同步（因为本质所有的线程在大尺度上是串行，小尺度上才是并行的。即硬件计算单元是有限）。所以没办法在一次kernel调用中得到最终解。再kernel外还有一次迭代的发射GPU任务的过程。网上的论述往往不介绍这一部分，本文会一起把代码列出来。

kernel代码的流程：先每个线程load一个元素到一个local mem中。这样整个group就可以高效的操作load进来的数据。因为我们的问题是二维求解，所以第一步现在X方向每一行用一次归约算法。所有的数据最终被加到了第一列里。在对第一列做一次归约算法最终整个group的求和结果就放在group内的0号线程里。我们让其把内存到结果里即可。

__kernel void reductionInner(__global  float *src,
					int blockH, // group height
					int blockW, // group width
					int srcStride,
					int dstStride,
					__global float *dst,
					__local float *aLocalArray /// local mem 
	int x_gid = get_global_id(0);
	int y_gid = get_global_id(1);
	int x_lid = get_local_id(0);
	int y_lid = get_local_id(1);
	// load data
	aLocalArray[y_lid*blockW + x_lid] = src[y_gid*srcStride + x_gid];
	barrier(CLK_LOCAL_MEM_FENCE);
	// reduction in x
	int step = blockW >> 1;
	for (; step > 0; step = (step >> 1))
		if (x_lid < step)
			aLocalArray[y_lid*blockW + x_lid] += aLocalArray[y_lid *blockW + x_lid + step];
		barrier(CLK_LOCAL_MEM_FENCE);
	// reduction in y
	step = blockH >> 1;
	for (; step > 0; step = (step >> 1))
		if (y_lid < step && x_lid == 0)
			aLocalArray[y_lid*blockW] += aLocalArray[(y_lid + step) *blockW];
		barrier(CLK_LOCAL_MEM_FENCE);
	// save
	if (x_lid == 0 && y_lid == 0)
		int x_bid = x_gid / get_local_size(0);
		int y_bid = y_gid / get_local_size(1);
		dst[y_bid*dstStride + x_bid] = aLocalArray[0];
	return;
比如目标大小为4096x4096, group 设置大小为8x8,则我们需要调用4次kernel。每轮问题规模都减小一半。每轮global size为输入的尺寸。注意和下面另一种方法做对比。 
4096 -> 512 - > 64 -> 8 -> 1
这个例子中我们要分配5块内存(4096x4096、512x512…1)。显然这样很没效率，每一次新的轮次里上一个轮次里的输入已经不需要保存了，我们直接复用这个内存作为本轮次的输出。这样我们只要分配前两块内存即可，用A B表示。若迭代次数为奇数则结果保存在B中，偶数就在A中。 
 A	     B      A     B    A
4096 -> 512 - > 64 -> 8 -> 1
以上就是host端代码需要注意的地方。下面为代码。 
		cl_kernel kernel = env.OEW_GetKernel("reductionInner");
		size_t curH = h;
		size_t curW = w;
		bNum = w / blockSize;
		bTurn = true;
		int iter = 0;
		timerWrite.Start();
		bufferGPU0.OIW_WriteBuffer(static_cast<void*>(pBufferHost0), w*h*sizeof(float));
		timerWrite.End();
		while (bNum > 0)
			//enTimer timer("testInnerReduction", 1);
			cl_mem cl_input, cl_output;
			size_t srcStride, dstStride;
			if (bTurn)
				cl_input = bufferGPU0.OIW_GetData();
				cl_output = bufferGPU1.OIW_GetData();
				srcStride = ori_srcStride;
				dstStride = ori_dstStride;
				bTurn = false;
				cl_output = bufferGPU0.OIW_GetData();
				cl_input = bufferGPU1.OIW_GetData();
				dstStride = ori_srcStride;
				srcStride = ori_dstStride;
				bTurn = true;
			size_t idx = 0;
			err |= clSetKernelArg(kernel, idx++, sizeof(cl_mem), &cl_input);
			err |= clSetKernelArg(kernel, idx++, sizeof(cl_int), &blockSize);
			err |= clSetKernelArg(kernel, idx++, sizeof(cl_int), &blockSize);
			err |= clSetKernelArg(kernel, idx++, sizeof(cl_int), &srcStride);
			err |= clSetKernelArg(kernel, idx++, sizeof(cl_int), &dstStride);
			err |= clSetKernelArg(kernel, idx++, sizeof(cl_mem), &cl_output);
			err |= clSetKernelArg(kernel, idx++, sizeof(float)*blockSize*blockSize, NULL);
			// run
			size_t global_offset_size[] = { 0, 0 };
			size_t global_work_size[] = { curW, curH };
			size_t global_local_size[] = { blockSize, blockSize };
			cl_event event;
			err |= clEnqueueNDRangeKernel(env.OEW_GetCommand(),
				kernel,
				global_offset_size,
				global_work_size,
				global_local_size,
				NULL,
				&event
			clWaitForEvents(1, &event);
			//update 
			iter++;
			LOG_DEBUG_INFO("iter NO.%2d [%3d,%3d]->", iter, curH, curW);
			curH = bNum;
			curW = bNum;
			bNum = curH / blockSize;
			LOG_DEBUG_INFO("[%3d,%3d] \n", curH, curW);
	timerTotal.End();
	auto total = timerTotal.Print();
	auto wtime = timerWrite.Print();
	LOG_WARRING_INFO("Run Time : %lf ms\n", total - wtime);
大致解释一下，bNum为每轮输出的尺寸最后一轮为1,变为0的时候迭代结束。每轮输入和输出的尺寸都要更新一次，每轮输入和输出内存的轮换用bTurn控制。整个工程会在文末给出链接，所以这里先只看主要代码即可。 
这个就是网上大部分入门文章所说的归约。我们来分析一下他们的效率。在kernel层面上的归约中，没轮迭代都会少一半的线程干活，空等(idle)别的线程干活。到最后一轮的时候只有一个线程干活，剩下的线程都在等待。感觉上是很没有效率啦。只有0号线程最勤劳，每轮迭代都在干活。。
  
cuda的教程里有优化手段，但是还是少不了空等。而且像最后的模板方法opencl也用不了。那我们试试让每个线程都不闲着吧！ 
归约算法2(reductionOuter)
 
既然分配那么多线程但是却闲着不干活，那就学学黑心老板吧。每个线程安排N多活。我们安排每个线程计算NxN个元素。没有local mem，没有同步。干就完事了。 
__kernel void reductionOut(__global  float *src,
	int blockH,
	int blockW,
	int srcStride,
	int dstStride,
	__global float *dst
	int x_gid = get_global_id(0);
	int y_gid = get_global_id(1);
	int x_start = x_gid*blockW;
	int y_start = y_gid*blockH;
	int x_end = x_start + blockW;
	int y_end = y_start + blockH;
	float fsum = 0.0f;
	for (; y_start < y_end; y_start++)
		__global  float *srcCur = src + y_start*srcStride;
		for (int x = x_start; x < x_end; x++)
			fsum += srcCur[x];
	dst[y_gid*dstStride + x_gid] = fsum;
	return;
与实现1不同，这里global size是每轮host端迭代的输出尺寸。有些卷积那个意思。外部代码调用代码类似，不贴了。具体看链接里的代码吧。 
 A	     B      A     B    A
4096 -> 512 - > 64 -> 8 -> 1
这份代码解决了偷懒线程的问题，但是会有一个并行度不够的问题。首先总线程数就没上个实现多。另外在问题规模变小的时候，很显然没有充分利用GPU的并行度。最极端的在最后一轮里只分配了一个线程来对64个值求和。不过此时问题规模也小，花费时间很小，关键性能消耗都在前几轮。 
GTX 1050实测。尺寸为4096，group size为16，所以有三个轮次。 
归约1的三个轮次性能如下： 
 round: ave 10.878376 ms in 1 times
 round: ave 0.177819 ms in 1 times
 round: ave 0.089080 ms in 1 times 
归约2的三个轮次性能如下： 
 round: ave 7.668755 ms in 1 times
 round: ave 0.138569 ms in 1 times
 round: ave 0.076111 ms in 1 times 
而在Adreno509上，归约2的会比归约1快几倍(几百ms降低到一百多ms)。
 可以看到在安排的线程数大于实际硬件并行度的时候，不让任何一个核闲着才是正途。所以在数组求和这个问题上，是否真的用归约1来做需要具体分析。 
完整代码：https://github.com/duzeyan/AOCL.git
                    任务分析矩阵内求和任务在GPU上实现并不直观。假设每个线程load一个数据并累加到最后的结果标量中，这样的处理过程会使得整个任务退化到串行处理。原因是同一时刻只有一个线程可以操作目的标量，其他线程只能等待。解决此类问题可以套用归约算法。若读者对这个问题感兴趣，请先对归约算法有个大概的了解。下文不会介绍基本概念，主要是讨论实现细节和优化点。归约算法1(reductionInner)归约可以分为两个层次。一个层次是kernel层面的。在下图中需要对一维的数据求和，假设整个数组是在同一个group中，自然
				2008年，苹果公司向KhronosGroup提交了一份关于跨平台计算框架的草案，该草案由苹果公司开发，并与AMD、IBM、Intel和NVIDIA公司合作逐步晚上。这个跨平台计算框架就是OpenCL。20088年12月8日，OpenCL1.0技术规范发布。2010年6月14日，OpenCL1.1发布，2011年11月19日，OpenCL1.2发布，2013年11月19日，OpenCL2.0分布。OpenCL是一个异构并行计算平台编写程序的工作标准，此异构计算可映射到CPU、GPU、DSP和FPGA等计算设备。OpenCL提供了底层硬件结构的抽象模型，旨在提供一个通用的开发的API。开发人员可
				初学OpenCL，写的不好的地方还请见谅。
本次说的是一个累加和的问题，给定一个一维数组A[n]，然后计算B[n]，其中，B[0]=A[0],B[1]=A[1]+B[0],B[2]=A[2]+B[1]...依此类推。
在我们的CPU中计算的时候使用一个循环来计算，程序是这样的：
B[0]=A[0];
for(int i=1;i<n;i++){
B[i]=B[i-1]+A[i];
5）在系统中正确的组件上按正确的顺序执行内核。
6）收集最终结果。
这些步骤通过OpenCL中的一系列API再加上一个面向内核的编程环境来完成。我们将采用一种“分而治之”的策略解释以上步骤的所有工作。我们把问题分解为以下模型：
1）平台模型 (platform model)：
				数值归并运算实质上就是求N项常数的总和的一种并行计算方法。
设数组的长度为N，若依次去运算数组中的每个元素(sum+=array[i])，显然时间复杂度就是O(N)；
对于数值归并运算还有一种更快捷的方法，可以通过将数组第i个元素与第（i+N/2）个元素相加，得到一个长度为N/2的数组；
得到的新数组中的每个元素，即为前一数组中元素两两相加所得（若原数组为奇数个，则保留最后一个元素至新数组）；
#include "device_functions.h"
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "stdlib.h"
#include <iostream>
using namespace ...
				原文地址:https://www.cnblogs.com/xudong-bupt/p/3586518.html
Reduction操作：规约操作就是由多个数生成一个数，如求最大值、最小值、向量点积、求和等操作，都属于这一类操作。
有大量数据的情况下，使用GPU进行任务并行与数据并行，可以收到可好的效果。
group同步：OpenCL只提供了工作组内的各线程之间的同步机制，并没有提供所有线程的...
				高斯消元是一种串行的算法，它的计算复杂度主要集中在矩阵上，因此它不太适合并行化。但是，可以通过一些方法对其进行优化，以提高其计算效率。
一种常用的方法是使用并行计算，将矩阵分成多个小块，然后在多个处理器上并行处理，以加快计算速度。这种方法在大规模的矩阵计算中非常有效，但也需要合适的硬件和软件架构来支持。
另一种方法是使用GPU加速计算，GPU拥有大量的计算核心和内存，可以同时处理多个数据和运算。可以使用CUDA或OpenCL等技术将高斯消元算法移植到GPU上进行加速计算，从而显著提升计算效率。
除了并行计算和GPU加速外，还可以使用一些优化技术来加速高斯消元算法。例如，使用矩阵重排列、矩阵分块、缓存优化等方法，可以减少计算过程中的数据访问次数和内存访问延迟，从而提高算法的计算效率。
总之，高斯消元算法并不容易并行化，但是可以通过一些优化技术和硬件支持来加速计算，以提高算法的实用性和适用性。