c++ - 高效布局和减少虚拟2D数据(抽象)

Question

我使用C++和CUDA/C，想为特定问题编写代码，但遇到了一个非常棘手的简化问题。

我在并行编程方面的经验不容忽视，但相当有限，我无法完全预见到此问题的特殊性。
我怀疑是否有一种方便甚至“轻松”的方式来处理我面临的问题，但也许我错了。
如果有任何涵盖此问题或类似问题的资源(例如文章，书籍，网络链接等)或关键字，请告知我。

我试图将整个案例尽可能地概括，并使其抽象化，而不是发布过多的代码。

布局 ...

我有一个由N个初始元素和N个结果元素组成的系统。 (例如，我将使用N = 8，但N可以是大于3的任何整数值。)

static size_t const N = 8;
double init_values[N], result[N];

我需要计算init值的几乎每一个(我怕不是全部)唯一的排列，而不会自我干扰。

这意味着计算 f(init_values[0],init_values[1])， f(init_values[0],init_values[2])，...， f(init_values[0],init_values[N-1])， f(init_values[1],init_values[2])，...， f(init_values[1],init_values[N-1])，...等等。

实际上，这是一个虚拟的三角形矩阵，其形状如下图所示。

 P     0    1    2    3    4    5    6    7
   |---------------------------------------
  0|   x 
   |
  1|   0    x
   |
  2|   1    2    x 
   |
  3|   3    4    5    x
   |
  4|   6    7    8    9    x
   |
  5|  10   11   12   13   14    x
   |
  6|  15   16   17   18   19   20    x
   |
  7|  21   22   23   24   25   26   27    x

每个元素都是 init_values中相应列和行元素的函数。

P[i] (= P[row(i)][col(i]) = f(init_values[col(i)], init_values[row(i)])

即

P[11] (= P[5][1]) = f(init_values[1], init_values[5])

使用示例 (N*N-N)/2 = 28，有 P[1][5]==P[5][1]可能是唯一的组合(注意: N = 8，所以我们只有一个较低(或较高)的三角矩阵)。

基本问题

从P计算结果数组，将其作为行元素的总和减去各个列元素的总和。
例如，位置3的结果将计算为第3行的总和减去第3列的总和。

result[3] = (P[3]+P[4]+P[5]) - (P[9]+P[13]+P[18]+P[24])
result[3] = sum_elements_row(3) - sum_elements_column(3)

我试图在N = 4的图片中进行说明。

因此，以下是正确的:

N-1操作(潜在的并发写入)将在每个result[i]上执行

result[i]将通过减法和N-(i+1)加法得到i写入

从每个P[i][j]中传出的内容将对r[j]进行减法运算，并对r[i]进行加法运算

这是主要问题出现的地方:

使用一个线程来计算每个P并直接更新结果将导致多个内核尝试写入相同的结果位置(每个N-1个线程)。

另一方面，存储整个矩阵P以便进行后续的缩减步骤在内存消耗方面非常昂贵，因此对于大型系统而言是不可能的。

对于每个线程块都具有唯一的，共享的结果 vector 的想法也是不可能的。
(N个50k构成25亿个P元素，因此[假定每个块最多有1024个线程]如果每个块都有其自己的结果数组(具有50k个double元素)，则最小的240万个块消耗900GiB的内存。)

我认为我可以处理减少操作以获得更静态的行为，但是就潜在的并发内存写访问而言，此问题相当动态。
(或者是否可以通过某种“基本”减少量来处理？)

增加一些并发症...

遗憾的是，取决于(任意用户)输入，该输入与初始值无关，需要跳过P的某些元素。
假设我们需要跳过排列P [6]，P [14]和P [18]。因此，我们还有24种组合需要计算。

如何告诉内核哪些值需要跳过？
我想出了三种方法，如果N非常大(如几万个元素)，则每种方法都有明显的缺点。

1.存储所有组合...

...及其各自的行和列索引 struct combo { size_t row,col; };，需要在 vector<combo>中计算并在此 vector 上进行操作。 (由当前实现使用)

std::vector<combo> elements;
// somehow fill
size_t const M = elements.size();
for (size_t i=0; i<M; ++i)
{
    // do the necessary computations using elements[i].row and elements[i].col  
}

由于仅“几个”(甚至可能是一万个元素，但这与数十亿个元素相比并没有太大关系)，所以该解决方案消耗了大量内存，但是它避免了

索引计算

查找已删除元素

对于P的每个元素，这是第二种方法的缺点。

2.对P的所有元素进行操作并查找已删除的元素

如果我想对P的每个元素进行操作并避免嵌套循环(在cuda中无法很好地重现)，则需要执行以下操作:

size_t M = (N*N-N)/2;
for (size_t i=0; i<M; ++i)
{
   // calculate row indices from `i`
   double tmp = sqrt(8.0*double(i+1))/2.0 + 0.5;
   double row_d = floor(tmp);
   size_t current_row = size_t(row_d);
   size_t current_col = size_t(floor(row_d*(ict-row_d)-0.5));
   // check whether the current combo of row and col is not to be removed
   if (!removes[current_row].exists(current_col))
   {
     // do the necessary computations using current_row and current_col
   }
}

与第一个示例中的 removes vector 相比， vector elements非常小，但是用于获得 current_row， current_col和if分支的额外计算效率非常低。
(请记住，我们仍在谈论数十亿次评估。)

3.操作P的所有元素，然后删除元素

我的另一个想法是独立计算所有有效和无效组合。
但是不幸的是，由于求和错误，以下语句是正确的:

calc_non_skipped() != calc_all() - calc_skipped()

是否有一种方便，已知的高性能方法来从初始值中获得所需结果？

我知道这个问题相当复杂，相关性可能有限。不过，我希望一些启发性的答案能够帮助我解决问题。

当前执行

当前，这是通过OpenMP实现为CPU代码。
我首先建立一个上述 combo的 vector ，该 vector 存储需要计算的每个P并将其传递给并行的for循环。
每个线程都具有私有(private)结果 vector ，并且并行区域末尾的关键部分用于适当求和。

Answer 1

首先，我有些困惑，为什么(N**2 - N)/2在N = 7时会产生27 ...但是对于索引0-7，N = 8来说，P中有28个元素。所以在今晚不要尝试回答这样的问题那天。 :-)

但是有一个潜在的解决方案:您是否需要将阵列P保留用于其他目的？如果没有，我认为您可以使用两个中间数组来获得所需的结果，每个中间数组的长度为N:一个用于行的总和，一个用于列的总和。

这是我想尝试做的一个简单示例(子例程direct_approach())，以及如何使用中间数组(子例程refined_approach())实现相同的结果:

#include <cstdlib>
#include <cstdio>

const int N = 7;
const float input_values[N] = { 3.0F, 5.0F, 7.0F, 11.0F, 13.0F, 17.0F, 23.0F };
float P[N][N];      // Yes, I'm wasting half the array.  This way I don't have to fuss with mapping the indices.
float result1[N] = { 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F };
float result2[N] = { 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F };

float f(float arg1, float arg2)
{
    // Arbitrary computation
    return (arg1 * arg2);
}

float compute_result(int index)
{
    float row_sum = 0.0F;
    float col_sum = 0.0F;
    int row;
    int col;

    // Compute the row sum
    for (col = (index + 1); col < N; col++)
    {
        row_sum += P[index][col];
    }

    // Compute the column sum
    for (row = 0; row < index; row++)
    {
        col_sum += P[row][index];
    }

    return (row_sum - col_sum);
}

void direct_approach()
{
    int row;
    int col;

    for (row = 0; row < N; row++)
    {
        for (col = (row + 1); col < N; col++)
        {
            P[row][col] = f(input_values[row], input_values[col]);
        }
    }

    int index;

    for (index = 0; index < N; index++)
    {
        result1[index] = compute_result(index);
    }
}

void refined_approach()
{
    float row_sums[N];
    float col_sums[N];
    int index;

    // Initialize intermediate arrays
    for (index = 0; index < N; index++)
    {
        row_sums[index] = 0.0F;
        col_sums[index] = 0.0F;
    }

    // Compute the row and column sums
    // This can be parallelized by computing row and column sums
    //  independently, instead of in nested loops.
    int row;
    int col;

    for (row = 0; row < N; row++)
    {
        for (col = (row + 1); col < N; col++)
        {
            float computed = f(input_values[row], input_values[col]);
            row_sums[row] += computed;
            col_sums[col] += computed;
        }
    }

    // Compute the result
    for (index = 0; index < N; index++)
    {
        result2[index] = row_sums[index] - col_sums[index];
    }
}

void print_result(int n, float * result)
{
    int index;

    for (index = 0; index < n; index++)
    {
        printf("  [%d]=%f\n", index, result[index]);
    }
}

int main(int argc, char * * argv)
{
    printf("Data reduction test\n");

    direct_approach();

    printf("Result 1:\n");
    print_result(N, result1);

    refined_approach();

    printf("Result 2:\n");
    print_result(N, result2);

    return (0);
}

并行化计算并不是那么容易，因为每个中间值都是大多数输入的函数。您可以单独计算总和，但这意味着多次执行f(...)。对于N的非常大的值，我可以想到的最佳建议是使用更多的中间数组，计算结果的子集，然后将部分数组求和以得出最终的总和。当我不那么累的时候，我不得不考虑那个。

要解决跳过问题:如果仅是“不要使用输入值x，y和z”，则可以将x，y和z存储在do_not_use数组中，并在循环计算时检查这些值总和。如果要跳过的值是行和列的某种功能，则可以将它们存储为对并检查对。

希望这能为您提供解决方案的想法!

更新，现在我已经清醒了:处理“跳过”在很大程度上取决于需要跳过哪些数据。第一种情况的另一种可能性-“不要使用输入值x，y和z”-对于大型数据集，更快的解决方案是添加一个间接级别:创建另一个数组，这个数组是整数索引，并只存储良好输入的索引。在第二个实例中，如果输入2和5中包含无效数据，则有效数组为:

int valid_indices[] = { 0, 1, 3, 4, 6 };

遍历数组 valid_indices，并使用这些索引从输入数组中检索数据以计算结果。另一方面，如果要跳过的值取决于P数组的两个索引，那么我看不到如何避免某种查找。

返回并行化-无论如何，您将要处理(N ** 2-N)/2个计算
的f()。一种可能性是只接受对总和的争执
数组，如果计算f()花费的时间长于
这两个增加。当您到达大量并行路径时，竞争将
仍然是一个问题，但是应该有一个“最佳点”来平衡并行数
计算f()所需时间的路径。

如果仍然存在争用，则可以采用几种方式对问题进行分区。一种方法是
一次计算一行或一列:对于一次一行，每个列的总和可以为
独立计算，并且可以为每行总和保留运行总计。

另一种方法是将数据空间划分为
子集，其中每个子集都有自己的行和列求和数组。每个块之后
计算后，可以对独立数组求和以生成您需要的值
计算结果。