[WIP]add mluop cholesky #1018
[WIP]add mluop cholesky #1018
Conversation
dglr
changed the title
|
添加验收计划： |
kernels/cholesky/cholesky_union1.mlu
Outdated
| factor=sqrt(diag[iter*POTF_NB+iter]); | ||
| factor = 1.0/factor; | ||
| for(int i = 0; i < span; i++) | ||
| { |
There was a problem hiding this comment.
这个for循环可以用bangc的向量化指令替换下
There was a problem hiding this comment.
这里的逻辑是将一列数据乘以factor，但是矩阵以行主序，所以在内存中数据不连续，无法直接使用bangc指令替换
There was a problem hiding this comment.
片上可以transpose成连续然后simd运算再transpose回来
kernels/cholesky/cholesky_union1.mlu
Outdated
| { | ||
| for(int h = 0; h < k; h++) | ||
| { | ||
| rC[i*span_b+j] += rA[i*NB+h] * rB[j*NB+h]; |
kernels/cholesky/cholesky_union1.mlu
Outdated
| { | ||
| if(j < i) | ||
| continue; | ||
| A[j * lda + i ] -= A[i*lda+iter] * A[j * lda + iter]; |
There was a problem hiding this comment.
可以用bang_fusion的FMA向量化指令替换
|  | ||
| 图7 最后一步划分 | ||
|
|
||
| 每个列块，仍然需要先计算该列块的外部依赖（该列块左侧的所有列块），然后对列块中的每一列分别计算内部依赖，对于这两个部分可以分别用两个kernel来实现。由于这一步骤是严重的串行瓶颈，因此在划分小块时需要尽量让计算的快更小，减少串行瓶颈对性能的影响 |
There was a problem hiding this comment.
算法中的TRSM，POTRF, SYRK等子函数缺少MLU上的具体的拆分逻辑，片上空间使用和MLU上具体的实现过程，伪代码，看完这个方案还不是不太明确MLU上是具体怎么实现这个算法的
|
当前进度：complex float类型完成64*64以下的规模，现在正在编写测试大规模的复数矩阵乘以及大规模TRSM算子，预计下周二能够完成 |
| test_param: { | ||
| error_func: DIFF1 | ||
| error_func: DIFF2 | ||
| error_threshold: 0.003 |
There was a problem hiding this comment.
补充 docs/bangc-docs/user_guide/9_operators/index.rst 算子说明，可参考 算子涉及文档
可参考 https://github.com/Cambricon/mlu-ops/pull/662/files#diff-7f0a558d8f985a4ebd89cd6674a4bf1a91549ddcc6e708a897f351cb2006f0e8
|
|
||
| ### 1.2 Cholesky分解 | ||
|
|
||
| 对正定厄密特矩阵$`A`$进行Cholesky分解，即求矩阵$`L`$使下式成立： |
There was a problem hiding this comment.
PR 冲突了，建议本次rebase 到最新的 cambricon/master 后，再push
|
添加测试文档： 本测试报告模板是希望能帮助算子开发者在完成算子开发后进行有效充分的自检，开发出功能、性能都满足要求的高质量算子。 1. Cholesky测试报告添加算子描述
1.1 精度验收标准输出结果的动态阈值 diff1， diff2，diff3_2 精度验收通过 1.2 算子方案CHECKLIST
1.3 新特性测试
1.4 参数检查提交新算子时，给出测试点，并说明测试结果。
2. 功能测试对于 New Feature Test 部分中使用的案例，此处记录了特征、案例数量和结果。当测试多个操作时，需要多个表来包含这些操作的详细信息。
3. 性能测试float类型单batch性能测试如下，表格中数字为运行时间，单位为微秒（us），最右侧一列为mlu的运行时间与pytorch在gpu上的运行时间的比值：
float类型多batch性能测试：
float类型的cholesky分解在mlu端运行时间在pytorch运行时间的10倍以内。
complex类型多batch性能测试：
对于mlu/pytorch>10的规模，例如batch为32，N为128时，使用cnperf-cli进行性能分析，如下图所示 4. 总结分析实现了Cholesky分解的功能，与pyTorch一致。大部分规模的性能测试时间小于v100中运行时间的10倍，小部分受限于复数矩阵乘法的性能。 |
|
性能验收标准上要求测试规模及格线是v100的10倍，性能不足部分还需要分析下原因做进一步优化。 |
|
函数覆盖率我看有一个没有被覆盖到，可以分析下，要做到100%的函数覆盖吧 分支覆盖率的数据也可以贴下哈 |
原因已经分析出结果，受限于没有原生的复数乘法，需要使用4次实数乘法来实现复数乘法，造成了性能瓶颈，使得矩阵乘法的耗时已经超过了v100运行时间的10倍。在测试报告中可以看到对应的截图和分析 |
已经修正，所有函数均被覆盖，修正后的代码已上传。分支覆盖率的数据是指什么呢，似乎没看到这个数据在哪里 |
可以把覆盖率生成的目录拷贝到本地，打开index.html可以看到 |
|
另外pytorch上的性能数据这个是怎么统计的呢？是用ncu吗？ |
使用torch.cuda.Event，torch.cuda.synchronize()和start.record()进行gpu上的时间记录 |
|
|
|
||
| for (int iter = 0; iter < k; iter += POTF_NB) { | ||
| __bang_move(rA, rp, POTF_NB * span * sizeof(float)); | ||
| __memcpy(rB, sB, POTF_NB * POTF_NB * sizeof(float), SRAM2NRAM); |
There was a problem hiding this comment.
可以使用异步接口__memcpy_asyn ? 后面也有很多同步memcpy看能否改成异步
There was a problem hiding this comment.
这里在memcpy之后后续会立即用到，其他地方也是memcpy之后会立即用到 同步改成异步没有显著收益
| if (if_execute) { | ||
| for (int i = iter + 1; i < iter_num; i++) { | ||
| for (int j = finish; j < finish + span; j++) { | ||
| if (j < i) continue; |
There was a problem hiding this comment.
设计文档中图3是gemm的计算。如果width按照32切分，则gemm可能的规模为[H, 64] * [64, 32] = [H, 32].这个规模使用bang_conv，性能应该远好于for循环计算。
这个gemm的计算是使用的for循环计算还是bang_conv呢？
| __sync(); | ||
| for (int i = iter + 1; i < width; i++) { | ||
| for (int j = 0; j < m; j++) { | ||
| dst[j * width + i] -= dst[i * width + iter] * dst[j * width + iter]; |
There was a problem hiding this comment.
__bang_transpose可以实现nram上的transpose。
另外，也可以使用memcpy_async或者bang_move搬数，将dst摆成片上连续的。计算完成后将数据按照对应的stride拷贝回去
| func_type = CNRT_FUNC_TYPE_UNION8; | ||
| carry_batch = batch < 8 ? 8 : batch; | ||
| } | ||
| dim.x = carry_batch * 4; |
There was a problem hiding this comment.
batch==31的时候，dim.x = 124 ? 有测试这种场景吗？
|
|
||
| temp_b = 0; | ||
| for (int j = 0; j < m - 1; j++) { | ||
| temp_b += rC[i * calc_length + j]; |
There was a problem hiding this comment.
此处float数不足32 ///// 数据个数不足32，将后面多余的数设置为0， 然后按照32个数进行计算
用了4次矩阵乘法这个也可以优化的，可以做成一个kernel，然后把4个输入都拷贝到片上做计算，这样可以节省2倍的io，当前每个输入的io会重复加载两次 |
|
|
| @@ -3835,6 +3835,10 @@ mluOpDynamicPointToVoxelForward(const mluOpHandle_t handle, | |||
| /*! | |||
| * @brief Gets extra space size that is needed in the GenerateProposalsV2 operation. | |||
| * | |||
| * @par Deprecated | |||
| * - ::mluOpGetGenerateProposalsV2WorkspaceSize is deprecated and will be removed in the future | |||
| * - None. | ||
| */ | ||
| mluOpStatus_t MLUOP_WIN_API | ||
| mluOpGetGenerateProposalsV2WorkspaceSize_v2(mluOpHandle_t handle, |
| MLUOP_CHECK(mluOpGetSizeOfDataType(dtype, &type_size)); | ||
| total_size = type_size * size_a * lda * ((uint64_t)batch_size); | ||
| PARAM_CHECK("mluOpCholesky", total_size < size_limit); | ||
| if (type_size == 8 && batch_size > 16 && size_a > 2000) { |
There was a problem hiding this comment.
这里的8,16,2000建议修改成有含义的变量
另外这里建议增加注释说明为啥有两个分支？
| calculate_body(handle, 16, input_desc, d_input, output_desc, d_output, | ||
| upper, (float*)workspace); | ||
| cnrtQueueSync(queue); | ||
| calculate_body(handle, ((uint64_t)batch_size) - 16, input_desc, |
There was a problem hiding this comment.
16 magic number建议修改成有含义的变量，提升可读性
| #ifndef __CHOLESKY_H | ||
| #define __CHOLESKY_H | ||
|
|
||
| #define DEBUG |
|
|
||
| #define CNB (32) | ||
| #define REC_NB (16) | ||
| #define POTF_NB ((REC_NB) / 4) |
There was a problem hiding this comment.
4这里是magic number，建议修改成有意义的变量
| } else if (batch <= 4) { | ||
| carry_batch = 4; | ||
| } else if (batch <= 8) { | ||
| carry_batch = 8; |
| cnrtFunctionType_t func_type = CNRT_FUNC_TYPE_UNION1; | ||
| dim.y = 1; | ||
| dim.z = 1; | ||
| if (batch < 8) { |
There was a problem hiding this comment.
这个逻辑可以根据板卡的cluster数来设置，当前写死了只能适用于8 cluster的板卡
MLU Cholesky 分解实现方案
1 Cholesky分解算法介绍
Cholesky分解是科学和数值领域中最重要的算法之一。Cholesky算法是指将一个厄密特矩阵分解成一个下三角矩阵与其共轭转置之乘积，这种分解方式可以提高代数运算效率。
1.1 厄密特矩阵
厄密特矩阵，又称自伴随矩阵，是共轭对称的方阵。厄密特矩阵中对角线元素均为实数，且每个第i行j列的元素都与第j行i列的元素互为共轭转置。例如：
对于一个矩阵$A$，如果其是厄密特矩阵，则可以对其进行Cholesky分解，如果其是正定矩阵（对于所有的非零实数$x$ ，都有$x^TAx>0$ ）则Cholesky分解的结果唯一，否则结果不唯一。
1.2 Cholesky分解
对正定厄密特矩阵$A$进行Cholesky分解，即求矩阵$L$使下式成立：
其中，$L$ 是一个下三角矩阵且对角元素均为正实数，$L^*$ 表示$L$的共轭转置，是一个上三角矩阵。当$A$ 是一个实数矩阵时，Cholesky分解可以改写为
下文中为表述方便，所有矩阵$A$均为实数矩阵。
对于一个$n\times n$ 的实矩阵$A$ ，Cholesky分解可以被写作如下过程：
根据上式不难看出，每个$a_{i,j}$ 等于由$l_{i,j}$ 和$L$ 矩阵的其它元素组成的多项式，例如$a_{32}=l_{21}l_{31}+l_{32}l_{22}$ ，并且多项式中只有一个项包含了$l_{i,j}$ （$a_{32}$ 等价的多项式中只有$l_{22}l_{32}$ 这一项），包含了$l_{i,j}$ 的项另一个因子都为对角线元素，因此为了计算$l_{i,j}$ ，可以由$a_{i,j}$ 减去不包含$l_{i,j}$ 的其它项然后除以对角线元素，这样就能算出每个$l_{i,j}$ 。
2 Cholesky分解实现
将输入矩阵进行分块，然后使用以下流程计算Cholesky分解：
上图中，假设矩阵$L$的左边两列块已经计算完毕（黄色部分的非对角元和红色的对角元），这个流程展示了计算中间列块的过程（蓝色部分和橙色部分），完整的Cholesky计算只需要对分块后的所有列重复执行此流程。
SYRK(HERK)、GEMM和TRSM均为标准BLAS库中的操作，POTRF为计算对角块（完整矩阵的对角元素所在的块）内部依赖的kernel。下面将按照计算顺序依次介绍。
2.1 SYRK(HERK)
SYRK是BLAS的标准操作（数据类型是复数时为HERK），定义为：
其中$C$ 为$n\times n$ 的方阵，$A$ 为$n\times m$ 的矩阵，$\alpha$ 和$\beta$ 是标量。
此处使用SYRK是为了计算橙色块的外部依赖，上式中的$C$ 代表橙色对角块（完整矩阵的对角元素所在的块），$A$ 代表橙色块左侧的所有黄色块，$\alpha$ 、$\beta$ 分别取-1和1。
2.2 GEMM
GEMM是BLAS的标准操作，定义为：
其中$C$ ，$A$ ，$B$ 分别是$m\times n$ ，$m\times k$ ，$k\times n$ 的矩阵，$\alpha$ 和$\beta$ 是标量。
这里使用GEMM计算蓝色非对角块的外部依赖，上式的$C$ 代表蓝色块，$A$ 和$B$ 分别代表橙色块左侧的黄色块和蓝色块左侧的黄色块。$\alpha$ 和$\beta$ 分别为-1和1。
2.3 TRSM
TRSM是BLAS的标准函数，定义为：
已知下三角矩阵$A$ 和矩阵$B$ ，TRSM解出矩阵$X$ ，$A$ 为$n\times n$ 方阵，$X$ 和$B$ 为$m\times n$ 的矩阵。
对角块在SYRK后需要经过POTRF完成后续计算，这里假设已经计算完毕，于是可以通过TRSM完成蓝色块的剩余计算，TRSM执行后蓝色部分计算完毕。上式中$A$ 为红色块，$X$ 和$B$ 均为蓝色块，计算结果覆盖原矩阵。
2.4 POTRF
POTRF这个函数名取自LAPACK中Cholesky分解的函数，POTRF的目的是计算橙色对角块的所有依赖，POTRF执行后对角块中的所有元素计算完毕。
对于POTRF计算的块边长的典型取值为512，这仍然是一个较大的规模，为了进一步分解，将其分成四个部分：
由于输入矩阵是对角块，因此右上角部分忽略不计，剩下三个部分分别称作P1、P2、P3。
对于P1，它和POTRF的输入矩阵（完整的橙色矩阵）结构完全一致，因此直接递归调用POTRF进行计算，当P1的规模小于设定值时停止递归开始计算，后文详细介绍计算方法。
对于P2，使用TRSM即可完成对P2部分的计算，使用方式和上文相同。
对于P3，使用syrk可以完成P3外部依赖的计算，剩下的内部依赖继续调用POTRF即可完成计算。
接下来介绍递归停止时计算POTRF的实现，此时输入矩阵的典型规模为128，将其分成若干8x8的小块，然后计算每个列块（由小块组成的列）
每个列块，仍然需要先计算该列块的外部依赖（该列块左侧的所有列块），然后对列块中的每一列分别计算内部依赖，对于这两个部分可以分别用两个kernel来实现。由于这一步骤是严重的串行瓶颈，因此在划分小块时需要尽量让计算的快更小，减少串行瓶颈对性能的影响
3 MLU层需求分析
3.1 算子需求分析
3.2 算子功能和应用场景描述
厄密特矩阵，又称自伴随矩阵，是共轭对称的方阵。
对正定厄密特矩阵$A$ 进行Cholesky分解，即求矩阵$L$使下式成立：
其中，$L$ 是一个下三角矩阵且对角元素均为正实数，$L^*$ 表示$L$ 的共轭转置，是一个上三角矩阵。当$A$ 是一个实数矩阵时，Cholesky分解可以改写为
3.3 算子输入输出参数要求
4 算子接口设计
接口为：
变量含义为上文所述。
5 总结
本文介绍了在MLU上实现Cholesky分解的方案和需求分析。Cholesky分解是一种分解正定厄密特矩阵为下三角矩阵及其共轭转置的算法，广泛应用于科学和数值计算。本文首先解释了厄密特矩阵和Cholesky分解的基本原理，随后通过将输入矩阵分块，并利用BLAS标准操作中的的SYRK、GEMM和TRSM函数，以及自定义POTRF函数，展示了如何逐步实现分解。然后本文详细描述了算子的需求，包括支持的数据类型、形状、布局，以及特定的计算需求，如原位操作和步长机制，并提供了算子的接口设计。