最近在学习LU的并行加速，从paper中得到了一些idea，就想着用GPU来实现一下。学习CUDA的过程中踩了不少坑，不过最终还是完成了测试。

一、LU基本算法

1、LU 分解是计算机做矩阵运算过程中重要的一步，通过将矩阵分解为一个上三角矩阵U和下三角矩阵L，能够有效的缩短计算时间。

LU分解的计算过程如下，采用高斯消元法。

2、基本算法

• 每一次循环都将A的第i行和第i列更新为L\U的一部分。
• 每一次循环分为三个部分。如果先计算左上角的元素，则下侧的L矩阵部分和右侧的U矩阵部分可以同时运算，没有干扰。比如当i==1时，因为更新L的第一列（l32、l42）需要先得到u22的值，而u22要通过更新U的第二行得到，所以如果先计算出u22，L和U的更新就可以并行了。

void lud_base(int *a, int size)
{int i,j,k;int sum;for (i=0; i<size; i++){//先计算左上角的U元素sum=a[i*size+i];for (k=0; k<i; k++) sum -= a[i*size+k]*a[k*size+i];a[i*size+i]=sum;//计算下侧的L矩阵部分for (j=i+1;j<size; j++){sum=a[j*size+i];for (k=0; k<i; k++) sum -=a[j*size+k]*a[k*size+i];a[j*size+i]=sum/a[i*size+i];}//计算右侧的U矩阵部分for (j=i+1; j<size; j++){sum=a[i*size+j];for (k=0; k<i; k++) sum -= a[i*size+k]*a[k*size+j];a[i*size+j]=sum;}}
}

上述代码其实存在一个问题，对左上角元素的更新是一个累计操作。CUDA进行累计操作有两种方法，但是都不太好用

1. 单thread执行
2. 归约的编程方式

所以我对代码再次进行了修改，还是以i==1时为例，如果我们直接对L和U进行更新呢，得到的结果是这样的：
    U（第1行）：U22、U23、U24
    L（第1列）：l32**U22、l42*U22
所以我们再加一个步骤将L的每一项都除以U22即可，而这一步可以很容易的并行。代码如下：

void lud_base(int *a, int size)
{int i,j,k;int sum;for (i=0; i<size; i++){//计算下侧的L矩阵部分for (j=i+1;j<size; j++){sum=a[j*size+i];for (k=0; k<i; k++) sum -=a[j*size+k]*a[k*size+i];a[j*size+i]=sum;}//计算右侧的U矩阵部分for (j=i; j<size; j++){sum=a[i*size+j];for (k=0; k<i; k++) sum -= a[i*size+k]*a[k*size+j];a[i*size+j]=sum;}//对下侧的L矩阵后续处理for (j=i+1;j<size; j++)a[j*size+i]=sum/a[i*size+i];}}
}

二、LU分解常用的并行化算法

    LU分解的并行化算法，可以分类为right-looking和left-looking。我采用right-looking算法进行测试。

Left-looking算法存在循环间的依赖关系，因为左侧的1，2，3列都要对第4列的数据进行更新，所以不能并行。
right-looking算法虽然并行度非常高，但是存在大量重复的访问和读取，因为要不断更新右侧的子矩阵。

left-looking算法如下：

right-looking算法如下：

三、CUDA加速

先记录一下花了好久才找出的Bug和解决的问题（学习还不深入，如果有误烦请指出）：

1. cuda支持累计操作，但是必须在单线程(sp)上，而且左值必须是声明的局部变量，如果左值是数组地址的话，那么每个线程都会去该地址寻址并写回，产生冲突导致结果错误。所以累加操作要用归约编程。
2. if-else内部不能有__syncthreads()，因为部分线程不会运行。
3. 线程束分化。一个warp中如果部分thread执行if块内的代码，另一部分线程condition不成立。则不成立的thread限制，等待其他thread执行完毕。会造成严重的性能下降。所以最好通过手动调整让一个warp中的线程都执行if或else。
4. 一些size变量要设置为全局变量才能被GPU识别，所有数据都要手动搬入显存。
5. 要根据显卡的硬件参数来设置block和grid，包括SM最大block数，最大线程数目等。
6. 出现《核心段错误》是由于GPU占用率接近100%了。

right-looking算法的CUDA代码

 __global__ void  lud_right_looking(float *m)
{__shared__ float shadow[BLOCK_SIZE][BLOCK_SIZE];for (int tix = threadIdx.x; tix<BLOCK_SIZE; tix += blockDim.x) shadow[tix][threadIdx.y] =  m[tix*BLOCK_SIZE+threadIdx.y];__syncthreads();for (int k=0; k< BLOCK_SIZE-1; k++){//这行语句不适合放进下面的for循环中，多执行一次就多浪费一次资源。if (threadIdx.y>k && threadIdx.x==0)shadow[threadIdx.y][k]=shadow[threadIdx.y][k]/shadow[k][k];__syncthreads();for (int tix = threadIdx.x; tix<BLOCK_SIZE; tix += blockDim.x) {if(tix>k && threadIdx.y>k)shadow[tix][threadIdx.y] -= shadow[tix][k]*shadow[k][threadIdx.y];__syncthreads();}}for (int tix = threadIdx.x; tix<BLOCK_SIZE; tix += blockDim.x) m[tix*BLOCK_SIZE+threadIdx.y]=shadow[tix][threadIdx.y];

基本算法的CUDA代码

__global__ void  lud_base_cuda(float *m)
{__shared__ float shadow[BLOCK_SIZE][BLOCK_SIZE];int i,j;for (int tix = threadIdx.x; tix<BLOCK_SIZE; tix += blockDim.x) shadow[tix][threadIdx.y] =  m[tix*BLOCK_SIZE+threadIdx.y];__syncthreads();for(i=0; i < BLOCK_SIZE-1; i++){if ( threadIdx.y>i && threadIdx.x==0 )//计算下侧的L部分{for(j=0; j < i; j++)shadow[threadIdx.y][i] -= shadow[threadIdx.y][j]*shadow[j][i];}else if(threadIdx.y>=i   && threadIdx.x==1 ) //计算上侧的U部分{for(j=0; j < i; j++)shadow[i][threadIdx.y] -= shadow[i][j]*shadow[j][threadIdx.y];}__syncthreads();if(threadIdx.y>i && threadIdx.x==0)//对L进行后续处理{shadow[threadIdx.y][i] /= shadow[i][i];}__syncthreads();}for (int tix = threadIdx.x; tix<BLOCK_SIZE; tix += blockDim.x) m[tix*BLOCK_SIZE+threadIdx.y]=shadow[tix][threadIdx.y];
}

输入矩阵大小为64*64，运行时间如下

lud_right_looking的并行度是最高的，所以分配了1024的线程， lud_base_cuda并行度不高所以只需要128线程（实际上还有可提取的部分，将累计操作用归约实现）。综上，虽然lud_right_looking会产生额外的大量访问，但是在足够资源下并行度是最高的，性能优于ud_base_cuda。因为数据量太小，所以性能不如CPU。

引用：
1、2005.LU-GPU Efficient Algorithms for Solving Dense Linear Systems on Graphics Hardware
2、2016.GPU-Accelerated Parallel Sparse LU Factorization
3、2017.GPU-Based Batch LU-Factorization Solver for Concurrent
4、2010.Blocking LU Decomposition for FPGAs
5、求大神帮忙，怎么在MATLAB上用LU分解法？ – 琦钒的回答 – 知乎