新人求线性方程组的优化策略

system · 2013 年11 月 7 日 06:52

自己编写了一个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速度反而慢了四分之一？

system · 2013 年11 月 7 日 07:05

LZ您好：

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

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

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

祝您编码顺利~

system · 2013 年11 月 7 日 08:54

可以用GPU/CPU混合的办法，比如你GPU要处理M列，可以改为GPU处理M-4m列，CPU处理4m列（开omp，一核算m列）
cudaEliminate<<>>(…);
parallel omp
{
cpuEliminate(…);
}
再主机和设备同步(不好意思函数名字我忘记怎么拼写了，抱歉:$)
假如你矩阵规模不大，不要开omp，开了反而慢，CPU直接串行算
这个方法可以应用在很多高等代数的优化上~但这个方法不是最优的，还一种look ahead的CPU/GPU异步模式，可以去看看volkov的论文，比较麻烦
2，3. 如ice哥哥所述，但是规约的时候你可以一个线程比较多个元素的最值，比如比较4个，8个，16个…，总有一个可以达到最最优化，这样你开适量的线程就OK的，号称世界上最快的MAGMA中的矩阵向量乘法就是这么尝试出来的

system · 2013 年11 月 8 日 01:11

[
果真如此！我将它展开后再调整了一下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：你说的方法对我来说就是天书，我目前不准备用混合编程，不过还是要谢谢你了。

system · 2013 年11 月 8 日 03:05

合并访问与显存位宽无关，需要相邻线程访问一段连续的(最好是满足对齐要求的)显存。

system · 2013 年11 月 12 日 10:22

接着问问题：
（1）解线性方程组的时候当维数低于1500时工作正常，超过后出现unknown error错误，问一下有经验的版主，一般情况下是哪里出了问题？
（2）维数为1000的双精度，与CPU的答案大概只有前4、5位数字是相同的，单精度的就更低了，正常么？

system · 2013 年11 月 12 日 10:30

LZ您好：

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

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

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

大致如此，祝您编码顺利~

system · 2013 年11 月 29 日 03:25

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