利用gsl库实现线性代数的并行化

Question

在我的c++脚本中，我有许多用于计算线性代数运算的for循环。我想知道什么是使循环平行的最佳方法？一个例子是下面的函数，它计算两个矩阵的克罗内克乘积。

void Kronecker(const gsl_matrix *K, const gsl_matrix *V, gsl_matrix *H) 
{
for (size_t i=0; i<K->size1; i++) {
for (size_t j=0; j<K->size2; j++) {
gsl_matrix_view H_sub=gsl_matrix_submatrix (H, i*V->size1, j*V->size2, V->size1, V->size2);
gsl_matrix_memcpy (&H_sub.matrix, V);
gsl_matrix_scale (&H_sub.matrix, gsl_matrix_get (K, i, j));
}
}
return;
}

当我有可以并行的for loops时，我如何提高代码的计算时间？

Answer 1

在不了解底层gsl调用中的内存布局、分配、系统调用和潜在副作用的情况下，通过OpenMP实现并行化是一种非常简单的方法。当然，这会引入依赖项并需要编译器支持，但它在像您这样的简单循环中特别有效。未经测试，可能需要更多来确保H正确写入，但类似于：

#pragma omp parallel for private(i, j)
for (size_t i=0; i<K->size1; i++) {
for (size_t j=0; j<K->size2; j++) {
gsl_matrix_view H_sub=gsl_matrix_submatrix (H, i*V->size1, j*V->size2, V->size1, V->size2);
gsl_matrix_memcpy (&H_sub.matrix, V);
gsl_matrix_scale (&H_sub.matrix, gsl_matrix_get (K, i, j));
}
}

请参阅https://curc.readthedocs.io/en/latest/programming/OpenMP-C.html了解更多详细信息。

如果你不想引入依赖项或有其他约束(例如，OpenMP在库代码中可能会有问题(，你总是可以自己做，在线程中设置内部for循环，在开始时启动N个线程，在结束时加入。当然，这是假设你有足够的工作，如果矩阵足够大的话，你可能会这样做。

Answer 2

不确定这是否会有多大帮助，但我有一个使用pthread.h库计算带有部分枢轴矩阵的高斯消去的老例子。

总之，亮点是：

• 创建线程数组pthread_t threads[N];
• 初始化线程运行到pthread_barrier_init(&barrier, NULL, numThreads);的停止屏障
• 在尝试多线程的函数中设置障碍，使其等待每个函数都具有继续执行所需的依赖项。在您的点数处添加pthread_barrier_wait(&barrier);
• 启动线程

for (i = 0; i < nthreads; i++)
{
pthread_create(&threads[i], NULL, functionWithThreading, (void *)i);
}

• 最后，等待所有线程完成并将它们连接起来

for (i = 0; i < nthreads; i++)
{
pthread_join(threads[i], NULL);
}

我知道这可能不是你想要的确切解决方案，但我希望这个例子能帮助