GPU BVH Construction

对100000面的Dragon生成的BVH

最近需要在GPU上实现BVH的并行高速构建方法,于是就照着一篇NVidia的文章《Thinking Parallel, Part III: Tree Construction on the GPU》做了。

此遍历的算法是完全并行的,因为在建树时早早已经将所有的片元(譬如三角形或是Vec3)根据Morton Code来进行排序。神奇之处也就在于在如此排序之后,对于某个片元,即可直接通过Morton Code来得到它在最终的表示BVH的二叉树中的位置,而这样所以能成立的原因则是这里的Code与Z-Curve是等价的,按数值大小从小到大排列就会呈现出以递归的方式填满空间的规律。

NVidia的文章在讲解此算法时引用了其文章作者本人在2012年的文章《Maximizing Parallelism in the Construction of BVHs, Octrees, and k-d Trees》( http://devblogs.nvidia.com/parallelforall/wp-content/uploads/2012/11/karras2012hpg_paper.pdf )。此篇文相比起先前的工作最大的进步就在于通过这个Morton Code的方法使得在GPU中每个线程之间完全没有任何的依赖关系,从而达成非常高的执行效率;相比起先前需要从上到下一层层扩张,其进步尤为明显,可以达到100%的核心利用率。不过这只是整个计算BVH的过程当中的一步而非全部;整体的过程应该是:

  1. 给所有的片元赋上一个Morton Code。如果是三角形,可以取其几何中心。再把Morton Code进行排序。如果有重复的Morton Code,还要去除重复。这也就意味着一个Morton Code可能对应多个片元。
  2. 用Morton Code进行排序并生成BVH。
  3. 把片元填充进BVH,再从叶子结点到根结点进行Reduce,写入包围盒的值。

在以上的第3步中,利用率还仍是会因为层层Reduce时树的结点数变少而变低,所以性能的提升主要来自于第2步。

在第二步中,需要在两个地方用到二分搜索的模式;第一个是找「覆盖的结点的下标的范围」,第二个是在第一步的范围中找「分界点」。写起来时和刷题中经常遇到的二分搜索是一样的。

此篇文章的算法只在所有的Morton Code都没有重复的情况下才有用。如果有重复,则会出现一个叶子结点对应多个父结点的情况,会导致Reduce的过程无法进行。这也是因为如果有重复的Morton Code时,第二步中的找「范围」会在不同次迭代时找到同一个分界点,这就导致了这个问题。这个问题也可以在BVH生成以后做一个简单的一致性检查来确认:只要看每个结点的parent和left/right是否能匹配就可以知道了。

在Reduction时,对每个结点的更新需要用加锁来确保原子操作,而对每个结点的访问亦要用atomicAdd来统计次数,以避免重复计算,不然在有锁的情况下,几千个线程从叶子结点涌向极结点造成的锁竞争可是要让程序永远都完不成了。

既然用了锁,就要设法回避死锁的问题。一种方法是加上 -G 打开 debug模式;更好方式是在 if{} else{} 之后加上一个随便什么语句,譬如 __asm(“pmevent 0;”) 。不过即使是打开 debug模式,性能会都比 CPU上 用最朴素的方法做要高出许多,更不用说打开 -O3 模式了。其中两个Kernel的CPU与GPU的性能对比如下:

CPU: i5-3210M
GPU: NVidia NVS 5200M

Morton Code生成与排序时间
CPU上排序:19 ms

构建BVH的时间(100000面,98848个不重复的Morton Code)
CPU,层次法 = 30 ms
CPU,Karras文中的方法,单线程 = 59 ms
GPU,Karras文中的方法,<<<120, 128>>> = 1.3 ms !!!!

计算包围盒的时间
CPU,从下往上 = 196 ms
GPU,并行 reduction = 7 ms !!!

所以可见这个方法是相当之有效的,而且适应不同大小的GPU。论文中用的是GTX480,前述网页的文章中用的是 GTX690,横跨了Fermi与Kepler这两代,说明在不同代的GPU之间也是应当能够通用的。(可能就除了一些细微的可能死锁的情况、或者可能发生非法内存访问但勉强蒙混过头的情况下不通用外,其它情况都应当是通用的吧)

发表在 Programming | 留下评论

HackerRank University CodeSprint 3

这是上周周末(9月29日~10月1日)的University CodeSprint 3。一共有7道题,但是只会做4道…其中第4道还是在网上搜索时搜出来的。

难度就用颜色表示:■■■■■■■■■■

  1. A Small Step Towards Calculators
    就是给定A{四则运算}B这么个字串,输出计算的结果。
    因为题目中不会有负数,所以可以用std::cin读入第一个数,再读运算符,再读第二个数,就完成了。
  2. Erupting Volcanoes
    就是格子中有若干火山,火山有随距离而减小的影响值,这值可互相叠加。那末要在给定所有火山时,求格子中最大的影响值。
    比较像LeetCode或是LintCode中的数邮局的题目。
    可以直接遍历每个格子,什么优化都不用加,也能完成。
  3. The Snake Vs The Wind
    贪吃蛇从格子的四角中的某角以某个方向出发,以贪婪的方式选定当前方向,直到遍历完所有格子,问蛇的轨迹。
    或者说,遇墙或遇到自己的身体就会往左或往右转,具体往哪里转依赖于当前风向:贪吃蛇偏向于顺风而行,垂直于风而行次之,逆风而行再次之。
    那就应该用模拟的方法,把蛇的方向模拟出来。用了一个Vec2 { int x, y; }来表示方向。在选定方向时以风向与方向的点积作为排序依据取最大者。
  4. Bob’s Game
    棋盘格上有墙和国王,国王可往左/上/左上移动,不可移至墙中;国王可以相互重叠。Bob与Alice依次行动,每次行动可以移动一个国王。两个玩家都是理性的。当前玩家若无法行动则输。求问给定棋盘Bob是否必胜,若必胜则第一步有多少种不同的走法。
    这篇文章所说,这个题目是一个「组合博弈」=「Combinational Game Theory」问题。是由「棋盘上仅有一个国王」这种简单的博弈问题重叠在一起而成的。只要对棋盘格求「S-G值」=「Sprague-Grundy」值,然后把所有国王所在的格子的S-G值进行XOR,就可知当前局势下,Bob是否能必胜。具体求取的方法是对所有出结点的S-G值的集合进行Mex=Minimum Excludant运算。在此题目中,该运算只可能得到0、1、2、3这四个值。
    题目中问的是Bob有多少种行动方式。那可以罗列不同的行动方式,再判定当前Alice是否必输。在罗列时欲求新的总体XOR的方法有点像快速求矩阵中「除了这个元素以外的所有元素的乘积」的方法。这样就能知道Bob的第一步有多少种行动方式了。
  5. Black White Tree 【事后知道了解法】
    无向树有黑色和白色结点。定义在某子树中白色结点数减去黑色结点数的差值的绝对值为「怪异指数」。求此树中最怪异的子树为何,并列出该子树中所有结点。
    如果从所有叶子结点开始向中心「爬山」,然后从最高处开始向外搜索,这样可以得到部分的分数…错误原因需要检查一下
    Update:解法类似打劫房屋;随意取某个点为树的根,然后罗列出从每个结点开始的子树中「怪异指数」最大者;取得之后再在回溯的过程中将解重建出来便是。在这里得到的收获就是「随意取一个结点作为根」并「罗列以每个结点开始的子树的情况」
  6. March of the King 【只得了部分分数】
    乍一看好似「Word Search」的变体。不过格子大小是固定的8×8,并且欲搜索的词也只有一个。
    直接做会有一半超时,如果将欲搜索的词的后四位存入std::vector中,则可以只有最后两个输入超时…(后4位是试出来的比较好的参数,太多或太少了都会更慢)还需要检查一下
  7. Simple Tree Counting 【只得了部分分数】
    无向树中的边有若干元素。为每个同样颜色的边组成的连通分量定义「Productivity」为N(N+1)/2,其中N为结点数。
    有许多Query,可能会改变边的颜色。需要频繁地求「某种颜色的所有连通分量的Productivity之和」与「包含某边的连通分量的Productivity」。直接做能拿到部分分数。
    猜想可能需要一些高级数据结构进行优化。

接下来…要对照一下正确解答,看一下后面三题怎么做。

话说在HackerRank上如果表现好是会有猎头来联系的!不过这回因为已经早早接了别家了所以就只能说「以后再保持联系」啦。

发表在 Programming | 留下评论

这是一间小房子

目前这个渣了几个月的体素显示程序(还不能称为是「引擎」吧)在加上了Greedy Meshing之后,可以在HD 4000上以比较高的帧率显示带AO但是不带阴影的场景,大小大概是300x300x80,并且包括数十个sprite。

在Greedy Meshing的过程中曾经想过关于u、v的延展方式的问题。在《0FPS》网站的链接中,这个Meshing的做法是先沿着u方向延展到最长可能的地方,然后再在v方向上延展的。起初还有点担心是否需要在延展u的过程当中试探不同的v变化量(这样想得到最佳解需要的就是DP而不是贪心了),不过再仔细想一想就可以知道其实两者是等价的,因为如果将当前正在处理中的平面旋转90度,贪心的解法就可以变成DP的解法了。而在「分割成尽可能少的矩形」这个目标中,中途尝试不同的v变化量未必就能使得最终的面数更少。所以就可以放心地用贪心法。

下一步的话需要将阴影部分好好修改一下。目前只有非常naive的阴影图,出现了漏光和锯齿边等不难看的视觉瑕疵。

螢幕快照_2017-07-20_10-13-09

螢幕快照_2017-07-20_10-13-34

发表在 Programming | 留下评论

二分法+范围最大值查询:Distant Pairs

Abstract:
近日看了 Hackerrank 上的一个题, Distant Pairs ,来自于 101 Hack 第 45 辑的。不会做,看了题目作者的解答才知道正确的做法。但是这道题中,涉及到了二分法与 Range Maximum Query 这两个有用的工具。其中,二分法是和先前遇到的 LintCode 617, Maximum Average Subarray 一样针对答案进行二分;而 Range Maximum Query 则是线段树的应用,将范围 aggregate 起来后通过树形结构快速查找的例子。

题干是这样的:

将 [0 … C-1] 这 C 个数,等距离围在一个周长为 C 的圆形上。举例来说,若是 C = 4 ,那就像是在钟面上有 0点、3点、6点 和 9点 这四个点。
将两点 a 和 b 所组成的点对记为 p(a, b) 。然后,我们在这个圆上选取 n 个这样的点对。每个点对中两点之间的距离记为 d(a, b),就是从 a 到 b 的比较近的那条弧的长度,即是 min(|a-b|, C-|a-b|)。
接下来我们考虑两个点对 p(a_i, b_i) 与 p(a_j, b_j) 。把这四个点两两相连就会有 6 个距离,定义两个点对之间的距离, d(p(a_i, b_i), p(a_j, b_j) ) 为这 6 个距离中的最小者。即是:
d(p(a_i, b_i), p(a_j, b_j)) = min( d(a_i, b_i), d(a_i, b_j), d(a_i, a_j), d(b_i, a_j), d(b_i, b_j), d(a_j, b_j) )。
那么想要问的就是:给定了 n 个点对之后,求所有 d(p(a_i, b_i), p(a_j, b_j) ) 中的最大值。

这题目如果用暴力做法就是枚举所有的 p(a_i, b_i) 与 p(a_j, b_j) 组合,这样需要 O(n^2) 次枚举。

二分的解法是针对可能的答案进行二分,这和 LintCode 617 中的 Maximum Average Subarray 是类似的思想。不过本题的答案必须介于 0 与周长 C 之间,所以二分的次数就是 log(C),不像 LintCode 617 需要对 double 类型的答案手动设定一个精度阈值。

在本题中,作者是通过这些观察来构造一个高效的答案的:

  1. p(a, b) 与 p(b, a) 没有区别。那就规定 a <= b。这样能方便很多,而且能去除重复。
  2. 对所有的 p(a, b),按先 a 后 b 的顺序升序排列。
  3. 在枚举 p(a_i, b_i) 与 p(a_j, b_j) 时,因为是组合,所以与顺序无关,可以规定 i 一定大于 j 或者 i 一定小于 j。这样可以对所需要枚举的空间进行更有效的裁剪。
  4. 在二分测试中,需要判定“是否可能 d(p(a_i, b_i), p(a_j, b_j))”大于或等于当前的阈值 thresh。这就意味着所罗列出的点对本身的弧长需要大于等于阈值。这样就又可以再裁剪一些。
  5. 在二分测试中,作者仍然是有两重的循环,即是类似暴力作法中的 i 循环与 j 循环。然后,对于每个 i 循环中的 p(a_i, b_i),满足阈值的 p(a_j, b_j) 若是存在,则 a_j 与 b_j 必须都离开 a_i, b_i 至少 thresh 的距离,因为只有这样,才能使得 d 中的六条弧都大于等于 thresh 。画成图,再展成一维数组,就是以下这样:
    1
    图例:
    红色 表示在 a_i 与 b_i 周围长度为 thresh 的禁止区域;
    绿色 表示 a_j 与 b_j 可以取值的区域;
    0 至 C-1 是圆上的点的一维坐标。在展成一维数组时有三种情况;如果规定 i 必须大于 j,那就只需考虑其中两种情况。
    Case 1:[0, a_i) 之间 没有绿色区域。(如果规定 i 必须大于 j,这种情况是不需要考虑的。因为绿色区域都)
    Case 2:[0, a_i) 之间有一些绿色区域,而且它横跨了 0 点,另外半边在 (b_i, C-1] 之间。
    Case 3:[0, a_i) 之间有绿色区域,并且没有横跨 0 点。

    这样,在 j 循环中的测试就变成了:“有没有一个点对 p(a_j, b_j),使得 a_j 和 b_j 都处于绿色的区域中。(如果处于相同的绿色区域中,那么 d(a_j, b_j) 也要大于等于 thresh。)”

  6. 以上的测试转变成区域查询,就是:给定绿色区间中的终点 b_j,它是否有对应的起点 a_j,也处于绿色区间中。这个内层循环中会向查询结构中不断增加有可能被选中的点对 p(a_j, b_j) 。内层循环是可以重复使用已有信息的,原因是:
    * 因为规定了 i 必须大于 j,所以内层循环滞后于外层循环;* 又因为 a_i 的左边 (a_i – thresh, a_i] 都是红色区域,所以内层循环滞后外层循环 thresh 个下标。
  7. 因为绿色区域会随着外层循环中的 i 的增加而向后挪动,所以后加入的点对 p(a_j, b_j) 可以替换掉先加入的点对。所以查询结构是 b => a 取最大值。
  8. 总的来说,在 i 一路增大的过程中,查询结构没有遗漏任何情况。查询结构进行的查询类型是 Range Maximum Query,其结构是给定 b_j,查询所有能到达此 b_j 的边中,a_j 最大的那条边。

总之,这个题目将二分查找法与 Range Maximum Query 结合到了一起,复杂度也从 O(n^2) 降低到了 O(n * log^2 C) 。n 是外层循环的次数, log 平方 C 是 Segment Tree 的查询(树的深度是 O(logC),每层中又需要将所查区间分割,一共会分割出 O(logC) 段。)

发表在 Programming | 留下评论

PA10之开始:由PA6的Parser想到的一些问题

最近打算着实地开始将两年前结束的 CPPGM 的 PA6 解析器重新用起来,用于解析(文曲星上的)GVMaker 文件。选择 GVMaker 主要还是想要完成以前的一个愿望,而利用 PA6 的解析器则是为了让 CPPGM 这个课的参与能稍微再完整一些(虽然已经结束两年多了,)能够在一个比较完整但是整体规模(不光是语法,还包括之后的代码生成和运行时之类的)也不是很大的语言上做一个相对完整的编译器流程。而取名为 PA10 ,只是因为它在 PA9 之后而已。

这次主要是想清楚了两件事情:

  1. 意识到一次解析的过程和结果都能用来表示,过程是「决策树」,结果是「语法树」,而结果是过程的一个子集;
  2. 用了一个类来维护在遍历和修改这棵树时的各种状态,减少代码累赘。

现在分别说以上两点。

※ 决策树与语法树 ※

决策树的建立方法是这样的:在解析到一个规则时,每根据这条规则右边的产生式产生一个新的非终结符,就给当前结点添加一个子结点。也就是说一个结点可以有N个子结点。这决策树是在解析器试探输入文件时的调用栈所扫过的痕迹。

而相比之下,语法树和决策树一模一样,但是只是只在成功的情况下才保留子结点,否则则删除子结点。另外语法树的结点中还包括终结符,不像决策树的结点中只包括非终结符。

所以如果给了以下输入:

int a;

这个文件的全部内容就是声明一个全局变量 int a。它最后解出来的语法树也反映了这一点:

Parse Tree of

树中,“int”是一个“simple_type_specifier”,“a”是一个“id_expression”,加上分号,三者构成了一个“simple_declaration”,最后组成了一个“translation_unit”。

这个语法树所对应的决策树是这样的(其中绿色的结点对应着语法树中的结点):

Decision tree generated while parsing

可见在解析时,解析器进行了两次“declaration”的试探。其中第一次试探发现了“int a;”的声明,而第二次试探完成后发现并没有更多的声明,就结束了。

整个 parsing 的过程,也就是根据每次解析的结果成功与否不断对决策树进行调整,直到找到符合语法规则且吃掉了所有 token 决策树时将树输出的过程。

※ 记录(bookkeeping)的事情 ※

记录的事情有以下这些:

  • 表示当前这个非终结符走到了第几个分岔(即是决策树中的分岔),例如
    noptr-declarator-root:
        declarator-id attribute-specifier*      
    
  • 每个分支成功与否
  • 当前的 Token 的位置。
  • 当前的递归深度。递归深度并不一定等于分岔的位置,就是因为在一个非终结符所对应的一条规则中,如果出现了加号或是星号,就可能进行多次尝试。
  • 当前所正在生成的语法树结点以及父结点。

这几点在 PA6 的解析器中,都直接以变量的形式存在。这样每次新增一条语法规则时都太容易出错了。为了减少出错,这回把这些信息都存在一个称为 ParseState 的表示解析状态的类中。它与其它部分的联系是:

  • 从「整个解析器」的角度看,在进入最外一层(translation_unit)时:
    第一种可能是刚进入,此时的「决策树」是空的;或者
    第二种可能是接续上一次失败或成功的解析,将「决策树」前进一格,然后再进行尝试;
  • 从「非终结符」的角度看,在进入一个非终结符的这条规则时,相应地也有两种可能:
    第一种可能、刚进入时,需要往「决策树」中新加入这次尝试所走的分岔;
    而在第二种可能、也就是因循上一次解析时,需要从「决策树」中读出这回应该尝试的岔路;而在「因循上一次的决策」结束时,又要回到第一种情况,向决策树中添加新的结点;
    所以,就有两个地方会对「决策树」进行修改:第一个地方是完成一次解析后的「进一格」,第二个地方就是在非终结符中新增新和结点。
  • 从「一条规则」的角度看,在需要解析一个非终结符时:
    如果是首次出现的非终结符,那就需要将语法树中深入一层;
    而如果不是首次出现的非终结符,那就保持语法树的层次。
    也就是类似于在深度优先搜索时需要记下每一次访问到了第几个子结点那样的。
    而这个地方的处理也要和分岔的进位一起考虑起来。

因为在解析的过程中 ParseState 会一直被改变,比如他的 Token 位置会一直前进和后退,分岔和递归深度需要一直更新或回退。所以在设计时就有了以下函数:

ParseState::OnFirstNT() 此条规则中,所遇到第一个非终结符时的拷贝构造
ParseState::OnNextNT() 这条规则中,所遇到的其后的非结结符时的拷贝构造
ParseState::NextChoiceFrom() 是在这一条分岔走完时用的拷贝构造函数

PA10Parser::ConsumeToken() 是搭配 OnFirstNT 所使用的
PA10Parser::ConsumeTokenToParent() 是搭配 OnNextNT 的。(*)

PA10Parser::OnEntry() 在进入一个非终结符时需要做的事情
PA10Parser::OnExit() 在离开一个非终结符时需要做的事情

(*):这两个 ConsumeToken 的存在是因为一条规则中所沿用的 ParseState 有数据依赖关系,在而这些相依赖的 ParseState 会因为某时调用了 OnFirstNT() 而使得 AST 向下进了一层。对于非终结符来说,一定要进一层;但是对于 Token 来说,并不一定要进一层。所以才会有「对于 Token,在进一层之前,调用 ConsumeToken ,而在进一层之后,调用 ConsumeTokenToParent 」的现象。

以上这些,还要搭配正确的赋值,并且将每个非终结符中的不同分岔以一个 switch 语句括起来,才能正确地处理解析中的其它情况,比如说带有星号、加号、问号的规则。

总结起来

这个 parser 目前尝试过的最大的输入是 yan 的编译器,算上所有的包含命令,一共有大约 81000 个 token 。每次 parse 用时大约 20 秒。

发表在 Uncategorized | 留下评论