高效布局和减少虚拟2d数据(摘要)

Question

我使用C++和CUDA/C并且想要为特定问题编写代码,我遇到了一个相当棘手的减少问题.

我在并行编程方面的经验不容忽视,但非常有限,我不能完全预见到这个问题的特殊性.我怀疑有一种方便甚至"简单"的方式来处理我所面临的问题,但也许我错了.如果有任何资源(即文章,书籍,网络链接......)或关键词覆盖此类或类似问题,请告诉我.

我试图尽可能地概括整个案例并保持抽象而不是发布太多代码.

布局 ...

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

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

我需要计算初始值的几乎所有(不是所有我害怕的)唯一置换而没有自干扰.

这意味着计算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在存储器消耗方面是非常昂贵的,因此对于非常大的系统是不可能的.

为每个线程块提供unqiue,共享结果向量的想法也是不可能的.(50k的N个产生25亿个P元素,因此[假设每个块最多1024个线程],如果每个块都有自己的具有50k双元素的结果数组,则最少240万个块消耗超过900GiB的内存.)

我认为我可以处理更多静态行为的减少,但就潜在的并发内存写访问而言,这个问题相当动态.(或者是否可以通过一些"基本"类型的减少来处理它？)

增加一些并发症......

不幸的是,根据(任意用户)输入,它与初始值无关,需要跳过P的一些元素.假设我们需要跳过排列P [6],P [14]和P [18].因此,我们剩下24个组合,需要进行计算.

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

1.存储所有组合......

...具有各自的行和列索引struct combo { size_t row,col; };,需要在a中计算vector<combo>并对此向量进行操作.(由当前实现使用)

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是在对比的非常小的elements在第一示例矢量但附加计算,以获得current_row,current_col并且如果分支是非常低效的.(请记住,我们仍然在讨论数十亿的评估.)

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

我的另一个想法是独立计算所有有效和无效的组合.但不幸的是,由于总和错误,以下陈述是正确的:

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

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

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

---

目前的实施

目前,这是作为带有OpenMP的CPU代码实现的.我首先建立一个上面提到的combos 的向量,存储每个需要计算的P并将其传递给并行for循环.每个线程都有一个私有结果向量,并行区域末端的临界区用于正确的求和.

Answer 1

首先,我感到困惑的是为什么(N**2 - N)/2N = 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数组的两个索引,我看不出如何避免某种查找.

回到并行化 - 无论如何,你将处理f()的(N**2 - N)/ 2计算.一种可能性就是接受对和数组的争用,如果计算f()花费的时间比两次加法长得多,这就不是一个大问题.当你到达非常大量的并行路径时,争用将再次成为一个问题,但应该有一个"最佳点"平衡并行路径的数量与计算f()所需的时间.

如果争用仍然存在问题,您可以通过多种方式对问题进行分区.一种方法是一次计算一行或一列:对于一次一行,每列总和可以独立计算,并且可以为每个行总和保持运行总和.

另一种方法是将数据空间划分,从而将计算划分为子集,其中每个子集具有其自己的行和列和数组.在计算每个块之后,可以对独立数组求和以生成计算结果所需的值.