使用cuda并联特征值求解器
Eigenvalue solver in parallel using CUDA
我一直在搜索和搜索网络,我似乎找不到正在寻找的答案。我有一个特殊的问题。
我正在编辑此问题,以简单地解决问题,并希望它更可读和可理解。
假设我有5000 20x20对称,密集的矩阵。我想在CUDA中创建一个内核,该内核将使每个线程负责计算每个对称矩阵的特征值。
如果可能的话,cuda内核的示例代码将很棒。
任何帮助/建议都将不胜感激!
谢谢,
Johnathan
我想在CUDA中创建一个内核,该内核将使每个线程负责计算每个对称矩阵的特征值。
对我来说,这是否是最快的方法是值得怀疑的,但这可能是很小的矩阵。即使在这种情况下,也可能会进行一些数据存储优化(跨线程的全局数据交织),但这会使事情复杂化。
如前所述,该请求可以映射到"令人尴尬的并行"算法中,其中每个线程都在完全独立的问题上工作。我们只需要找到合适的单线"供体代码"。经过快速的Google搜索,我遇到了这个。修改该代码以这种无关的方式运行非常简单。我们只需要借3个例程(jacobi_eigenvalue
,r8mat_diag_get_vector
和r8mat_identity
),并且可以在GPU上使用__host__ __device__
适当地装饰这些例程,同时进行没有其他更改。
所讨论的代码似乎是佛罗里达州立大学J Burkardt许可的GNU LGPL。因此,考虑到这一点,并且遵循传统的观点,我没有在此答案中包含任何大量该代码。但是您应该能够使用我提供的说明在实验中重建我的结果。
注意:我不确定使用此代码的法律后果是什么,该代码声称已获得GNU LGPL许可。如果您选择使用此代码或部分代码,则应确保遵守任何必要的要求。我在这里使用它的主要目的是演示单线程问题求解器的相对微不足道的"尴尬并行"扩展的概念。
通过转到此处并将3个指示功能复制到其余代码 - 骨骼中指示的位置中,重建我的完整代码应该是微不足道的。但这不会改变前面提到的任何通知/免责声明。
使用它。再次,从性能的角度来看,没有其他更改可能不是最好的想法,但是它会导致琐碎的努力,并且可以作为可能有用的起点。一些可能的优化可能是:
- 寻找数据交织策略,以便相邻的线程更有可能正在读取相邻的数据
- 从线程代码中消除
new
和delete
功能,并用固定的分配替换(这很容易做到) - 删除不必要的代码 - 例如,如果该数据不需要,则计算和分类特征向量的代码
无论如何,使用上述装饰的供体代码,我们只需要围绕它包裹一个琐碎的内核(je
),即可启动每个在单独的数据集(即矩阵)上操作的每个线程,并且每个线程都会产生自己的特征值集合(和特定代码库)。
我将其精心设计以使用3个线程和3个4x4矩阵用于测试目的,但是将其扩展到您想要的许多矩阵/线程应该很微不足道。
为了简短的演示,我已经分配了通常的错误检查,但是我建议您使用它或至少使用cuda-memcheck
运行代码。
我还构建了代码以向上调整设备堆的大小,以适应内核new
操作,具体取决于矩阵数(即螺纹)和矩阵尺寸。如果您在上面提到的第二个优化方面工作,则可以删除此。
t1177.cu:
#include <stdio.h>
#include <iostream>
const int num_mat = 3; // total number of matrices = total number of threads
const int N = 4; // square symmetric matrix dimension
const int nTPB = 256; // threads per block
// test symmetric matrices
double a1[N*N] = {
4.0, -30.0, 60.0, -35.0,
-30.0, 300.0, -675.0, 420.0,
60.0, -675.0, 1620.0, -1050.0,
-35.0, 420.0, -1050.0, 700.0 };
double a2[N*N] = {
4.0, 0.0, 0.0, 0.0,
0.0, 1.0, 0.0, 0.0,
0.0, 0.0, 3.0, 0.0,
0.0, 0.0, 0.0, 2.0 };
double a3[N*N] = {
-2.0, 1.0, 0.0, 0.0,
1.0, -2.0, 1.0, 0.0,
0.0, 1.0, -2.0, 1.0,
0.0, 0.0, 1.0, -2.0 };
/* ---------------------------------------------------------------- */
//
// the following functions come from here:
//
// https://people.sc.fsu.edu/~jburkardt/cpp_src/jacobi_eigenvalue/jacobi_eigenvalue.cpp
//
// attributed to j. burkardt, FSU
// they are unmodified except to add __host__ __device__ decorations
//
//****************************************************************************80
__host__ __device__
void r8mat_diag_get_vector ( int n, double a[], double v[] )
/* PASTE IN THE CODE HERE, FROM THE ABOVE LINK, FOR THIS FUNCTION */
//****************************************************************************80
__host__ __device__
void r8mat_identity ( int n, double a[] )
/* PASTE IN THE CODE HERE, FROM THE ABOVE LINK, FOR THIS FUNCTION */
//****************************************************************************80
__host__ __device__
void jacobi_eigenvalue ( int n, double a[], int it_max, double v[],
double d[], int &it_num, int &rot_num )
/* PASTE IN THE CODE HERE, FROM THE ABOVE LINK, FOR THIS FUNCTION */
// end of FSU code
/* ---------------------------------------------------------------- */
__global__ void je(int num_matr, int n, double *a, int it_max, double *v, double *d){
int idx = threadIdx.x+blockDim.x*blockIdx.x;
int it_num;
int rot_num;
if (idx < num_matr){
jacobi_eigenvalue(n, a+(idx*n*n), it_max, v+(idx*n*n), d+(idx*n), it_num, rot_num);
}
}
void initialize_matrix(int mat_id, int n, double *mat, double *v){
for (int i = 0; i < n*n; i++) *(v+(mat_id*n*n)+i) = mat[i];
}
void print_vec(int vec_id, int n, double *d){
std::cout << "matrix " << vec_id << " eigenvalues: " << std::endl;
for (int i = 0; i < n; i++) std::cout << i << ": " << *(d+(n*vec_id)+i) << std::endl;
std::cout << std::endl;
}
int main(){
// make sure device heap has enough space for in-kernel new allocations
const int heapsize = num_mat*N*sizeof(double)*2;
const int chunks = heapsize/(8192*1024) + 1;
cudaError_t cudaStatus = cudaDeviceSetLimit(cudaLimitMallocHeapSize, (8192*1024) * chunks);
if (cudaStatus != cudaSuccess) {
fprintf(stderr, "set device heap limit failed!");
}
const int max_iter = 1000;
double *h_a, *d_a, *h_v, *d_v, *h_d, *d_d;
h_a = (double *)malloc(num_mat*N*N*sizeof(double));
h_v = (double *)malloc(num_mat*N*N*sizeof(double));
h_d = (double *)malloc(num_mat* N*sizeof(double));
cudaMalloc(&d_a, num_mat*N*N*sizeof(double));
cudaMalloc(&d_v, num_mat*N*N*sizeof(double));
cudaMalloc(&d_d, num_mat* N*sizeof(double));
memset(h_a, 0, num_mat*N*N*sizeof(double));
memset(h_v, 0, num_mat*N*N*sizeof(double));
memset(h_d, 0, num_mat* N*sizeof(double));
initialize_matrix(0, N, a1, h_a);
initialize_matrix(1, N, a2, h_a);
initialize_matrix(2, N, a3, h_a);
cudaMemcpy(d_a, h_a, num_mat*N*N*sizeof(double), cudaMemcpyHostToDevice);
cudaMemcpy(d_v, h_v, num_mat*N*N*sizeof(double), cudaMemcpyHostToDevice);
cudaMemcpy(d_d, h_d, num_mat* N*sizeof(double), cudaMemcpyHostToDevice);
je<<<(num_mat+nTPB-1)/nTPB, nTPB>>>(num_mat, N, d_a, max_iter, d_v, d_d);
cudaMemcpy(h_d, d_d, num_mat*N*sizeof(double), cudaMemcpyDeviceToHost);
print_vec(0, N, h_d);
print_vec(1, N, h_d);
print_vec(2, N, h_d);
return 0;
}
编译和样本运行:
$ nvcc -o t1177 t1177.cu
$ cuda-memcheck ./t1177
========= CUDA-MEMCHECK
matrix 0 eigenvalues:
0: 0.166643
1: 1.47805
2: 37.1015
3: 2585.25
matrix 1 eigenvalues:
0: 1
1: 2
2: 3
3: 4
matrix 2 eigenvalues:
0: -3.61803
1: -2.61803
2: -1.38197
3: -0.381966
========= ERROR SUMMARY: 0 errors
$
我的输出对我来说似乎是合理的,主要与此处的输出匹配。
- 编译时未启用intel oneApi CUDA支持
- 在cuda线程之间共享大量常量数据
- 为什么即使使用-cudart-static进行编译,库用户仍然需要链接到cuda运行时
- Cuda C++:设备上的Malloc类,并用来自主机的数据填充它
- CUDA内核和数学函数的显式命名空间
- CUDA:统一内存和指针地址的更改
- 调试 CUDA MMU 故障
- 使用 CUDA 和纹理进行图像减法
- 将 2D 推力::d evice_vector 复矩阵传递给 CUDA 内核函数
- 编译 CUDA 与数学函数的叮当
- 为什么 CUDA 不会导致C++代码加速?
- 如何防止 CUDA-GDB 中的<优化输出>值
- 通过Python Distutils(用于Python C扩展)使用可重定位的设备代码编译CUDA代码
- CUDA三角函数中的数学保证
- CUDA 使用共享内存平铺 3D 卷积实现
- CUDA:cudaMallocManage处理退出吗?
- Opencv 加速与 CUDA 在 C++.
- Cuda:具有位集数组的 XOR 单位集
- 用于构建 cuda .so 文件(共享库)的生成文件
- 使用cuda并联特征值求解器