c++ - 我怎样才能最有效地映射 OpenCL 中厄米特(对称)矩阵的内核范围？

Question

我正在开展一个 OpenCL 项目以生成非常大的厄尔米特(对称)矩阵，并且我正在尝试确定生成工作 ID 的最佳方式。

厄密矩阵沿对角线对称，因此 M(i,j) = M*(j,i)。

在暴力方式下，for 循环看起来像:

for(int i = 0; i < N; i++)
{
    for(int j = 0; j < N; j++)
    {
      complex<float> result = doSomeCalculation();
      M(i,j) = result;
    }
}

但是，利用厄尔米特特性，只需计算矩阵的上三角部分并将结果复制到下三角部分，循环的效率可以提高一倍:

for(int i = 0; i < N; i++)
{
    for(int j = i; j < N; j++)
    {
      complex<float> result = doSomeCalculation();
      M(i,j) = result;
      M(j,i) = conj(result);
    }
}

在这两个循环中，doSomeCalculation() 都是一项昂贵的操作，矩阵中的每个条目与其他每个条目完全不相关(即问题是愚蠢的并行)。

我的问题是:

我如何将带有 doSomeCalculation 的第二个循环实现为 OpenCL 内核，以便最有效地使用线程 ID(即，以便线程计算 M(i,j) 和 M(j,i) 而无需必须调用 doSomeCalculation() 两次)？

Answer 1

您需要使用线性索引，例如您可以用这种方式索引矩阵的每个元素:

0  1  2   ...  N-1
* N-2     ...  2N-2
        ....
*  * 2N-1 ...  N(N+1)/2 -1

也就是说，索引 K 由下式给出:

 k=iN-i*(i+1)/2+j

其中 N 是矩阵的大小，(i,j) 分别是行和列的从 0 开始的索引。

这种关系是可以颠倒的；查看this的答案问题，我在这里报告的完整性:

i = floor( ( 2*N+1 - sqrt( (2N+1)*(2N+1) - 8*k ) ) / 2 ) ;
j = k - N*i + i*(i+1)/2 ;

因此您需要将一个具有 N(N+1)/2 个工作项的 1D 内核入队，并且您可以自行决定工作组的大小(通常每个工作组 64 个项目是一个不错的选择)。

然后在 OpenCL 代码中，您可以使用以下方法检索索引 K:

int k = get_group_id(0)*64 + get_local_id(0);

然后使用你需要计算的矩阵元素索引上面的两个关系。

此外，请注意，您还可以通过将厄密矩阵表示为具有 N(N+1)/2 个元素的线性 vector 来节省空间。