CUDA矩阵求逆的问题

weizhuoqun 2009-09-22 09:10:30
加精
我想用cuda实现矩阵的求逆，用的方法是高斯消元法，依次令每一行作为主行去与其他行进行消元操作，这个过程是并行的。当我求100*100的矩阵逆的时候，kernal的执行时间很短，远远大于串行的CPU程序，但是换做是500*500的时候，GPU的时间却是CPU串行的3倍！用了n个线程，请教这是怎么回事！是我的程序有问题还是其他原因，下面我将我的核函数部分贴出来，用的是高斯消元法


__global__ void matrixinver(float *A,float *B,int i)

{

	const int tid=threadIdx.x;

	float temp=0.0;

 	if (tid!=i)

 	{

 		temp=A[tid*m+i]/A[i*m+i];

 		for (int k=0;k<m;k++)

 		{

 			A[tid*m+k]=A[tid*m+k]-A[i*m+k]*temp;

 		    B[tid*m+k]=B[tid*m+k]-B[i*m+k]*temp;

 		}

 	}

 	else{}

}



	cudaMalloc((void**) &tt_gpu,sizeof(float)*m*m);

	cudaMalloc((void**) &B_gpu,sizeof(float)*m*m);

	

	cudaMemcpy2D(tt_gpu,sizeof(float)*m,tt,sizeof(float)*m,sizeof(float)*m,m,cudaMemcpyHostToDevice);

	cudaMemcpy2D(B_gpu,sizeof(float)*m,B,sizeof(float)*m,sizeof(float)*m,m,cudaMemcpyHostToDevice);



	now=clock();

	for (int i=0;i<m;i++)

	{

		matrixinver<<<1,m,0>>>(tt_gpu,B_gpu,i);

	}

	cudaMemcpy2D(tt,sizeof(float)*m,tt_gpu,sizeof(float)*m,sizeof(float)*m,m,cudaMemcpyDeviceToHost);

	cudaMemcpy2D(B,sizeof(float)*m,B_gpu,sizeof(float)*m,sizeof(float)*m,m,cudaMemcpyDeviceToHost);

	then=clock();
...全文
4854 47 打赏 收藏 转发到动态 举报
写回复
用AI写文章
47 条回复
切换为时间正序
请发表友善的回复…
发表回复
jolynne1202 2011-09-19
  • 打赏
  • 举报
 回复
诚心学习您的LU分解方法和cholesky，能将你的代码发给我吗？先谢谢你了！邮箱：mumubenmao110@163.com
[Quote=引用 39 楼 cyrosly 的回复:]
标准的高斯约旦矩阵求逆和LU分解过程都是一个高斯消元过程，如果并行的话需要按行或按列（分解的话每步一行一列）求出然后更新剩余的值。所以N阶矩阵需要N个遍。除非找到可以将遍数降到log(N)复杂度的算法（事实上在GPU上实现这些已经变成线性复杂度了）.另外LL分解的程序不是已经发你了吗，虽然和LU分解不一样，但是还是有很多类似的地方，如果需要LU分解再发你，逆的话没做，有时间再做
[/Quote]
cicelyxq1987 2010-05-11
  • 打赏
  • 举报
 回复
[Quote=引用 23 楼 zwang4000 的回复:]
if ((xIndex < width ) && (yIndex < height))
{
unsigned int index_in = yIndex * width + xIndex;
block[threadIdx.y][threadIdx.x] = idata[index_in];
}
[/Quote]
这段是进行计算吗？我没理解，讲解下好吗？
yangxiaoxiaogang 2009-12-01
  • 打赏
  • 举报
 回复
[Quote=引用 1 楼 cyrosly 的回复:]
方法不对，可以参考下我做的cholesky分解（http://www.nvidia.com/object/cuda_home.html#state=detailsOpen;aid=d86e78ef-8ada-44e6-ad42-fc6386e55cc0）或者LU分解（需要邮箱），要尽量使用shared memory优化
[/Quote]

我的邮箱是 yangxiaoxiaogang@163.com ，我想学学你的LU分解那个cuda程序，能发给我一下吗？诚心谢过。
wjg_xaut 2009-11-01
  • 打赏
  • 举报
 回复
学习下
ff8ff7 2009-10-31
  • 打赏
  • 举报
 回复
学习了，
前段时间写了个SOR法的电场模拟的CUDA的计算程序，速度很慢，学到点东西
无心人_过过小日子 2009-10-26
  • 打赏
  • 举报
 回复
[Quote=引用 41 楼 weizhuoqun 的回复:]
引用 39 楼 cyrosly 的回复:
还有就是，kernal中对内存的读取很多，这是影响速度的一个主要原因，我用的是一个内核，就是直接进行普通的消元，不知道有什么解决方案。用了分块的做法，但是也遇到了不会找主行的问题
[/Quote]

分块是必须的吧.注意计算上的同步.
weizhuoqun 2009-10-26
  • 打赏
  • 举报
 回复
[Quote=引用 39 楼 cyrosly 的回复:]
标准的高斯约旦矩阵求逆和LU分解过程都是一个高斯消元过程，如果并行的话需要按行或按列（分解的话每步一行一列）求出然后更新剩余的值。所以N阶矩阵需要N个遍。除非找到可以将遍数降到log(N)复杂度的算法（事实上在GPU上实现这些已经变成线性复杂度了）.另外LL分解的程序不是已经发你了吗，虽然和LU分解不一样，但是还是有很多类似的地方，如果需要LU分解再发你，逆的话没做，有时间再做
[/Quote]
还有就是，kernal中对内存的读取很多，这是影响速度的一个主要原因，我用的是一个内核，就是直接进行普通的消元，不知道有什么解决方案。用了分块的做法，但是也遇到了不会找主行的问题
weizhuoqun 2009-10-25
  • 打赏
  • 举报
 回复
[Quote=引用 39 楼 cyrosly 的回复:]
如果需要LU分解再发你，逆的话没做，有时间再做
[/Quote]
恩，正在看您给我的cholesky分解的代码，现在也是在找关于高斯消元的突破点，觉得ＬＵ分解的高斯消元过程，与我现在正在想解决的问题比较接近，另外因为老师给的任务比较紧，想尽快解决矩阵求逆的问题，但是基础不是很好，所以才急切的想向您请教ＬＵ分解的例子，还是真的很需要您的例子的，谢谢你！
Cyrosly 2009-10-24
  • 打赏
  • 举报
 回复
标准的高斯约旦矩阵求逆和LU分解过程都是一个高斯消元过程，如果并行的话需要按行或按列（分解的话每步一行一列）求出然后更新剩余的值。所以N阶矩阵需要N个遍。除非找到可以将遍数降到log(N)复杂度的算法（事实上在GPU上实现这些已经变成线性复杂度了）.另外LL分解的程序不是已经发你了吗，虽然和LU分解不一样，但是还是有很多类似的地方，如果需要LU分解再发你，逆的话没做，有时间再做
weizhuoqun 2009-10-24
  • 打赏
  • 举报
 回复
[Quote=引用 35 楼 cyrosly 的回复:]
这样N阶的矩阵需要N个循环遍，对于数千阶数万阶的矩阵虽然需要的循环数看似比较多，但是这期间不存在host memory和device memory之间的存储器通信，所以影响微乎其微。事实一定尺寸的情况下可以达到2个数量级的加速比（CPU考虑保守的SSE指令优化（亦即假设SSE ISA达到理想的4倍速）。情况和列选主元LU分解差不多。
[/Quote]
可不可以问一下，能不能用具体的例子解释一下这种方法，小弟初涉CUDA还不是很明白，谢谢了！
weizhuoqun 2009-10-24
  • 打赏
  • 举报
 回复
[Quote=引用 36 楼 l7331014 的回复:]
同意.
[/Quote]
不是很理解楼上的问题，请用代码解释一下可以吗？
Cyrosly 2009-10-23
  • 打赏
  • 举报
 回复
这样N阶的矩阵需要N个循环遍，对于数千阶数万阶的矩阵虽然需要的循环数看似比较多，但是这期间不存在host memory和device memory之间的存储器通信，所以影响微乎其微。事实一定尺寸的情况下可以达到2个数量级的加速比（CPU考虑保守的SSE指令优化（亦即假设SSE ISA达到理想的4倍速）。情况和列选主元LU分解差不多。
Cyrosly 2009-10-23
  • 打赏
  • 举报
 回复
[Quote=引用 33 楼 leightonlz 的回复:]
引用 29 楼 zwang4000 的回复:
兄弟，上面不是有Kernel部分的code吗，你照着改。（不懂你为什么要在Kernel中用循环）


你好，谢谢你的方法，已经取得了很大的进步，还有一个知识点，就是因为我的是用来做矩阵求逆的，需要对矩阵的每一行作为主行依次去与其他行进行消元，我还需要在计算的时候把这一行单独取出来，求解决方案，谢谢！
[/Quote]

第一个内核选择主元，并记录下是否进行过交换（可通过行索引与自然顺序是否相同判断）
第二个内核如果发生过交换则交换右端向量
第三个内核计算一列
第四个内核传播，用上个内核的结果对其它incompleted的待求列列进行更新
无心人_过过小日子 2009-10-23
  • 打赏
  • 举报
 回复
[Quote=引用 35 楼 cyrosly 的回复:]
这样N阶的矩阵需要N个循环遍，对于数千阶数万阶的矩阵虽然需要的循环数看似比较多，但是这期间不存在host memory和device memory之间的存储器通信，所以影响微乎其微。事实一定尺寸的情况下可以达到2个数量级的加速比（CPU考虑保守的SSE指令优化（亦即假设SSE ISA达到理想的4倍速）。情况和列选主元LU分解差不多。
[/Quote]

同意.
weizhuoqun 2009-10-20
  • 打赏
  • 举报
 回复
[Quote=引用 29 楼 zwang4000 的回复:]
兄弟，上面不是有Kernel部分的code吗，你照着改。（不懂你为什么要在Kernel中用循环）
[/Quote]
我的程序是一个线程控制矩阵的一行进行消元，循环的意思是这一行的每一个元素，刚知道要在kernal中减少控制流程，最好用分块的解法解决大矩阵，正在学习你的code，我分的是16X16的小矩阵，我在试一试。谢谢你，朋友
Leightonlz 2009-10-20
  • 打赏
  • 举报
 回复
[Quote=引用 29 楼 zwang4000 的回复:]
兄弟，上面不是有Kernel部分的code吗，你照着改。（不懂你为什么要在Kernel中用循环）
[/Quote]

你好，谢谢你的方法，已经取得了很大的进步，还有一个知识点，就是因为我的是用来做矩阵求逆的，需要对矩阵的每一行作为主行依次去与其他行进行消元，我还需要在计算的时候把这一行单独取出来，求解决方案，谢谢！
weizhuoqun 2009-10-20
  • 打赏
  • 举报
 回复
楼上的问题解决了，每个block的thread过多，超过了512
weizhuoqun 2009-10-20
  • 打赏
  • 举报
 回复
[Quote=引用 29 楼 zwang4000 的回复:]
兄弟，上面不是有Kernel部分的code吗，你照着改。（不懂你为什么要在Kernel中用循环）
[/Quote]
我自己写了一下，不知道错在哪，就是简单做了个测试（矩阵每个元素的数值变为100），然后输出看结果，但是没有变化，下面是我的代码，请求修改指正，谢谢


//只有调用kernal，和kernal部分

#define N 1000

#define THREAD_NUM 100

#define BLOCK_NUM 10

#define BLOCK_SIZE 32

#define BLOCK_DIM 32

clock_t matinvCUDA(float* a, int lda, float* b, int ldb)

{

 clock_t start,finish;

 float *ac, *bc;

 start=clock();

 size_t pitch_a,pitch_b;

 cudaMallocPitch((void**) &ac, &pitch_a, sizeof(float) * N, N);

 cudaMallocPitch((void**) &bc, &pitch_b, sizeof(float) * N, N);

 cudaMemcpy2D(ac, pitch_a, a, sizeof(float) * lda,sizeof(float) * N, N, cudaMemcpyHostToDevice);

 cudaMemcpy2D(bc, pitch_b, b, sizeof(float) * ldb,sizeof(float) * N, N, cudaMemcpyHostToDevice);

 

 int bx = (N + BLOCK_SIZE - 1) / BLOCK_SIZE;

 dim3 blocks(bx, bx);

 dim3 threads(BLOCK_SIZE, BLOCK_SIZE);

 //for (int i=0;i<N;i++)

 //{

  matrixinv<<<blocks,threads>>>(ac, pitch_a/sizeof(float), bc, pitch_b/sizeof(float), N, 0);

 //}

 cudaMemcpy2D(a, sizeof(float) * N, ac, pitch_a,sizeof(float) * N, N, cudaMemcpyDeviceToHost);

 cudaMemcpy2D(b, sizeof(float) * N, bc, pitch_b,sizeof(float) * N, N, cudaMemcpyDeviceToHost);

 cudaFree(ac);

 cudaFree(bc);

 finish=clock();

 return finish-start;

}

__global__ static void matrixinv(float* a, size_t lda, float* b, size_t ldb, int n, int k)

{

 __shared__ float matA[BLOCK_SIZE][BLOCK_SIZE];

 __shared__ float matB[BLOCK_SIZE][BLOCK_SIZE];

 

 const int tidc = threadIdx.x;

 const int tidr = threadIdx.y;

 

 const int bidc = blockIdx.x * BLOCK_SIZE;

 const int bidr = blockIdx.y * BLOCK_SIZE;

 

 unsigned int xIndex = blockIdx.x * blockDim.x + threadIdx.x; 

 unsigned int yIndex = blockIdx.y * blockDim.y + threadIdx.y;

 

 if ((xIndex<N)&&(yIndex<N))

 {

  unsigned int index_ain = yIndex * lda + xIndex;

  unsigned int index_bin = yIndex * ldb + xIndex;

  matA[tidr][tidc] = a[index_ain];

  matB[tidr][tidc] = b[index_bin];

 }

 //测试

 matA[tidr][tidc] += 100;

 matB[tidr][tidc] += 100;

 __syncthreads();

 if ((xIndex<N)&&(yIndex<N))

 {

  unsigned int index_aout = yIndex * lda + xIndex;

  unsigned int index_bout = yIndex * ldb + xIndex;

  a[index_aout] = matA[tidr][tidc];

  b[index_bout] = matB[tidr][tidc];

 }

}
zwang4000 2009-10-19
  • 打赏
  • 举报
 回复
兄弟，上面不是有Kernel部分的code吗，你照着改。（不懂你为什么要在Kernel中用循环）
weizhuoqun 2009-10-18
  • 打赏
  • 举报
 回复
[Quote=引用 25 楼 zwang4000 的回复:]
基本思路
1、把矩阵分成32X32的小矩阵
2、每个小矩阵都在一个Block里执行
3、把计算结果先存入share memory里的block变量
__shared__ float block[BLOCK_DIM][BLOCK_DIM+1];
4、在把block变量里的值拷贝到odata里
[/Quote]
理解起来 感觉不是很难 但是还是没有实现，可不可以麻烦这位朋友帮我实现以下kernal部分啊，我程序现在改成了这样，请过目。怎么进行修改一下，看CUDA时间不长，实在不知道如何解决，谢谢了！


//调用核函数部分，外面的循环意思就是从第一行开始作为主行，进行消元，i是获得主行的行号，但是这么实现起来速度还是不好，之比串行快4s，原因还是内存读取太多

	for (int i=0;i<N;i++)

	{

		matrixinv<<<BLOCK_NUM,THREAD_NUM,sizeof(float)*N*2>>> 

               (ac,pitch_a/sizeof(float), bc, pitch_b/sizeof(float), i); 

		cudaThreadSynchronize();

	}

__global__ static void matrixinv(float* a, size_t lda, float* b, size_t ldb, int i)

{

	const int bid=blockIdx.x;

	const int tid=threadIdx.x;

	const int row=bid*THREAD_NUM+tid;

	extern __shared__ float shared[];

	for (int k=0;k<N;k++)

	{

		shared[k]=a[i*lda+k];

		shared[N+k]=b[i*ldb+k];

	}

	__syncthreads();

	__shared__ float ii;

	ii=a[i*lda+i];

	float temp=0.0;

	if (row!=i)

	{

		temp=a[row*lda+i]/ii;

		for (int k=0;k<N;k++)

		{

			a[row*lda+k]=a[row*lda+k]-shared[k]*temp;

			b[row*ldb+k]=b[row*ldb+k]-shared[N+k]*temp;

		}

	}

	else{}

}
加载更多回复（23）
CUDA安装配置和求逆举例
本文主要介绍了基于CUDA的CULA的安装配置，并附带了能够调试通过并有结果的求逆程序代码
cuda_code.zip_CUDA求逆_GPU_cuda 矩阵_矩阵求逆cuda_逆矩阵
利用GPU的多线程，借助cuda实现给定矩阵的求逆代码
FFmpeg系列之35：FFmpeg+CUDA硬件加速原理与案例
FFmpeg+CUDA硬件加速原理与案例实战FFmpeg系列之35FFmpeg第2季编解码专题之5：FFmpeg+CUDA硬件加速原理与案例实战本课程主要讲解的知识点包括：GPU高性能编程CUDA入门、CUDA编程模型的原理解析、CUDA编程小白案例入门实战、CUDA官方经典案例实战与解析、FFmpeg+CUVID硬件加速原理与命令行、Qt+FFmpeg实现cuda硬解码加速代码实战等。---------------------------------------------------------------目前FFmpeg第二季【编解码专题】，主要包括以下课程：FFmpeg4.3系列之31：图像格式与压缩FFmpeg4.3系列之32：MPEG-1/2视频编解码FFmpeg4.3系列之33：H.264/H.265视频编码FFmpeg4.3系列之34：FFmpeg+DXVA2+D3D硬件加速FFmpeg4.3系列之35：FFmpeg+CUDA硬件加速
矩阵求逆测试
程序中比较了LU分解和快算求逆算法的耗时的比较，对于需要这方面工作的同学有学习的必要，最大的优点是我写的一个矩阵类，非常好用。
基于GPU的矩阵求逆性能测试和分析 (2010年)
在CPU串行运算模式下实现大规模矩阵求逆是一个非常耗时的过程。为了解决这一问题,基于NVIDIA公司专为GPU(图形处理器)提供的CUDA(计算统一设备架构),从新的编程角度出发,利用GPU多线程并行处理技术,将矩阵求逆过程中大量的数据实现并行运算,从而获得了较大的加速比。同时,根据程序的执行结果,分析了GPU的单精度与双精度的浮点运算能力及其优、劣势。最后,通过分析数据传输时间对 GPU性能的影响,总结出适合GPU的算法特征。
CUDA高性能计算讨论

353

社区成员

615

社区内容

发帖
与我相关
我的任务
社区描述
CUDA高性能计算讨论
社区管理员
  • CUDA高性能计算讨论社区
加入社区
  • 近7日
  • 近30日
  • 至今

加载中

查看更多榜单
社区公告
暂无公告

试试用AI创作助手写篇文章吧

+ 用AI写文章