在HPC (high performance computing)领域, 大规模计算问题一直是科学和工程计算研究的热点. 然而, 在机器学习^[1], 数据挖掘^[2], 图像及信号处理^[3-5], 计算流体力学^[6], 天体物理学^[7] 和量子化学^[8]等许多科学计算应用中存在着大量小规模矩阵计算问题, 目前对同时处理大量小规模矩阵计算问题的研究得到越来越多的关注. 此外, 在一些超大规模计算问题的算法中, 也存在着对通过区域分解产生的大量小区域同时进行局部处理的问题. 这些领域需要处理的矩阵的行、列数通常在几十到几百之间, 但矩阵数量可能多至上万. 此类问题需要采用适合同时处理大量小规模矩阵计算的高效算法. 我们称这类同时处理大量小规模矩阵的问题为批量矩阵计算问题.

在解决Cholesky分解、LU分解和QR分解等常见的矩阵计算问题时, 为了更好地实现数据的复用, 通常会采用将矩阵分成一系列列条块, 然后逐个处理这些列条块的块算法. 这些基于块算法的矩阵分解通常分为两步: 条块内分解和尾部矩阵更新. 对于在GPU (graphics processing unit)架构上实现大规模矩阵分解问题, 文献[9-11]已经给出了一个异构并行算法及其实现方法, 由于条块内分解计算为BLAS-1 (向量和向量运算)和BLAS-2 (矩阵和向量运算), 这并不能充分利用GPU的计算能力, 而尾部矩阵的更新是一个可以充分发挥GPU计算能力的BLAS-3 (矩阵和矩阵运算)计算, 因此可将条块内分解放在CPU上计算, 而尾部矩阵更新放在GPU上进行计算, 该异构算法可将与条块内分解相关的CPU-GPU间数据传输和GPU上计算进行重叠. 对于大规模问题, 这种CPU-GPU异构算法能够得到较好的性能, 但对于小规模矩阵问题, 由于GPU计算时间相比于CPU-GPU之间的数据传输耗时较少, 因此对批量小规模矩阵, 这种异构算法带来的性能提升是有限的^[12]. 因此对于批量小规模矩阵的LU分解和求逆在GPU上的实现, 通常采用将条块内分解和尾部矩阵更新都放在GPU上处理的方法^[12-17].

GPU是一种由大量运算单元组成的大规模并行计算架构, 可以支撑大量数据的并行计算, 适合批量处理具有相同计算任务的工作. 近年来已有许多工作关注于在GPU内批量矩阵处理运算, 文献[18-21]研究了在GPU内关于批处理小规模矩阵-矩阵乘法GEMM问题的实现. 文献[12]描述了批量处理矩阵常见的Cholesky分解、LU分解和QR分解在GPU上的实现. 批量LU分解和批量矩阵求逆是许多科学计算和应用领域的关键计算问题^{[11, 22-25]}. 针对GPU架构, 基于批量LU分解的矩阵求逆算法实现方法得到了深入研究, 著名GPU生产商NVIDIA公司研发的CUBLAS (CUDA basic linear algebra subroutine)库^[26]提供了批量LU分解与批量求逆的实现. GPU数值代数数学库MAGMA (matrix algebra for GPU and multicore architectures)中也提供了批量LU分解与批量求逆在GPU上的实现^[27]. Villa等人^[14,15]给出了部分选主元与完全选主元的方法批量LU分解算法. Dong和Haidar等人在文献[12,13]中给出了在GPU上批量多级Right-looking LU分解算法. Abdelfattah等人在文献[17]中给出了极小规模的矩阵(<32)的批量 LU分解与批量求逆算法.

本文第1节给出在GPU上实现批量LU分解与求逆的相关工作, 以及我们实现方法的优势与取得的效果. 第2节描述不同的求解LU分解和求逆的算法. 第3节给出在GPU上批量LU分解算法与批量求逆算法的设计实现和优化. 第4节给出我们实现的算法的实验结果与分析. 最后给出本文工作的总结.

1 相关工作

对于批量LU分解算法在GPU上的实现方法, Villa等人在文献[14]中讨论了针对GPU架构的3种不同级别的并行方法. 第1种基于CUBLAS的Warp级并行(一个Warp对应一个矩阵), 批处理函数为每个矩阵分配一个Warp (32个线程), 因此矩阵的规模不能超过32×32. 此实现严重依赖共享内存, 但是共享内存的内容不会在不同的内核调用之间保留, 因此每个内核都必须进行额外的工作才能将数据重新缓存在共享内存中. 第2种Thread-Block级并行(一个Block对应一个矩阵). 出于性能方面的考虑, 将矩阵完整地加载到共享内存中, 这意味着可以处理的矩阵的最大尺寸受到可用共享内存的限制. 对于Fermi和Kepler之前的GPU架构, 此实现可以处理尺寸为76×76 (双精度)的系统. 第3种针对矩阵规模小于128时, Thread级并行(一个线程对应一个矩阵), 此种方法每个线程在其对应的矩阵上进行操作, 在选主元的过程中, 由于每个矩阵的主元是不同的, 则每个线程寻找主元的过程是不同的, 则同一个Warp中的不同线程被分配了不同的任务, 这可能会由于Warp中的线程发散导致性能不佳. 第1种和第2种方法仅支持部分选主元, 第3种方法支持部分选主元和完全选主元中方法. 特别地, 第3种Thread级并行方法用于求解一个线性方程组时性能不如CUBLAS中的实现方法^[13]. 该方法已用于地下运输模拟, 对流动路径中的许多化学和微生物反应进行了并行模拟^[15].

Dong和Haidar等人在文献[12,13]中实现了在GPU上批量Right-looking LU分解块算法. 提出了多级Right-looking LU分解块算法. 将LU分解分为4个过程, 在GPU上启动4个Kernel完成. 第1步是当前列条块的分解, 在分解过程中并没有将整个列条块从全局内存写到共享内存, 仅将一列写到了共享内存. 第2步是行交换, 给出了NB次行交换并行交换方法, 这里的NB表示块算法中块的大小. 第3步求解 $ \boldsymbol{U} $ 的子矩阵块 $ {\boldsymbol{U}}_{12} $ . 第4步修正尾部矩阵, 这也是多级Right-looking LU分解算法与普通Right-looking块算法的不同之处. 通常的Right-looking块算法对当前列条块分解完成后, 对整个右边未求解的尾部矩阵进行修正, 而在多级块算法对尾部矩阵进行修正过程中, 仅对下一个列条块部分进行修正, 推迟对剩余部分的修正, 这是软件包MAGMA采用的方法.

Abdelfattah等人在文献[17]中实现了部分选主元的批量矩阵LU分解, 也实现了求批量矩阵的逆矩阵, 但他们考虑的矩阵大小不超过32. 在实现时因为矩阵非常小, 为了加速数据的访问, 因此将整个矩阵保存在寄存器中. Block中的线程采用了一维结构, 因此实现时不需要同步点. 另外, 由于矩阵规模小, 若一个Block对应一个矩阵, 则一个Block中的线程的数目太少, 为了提高占用率, 采用一个Block同时分解多个矩阵的方法. 对于LU分解中的行交换步骤, 采取了lazy swap技术, 每个线程记录它对应行是否是最大主元所在行, 只有不是主元所在行对应的线程对尾部矩阵进行修正, 并且在下一次分解时选取最大主元. 分解完成后, 每个线程都知道其行的最终位置, 并将其直接写入全局内存, 此文采用的是Right-looking LU分解块算法.

国内对于在GPU上实现批量LU分解与求逆的研究比较少. 对于单个矩阵的LU分解在GPU上的实现, 文献[28]实现了部分选主元的Right-looking LU分解算法. 其中, 选取主元的部分在CPU上实现, CPU每次找到主元并进行行交换后, 再把数据传到 GPU上, GPU负责处理主元列的并行计算和子矩阵元素的并行计算. 在GPU上计算时, 由于共享内存的大小限制, 当矩阵规模过大时, 不能将整个矩阵放到共享内存中, 因此采取了矩阵分段计算: 每次将两列元素放入共享内存计算. 将LU分解的过程分成在CPU和GPU上进行, 则每次迭代都需要CPU和GPU之间数据传输, 频繁的数据传输将会影响程序的性能. 此外, 该论文中实现的是非块的Right-looking LU分解算法, 大多数计算都是BLAS-1操作, 并不能充分利用GPU的计算性能. 此算法并不适用于批量的LU分解. 文献[29]中讨论了电力系统潮流计算当中对稀疏矩阵的基于GPU并行加速的批量 Left-looking LU 分解算法.

在GPU上实现批量LU分解算法时, 我们采用了一个Block对应一个矩阵的Thread-Block级并行方法, 不同于Villa等人在文献[14]中Thread-Block级并行方法, 我们通过使用块算法, 只需将要处理的当前列条块放在共享内存中, 而不是将整个矩阵放在共享内存中, 这样处理的矩阵规模就不受共享内存大小的限制. 在GPU内实现LU分解之后的批量求逆算法时, 我们同样采取了Thread-Block级并行块方法. 由于一个线程对应矩阵的一行, 并且列条块内每行的求解相互独立, 因此可将当前列条块存储在线程的私有寄存器中, 加快数据的读写速度. 我们的算法通过更好地利用GPU的各级内存, 提升了算法的性能. Dong等人^[13]在GPU上实现批量多级Right-looking LU分解算法时, 将LU分解分为4个过程, 在GPU上启动4个Kernel完成. 多次启动Kernel, 共享内存的内容不会在不同的Kernel调用之间保留, 因此每个Kernel都必须将数据重新缓存在共享内存中, 这存在数据的重复读取, 对性能造成一定的影响, 当矩阵规模较小时造成的影响更大. 因为在GPU上影响批量小矩阵处理性能的主要因素是数据的读写时间^[30]. 在GPU上实现批量LU分解时, 由于Left-looking算法比Right-looking算法具有更少的对全局内存数据访问量, 因此我们采用了基于Left-looking方法的批量LU分解算法. 并且将整个LU分解过程放在一个Kernel中完成, 这样只需要将要处理的当前列条块从全局内存中读取一次, 求解完成后再写回即可, 避免了数据的重复读取.

在GPU上实现批量LU分解算法时, 由于不选主元的LU分解算法数值不稳定, 因此我们采用了列选主元的LU分解算法. 以列优先存储的矩阵在GPU上进行LU分解时, 由于选主元存在矩阵的行交换过程, 而同一行的矩阵元素地址是不连续的, 造成线程的访存地址不连续(non-coalesced memory accesses), 大大降低了在GPU中实现的LU分解算法的性能. 为此, 我们提出了Warp分组行交换和行交换延迟2个优化方案降低了行交换过程中访存地址不连续带来的性能影响. 在TITAN V GPU上对10000个规模在33–190之间的随机矩阵进行了测试, 测试的数据类型为单精度复数、双精度复数、单精度实数和双精度实数. 我们实现的批量LU分解算法的浮点计算性能分别可达到约2 TFLOPS、1.2 TFLOPS、1 TFLOPS、0.67 TFLOPS, 并同CUBLAS和MAGMA中的批量LU分解算法的实现进行了比较, 与CUBLAS中的实现相比加速比分别最高达到了约9×、8×、12×、13×, 与MAGMA中的实现相比加速比分别达到了约1.2×–2.5×、1.2×–3.2×、1.1×–3×、1.1×–2.7×.

LU分解之后的求逆块算法包括列条块求逆和尾部矩阵修正两个阶段. 为了减少修正过程对全局内存的读写数据量, 我们在批量求逆算法的GPU实现中提出了延迟修正的矩阵求逆块算法. 此外, 为了减少对全局内存的访问和加快数据的读写速度, 求逆算法采用了更多利用寄存器和共享内存的优化方法, 并设计了减少访存数据量的列交换方法. 另外, 为了避免采用一次性分配GPU资源造成的线程闲置和共享内存浪费问题, 本文设计了运行时按需分配线程、共享内存等GPU资源的动态资源分配方法, 相比一次性分配的静态资源分配方法, 性能得到了一定的提升. 在TITAN V GPU上对10000个规模在33–190之间的随机矩阵进行了测试, 测试的数据类型为单精度复数、双精度复数、单精度实数和双精度实数. 我们实现的批量求逆块算法浮点计算性能分别可达到约4 TFLOPS、2 TFLOPS、2.2 TFLOPS、1.2 TFLOPS, 与CUBLAS中的实现相比加速比分别最高达到了约5×、4×、7×、7×, 与MAGMA中的实现相比加速比分别达到了约2×–3×、2×–3×、2.8×–3.4×、1.6×–2×.

2 基于LU分解的矩阵求逆算法

通过矩阵A的LU分解求解A的逆比直接对矩阵A求逆更加高效. 基于LU分解的矩阵求逆算法主要分为两步, 首先对矩阵A进行选主元的LU分解得到 $ {\boldsymbol{A}} = {\boldsymbol{PLU}} $ , 其中L为单位下三角矩阵, U为上三角矩阵, P为置换矩阵. 然后对三角矩阵L和U分别求逆, 最后由置换矩阵P得到原始矩阵A的逆.

2.1 选主元的LU分解算法

选主元的LU分解就是将一个 $ n \times n $ 的可逆方阵A分解成 $ n \times n $ 置换矩阵P、 $ n \times n $ 的单位下三角矩阵L和 $ n \times n $ 的上三角矩阵U的乘积, 即 $ {\boldsymbol{A}} = {\boldsymbol{PLU}} $ . 在实现时, 置换矩阵P可以压缩存储在一个向量中. 算法1给出了对矩阵A进行LU分解的非块算法, 首先寻找当前第i列中最大主元所在的行号 $ maxID $ , 并存储在 $ {\boldsymbol{P}}[i] $ 中(第2行), 然后交换第i和第 $ maxID $ 行(第3行), 交换完成后对第i列元素进行修正(第4行), 最后对尾部矩阵进行秩–1修正(第5行) , 这里的计算为BLAS-1运算.

算法1. 部分选主元的LU分解非块算法.

1　for $ i = 1 $ to $ n $

2　　 ${\boldsymbol{P}}[i] = maxID = \max \left( {{\rm abs} ({\boldsymbol{A}}[i:n, i])} \right)$ ; /*寻找最大主元所在的行*/

3　　 $ {\boldsymbol{A}}[i, 1:n] \leftrightarrow {\boldsymbol{{ A}}}[maxID, 1:n] $ ; /*行交换*/

4　　 $ {\boldsymbol{A}}[i + 1:n, i] = {\boldsymbol{A}}[i + 1:n, i]/{\boldsymbol{A}}[i, i] $ ;

5　　 $ {\boldsymbol{A}}[i + 1:n, i + 1:n] - = {\boldsymbol{A}}[i + 1:n, i] \times {\boldsymbol{A}}[i, i + 1:n] $ ;

6　endfor

为了更好地实现数据的复用, 通常LU分解采用具有更多BLAS-3运算的块算法, LU分解的块算法将矩阵分成一系列宽度为 $ \mathit{NB} $ 列条块, 然后逐个处理这 $ bn\;\left( {bn = \left\lceil {n/\mathit{NB}} \right\rceil } \right) $ 个列条块. 常用的LU分解的块算法分为Right-looking LU和Left-looking LU两种^[20]. 算法2和算法3分别给出了这两种LU分解块算法, 在算法中求得的单位下三角矩阵L和上三角矩阵U分别覆盖矩阵A的上下三角部分, 矩阵A对角线存储的是矩阵U的对角线元素.

算法2给出了Right-looking LU分解的块算法, 首先对如图1所示的当前第i个列条块 $\left[ {\begin{array}{*{20}{c}} {{ {\boldsymbol{A}}_{11}} } \\ {{ {\boldsymbol{A}}_{21}} } \end{array}} \right]$ 进行部分选主元的LU分解得到 $ {{\boldsymbol{L}}_{11}} $ 、 $ {{\boldsymbol{L}}_{21}} $ 和 $ {{\boldsymbol{U}}_{11}} $ 并覆盖A的对应部分(第3行), 分解完成后对当前列条块外的两侧进行行交换(第5行), 行交换完成后立即求解矩阵U的子块 $ {{\boldsymbol{U}}_{12}} $ (第7行), 最后用 $ {{\boldsymbol{A}}_{21}}\left( {{{\boldsymbol{L}}_{21}}} \right) $ 和 $ {{\boldsymbol{A}}_{12}}\left( {{{\boldsymbol{U}}_{12}}} \right) $ 来修正如图1所示的当前列条块右边未求解的尾部矩阵 $ {{\boldsymbol{A}}_{22}} $ (第8行), 修正完成后向右继续处理下一个列条块, 直到最后一个列条块. 如图1所示, 由于算法2只关注了当前列条块及其右边部分, 因此称为Right-looking LU分解块算法.

算法3给出了Left-looking LU分解算法, 在求解如图2所示当前第i个列条块之前, 先用其左边前 $ i - 1 $ 个列条块已经求解出的 $ {{\boldsymbol{L}}_{00}} $ 和 $ {{\boldsymbol{A}}_{01}} $ 来求解 $ {{\boldsymbol{U}}_{01}} $ (第3行), 然后用左边前 $ i - 1 $ 个列条块已经求解出的 $\left[ {\begin{array}{*{20}{c}} { {{\boldsymbol{A}}_{10}} } \\ { {{\boldsymbol{A}}_{20}} } \end{array}} \right]\left( {\left[ {\begin{array}{*{20}{c}} { {{\boldsymbol{L}}_{10}} } \\ { {{\boldsymbol{L}}_{20}} } \end{array}} \right]} \right)$ 和 $ {{\boldsymbol{A}}_{01}}\left( {{{\boldsymbol{U}}_{01}}} \right) $ 对当前第i个列条块进行修正(第4行), 修正完成后对当前列条块进行部分选主元的LU分解算法(第6行), 最后对列条块外的两侧进行行交换(第8行), 然后向右处理下一个列条块, 直到最后一个列条块. 如图2所示, 由于算法3只关注了当前列条块及其左边的部分, 因此可以称为Left-looking LU分解块算法.

算法2. Right-looking LU分解块算法.

1　for $ i = 1 $ to $ bn $

2　　/*当前列条块的选主元LU分解*/

3　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{11}}} \\ {{{\boldsymbol{A}}_{21}}} \end{array}} \right] = {\boldsymbol{P}}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{L}}_{11}}} \\ {{{\boldsymbol{L}}_{21}}} \end{array}} \right] \times {{\boldsymbol{U}}_{11}}$ ;

4 　　/*列条块外部分的行交换*/

5　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{10}}}&{{{\boldsymbol{A}}_{12}}} \\ {{{\boldsymbol{A}}_{20}}}&{{{\boldsymbol{A}}_{22}}} \end{array}} \right] = {\boldsymbol{P}}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{10}}}&{{{\boldsymbol{A}}_{12}}} \\ {{{\boldsymbol{A}}_{20}}}&{{{\boldsymbol{A}}_{22}}} \end{array}} \right]$ ;

6　　/*修正尾部矩阵*/

7　　 ${{\boldsymbol{A}}_{12}}: = {{\boldsymbol{U}}_{12}} = {\boldsymbol{L}}_{11}^{ - 1} \times {{\boldsymbol{A}}_{12}}$ ;

8　　 ${{\boldsymbol{A}}_{22}} - = {{\boldsymbol{A}}_{21}} \times {{\boldsymbol{A}}_{12}}$ ;

9　endfor

算法3. Left-looking LU分解块算法.

1　for $ i = 1 $ to $ bn $

2　　/*修正当前列条块*/

3　　 ${{\boldsymbol{A}}_{01}}: = {{\boldsymbol{U}}_{01}} = {\boldsymbol{L}}_{00}^{ - 1} \times {{\boldsymbol{A}}_{01}}$ ;

4　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{11}}} \\ {{{\boldsymbol{A}}_{21}}} \end{array}} \right] - = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{10}}} \\ {{{\boldsymbol{A}}_{20}}} \end{array}} \right] \times {{\boldsymbol{{\boldsymbol A}}}_{01}}$ ;

5　　/*当前列条块的选主元LU分解*/

6　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{11}}} \\ {{{\boldsymbol{A}}_{21}}} \end{array}} \right] = {\boldsymbol{P}}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{L}}_{11}}} \\ {{{\boldsymbol{L}}_{21}}} \end{array}} \right] \times {{\boldsymbol{U}}_{11}}$ ;

7　　/*列条块外部分的行交换*/

8　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{10}}}&{{{\boldsymbol{A}}_{12}}} \\ {{{\boldsymbol{A}}_{20}}}&{{{\boldsymbol{A}}_{22}}} \end{array}} \right] = {\boldsymbol{P}}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{10}}}&{{{\boldsymbol{A}}_{12}}} \\ {{{\boldsymbol{A}}_{20}}}&{{{\boldsymbol{A}}_{22}}} \end{array}} \right]$ ;

9　endfor

2.2 LU分解后矩阵求逆算法

在矩阵A的LU分解之后得到 ${\boldsymbol{A}} = {\boldsymbol{PLU}}$ , 则LU分解之后的求逆过程可以分为如下3步.

(1) 求上三角矩阵 ${\boldsymbol{U}}$ 的逆 ${{\boldsymbol{U}}^{ - 1}}$ ;

(2) 求解下三角方程组 ${\boldsymbol{XL}} = {{\boldsymbol{U}}^{ - 1}}$ , 得到 ${\boldsymbol{X}}$ ;

(3) 对X进行列变换P, 得到 ${{\boldsymbol{A}}^{ - 1}}$ .

上述第(1)步求解矩阵 ${\boldsymbol{U}}$ 的逆可以通过求解方程组 ${\boldsymbol{XU}} = {\boldsymbol{I}}$ 或者 ${\boldsymbol{UX}} = {\boldsymbol{I}}$ 实现, 这里I为单位矩阵. 其中方程组 ${\boldsymbol{XU}} = {\boldsymbol{I}}$ 适合按列优先存储的矩阵, 方程组 ${\boldsymbol{UX}} = {\boldsymbol{I}}$ 适合按行优先存储的矩阵. 由于本文矩阵采用了按列优先的存储方式, 因此 ${\boldsymbol{U}}$ 的逆通过求解上三角方程组 ${\boldsymbol{XU}} = {\boldsymbol{I}}$ 获得.

与LU分解算法一样求逆过程也通过使用块算法使其具有更多BLAS-3运算. 求逆块算法将矩阵分成一系列大小NB的列条块, 然后逐次对这 $bn\;\left( {bn = \left\lceil {n/NB} \right\rceil } \right)$ 个列条块进行求解. 求逆块算法包括列条块求逆和尾部矩阵修正2个模块. 根据对尾部矩阵是否及时修正, 可以分为及时修正与延迟修正2种块算法. 算法4和算法5分别给出了这2种求逆块算法. 这里假设LU分解得到的上三角矩阵为 ${\boldsymbol{U}}$ , ${{\boldsymbol{U}}^{ - 1}}$ 通过求解上三角方程组 ${\boldsymbol{XU}} = {\boldsymbol{W}}$ 获得, 且求得的逆矩阵的结果 ${{\boldsymbol{U}}^{ - 1}}$ 覆盖矩阵 ${\boldsymbol{W}}$ , 初始时 ${\boldsymbol{W}}$ 为单位矩.

算法4给出了及时修正的求解方程组 ${\boldsymbol{XU}} = {\boldsymbol{W}}$ 的块算法, 求解过程如图3, 初始时 ${\boldsymbol{W}} = {\boldsymbol{I}}$ , 且结果覆盖 ${\boldsymbol{W}}$ . 首先如图3(a)中当前第 $ i $ 个列条块(第2行), 当前列条块的结果可通过求解上三角方程组 $\left[ {\begin{array}{*{20}{c}} { {{\boldsymbol{X}}_{01}} } \\ { {{\boldsymbol{X}}_{11}} } \end{array}} \right]{{\boldsymbol{U}}_{11}} = \left[ {\begin{array}{*{20}{c}} { {{\boldsymbol{W}}_{01}} } \\ { {{\boldsymbol{W}}_{11}} } \end{array}} \right]$ 得到, 且求解的结果覆盖矩阵 ${\boldsymbol{W}}$ 的列条块 $\left[ {\begin{array}{*{20}{c}} { {{\boldsymbol{W}}_{01}} } \\ { {{\boldsymbol{W}}_{11}} } \end{array}} \right]$ , 求解完成后与 ${\boldsymbol{U}}$ 的子矩阵 ${{\boldsymbol{U}}_{12}}$ 修正列条块右边剩余未求解的尾部矩阵 (第3行), 然后向右进行下一个列条块的求解, 直到最后一个列条块.

算法5给出了延迟修正的求 ${\boldsymbol{XU}} = {\boldsymbol{W}}$ 的块算法, 在求解如图3(b)中当前第i个列条块之前, 首先用当前列条块左边已经求解出前 $ i - 1 $ 个列条块的结果 ${{\boldsymbol{W}}_{00}}$ 与 $ {\boldsymbol{U}} $ 的子矩阵 ${{\boldsymbol{U}}_{01}}$ 来修正当前第i个列条块 ${{\boldsymbol{W}}_{01}}$ (第2行), 修正完成后再求解当前列条块(第3行), 然后向右进行下一个列条块的求解, 直到最后一个列条块.

算法4. 及时修正的求 ${\boldsymbol{XU}} = {\boldsymbol{W}}$ 的块算法.

1　for i=1 to bn

2　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{01}}} \\ {{{\boldsymbol{W}}_{11}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{01}}} \\ {{{\boldsymbol{W}}_{11}}} \end{array}} \right] \times {\boldsymbol{U}}_{11}^{ - 1}$ ; /*求解当前列条块*/

3　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{02}}} \\ {{{\boldsymbol{W}}_{12}}} \end{array}} \right] - = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{01}}} \\ {{{\boldsymbol{W}}_{11}}} \end{array}} \right] \times {{\boldsymbol{U}}_{12}}$ ; /*修正剩余尾部矩阵*/

4　endfor

算法5. 延迟修正的求 ${\boldsymbol{XU}} = {\boldsymbol{W}}$ 的块算法.

1　for $ i = 1 $ to $ bn $

2　　 ${{\boldsymbol{W}}_{01}} - = {{\boldsymbol{W}}_{00}} \times {{\boldsymbol{U}}_{01}}$ ; /*修正当前列条块*/

3　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{01}}} \\ {{{\boldsymbol{W}}_{11}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{01}}} \\ {{{\boldsymbol{W}}_{11}}} \end{array}} \right] \times {\boldsymbol{U}}_{11}^{ - 1}$ ; /*求解当前列条块*/

4　endfor

前文所述求逆过程的第(2)步中, 通过求解一个右端项为 ${{\boldsymbol{U}}^{ - 1}}$ 的下三角方程组 ${\boldsymbol{XL}} = {{\boldsymbol{U}}^{ - 1}}$ 得到 ${\boldsymbol{X}}$ , 因为第(1)步得到 ${{\boldsymbol{U}}^{ - 1}}$ 的结果存储在矩阵 ${\boldsymbol{W}}$ 中, 因此可以通过求解方程组 ${\boldsymbol{XL}} = {\boldsymbol{W}}$ 实现, 并且求得的结果覆盖矩阵 ${\boldsymbol{W}}$ . 第(2)步求解下三角方程组的过程与第(1)步求解上三角方程组的过程类似, 块算法包括列条块求解和尾部矩阵修正两个模块. 根据对尾部矩阵是否及时修正, 可以分为及时修正与延迟修正2种块算法. 算法6和算法7具体描述了这2种方法, 对应的块算法示意图如图4所示, 初始时 ${\boldsymbol{W}} = {{\boldsymbol{U}}^{ - 1}}$ , 且结果覆盖 ${\boldsymbol{W}}$ . 求解过程从右到左. 这里假设LU分解得到的单位下三角矩阵为 ${\boldsymbol{L}}$ .

算法6. 及时修正的求 ${\boldsymbol{XL}} = {\boldsymbol{W}}$ 的块算法.

1　for $ i = 1 $ to $ bn $

2　　 ${{\boldsymbol{W}}_1} = {{\boldsymbol{W}}_1} \times {\boldsymbol{L}}_{11}^{ - 1}$ ; /*求解当前列条块*/

3　　 ${{\boldsymbol{W}}_0} - = {{\boldsymbol{W}}_1} \times {{\boldsymbol{L}}_{10}}$ ; /*修正剩余尾部矩阵*/

4　endfor

算法7. 延迟修正的求 $ {\boldsymbol{XL}} = {\boldsymbol{W}} $ 的块算法.

1　for $ i = 1 $ to $ bn $

2　　 ${{\boldsymbol{W}}_1} - = {{\boldsymbol{W}}_2} \times {{\boldsymbol{L}}_{21}}$ ; /*修正当前列条块*/

3　　 ${{\boldsymbol{W}}_1} = {{\boldsymbol{W}}_1} \times {\boldsymbol{L}}_{11}^{ - 1}$ ; /*求解当前列条块*/

4　endfor

3 基于批量LU分解的矩阵求逆算法在GPU上的设计实现和优化

在GPU上实现批量算法时, 如图5所示一个Block对应一个矩阵. 而对于Block内线程的配置, 可以是1D的也可以是2D的. 由于2D线程结构需要使用 $ n \times n $ 个线程对应矩阵的 $ n \times n $ 个元素, 限制了CUDA运行时每个SM (streaming multiprocessor)同时调度多个Block的能力. 另外对于大多数矩阵规模而言, 采用2D结构的线程个数通常会超出每个Warp的线程数目(32个), 这会导致在线程需要共享数据时2D会比1D结构更多必要的同步点. 并且, 当采用2D结构对当前列条块进行处理时, 由于一个线程对应矩阵的一个元素, 这会造成条块外的线程处于空闲状态, 导致线程利用率较低. 另外, 1D结构中一个线程对应矩阵的一行, 这使得每个线程被分配更多的任务, 指令级并行性更好, 因此我们采用了如图5所示的1D结构. 假设批量矩阵已经存储于GPU的全局内存上, 并且矩阵元素是按照列主序存储, 共有bs个矩阵.

3.1 批量LU分解块算法在GPU上的有效实现

文献[15]指出, 影响批量小矩阵在GPU上性能的问题的主要因素是对全局内存中数据的读写时间, 而不是计算时间, 因此, 在选择算法时, 应选择具有更少数据读写次数的算法. 由于Left-looking LU分解块算法比Right-looking LU分解块算法具有更少的对全局内存的数据读写次数, 因此, 在GPU上的实现批量LU分解时, 我们选取了Left-looking LU分解块算法. 算法8和算法9给出了这两种算法在GPU上的实现, 并通过算法8和算法9分析说明了Left-looking对全局内存的数据读取量少于Right-looking对全局内存的数据读取次数, 算法中下标 “S”的变量表示存储在共享内存中的变量.

如算法8描述, 在GPU上实现Right-looking LU分解块算法时, 由于对当前列条块中的数据需要进行多次读写, 则首先将列条块从全局内存读到共享内存(第3行), 然后对当前列条块进行部分选主元的LU分解(第7行), 分解完成后再将列条块写回到全局内存(第10行), 此时对全局内存的数据读写数据量为 $ 2 \times \left( {n - \left( {i - 1} \right) \times NB} \right) \times NB $ . 然后, 对除当前列条块外的两侧进行行交换(第12行), 这个过程数据读写次数为 $ 4 \times \left( {n \times NB - NB \times NB} \right) $ . 之后将列条块右边的矩阵 $ {\boldsymbol{U}} $ 的子块 $ {{\boldsymbol{U}}_{12}} $ 读到共享内存(第4行)并进行计算(第15行), 计算完成后再写回到全局内存(第17行), 此时数据读写次数为 $ 2 \times NB \times \left( {n - i \times NB} \right) $ , 最后用列条块的求解结果与 ${{\boldsymbol{U}}_{{\text{S, 12}}}}\left( {{{\boldsymbol{A}}_{{\text{S, 12}}}}} \right)$ 修正尾部矩阵 $ {{\boldsymbol{A}}_{22}} $ (第18行), 数据读写次数为 $ 2 \times \left( {n - i \times NB} \right) \times \left( {n - i \times NB} \right) $ . 则完成整个Right-looking LU分解块算法对全局内存的读写总次数约为: ${V_{\text{R}}} = \dfrac{2}{{3NB}}{n^3} + 5{n^2} - \dfrac{{11}}{3} \times NB \times n$ .

算法8. GPU上的Right-looking LU分解块算法.

1　for $ i = 1 $ to $ bn $

2　　/*将当前列条块读到共享内存*/

3　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{{\rm{S},11}}}} \\ {{{\boldsymbol{A}}_{{\rm{S},21}}}} \end{array}} \right] \leftarrow \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{11}}} \\ {{{\boldsymbol{A}}_{21}}} \end{array}} \right]$ ;

4　　 ${{\boldsymbol{A}}_{{\rm{S},12}}} \leftarrow {{\boldsymbol{A}}_{12}}$ ;

5　　__syncthreads();

6　　/*当前列条块的选主元LU分解*/

7　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{{\rm{S},11}}}} \\ {{{\boldsymbol{A}}_{{\rm{S},21}}}} \end{array}} \right] = {\boldsymbol{P}}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{L}}_{{\rm{S},11}}}} \\ {{{\boldsymbol{L}}_{{\rm{S},21}}}} \end{array}} \right] \times {{\boldsymbol{U}}_{{\rm{S, 11}}}}$ ;

8　　__syncthreads();

9　　/*将当前列条块写回到全局内存*/

10　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{11}}} \\ {{{\boldsymbol{A}}_{21}}} \end{array}} \right] \leftarrow \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{{\rm{S},11}}}} \\ {{{\boldsymbol{A}}_{{\rm{S},21}}}} \end{array}} \right]$ ;

11　　/*列条块外部分行交换*/

12　　 $ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{10}}}&{{{\boldsymbol{A}}_{12}}} \\ {{{\boldsymbol{A}}_{20}}}&{{{\boldsymbol{A}}_{22}}} \end{array}} \right] = {\boldsymbol{P}}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{10}}}&{{{\boldsymbol{A}}_{12}}} \\ {{{\boldsymbol{A}}_{20}}}&{{{\boldsymbol{A}}_{22}}} \end{array}} \right] $ ;

13　　__syncthreads();

14　　/*修正尾部矩阵*/

15　　 ${{\boldsymbol{A}}_{\rm{S},12}} = {\boldsymbol{L}}_{\rm{S},11}^{ - 1} \times {{\boldsymbol{A}}_{\rm{S},12}}$ ;

16　　__syncthreads();

17　　 ${{\boldsymbol{A}}_{12}} \leftarrow {{\boldsymbol{A}}_{{\rm{S},12}}}$ ; /*将 $ {{\boldsymbol{A}}_{12}} $ 写回到全局内存*/

18　　 ${{\boldsymbol{A}}_{22}} - = {{\boldsymbol{A}}_{{\rm{S},21}}} \times {{\boldsymbol{A}}_{{\rm{S}, 12}}}$ ;

19　　__syncthreads();

20　endfor

在GPU上实现如算法9所描述的Left-looking LU分解块算法时, 首先将当前列条块从全局内存读到共享内存(第3行), 然后用当前列条块左边已经求解出的 $ {{\boldsymbol{L}}_{00}} $ 来求解矩阵 $ {\boldsymbol{U}} $ 的子块 ${{\boldsymbol{U}}_{{\rm{S}, 01}}}\left( {{{\boldsymbol{A}}_{{\rm{S}, 01}}}} \right)$ (第6行), 然后修正列条块 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{{\rm{S, 11}}}}} \\ {{{\boldsymbol{A}}_{{\rm{S, 21}}}}} \end{array}} \right]$ (第8行), 此时对全局内存的数据读取量为 $ n \times NB + n \times \left( {i - 1} \right) \times NB $ . 修正完成后, 对列条进行部分选主元的LU分解(第11行), 分解完成后将列条块写回到全局内存(第14行), 此时读写数据量为 $ n \times NB $ . 最后对列条块外的部分进行行交换(第16行), 这个过程数据读写次数为 $ 4 \times \left( {n \times NB - NB \times NB} \right) $ . 则完成整个Left-looking LU分解块算法对全局内存的读写总量约为: ${V_{\rm{L}}} = \dfrac{1}{{2NB}}{n^3} + \dfrac{{11}}{2}{n^2} - 4 \times NB \times n$ .

算法9. GPU上的Left-looking LU分解块算法.

1　for $ i = 1 $ to $ bn $

2　　/*将当前列条块读到共享内存*/

3　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{{\rm{S, 01}}}}} \\ {{{\boldsymbol{A}}_{{\rm{S, 11}}}}} \\ {{{\boldsymbol{A}}_{{\rm{S, 21}}}}} \end{array}} \right] \leftarrow \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{01}}} \\ {{{\boldsymbol{A}}_{11}}} \\ {{{\boldsymbol{A}}_{21}}} \end{array}} \right]$ ;

4　　__syncthreads();

5　　/*修正当前列条块*/

6　　 ${{\boldsymbol{A}}_{{\text{S, 01}}}} = {\boldsymbol{L}}_{00}^{ - 1} \times {{\boldsymbol{A}}_{{\text{S, 01}}}}$ ;

7　　__syncthreads();

8　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{{\rm{S, 11}}}}} \\ {{{\boldsymbol{A}}_{{\rm{S, 21}}}}} \end{array}} \right] - = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{10}}} \\ {{{\boldsymbol{A}}_{20}}} \end{array}} \right] \times {{\boldsymbol{A}}_{{\rm{S, 01}}}}$ ;

9　　__syncthreads();

10　　/*当前列条块的选主元LU分解*/

11　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{{\rm{S, 11}}}}} \\ {{{\boldsymbol{A}}_{{\rm{S, 21}}}}} \end{array}} \right] = {\boldsymbol{P}}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{L}}_{{\rm{S, 11}}}}} \\ {{{\boldsymbol{L}}_{{\rm{S, 21}}}}} \end{array}} \right] \times {{\boldsymbol{U}}_{{\rm{S, 11}}}}$ ;

12　　__syncthreads();

13　　/*将当前列条块写回到全局内存*/

14　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{01}}} \\ {{{\boldsymbol{A}}_{11}}} \\ {{{\boldsymbol{A}}_{21}}} \end{array}} \right] \leftarrow \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{{\rm{S, 01}}}}} \\ {{{\boldsymbol{A}}_{{\rm{S, 11}}}}} \\ {{{\boldsymbol{A}}_{{\rm{S, 21}}}}} \end{array}} \right]$ ;

15　　/*列条块外部分行交换*/

16　　 $ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{10}}}&{{{\boldsymbol{A}}_{12}}} \\ {{{\boldsymbol{A}}_{20}}}&{{{\boldsymbol{A}}_{22}}} \end{array}} \right] = {\boldsymbol{P}}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_{10}}}&{{{\boldsymbol{A}}_{12}}} \\ {{{\boldsymbol{A}}_{20}}}&{{{\boldsymbol{A}}_{22}}} \end{array}} \right] $ ;

17　　__syncthreads();

18　endfor

在GPU上实现Left-looking LU分解块算法时, 除了将当前列条块读到共享内存, 计算完成后再将其写回到全局内存外, 另外有3个关键步骤: 修正列条块、对列条块进行部分选主元的LU分解和列条块外部分的行交换. 接下来我们将对这3个步骤在GPU上的实现设计与优化进行描述.

(1) 修正当前列条块

修正当前第 $ i $ 列条块(已存储于共享内存中)的每列数据需要列条块左边已经求解出的前 $ i-1 $ 个列条块(存储于全局内存)中的每列数据, 为了减少从读取效率低的全局内存读取数据量, 所有线程同时从当前列条块的左边读取一列数据到寄存器, 然后修正当前列条块的每列数据. 修正当前列条块所有列数据之后, 所有线程再同时读取下一列数据到寄存器, 继续修正当前列条块的所有列. 此设计有2个优点, 首先从全局内存读取数据没有Warp发散. 其次, 只对存储在全局内存中的当前列条块左边部分进行了一次读取就更新了当前列条块中的所有数据. 如图6所示.

(2) 对列条块进行部分选主元的LU分解

对列条块进行部分选主元的LU分解时, 需要找出绝对值最大的元素所在行, 为了让更多的线程参与搜索过程, 减少搜索时间, 我们采用线程访存连续的(coalesced)并行二叉树搜索算法进行, 每一个线程对应一列中两个元素的比较. 如图7所示, 假设当前列需要比较的元素个数为8, 则第1层共需要4个线程并行地对 ${{a}}_{1}$ 和 $ {{a}}_{5} $ (黑色)、 $ {{a}}_{2} $ 和 ${{a}}_{6}$ (蓝色)、 $ {{a}}_{3} $ 和 $ {{a}}_{7} $ (红色)、 $ {a_4} $ 和 $ {a_8} $ (橙色)找到最大值, 显然这4个线程对应的元素 ${a_1} - {a_4}$ 、 ${a_5} - {a_8}$ 是连续的(coalesced)的, 因此二叉树搜索算法的在这里效率很高. 然后第2层再对找到的4个最大值执行相同的过程, 直到找到本列的最大主元.

(3) 列条块外两侧的行交换

对列条外的部分进行行交换时, 通常的做法是NB次行交行逐次进行, 每个线程对应一列的2个元素交换, 这种方法虽然线程利用率高, 但是由于矩阵是以列主序存储在全局内存中, 则在行交换过程线程访问全局内存的地址是不连续的(non-coalescent), 同一个Warp中的线程对应的数据地址跨步大, 这会使得行交换过程耗时较长, 对性能造成影响. 为了降低线程访问全局内存地址不连续带来的影响, 我们提出了Warp分组行交换技术, 将线程按照Warp进行分组, 组内重新设计线程与矩阵元素的对应关系, 再将线程细分成多组, 每组包含NB个线程. 由于行交换过程至多进行NB次行交换, 因此NB个线程刚好对应NB次行交换. 注意为了避免线程同步影响性能, 对矩阵同一列进行行交换的NB个线程必须属于同一个Warp. 如图8所示, 假设块的大小NB=8, 线程按照Warp分成了红、蓝两组, 组内再将线程分成4组, 每组8个连续线程, 同时考虑8次行交换, 也就是将线程重新设计成一个[8, 4]的二维结构. 红蓝两组同时进行, 完成这两部分的行交换后, 往右推进, 直到所有的行交换完成.

此外, 由于列条块右边的部分, 随着算法的进行会被取到共享内存中, 因此可将这部分的行交换推迟到后边取到共享内存中时再进行, 加快这部分的行交换的时间, 如图9所示. 此时Left-looking LU分解块算法对全局内存的数据读取量减少 $ 2{n^2} - 2 \times NB \times n $ .

3.2 批量求逆算法在GPU上的有效实现

在GPU上实现LU分解之后的求逆算法与在CPU上实现一样, 也分为及时修正与延迟修正两种块算法, 但在GPU上实现时, 这两种方法对全局内存的数据读写次数不同. 因为影响批量小规模矩阵问题在GPU上实现性能的主要原因是对全局内存的数据读写时间, 由于延迟修正求逆块算法比及时修正求逆块算法具有更少的对全局内存的数据读写次数, 因此, 在GPU上的实现批量求逆算法时, 我们选取了延迟修正块算法.

算法10–算法13给出了第2.2节求逆过程中第(1)步和第(2)步及时修正与延迟修正的块算法在GPU上的实现, 并通过算法10–算法13分析说明了延迟修正求逆块算法对全局内存的数据读写次数少于及时修正求逆块算法对全局内存的数据读取量. 在条块内求解时, 由于矩阵 $ {\boldsymbol{W}} $ 中的当前列条块中的数据需要多次读写, 并且每个线程对应的每行数据的求解、修正相互独立, 因此可以将其从全局内存读到线程私有的寄存器中, 加快数据的读写速度, 求解完成后再写回到全局内存即可. 对于使用到的矩阵 $ {\boldsymbol{U}} $ 和矩阵L中子矩阵行/列条块, 由于对其进行了多次读取操作, 因此可以将其从全局内存取到共享内存中, 加快数据的读取速度. 在算法10–算法13中, 下标“R”表示存储在寄存器中的变量, 下标“S”表示存储在共享内存中的变量.

算法10给出了求解方程组 $ {\boldsymbol{XU}} = {\boldsymbol{W}} $ 及时修正块算法在GPU上的实现, 首先从全局内存中将当前列条块 $ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{01}}} \\ {{{\boldsymbol{W}}_{11}}} \end{array}} \right] $ 读到寄存器(第3行), 然后进行当前列条块内的计算(第9行), 计算完成后再将列条块写回到全局内存(第12行), 这时对全局内存的数据读写次数为 $ 2 \times \left( {i \times NB} \right) \times NB $ . 然后与 $ {{\boldsymbol{U}}_{12}} $ 来修正列条块右边的部分 $ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{02}}} \\ {{{\boldsymbol{W}}_{12}}} \end{array}} \right] $ (第14行), 此时, 对全局内存中矩阵 $ {\boldsymbol{W}} $ 的数据读写次数为 $ 2 \times \left( {i \times NB} \right) \times \left( {n - i \times NB} \right) $ . 整个修正尾部矩阵的过程共对矩阵 $ {\boldsymbol{U}} $ 的上三角部分进行了一次访问, 对全局内存中的矩阵 $ {\boldsymbol{U}} $ 的数据读取次数约 $ \dfrac{1}{2}{n^2} $ , 则整个及时修正的求 $ {{\boldsymbol{U}}^{ - 1}} $ 的块算法对全局内存的读写次数约为: $ {V_{{\text{T1}}}} = \dfrac{1}{{3NB}}{n^3} + \dfrac{3}{2}{n^2} + \dfrac{{2NB}}{3}n $ .

算法10. GPU内及时修正的求 $ {\boldsymbol{XU}} = {\boldsymbol{W}} $ 的块算法.

1　for $ i = 1 $ to $ bn $

2　　/*将当前列条块从全局内存读到寄存器*/

3　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{{\rm{R, 01}}}}} \\ {{{\boldsymbol{W}}_{{\rm{R, 11}}}}} \end{array}} \right] \leftarrow \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{{\rm{01}}}}} \\ {{{\boldsymbol{W}}_{{\rm{11}}}}} \end{array}} \right]$ ;

4　　__syncthreads();

5　　/*将 $ {{\boldsymbol{U}}_{{\text{11}}}} $ 和 $ {{\boldsymbol{U}}_{{\text{12}}}} $ 从全局内存读到共享内存*/

6　　 $\left[ {{{\boldsymbol{U}}_{{\rm{S, 11}}}}\begin{array}{*{20}{c}} {{{\boldsymbol{U}}_{{\rm{S, 12}}}}} \end{array}} \right] \leftarrow \left[ {{{\boldsymbol{U}}_{{\rm{11}}}}\begin{array}{*{20}{c}} {{{\boldsymbol{U}}_{{\rm{12}}}}} \end{array}} \right]$ ;

7　　__syncthreads();

8　　/*求解当前列条块*/

9　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{{\rm{R, 01}}}}} \\ {{{\boldsymbol{W}}_{{\rm{R, 11}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{{\rm{R, 01}}}}} \\ {{{\boldsymbol{W}}_{{\rm{R, 11}}}}} \end{array}} \right] \times {\boldsymbol{U}}_{{\rm{S, 1}}1}^{ - 1}$ ;

10　　__syncthreads();

11　　/*将当前列条块写回全局内存*/

12　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{{\rm{01}}}}} \\ {{{\boldsymbol{W}}_{{\rm{11}}}}} \end{array}} \right] \leftarrow \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{{\rm{R, 01}}}}} \\ {{{\boldsymbol{W}}_{{\rm{R, 11}}}}} \end{array}} \right]$ ;

13　　/*修正剩余尾部矩阵*/

14　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{{\rm{02}}}}} \\ {{{\boldsymbol{W}}_{{\rm{12}}}}} \end{array}} \right] - = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{{\rm{R, 01}}}}} \\ {{{\boldsymbol{W}}_{{\rm{R, 11}}}}} \end{array}} \right] \times {{\boldsymbol{U}}_{{\rm{S, 12}}}}$ ;

15　　__syncthreads();

16　endfor

算法11给出了求解方程组 $ {\boldsymbol{XU}} = {\boldsymbol{W}} $ 延迟修正块算法在GPU上的实现, 首先对寄存器中的当前列条块进行初始化(第3行), 这是由于 $ {\boldsymbol{W}} $ 初始时为单位矩阵, 所以这里不需要从全局内存中读入数据到寄存器. 然后用当前列条块左边已经求解出的部分 $ {{\boldsymbol{W}}_{00}} $ 与 $ {{\boldsymbol{U}}_{{\rm{S, 01}}}} $ 来修正当前列条块 $ {{\boldsymbol{W}}_{{\rm{R, 01}}}} $ (第9行), 此时, 对全局内存中矩阵 $ {\boldsymbol{W}} $ 的数据读取量约为 $\dfrac{1}{2} \times \left( {i - 1} \right) \times NB \times \left( {i - 1} \right) \times NB$ , 修正当前列条块的方法与Left-looking LU分解算法中修正列条块的方法类似, 所有线程同时从已经求解出的部分 $ {{\boldsymbol{W}}_{00}} $ 中读取一列数据到寄存器, 然后与共享内存中的 $ {{\boldsymbol{U}}_{{\rm{S, 01}}}} $ 修正存储在寄存器中的当前列条块的每列数据, 修正当前列条块的所有列数据之后, 所有线程再同时读取下一列数据到寄存器, 继续修正当前列条块的所有列. 修正完成后, 然后求解当前列条块 (第12行), 完成后写回到全局内存(第15行), 此时数据写量为 $ \left( {i \times NB} \right) \times NB $ . 另外, 整个修正过程共对全局内存中的矩阵 $ {\boldsymbol{U}} $ 进行了一次访问, 对全局内存中矩阵 $ {\boldsymbol{U}} $ 的数据读取量约 $ \dfrac{1}{2}{n^2} $ . 则整个延迟修正的求 $ {{\boldsymbol{U}}^{ - 1}} $ 的块算法对全局内存的读写量约为: $ {V_{{\rm{D1}}}} = \dfrac{1}{{6NB}}{n^3} + \dfrac{3}{4}{n^2} + \dfrac{{7NB}}{{12}}n $ .

算法11. GPU内延迟修正的求 $ {\boldsymbol{XU}} = {\boldsymbol{W}} $ 的块算法.

1　for $ i = 1 $ to $ bn $

2　　/*初始化寄存器中的当前列条块*/

3　　 ${{\boldsymbol{W}}_{{\rm{R, 01}}}} = 0, {{\boldsymbol{W}}_{{\rm{R, 11}}}} = {\boldsymbol{I}}$ ;

4　　__syncthreads();

5　　/*将 $ {{\boldsymbol{U}}_{{\text{01}}}} $ 和 $ {{\boldsymbol{U}}_{{\text{11}}}} $ 读到共享内存*/

6　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{U}}_{{\rm{S, 01}}}}} \\ {{{\boldsymbol{U}}_{{\rm{S, 11}}}}} \end{array}} \right] \leftarrow \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{U}}_{{\rm{01}}}}} \\ {{{\boldsymbol{U}}_{{\rm{11}}}}} \end{array}} \right]$ ;

7　　__syncthreads();

8　　/*修当前列条块*/

9　　 ${{\boldsymbol{W}}_{{\rm{R, 01}}}} - = {{\boldsymbol{W}}_{{\rm{00}}}} \times {{\boldsymbol{U}}_{{\rm{S, 01}}}}$ ;

10　　__syncthreads();

11　　/*求解当前列条块*/

12　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{{\rm{R, 01}}}}} \\ {{{\boldsymbol{W}}_{{\rm{R, 11}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{{\rm{R, 01}}}}} \\ {{{\boldsymbol{W}}_{{\rm{R, 11}}}}} \end{array}} \right] \times {\boldsymbol{U}}_{{\rm{S, 1}}1}^{ - 1}$ ;

13　　__syncthreads();

14　　/*将当前列条块写回全局内存*/

15　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{{\rm{01}}}}} \\ {{{\boldsymbol{W}}_{{\rm{11}}}}} \end{array}} \right] \leftarrow \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{W}}_{{\rm{R, 01}}}}} \\ {{{\boldsymbol{W}}_{{\rm{R, 11}}}}} \end{array}} \right]$ ;

16　　__syncthreads();

17　endfor

求逆过程的第(2)步求解方程组 $ {\boldsymbol{XL}} = {\boldsymbol{W}} $ 的过程在GPU上的实现与第(1)步的在GPU上的实现类似, 也分为及时修正与延迟修正2种块算法, 算法12和算法13分别给出了这2种求解方程组 $ {\boldsymbol{XL}} = {\boldsymbol{W}} $ 的块算法在GPU上的实现. 其中, 及时修正算法对全局内存的读写量约为: $ {{\boldsymbol{V}}_{{\text{T2}}}} = \dfrac{1}{{NB}}{n^3} + \dfrac{3}{2}{n^2} $ . 延迟修正算法对全局内存的读写量为: $ {V_{{\text{D2}}}} = \dfrac{1}{{2NB}}{n^3} + 2{n^2} $ .

算法12. GPU内及时修正的求 $ {\boldsymbol{XL}} = {\boldsymbol{W}} $ 的块算法.

1　for $ i = 1 $ to $ bn $

2　　/*将当前列条块从全局内存取到寄存器*/

3　　 ${{\boldsymbol{W}}_{{\rm{R, 1}}}} \leftarrow {{\boldsymbol{W}}_1}$ ;

4　　__syncthreads();

5　　/*将 $ {{\boldsymbol{L}}_{10}} $ 和 $ {{\boldsymbol{L}}_{11}} $ 从全局内存读到共享内存*/

6　　 $\left[ {{{\boldsymbol{L}}_{{\rm{S, 10}}}}\begin{array}{*{20}{c}} {{{\boldsymbol{L}}_{{\rm{S, 11}}}}} \end{array}} \right] \leftarrow \left[ {{{\boldsymbol{L}}_{{\rm{10}}}}\begin{array}{*{20}{c}} {{{\boldsymbol{L}}_{{\rm{11}}}}} \end{array}} \right]$ ;

7　　__syncthreads();

8　　/*求解列条块*/

9　　 ${{\boldsymbol{W}}_{{\rm{R}, 1}}} = {{\boldsymbol{W}}_{{\rm{R, 1}}}} \times {\boldsymbol{L}}_{{\rm{S, 11}}}^{ - 1}$ ;

10　　__syncthreads();

11　　/*将列条块的求解结果写回到全局内存*/

12　　 ${{\boldsymbol{W}}_1} \leftarrow {{\boldsymbol{W}}_{{\rm{R, 1}}}}$ ;

13　　/*修正剩余尾部矩阵*/

14　　 ${{\boldsymbol{W}}_0} - = {{\boldsymbol{W}}_{{\rm{R, 1}}}} \times {{\boldsymbol{L}}_{{\rm{S, 10}}}}$ ;

15　endfor

算法13. GPU内延迟修正的求 $ {\boldsymbol{XL}} = {\boldsymbol{W}} $ 的块算法.

1　for $ i = 1 $ to $ bn $

2　　/*将当前列条块从全局内存取到寄存器*/

3　　 ${{\boldsymbol{W}}_{{\rm{R, 1}}}} \leftarrow {{\boldsymbol{W}}_1}$ ;

4　　__syncthreads();

5　　/*将 $ {{\boldsymbol{L}}_{11}} $ 和 $ {{\boldsymbol{L}}_{21}} $ 读到共享内存*/

6　　 $\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{L}}_{{\rm{S, 11}}}}} \\ {{{\boldsymbol{L}}_{{\rm{S, 21}}}}} \end{array}} \right] \leftarrow \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{L}}_{{\rm{11}}}}} \\ {{{\boldsymbol{L}}_{{\rm{21}}}}} \end{array}} \right]$ ;

7　　__syncthreads();

8　　/*修正列条块*/

9　　 ${{\boldsymbol{W}}_{{\rm{R, 1}}}} - = {{\boldsymbol{W}}_2} \times {{\boldsymbol{L}}_{{\rm{S, 21}}}}$ ;

10　　__syncthreads();

11　　/*求解列条块*/

12　　 ${{\boldsymbol{W}}_{{\rm{R, 1}}}} = {{\boldsymbol{W}}_{{\rm{R, 1}}}} \times {\boldsymbol{L}}_{{\rm{S, 11}}}^{ - 1}$ ;

13　　__syncthreads();

14　　/*将列条块的求解结果写回到全局内存*/

15　　 ${{\boldsymbol{W}}_1} \leftarrow {{\boldsymbol{W}}_{{\rm{R, 1}}}}$ ;

16　　__syncthreads();

17　endfor

显然, 求逆的延迟修正块算法对全局内存的数据读写次数小于及时修正块算法: $ {V_{{\text{D1}}}} + {V_{{\text{D2}}}} \lt {V_{{\text{T1}}}} + {V_{{\text{T2}}}} $ . 我们对采用这2种修正方法的批量求逆算法在GPU上进行了实现, 实验结果详见第4.4节.

在GPU上实现的批量算法中, 一个Block对应一个矩阵, Block中一个线程对应矩阵的一行. 若上述求逆步骤在一个Kernel中完成, 为了保证算法的正确性, 则申请的Block中线程的个数只能是 $ n $ , 共享内存的大小只能为最大的 $ n \times NB $ . 我们称这种一次性分配的方法为静态资源分配法. 如算法11所描述的GPU内延迟修正的求 $ {{\boldsymbol{U}}^{ - 1}} $ 块算法中, 所需的线程数目和共享内存的大小与当前列条块的大小有关, 而列条块的大小是动态变化的, 每一次列条块的求解需要的线程数目为 $ i \times NB $ , 需要的共享内存的大小为 $ \left( {i \times NB} \right) \times NB $ . 若采用静态资源分配法, 显然会造成线程的闲置和共享内存资源的浪费. 因此我们设计了运行时按需申请线程和共享内存资源的动态资源分配方法, 每一个列条块求解启动一次Kernel, 根据列条块的大小申请线程和共享内存, 这种方法避免了GPU资源的浪费, 进一步提升批量求逆算法的性能. 实验结果详见第4.4节.

3.3 减少访存数据量的列交换优化方法

以上过程结束后, 根据P交换相应的列即可得到 $ {{\boldsymbol{A}}^{ - 1}} $ . 在GPU的实现中, 若要进行两列间交换, 首先需要从全局内存中将这两列读出, 然后将交换后两列写回到全局内存中对应的列, 则列交换需要访问全局内存的数据量为 $ 4{n^2} $ . 为了进一步降低全局内存的访问量, 我们设计了减少访存数据量的列交换优化方法, 该方法将上述求解的结果直接存储到原始矩阵 $ {\boldsymbol{A}} $ , 而不是存储在 $ {{\boldsymbol{A}}^{ - 1}} $ 中. 则可根据 $ {\boldsymbol{P}} $ 直接将 $ {\boldsymbol{A}} $ 中对应的列从全局内存中读出, 然后写入全局内存中的逆矩阵 $ {{\boldsymbol{A}}^{ - 1}} $ 中去, 这样避免了显式的矩阵列交换, 对全局内存的访问数据量减少至 $ 2{n^2} $ .

4 实验分析

为了测试算法的性能和准确性, 我们在TITAN V GPU上实现并测试了本文描述的批量LU分解与批量求逆算法, TITAN V GPU的显存容量为12 GB, 峰值双精度浮点性能7.5 TFLOPS. 我们对10000个规模在33–190之间的随机矩阵进行了测试. 并在同一平台上用CUBLAS库和MAGAMA库中的算法进行了性能测试和对比分析. 测试中使用的CUBLAS版本是9.1.85, MAGMA版本是2.5.1.

4.1 块大小对算法实现性能的影响

(1) 块大小对批量LU分解块算法实现性能的影响

批量LU分解块算法中块的大小NB将会影响其在GPU上的实现性能. 这是因为, 块的大小NB决定着对全局内存的访问数据次数、每个列条块使用共享内存的大小以及条块外Warp分组行交换中线程的组织形式等. 因此, 在使用时需要衡量块的大小, 使程序达到最好的性能. 图10描述的是当矩阵类型为单精度复数、块大小NB=8, 10, 14, 18时批量LU分解块算法的浮点计算性能曲线示意图. 从图10中可以看出, 开始时随着NB的增大, 批量LU分解块算法的性能越来越好, 当NB=14时达到的性能相对较好, 当NB再增加时性能反而下降, 这是因为随着NB的增大, 程序所需的共享内存等GPU资源越多, 会造成资源的竞争, 程序的性能反而降低了.

(2) 块大小对批量求逆块算法实现性能的影响

不同的块的大小NB也会影响批量求逆算法的实现性能. 这是因为, 在批量求逆块算法中块的大小NB决定着对全局内存的访问数据量、每个线程使用寄存器的数量, 以及每个列条块使用共享内存的大小. 通常来说使用越多的寄存器和共享内存越能加快数据的读写速度, 但是GPU内的资源是有限的, 若使用过多的寄存器和共享内存会造成资源的竞争, 反而会降低算法实现的性能. 因此, 在使用时需要衡量块的大小, 使算法实现达到最好的性能. 图11描述的是当矩阵类型为单精度复数、块大小NB=8, 14, 20, 24, 28时批量求逆块算法的浮点计算性能曲线示意图. 从图11中可以看出, 同批量LU分解块算法类似, 批量求逆块算法的性能开始时随着NB的增大性能不断提升, 当NB=16–20时取得的性能相对较好, 当NB再增加时性能反而下降了.

4.2 本文批量算法与CUBLAS和MAGMA中算法实现的浮点计算性能比较

(1) 批量LU分解块算法与CUBLAS和MAGMA中实现的浮点计算性能比较

图12给出了矩阵类型为单精度复数(float2)、双精度复数(double2)、单精度实数(float)和双精度实数(double)时, 本文实现的批量LU算法(OUR_LU)与CUBLAS (CUBLAS_LU)和MAGMA中实现(MAGMA_LU)的批量LU算法的浮点计算性能曲线示意图. 从图12中可以看出, 不同矩阵类型的测试结果MAGMA_LU的性能高于CUBLAS_LU的性能, 而OUR_LU的性能又高于MAGMA_LU的性能. 此外, CUBLAS_LU的浮点计算性能随着矩阵规模的增加比较平稳. 而OUR_LU与MAGMA_LU的浮点计算性能随着矩阵规模的增大不断提升, 不同矩阵类型的OUR_LU的浮点计算性能分别可达约2 TFLOPS (float2)、1.2 TFLOPS (double2)、1 TFLOPS (float)、0.67 TFLOPS (double). 图13描述了本文OUR_LU算法与CUBLAS_LU和MAGMA_LU的不同矩阵类型的加速比曲线示意图. 从图13中可以看出, 开始时随着矩阵规模的增大, OUR_LU与CUBLAS_LU的加速比越大, 不同矩阵类型分别最高达到约9× (float2)、8× (double2)、12× (float)、13× (double), 当矩阵规模增大到一定程度后, 加速比趋于平稳, 分别约在7×–8× (float2)、5×–6× (double2)、10×–12× (float)、9×–12× (double)之间. OUR_LU与MAGMA_LU的加速比相对平稳, 不同矩阵类型分别约在1.2×–2.5× (float2)、1.2×–3.2×(double2)、1.1×–3× (float)、1.1×–2.7× (double)之间.

(2) 批量求逆块算法与CUBLAS和MAGMA中实现的浮点计算性能比较

图14给出了矩阵类型为单精度复数(float2)、双精度复数(double2)、单精度实数(float)和双精度实数(double)时本文实现的批量求逆算法(OUR_INV)与CUBLAS (CUBLAS_INV)和MAGMA中实现(MAGMA_INV)的批量求逆算法的浮点计算性能曲线示意图. 从图14中可以看出, 当矩阵规模较小时, 不同矩阵类型的测试结果中, CUBLAS_INV的性能要高于MAGMA_INV的性能, 随着矩阵规模的增大, MAGMA_INV的性能反超CUBLAS_INV的性能. 本文的OUR_INV的浮点计算性能一直高于CUBLAS_INV和MAGMA_INV的性能. 此外, CUBLAS_INV的性能比较平稳, OUR_INV与MAGMA_INV的浮点计算性能随着矩阵规模的增大不断提升. 不同矩阵类型的OUR_INV的浮点计算性能大约分别可达到4 TFLOPS (float2)、2 TFLOPS (double2)、2.2 TFLOPS (float)、1.2 TFLOPS (double). 后文图15显示了OUR_INV与CUBLAS_INV和MAGMA_INV的加速比曲线示意图. 从图15中可以看出, 开始时随着矩阵规模越大, OUR_INV与CUBLAS_INV的加速比越大, 最高分别达到约5× (float2)、4× (double2)、7× (float)、7× (double). 当矩阵规模达到一定程度后, 加速比趋于平稳, 分别约在4×–5× (float2)、3×–4× (double2)、5×–7× (float)、5×–6× (double)之间. 而OUR_INV与MAGMA_INV的加速比, 当矩阵规模较小时比较大, 随着矩阵规模的增加趋于稳定, 不同矩阵类型分别约在2×–3× (float2)、2×–3× (double2)、2.8×–3.4× (float)、1.6×–2× (double)之间.

4.3 批量LU分解算法优化分析

图16显示了当矩阵类型为单精度复数、块大小NB=14时, 批量LU分解块算法中列条块外的两侧进行行交换采用Warp分组行交换(OUR_LU_v2)以及行推迟优化技术(OUR_LU_v3)与优化前(OUR_LU_v1)的浮点计算性能曲线示意图. 从图16中可以看出, 对列条块外的两侧进行行交换进行优化后, OUR_LU_v3的浮点计算性能相对于OUR_LU_v1的性能得到2–5倍的提升, 这是因为列条块外两侧的行交换在批量LU分解中占据了60%以上的时间^[13], 因此, 对列条块外两侧的行交换进行优化, 将会明显提升批量LU算法的性能. 从图16还可以看出随着矩阵规模的增大提升幅度同样在增大, 这是因为随着矩阵规模的增大, 列条块外部分行交换对全局内存访问的数据量越多, 对列条块外在全局内存部分的行交换进行优化的效果越明显.

4.4 批量求逆算法优化分析

(1) 及时修正与延迟修正两种算法实现的浮点计算性能比较

图17显示了当矩阵类型为单精度复数、块的大小NB=8, 16时, 采用及时修正(timely)与延迟修正(delay)两种批量矩阵求逆算法方法在GPU内实现的浮点计算性能曲线示意图. 从图17中可以看出, 延迟修正算法的浮点计算性能优于及时修正算法的性能. 这是因为在GPU内实现及时修正的算法对于全局内存数据的读写数量多于延迟修正的算法. 并且NB越小, 提升越明显, 这是因为NB越小对全局内存的数据读写次数越多, 因此延迟修正算法的性能相较于及时修正算法的性能提升越明显. 这里及时修正与延迟修正算法使用的都是采用一次性分配资源的静态法.

(2) 静态资源分配法与动态资源分配方法实现的浮点计算性能比较

图18显示的是当矩阵类型为单精度复数、块的大小NB=8, 16时, 延迟修正算法采用一次性分配资源的静态资源分配方法(static)与运行时按需分配资源的动态资源分配方法(dynamic)的浮点计算性能曲线示意图. 从图18中可以看出, NB取不同的值时动态资源分配方法的浮点计算性能相较于静态资源分配方法都有所提升. 并且随着矩阵规模的增大, 提升幅度也在不断增大, 最大约30%. 这是因为随着矩阵规模的增大, 采用静态资源分配方法时线程闲置的情况和共享内存资源的浪费更加严重, 这时动态资源分配方法的优势就更加地显示出来, 因此提升的幅度越大.

5 总　结

本文给出了批量矩阵的LU分解和求逆算法及其在GPU上实现及优化方法. 针对批量LU分解问题, 分析了Left-looking和Right-looking等常用LU分解块算法在GPU上实现时对全局内存的数据读写次数, 由于Left-looking算法具有更少访存次数, 能更好地适应GPU架构下批量小矩阵的LU分解, 因此我们选择了具有较少访存的Left-looking块算法. 在LU分解的选主元过程, 采用了适合GPU架构的并行二叉树搜索算法. 此外, 为了降低选主元引起的行交换过程对算法性能的影响, 本文提出了Warp分组行交换和行交换延迟两个优化技术, 大大提升了算法的性能. 针对LU分解后的批量求逆问题, 分析了矩阵求逆过程中修正方法, 为了减少修正过程对全局内存的数据读写次数, 在批量求逆的GPU实现中采用了延迟修正的矩阵求逆块算法. 同时, 为了加快对全局内存的数据读写速度, 采用了更多利用寄存器和共享内存的优化方法和减少访存数据量的列交换优化方法. 另外, 为了避免线程的闲置和共享内存等GPU资源浪费, 提出了运行时按需申请GPU资源的动态资源分配方法. 最终, 在TITAN V GPU上, 我们对10000个规模在33–190之间的随机矩阵进行了测试. 当矩阵类型为单精度复数、双精度复数、单精度实数和双精度实数时, 我们实现的批量LU分解算法的浮点计算性能分别可达到约2 TFLOPS、1.2 TFLOPS、1 TFLOPS、0.67 TFLOPS, 与CUBLAS中的实现相比加速比分别最高达到了约9×、8×、12×、13×, 与MAGMA中的实现相比加速比分别达到了约1.2×–2.5×、1.2×–3.2×、1.1×–3×、1.1×–2.7×. 不同矩阵类型的批量求逆算法的浮点计算性能分别可达到约4 TFLOPS、2 TFLOPS、2.2 TFLOPS、1.2 TFLOPS, 与CUBLAS中的实现相比加速比分别最高达到了约5×、4×、7×、7×, 与MAGMA中的实现相比加速比分别达到了约2×–3×、2×–3×、2.8×–3.4×、1.6×–2×. 此外, 不同的条块的大小将会影响程序的性能, 在使用时需要权衡. 本文批量LU分解算法与批量求逆算法在GPU上的实现与优化方法同样适用于Cholesky分解和QR分解等常见的矩阵分解算法在GPU上的批量实现.