如何最有效地为 OpenCL 中的厄米特（对称）矩阵映射内核范围？答案

【问题标题】：How can I most efficiently map a kernel range for a hermitian (symmetric) matrix in OpenCL?如何最有效地为 OpenCL 中的厄米特（对称）矩阵映射内核范围？
【发布时间】：2014-07-31 18:06:27
【问题描述】：

我正在开发一个 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() 两次）？

【问题讨论】：

i==j 的对角线使用什么公式？是 conj（结果）吗？
在我的例子中，对于 i == j 它是一个特殊的缩放身份（即一个标量）。然而，由于恒等对角线在厄米矩阵中是纯实数，所以它是它自己的复共轭。

标签： c++ matrix opencl

【解决方案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 ;

所以你需要将一个1D内核加入N(N+1)/2个工作项，你可以自己决定工作组的大小（通常每个工作组64个项目是一个不错的选择）。

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

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

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

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

【讨论】：

谢谢，但似乎有一个错误。在计算 J 之后，需要将 I 添加到其中以获得 J 的“正确”值。
@stix。谢谢，已修复。

【解决方案2】：

如果您的矩阵真的很大，那么您可以将您的 NxN 矩阵切成 (N/k)x(N/k) 块，每个块的大小为 kxk。只要您只需要一半的数据，您就可以大致创建大小为local_group_size * (N/k)x(N/k)/2 的1D NDRange。

矩阵的每一块都由一个 LocalGroup 处理（LocalGroup 的大小由您选择）。这个想法是您在主机端创建一个数组，其中包含矩阵中每个工作组的位置。内核存根应如下所示：

void __kernel myKernel(
    __global int* coords,
    ....)
{
    int2 WorkGroupPositionInMatrix = vload2(get_group_id(0), coords);
    ...
    DoCalculation();
    ...
    WriteResultTwice();
    ...
    return;
}

您需要手动处理 - 处理将放置在矩阵对角线上的 thouse WorkGroups。如果矩阵大小很大，那么放置在对角线上的 LocalGroups 的开销可以忽略不计。

【讨论】：

【解决方案3】：

可以将直角三角形垂直切成两半，然后将较小的部分旋转以适应较大的部分，从而形成一个面积相等的矩形。因此，很容易将您的三角形全局工作区变成一个适合 OpenCL 的矩形。

在这里查看我的答案：OpenCL efficient way to group a lower triangular matrix

【讨论】：