c++ - guide - cuda9 documentation



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

我使用C ++和CUDA / C,并想写一个特定的问题的代码,我遇到了一个非常棘手的减少问题。

我在并行编程方面的经验不可忽略,但相当有限,我不能完全预见这个问题的特殊性。 我怀疑是否有一种方便或甚至“简单”的方法来处理我所面对的问题,但也许我错了。 如果有任何资源(如文章,书籍,网站链接等)或关键词涉及这个或类似的问题,请让我知道。

我尽量推广整个案例,并保持抽象而不是张贴太多的代码。

布局 ...

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

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相应列和行元素的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*NN)/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的图片来说明它。

结果如下:

  • 将对每个result[i]执行N-1操作(潜在的并发写入)
  • result[i]将会有减法和i加法的N-(i+1)写入
  • 从每个P[i][j]将会有r[j]的相减,并且r[i]

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

  • 使用一个线程计算每个P并直接更新结果将导致多个内核尝试写入相同的结果位置(每个线程有N-1个线程)。
  • 另一方面,存储整个矩阵P用于随后的减少步骤在内存消耗方面非常昂贵,因此对于非常大的系统是不可能的。

为每个线程块分配一个共享结果向量的想法也是不可能的。 (如果每个块有50k个双元素的结果数组,则50k的N中有25亿P元素,因此[假设每块最多1024个线程]最少数量为240万块,消耗900GiB的内存。

我认为我可以处理更多的静态行为,但是这个问题在潜在的并发内存写入访问方面是相当动态的。 (还是有可能通过一些“基本”减少来处理?)

添加一些复杂的...

不幸的是,取决于(独立用户)输入,这是独立于初始值,P的一些元素需要被跳过。 假设我们需要跳过排列P [6],P [14]和P [18]。 因此,我们有剩下的24个组合,需要计算。

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

1.存储所有组合...

...以及它们各自的行和列索引struct combo { size_t row,col; }; struct combo { size_t row,col; }; ,需要在一个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
   }
}

与第一个示例中的elements矢量相比,矢量removes非常小,但是获取current_rowcurrent_colcurrent_col分支的额外计算效率非常低。 (请记住,我们仍在谈论数十亿的评估。)

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

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

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

是否有一种方便的,已知的,高性能的方法从初始值中获得预期的结果?

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

目前的实施

目前这是用OpenMP的CPU代码实现的。 我首先建立了一个上述combo的向量,存储每个需要计算的P,并将其传递给并行for循环。 每个线程都提供了一个私有结果向量,并行区域末尾的关键部分用于进行适当的求和。

https://src-bin.com


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数组的两个索引,我不知道如何避免某种查找。

回到并行化 - 无论如何,你将会处理f()的(N ** 2 - N)/ 2个计算。 一种可能性就是只接受和数组的争用,如果计算f()比两个增加要花费更多的时间,这将不是一个大问题。 当你得到大量的并行路径时,争用将再次成为一个问题,但是应该有一个平衡路径数量和计算f()所需时间的“最佳位置”。

如果争用仍然是一个问题,您可以通过几种方法来分割问题。 一种方法是一次计算一行或一列:对于一次一行,每列的总和可以独立计算,并且可以为每一行的总和保留一个运行的总和。

另一种方法是将数据空间划分为子集,其中每个子集都有自己的行和列和数组。 在计算完每个块之后,可以将独立数组相加以产生计算结果所需的值。


Answer #2

这可能会是那些天真无用的答案之一,但也可能有帮助。 随意告诉我,我是完全错误的,我误解了整个事情。

所以...我们走吧!

基本问题

在我看来,你可以定义你的结果函数有点不同,它会至少解除你的中间值的一些争论。 假设你的P矩阵是下三角形的。 如果你(实际上)用较低值的负数(和全零的主对角线)填充上三角形,则可以将结果的每个元素重新定义为单行的和:(这里显示为N = 4 ,其中-i表示单元格中标记为i的值的负值)

 P     0    1    2    3 
   |--------------------
  0|   x   -0   -1   -3
   |
  1|   0    x   -2   -4
   |
  2|   1    2    x   -5
   |
  3|   3    4    5    x

如果启动独立线程(执行相同内核)来计算此矩阵的每一行的总和,则每个线程将写入一个结果元素。 看来你的问题的规模足够大,足以饱和你的硬件线程,并保持忙碌。

需要注意的是,你将会计算两次f(x, y) 。 我不知道这是多么昂贵,或者以前有多少内存争夺费用,所以我不能判断这是否是一个值得做的权衡。 但是,除非f真的很贵,我想这可能是。

跳过价值

你提到你可能有数以万计的P矩阵元素,你需要忽略你的计算(有效地跳过它们)。

为了处理上面提出的方案,我相信你应该将跳过的元素存储为(row, col)对,并且还必须添加每个坐标对的转置(所以你将有两倍的跳过数因此你的P[6], P[14] and P[18]跳过列表变成P(4,0), P(5,4), P(6,3) ,然后变成P(4,0), P(5,4), P(6,3), P(0,4), P(4,5), P(3,6)

然后你根据行,然后列排序这个列表。 这使得我们的列表是P(0,4), P(3,6), P(4,0), P(4,5), P(5,4), P(6,3)

如果你的虚拟P矩阵的每一行都是由一个线程(或者你的内核的一个单独的实例或者其他的)来处理的话,你可以把它传递给它所需要的值。 就个人而言,我将所有这些存储在一个大的一维数组中,只是传递每个线程需要查看的第一个和最后一个索引(我也不会将行索引存储在我传入的最终数组中,因为它可以)在上面的例子中,对于N = 8,传递给每个线程的开始和结束对将是:(注意,结尾是一个超过需要处理的最终值,只是像STL,所以一个空列表由begin == end表示)

Thread 0: 0..1
Thread 1: 1..1 (or 0..0 or whatever)
Thread 2: 1..1
Thread 3: 1..2
Thread 4: 2..4
Thread 5: 4..5
Thread 6: 5..6
Thread 7: 6..6

现在,每个线程继续计算并将所有中间值相加在一行中。 当它遍历列的索引时,它也逐步跳过这个跳过的值列表,并跳过列表中出现的任何列号。 这显然是一个高效而简单的操作(因为列表也是按列排序的,就像合并)。

伪实施

我不知道CUDA,但是我有一些使用OpenCL的经验,我想接口是相似的(因为它们所针对的硬件是相同的)。下面是内核的一个实现,它处理一行(即计算result一个条目)在伪C ++中:

double calc_one_result (
    unsigned my_id, unsigned N, double const init_values [],
    unsigned skip_indices [], unsigned skip_begin, unsigned skip_end
)
{
    double res = 0;
    for (unsigned col = 0; col < my_id; ++col)
        // "f" seems to take init_values[column] as its first arg
        res += f (init_values[col], init_values[my_id]);
    for (unsigned row = my_id + 1; row < N; ++row)
        res -= f (init_values[my_id], init_values[row]);
    // At this point, "res" is holding "result[my_id]",
    // including the values that should have been skipped

    unsigned i = skip_begin;
    // The second condition is to check whether we have reached the
    // middle of the virtual matrix or not
    for (; i < skip_end && skip_indices[i] < my_id; ++i)
    {
        unsigned col = skip_indices[i];
        res -= f (init_values[col], init_values[my_id]);
    }
    for (; i < skip_end; ++i)
    {
        unsigned row = skip_indices[i];
        res += f (init_values[my_id], init_values[row]);
    }

    return res;
}

请注意以下几点:

  1. init_values和函数f的语义如问题所述。
  2. 该函数计算result数组中的一个条目; 具体来说,它计算result[my_id] ,所以你应该启动这个N实例。
  3. 它写入的唯一共享变量是result[my_id] 。 那么,上面的函数不会写任何东西,但是如果你把它翻译成CUDA的话,我想你最后还得写下这个函数。 但是,没有其他人写入result特定元素,所以这个写入不会引起数据竞争的任何争用。
  4. 两个输入数组init_valuesskipped_indices在该函数的所有正在运行的实例中共享。
  5. 所有对数据的访问都是线性和顺序的,除了跳过的值,我认为这是不可避免的。
  6. skipped_indices包含应该在每一行中跳过的索引列表。 它的内容和结构如上所述,只有一个小的优化。 由于没有必要,我已经删除了行号,只留下了列。 行号将作为my_id传递到函数中,每个调用应该使用的skipped_indices数组的切片使用skip_beginskip_end来确定。

    对于上面的示例,传递给calc_one_result所有调用的calc_one_result将如下所示: [4, 6, 0, 5, 4, 3] calc_one_result [4, 6, 0, 5, 4, 3]

  7. 正如你所看到的,除了循环,这个代码中唯一的条件分支是第三个for循环中的skip_indices[i] < my_id 。 虽然我相信这是无害的,完全可以预测的,但是这个分支在代码中很容易避免。 我们只需要传入另一个名为skip_middle参数,告诉我们跳过的项目在主对角线上的位置(即,对于行my_id ,跳过的索引skipped_indices[skip_middle]是比my_id大的第一个参数)。

结论是

我绝不是CUDA和HPC的专家。 但是如果我正确地理解了你的问题,我认为这种方法可能会消除任何和所有的内存争用。 此外,我不认为这会导致任何(更多)数值稳定性问题。

实施这个的成本是:

  • 总共调用f两次(并记录调用row < col以便将结果乘以-1
  • 跳过的值列表中存储两倍的项目。 由于这个列表的大小是在成千上万(而不是数十亿),它应该不是什么大问题。
  • 对跳过的值列表进行排序; 这又由于其规模,应该没有问题。

更新 :添加了伪执行部分。)





cuda