并行与分布式计算导论——HW04

碎碎念：别问我为啥没有03orz 这学期网课太容易松懈了…等考试周过去看看能不能补一补orz（flag*1）
题目：
⼩明⽤OpenCL编写⼀份计算矩阵-向量乘法的代码，AX=y, 其中A是N*N的矩阵，X是⻓度为N的向量。 他的做法是对矩阵按⾏分组，每个work group计算每⾏的结果。在每个work group内，他⼜设置了16 个work item，其中第i个work item计算下标模16余i的元素的乘法与加法，最后使⽤规约 （reduction）来计算最终结果。
下⾯给出了部分主程序和kernel代码，未列出的部分（初始化或后⾯的代码）都可以假设是正确的。
1 . 在空缺部分补全代码，完成代码功能。
2 . 分析这样设计下的性能

vmul_kernel.cl:
=========================================================================
 __kernel void vmnul_kernel(int N, float *A, __global float *X, __global float *y){
    int i = get_group_id(0);
    int tid = get_local_id(0);  
    int tdim = get_local_size(0);  
    float l_val = 0.0;  
    for (int j = tid; j < N; j += tdim){     
    	_____________________________ //1  
    }  
    __local float tmp[16];  
    tmp[tid] = l_val;  
    barrier(CLK_LOCAL_MEM_FENCE);  
    for(int s = 8; s > 0; s >>= 1){
        if(tid < s){
              __________________________ //2      
              barrier(CLK_LOCAL_MEM_FENCE);    
        }  
    }  
    if(tid == 0) y[i] = tmp[tid]; 
}
main.cpp
=============================================================================
  cl_mem bufA = clCreateBuffer(    
  	context, CL_MEM_READ_ONLY, N*N*sizeof(float), NULL, &status);  
  cl_mem bufX = clCreateBuffer(    
  	context, CL_MEM_READ_ONLY, N*sizeof(float), NULL, &status);  
  cl_mem bufY = clCreateBuffer(    
  	context, CL_MEM_WRITE_ONLY,N*sizeof(float), NULL, &status);
 
  status = clSetKernelArg(vmul_kernel, 0, sizeof(int), N);  
  status = clSetKernelArg(vmul_kernel, 1, sizeof(cl_mem), bufA);  
  status = clSetKernelArg(vmul_kernel, 2, sizeof(cl_mem), bufX);  
  status = clSetKernelArg(vmul_kernel, 3, sizeof(cl_mem), bufY);
  status = clEnqueueWriteBuffer(vmul_kernel, bufA, CL_FALSE, 0,
      N*N*sizeof(float), A, 0, NULL, NULL);   
  size_t globalWorkSize[1];  
  size_t localWorkSize[1];  
  globalWorkSize[0] = ____________________; //3  
  localWorkSize[0] = ____________________;  //4  
  status = clEnqueueNDRangeKernel(cmdQueue, vmul_kernel,
                 1, NULL, globalWorkSize, localWorkSize, 0, NULL, NULL);
  clEnqueueReadBuffer(cmdQueue, bufY, CL_TRUE,
                 0, N*sizeof(float), Y, 0, NULL, NULL);

填空题。这次罗老师非常和蔼地没让run一下算什么程序实验报告里让写的speedup等等，只用写填空题&时间复杂度就ok了，（俗称嘴炮时间x）
先上答案：

一、补全代码：
1.l_val+=A[i*N+j]*X[j];
2.tmp[id]+=tmp[id+s];
3.N*16
4.16

二、性能分析（时间复杂度）
上述题目的串行时间复杂度应为O(N2)，而该并行程序时间复杂度： 每个work group的工作量为O(N)，每个item工作量为O(N/16)，故并行程序时间复杂度为O(N/16)=O(N)。可以看出，相比于串行程序的O(N2)，并行程序的时间复杂度有了较大改进。

分析过程：

感觉就是整个要算N^2次乘法，然后分了N个work group，每个group算y[i]=A[I]*X；然后又细分为item，每个item又隔16个算一个A[i][j]*X[j]；第一个框里的部分就是每个item算的那个，但是因为A不是以A[][]形式传过来的，而是把它从一个矩阵弄成了一条(N*N，1)大小的向量（如下图：可看出A是N*N float）

所以这里A[i][j]相对应成向量里的就应该是A[i*N+j]

第二个框框就是规约，就是算y[i]=之前算的一堆数的加和，例如s=8如下图，所有9~16的数被加到1~8上，s=4把5~8加到1~4上，然后直到s=1的时候tmp[0]=之前的所有的和，这就是y[i]了，所以底下令y[i]=tmp[0]

第三第四个填空的global work size和local work size分别是整个程序里的work item数和一个group里的work item数，所以分别是N*16和16。

分析复杂度就更简单了，答案里写的应该挺清楚了，就不做赘述了。