最近在研究GPU求解下三角矩阵方程Ax=b，
计算步骤：
for n=1:N
1对A的第n行的前n-1个元素和x的前n-1个元素做一次向量点乘；
2将b的第n个元素减去前一步内积的结果；
3最后将第2步中的结果除以A的第n行的第n个元素。
end
实现方法：
由于A为稀疏矩阵，存储时采用行稀疏存储，在global memory的三个矩阵value，index，count中分别记录了A每一行中非0元素的值、列坐标和第一个非0元素在value或index矩阵中的位置。
为了避免线程空置，循环在内核函数之外实现，每次循环内部的计算由内核函数完成，每次循环开始时将块（block）数设置为需要求解的方程个数，块中线程个数设置为A的第n行中非0元素的个数。第一步中先多线程并行完成两个向量的元素相乘，并将结果写入shared memory，然后将进行规约求和。最后再让第一个线程完成第二步和第三步。编程中避免了bank conflicts并实现了global memory的coalesces访问（即访问shared memory和global memory时让连续的线程访问连续的存储空间）；尽量减少了对global memory的访问；整个计算过程中内存和GPU之间无数据交换。
内核函数
global void HH_device(float* b,float* v_gpu,int* i_gpu,int* worder_gpu,int* wst_gpu)
{
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i,k,l,n,m,count1,count2,count3,ii,jj,windex;
float temp1r,temp2r,temp3r,temp4r;
windex=worder_gpu[bid+wst_gpu[0]];
i=i_gpu[0];
count1=blockDim.x;
count2=2windexnodecount_gpu[0];
count3=windex*(nodecount_gpu[0]+1);
shared float Mv[2050];
ii=tex1Dfetch(rT_Hcount2,i+count3);
jj=tex1Dfetch(rT_Hcount2,i+1+count3)-ii-1;
if (tid==0)
{
Mv[0]=0;
Mv[count1]=0;
Mv[2050-2]=tex1Dfetch(rT_HMre2,ii+jj);
Mv[2050-1]=tex1Dfetch(rT_HMim2,ii+jj);
}
// 向量元素相乘
if (tid<jj)
{
temp1r=tex1Dfetch(rT_HMre2,ii+tid);
temp2r=tex1Dfetch(rT_HMim2,ii+tid);
k=tex1Dfetch(rT_Hindex2,ii+tid);
temp3r=v_gpu[k+count2];
temp4r=v_gpu[k+count2+nodecount_gpu[0]];
Mv[tid]=temp1rtemp3r-temp2rtemp4r;
Mv[tid+count1]=temp1rtemp4r+temp2rtemp3r;
}
__syncthreads();
//规约求和
k=ceilf(log2f(jj*1.0));
m=jj;
n=jj;
for(l=0;l<k;l++)
{
m=floorf(n/2.0);
if(tid<m)
{
Mv[tid]=Mv[tid]+Mv[n-tid-1];
Mv[tid+count1]=Mv[tid+count1]+Mv[count1+n-tid-1];
}
n=n-m;
__syncthreads();
}
//第二步和第三步中的减法和除法
if (tid==0)
{
temp1r=Mv[2050-2]*Mv[2050-2]+Mv[2050-1]*Mv[2050-1];
temp3r=b[i+count2]-Mv[0];
temp4r=b[i+count2+nodecount_gpu[0]]-Mv[count1];
v_gpu[i+count2]=(Mv[2050-2]*temp3r+Mv[2050-1]*temp4r)/temp1r;
v_gpu[i+count2+nodecount_gpu[0]]=(Mv[2050-2]*temp4r-Mv[2050-1]*temp3r)/temp1r;
}
//行数加1
i=i+1;
if (tid==0)
{
i_gpu[0]=i;
}
}
主函数
tempsize[0]=0;
cudaMemcpy(i_gpu,tempsize, sizeof(int) *1,cudaMemcpyHostToDevice);
grid.x=wcount-wst[0];
for (i=0;i<nodecount;i++)
{
block.x=Hblock2;
HH_device<<<grid, block>>>(r_gpu,p_gpu,i_gpu,worder_gpu,wst_gpu);
}
//r_gpu为b，p_gpu为x。HMre为A矩阵非0元素实部，HMim为虚部，Hindex为非0元素列坐标，Hcount为每行非0元素起始位置，后面4个矩阵已事先拷入显存并绑定为纹理。
运行结果：
GTX580 3G 对比 I7 3770 RAM 32G，当只求解一个方程时，GPU完全处于劣势，求解方程个数增加到20个左右时GPU并行求解速度和CPU串行求解速度基本相当，当求解方程个数为54时，加速比还不到2，方程数为2592时加速比才5左右，似乎并行计算的效率不高啊。
PS：其实实际计算时间很短，每个方程的未知数个数为15000左右，求解2592个方程总的计算时间为2 秒，但不知是不是I7性能变态，搞得加速比很难看。
不知各位对这个问题有什么好办法没有？