【问题标题】:If C is row-major order, why does ARM intrinsic code assume column-major order?如果 C 是行优先顺序,为什么 ARM 内部代码假定列优先顺序?
【发布时间】:2021-08-18 04:01:54
【问题描述】:

我不确定在哪里问这个问题的最佳地点,但我目前正在使用 ARM 内部函数并遵循本指南:https://developer.arm.com/documentation/102467/0100/Matrix-multiplication-example

但是,那里的代码是假设数组以列优先顺序存储的。我一直认为 C 数组是按行存储的。他们为什么这么认为?

编辑: 例如,如果不是这样:

void matrix_multiply_c(float32_t *A, float32_t *B, float32_t *C, uint32_t n, uint32_t m, uint32_t k) {
    for (int i_idx=0; i_idx < n; i_idx++) {
        for (int j_idx=0; j_idx < m; j_idx++) {
            for (int k_idx=0; k_idx < k; k_idx++) {
                C[n*j_idx + i_idx] += A[n*k_idx + i_idx]*B[k*j_idx + k_idx];
            }
        }
    }
}

他们这样做了:

void matrix_multiply_c(float32_t *A, float32_t *B, float32_t *C, uint32_t n, uint32_t m, uint32_t k) {
    for (int i_idx=0; i_idx < n; i_idx++) {
        for (int k_idx=0; k_idx < k; k_idx++) {
            for (int j_idx=0; j_idx < m; j_idx++) {
                C[n*j_idx + i_idx] += A[n*k_idx + i_idx]*B[k*j_idx + k_idx];
            }
        }
    }
}

由于按 C[0]、C[1]、C[2]、C[3] 的顺序访问 C 的空间局部性,而不是按 C[0]、C[2 的顺序访问 C,代码将运行得更快]、C[1]、C[3](其中 C[0]、C[1]、C[2]、C[3] 在内存中是连续的)。

【问题讨论】:

  • 地址空间是线性的,不是二维的,没有“行”和“列”——有些元素就像更远和更近。 Why did they assume this? We have assumed a column-major layout of the matrices in memory. That is, an n x m matrix M is represented as an array M_array where Mij = M_array[n*j + i]. 他们可以做到n*j_idx + i_idxj_idx + n*i_idx。什么是列或行,只是您/他们的解释,您可以翻转术语并没问题。
  • @KamilCuk 好的,我想我明白你的意思了。所以我想我的问题应该是为什么他们决定低效地实现矩阵乘法而不是连续访问数组?
  • matrix multiply inefficiently 我没有看到更有效的方法......你是什么意思? not accessing the arrays contiguously 因为它是一个“矩阵乘法示例”——所以他们访问数组就好像它是一个矩阵一样,并且他们以特定的方式对其进行特定的操作。编程语言是一种抽象——总的来说,它都是机器指令。这只是一个抽象——当Mij = ...
  • @KamilCuk 通过在编辑中在 j 和 k 之间进行循环交换更有效
  • 您发布的代码无效 - C[n*j_idx + i_idx] 未初始化时使用,j_idx = 1 时使用,然后它会被归零。我首先会做memset(C, 0, ...)。不管怎样,页面上写着“矩阵乘法示例”——当然这并不是为了高效,而是展示它是如何工作的。

标签: c optimization matrix-multiplication neon row-major-order


【解决方案1】:

您没有使用像 C[i][j] 这样的 C 2D 数组,所以这不是 C 如何存储任何内容的问题,而是如何在这段代码中使用 n * idx_1 + idx_2 手动完成 2D 索引,可以选择在内部循环和外部循环中循环。

但是两个矩阵都未转置的 matmul 的难点在于,您需要为两个输入矩阵做出相反的选择: 天真的 matmul 必须跨越其中一个矩阵的遥远元素输入矩阵,所以它本质上是错误的。这就是为什么仔细的缓存阻塞/循环平铺对于矩阵乘法很重要的一个主要部分。 (O(n^3) 对 O(n^2) 数据进行处理 - 每次将其放入 L1d 缓存和/或寄存器时,您都希望充分利用它。)

循环交换可以加快速度,以利用最内层循环中的空间局部性,如果你做得对的话。

请参阅What Every Programmer Should Know About Memory? 中的缓存阻塞 matmul 示例,该示例遍历内部几个循环中所有 3 个输入中的连续内存,选择未在 3 个矩阵中的任何一个中缩放的索引作为内部矩阵。看起来像这样:

  for (j_idx)
    for (k_idx)
      for (i_idx)
          C[n*j_idx + i_idx] += A[n*k_idx + i_idx]*B[k*j_idx + k_idx];

请注意,B[k * j_idx + k_idx] 在循环内部循环上是不变的,并且您正在对连续内存进行简单的dst[0..n] += const * src[0..n] 操作(这很容易进行 SIMD 矢量化),尽管您仍在执行 2 次加载 + 1为每个 FMA 存储,所以这不会最大化您的 FP 吞吐量。

与缓存访问模式分开,这也避免了长依赖链到单个累加器(C 的元素)。但这对于优化的实现来说并不是真正的问题:您当然可以使用多个累加器。由于舍入误差,FP 数学不是严格关联的,但多个累加器是 closer to pairwise summation 并且可能比连续添加行 x 列点积的每个元素具有更小的 FP 舍入误差。 按标准简单 C 循环的顺序添加会有不同的结果,但通常更接近准确的答案。


您建议的循环顺序 i,k,j 是最差的。

您正在通过内部循环中 3 个矩阵中的 2 个的遥远元素,包括对 C[] 的不连续访问,这与您在上一段中所说的相反。

使用j 作为最内层循环,您将在第一次外部迭代中访问C[0]C[n]C[2n] 等。 B[] 也是如此,这真的很糟糕。

交换 ij 循环将使您可以在中间循环中连续访问 C[],而不是在最里面的循环。所以这将是一个严格的改进:是的,你是对的,这个天真的例子比它需要的更糟糕。

但关键问题是对 inner 循环中某些内容的跨步访问:这是性能灾难;这就是为什么仔细的缓存阻塞/循环平铺对于矩阵乘法很重要的一个主要部分。唯一从未与比例因子一起使用的索引是i

【讨论】:

    【解决方案2】:

    C 本质上不是行优先或列优先。

    在编写a[i][j] 时,由您决定i 是行索引还是列索引。

    虽然首先编写行索引(使数组以行为主)是一种常见的约定,但没有什么能阻止您做相反的事情。

    另外,请记住A × B = C 等价于Bt × At = Ctt 表示转置矩阵),并且读取行优先矩阵就像它是列优先(反之亦然)转置它,这意味着如果你想保持你的矩阵行优先,你可以颠倒操作数的顺序。

    【讨论】:

    • 由于数组存储在连续内存中的方式,通过增加第二个索引而不是第一个索引来扫描二维数组不是更有效吗?
    • @VeryConfused 是的,它可能更有效。但这与选择哪个索引是行哪个是列无关。
    • 我明白了,所以我想我的问题应该是为什么他们决定通过不连续访问数组来做一种效率较低的矩阵乘法方式
    • @VeryConfused 这不是真的。您可以增加内部循环中的第二个索引,无论它是指行还是列。没有人说你必须逐行而不是逐列遍历矩阵。
    • 问题编辑中的代码呢,如果进行循环交换,代码不会运行得更快吗?
    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2014-11-12
    • 1970-01-01
    相关资源
    最近更新 更多