用户
 找回密码
 立即注册
acone 该用户已被删除
发表于 2013-11-7 14:52:06
110377
自己编写了一个CUDA求解线性方程组的程序,能够正确实现,但是速度很不理想。我用的是联想T430笔记本自带的显卡NVS5400,只比CPU快4倍多一点,考虑到CPU是双核四线程,如果开四个线程粗粒度并行的话,GPU的速度只和CPU大致相当,比较失望。现在把最耗时的函数贴出来,求各位大佬给出一些优化策略来,不胜感激。

typedef double floatype;

//对 Ax=Y线性方程式进行变换,将系数矩阵变成三角矩阵
//pfDeviceA:系数矩阵;pfDeviceY:等式右边的Y向量;nVectorSize:方程组维数;nMainElem:对角线主元;
//nBlockThrdNum:每个块的线程数量;nGridThrdNum:每个Grid的线程数量
__global__ void cudaEliminate( floatype * pfDeviceA, floatype * pfDeviceY, int nVectorSize,
                                                        int nMainElem, int nBlockThrdNum, int nGridThrdNum )
{
        int nTid = blockIdx.x * nBlockThrdNum + threadIdx.x;
        while( nMainElem + nTid < nVectorSize - 1 )
        {
                int nOriginRowStart = nMainElem * nVectorSize;
                int nEliminatedRowStart = (nMainElem+nTid+1)*nVectorSize;
                floatype fFactor = pfDeviceA[ nEliminatedRowStart + nMainElem ]/pfDeviceA[ nOriginRowStart + nMainElem ];
                for( int i=nMainElem; i<nVectorSize; i++ )
                {
                        pfDeviceA[ nEliminatedRowStart + i ] -= fFactor * pfDeviceA[ nOriginRowStart + i ];
                }
                pfDeviceY[ nMainElem+nTid+1 ] -= fFactor * pfDeviceY[ nMainElem ];

                nTid += nGridThrdNum;
        }
}
(1)求上面函数的优化方法。
(2)上面的for循环如果用“动态并行”是不是很适合?速度大概能提高多少?
(3)求一个寻找矩阵中一列的最大值及其所在位置的实现方法,考虑到耗时占比比较少,可以不用归约,能并行实现就行。
(4)为什么我将floatype定义为float速度反而慢了四分之一?
使用道具 举报 回复
发表于 2013-11-7 15:05:20
LZ您好:

因为我并不熟悉线性代数相关的数值方法,因此无法告诉您您的方法是否合适而得当,亦无法给出通盘建议,下面仅就局部实现细节简要说下:

1:(本条略过,原因如上)。
2:目测您的for循环仅仅是读取,相减,写入的工作,这个直接用多个线程展开实现即可。
3:只考虑一列的话,你基本上只能用规约了,但是如果要处理很多列的话,你可以一个线程循环解决一列,相邻的一组线程解决相邻的多列,这样如果是行优先存储的矩阵,同时也是合并访问的,效果较好,实现也较为方便。
4:可能和您实现有关,尚不清楚。

大致如上,未及指出,请其他熟悉线性代数数值方法的网友/斑竹/总版主/NV原厂支持予以补充和讨论。

祝您编码顺利~
使用道具 举报 回复 支持 反对
发表于 2013-11-7 16:54:50
本帖最后由 seumx 于 2013-11-7 17:00 编辑

1. 可以用GPU/CPU混合的办法,比如你GPU要处理M列,可以改为GPU处理M-4m列,CPU处理4m列(开omp,一核算m列)
             cudaEliminate<<<grid block>>>(.....);
             parallel omp
             {
                 cpuEliminate(.....);
             }
             再主机和设备同步(不好意思函数名字我忘记怎么拼写了,抱歉)
假如你矩阵规模不大,不要开omp,开了反而慢,CPU直接串行算
这个方法可以应用在很多高等代数的优化上~但这个方法不是最优的,还一种look ahead的CPU/GPU异步模式,可以去看看volkov的论文,比较麻烦
2,3. 如ice哥哥所述,但是规约的时候你可以一个线程比较多个元素的最值,比如比较4个,8个,16个.....,总有一个可以达到最最优化,这样你开适量的线程就OK的,号称世界上最快的MAGMA中的矩阵向量乘法就是这么尝试出来的
使用道具 举报 回复 支持 反对
发表于 2013-11-8 09:11:28
2:目测您的for循环仅仅是读取,相减,写入的工作,这个直接用多个线程展开实现即可。

果真如此!我将它展开后再调整了一下block size,比CPU的速度由4倍多提升到了14倍!现将展开后的代码以及调用它的代码粘贴如下:
__global__ void cudaEliminate( floatype * pfDeviceA, floatype * pfDeviceY, int nVectorSize,
                                                        int nMainElem, int nBlockThrdNum, int nGridThrdNum )
{
        int nTid = blockIdx_x * nBlockThrdNum + threadIdx_x;
        while( nTid < (nVectorSize-nMainElem)*(nVectorSize-nMainElem-1) )
        {
                int nSubBlockWidth = nVectorSize - nMainElem;
                int nRowOfThisThrd = nMainElem + nTid/nSubBlockWidth + 1;
                int nColOfThisThrd = nMainElem + nTid%nSubBlockWidth;//nTid-(nTid/nSubBlockWidth)*nSubBlockWidth + 1;//
                floatype fFactor = pfDeviceA[ nRowOfThisThrd * nVectorSize + nMainElem ]
                                                                        / pfDeviceA[ nMainElem * nVectorSize + nMainElem ];
                       
                if( 0 == nTid % nSubBlockWidth )//if( nTid == (nTid/nSubBlockWidth)*nSubBlockWidth )//
                {
                        pfDeviceY[ nRowOfThisThrd ] -= fFactor * pfDeviceY[ nMainElem ];
                }
                else
                {
                        pfDeviceA[ nRowOfThisThrd * nVectorSize + nColOfThisThrd ]
                                -= fFactor * pfDeviceA[ nMainElem * nVectorSize + nColOfThisThrd ];
                }

                nTid += nGridThrdNum;
        }
}//*/
调用部分:
bool CudaSolveEquations( floatype * pfDeviceA, floatype * pfDeviceY, int nVectorSize )
{
        int i, nBlockThrdNum, nBlockNum;
        for( i=0; i<nVectorSize-1; i++ )
        {
                //nBlockNum = GetBlockNumAndBlockThrdNum( nVectorSize-1, & nBlockThrdNum );
                nBlockNum = GetBlockNumAndBlockThrdNum( (nVectorSize-i)*(nVectorSize-i-1), & nBlockThrdNum );
                cudaEliminate<<<nBlockNum, nBlockThrdNum>>>( pfDeviceA, pfDeviceY, nVectorSize, i, nBlockThrdNum, nBlockNum * nBlockThrdNum );
                cudaDeviceSynchronize();       
        }
......

(1)继续求优化措施。
(2)内存合并访问是否和显卡的位宽有关?比如显卡位宽为64位,如果顺序读入8个双精度数,是否就可以一次读入?
(3)怎样去掉那恼人的C4819警告?
(4)回seumx:你说的方法对我来说就是天书,我目前不准备用混合编程,不过还是要谢谢你了。
使用道具 举报 回复 支持 反对
发表于 2013-11-8 11:05:17
acone 发表于 2013-11-8 09:11
果真如此!我将它展开后再调整了一下block size,比CPU的速度由4倍多提升到了14倍!现将展开后的代码以及 ...

合并访问与显存位宽无关,需要相邻线程访问一段连续的(最好是满足对齐要求的)显存。
使用道具 举报 回复 支持 反对
发表于 2013-11-12 18:22:27
接着问问题:
(1)解线性方程组的时候当维数低于1500时工作正常,超过后出现unknown error错误,问一下有经验的版主,一般情况下是哪里出了问题?
(2)维数为1000的双精度,与CPU的答案大概只有前4、5位数字是相同的,单精度的就更低了,正常么?
使用道具 举报 回复 支持 反对
发表于 2013-11-12 18:30:34
acone 发表于 2013-11-12 18:22
接着问问题:
(1)解线性方程组的时候当维数低于1500时工作正常,超过后出现unknown error错误,问一下有 ...

LZ您好:

1:这个一般是访存越界了,您可以上一下nsight或者cuda memchecker以便为您快速指出越界位置。

2:这个不好判定,因为浮点数是不满足实数的交换律和结合律的,不同的计算过程(比如CPU上的循环VS GPU上的规约)可能会产生不同的结果,并且这个也不好说哪些是更为准确的,因为CPU结果也只是“一种浮点运算顺序”下的结果,同时,还可能因为不同的中间过程精度,比如CPU上如果使用X87指令,那么会有80bit的精度,到最后才截断,但是如果SIMD指令则没有这么好的中间精度;对于乘加运算,GPU有FMA计算,保留了较好的中间精度,而CPU上一般只能分开计算(最新的CPU才支持FMA的指令集)。

当然如果是您代码实现考虑不周也是有可能的。

大致如此,祝您编码顺利~
使用道具 举报 回复 支持 反对
发表于 2013-11-29 11:25:54
果真如此!我将它展开后再调整了一下block size,比CPU的速度由4倍多提升到了14倍!

汇报一下:收到别的帖子启发,我用每个block处理一行数据,那个乘数因子fFactor申请为共享内存,速度进一步提升了70%-80%,比CPU快了20倍多。
这个核函数就是一个高斯消元过程:一个n*n的方阵,将第二行减去乘了乘数因子的第一行,使得第二行的第一个元素为0;将第三行减去乘了乘数因子的第一行,使得第三行的第一个元素为0;以此类推,最终使得第一列除了第一行其他元素都为0.
由于这个核函数是个主耗时函数,且高斯消元的算法复杂度为N的立方,还是请各位大侠继续提出优化措施,我将感激不尽。将第一行映射为纹理内存是否可行?怎么做?
使用道具 举报 回复 支持 反对
发新帖
您需要登录后才可以回帖 登录 | 立即注册