超标量编程 101(矩阵相乘)第 2 部分(共 5 部分)

作者:Jim Dempsey

在我的上一篇文章中,我们以下列图表结束了讨论


在上图中,我们能够很明显地看出两幅图表有“本质的不同”。这两幅图表是将非缓存敏感的串行技巧与缓存敏感的并行技巧相比较。它们是宣传并行编程的不错材料,当然也能激励开发人员了解如何对程序进行并行化,但是却不适合作为经验丰富的并行程序员分析使用。

任何并行程序员首先要知道的一点是:先改进串行算法,然后通过并行化解决问题(将内部小型向量函数留在最后使用)。

如果您研究串行函数的最内层循环,我们将发现:

         for (int k = 0; k < size; k++)
         {
            temp += m1[i][k] * m2[k][j];
         }

我们利用传统的 C/C++ 编程实践使用指向行的指针数组解决这一问题,通过实施这种实践获得了 6 倍到 12 倍的性能提升。我们还可以做些什么?

暂时将 TLB 和行表问题放在一边。在检查主要计算语句的过程中,我们发现 m1[i][k] 中的 k 索引会顺序访问,但 m2[k][j] 的 k 索引却不是。这意味着,编译器无法对这个循环进行向量化。可能更糟的是,在以特定间距处理行时会造成高速缓存移出(例如在酷睿 i7 920 中看到的性能曲线的明显下降)。对 double 型数据进行向量化能够获得 2 倍性能提升,非常值得一试。因此,让我们提高该程序串行版本的向量化程度(vector-ability)。(当然也包括串行版本的并行变体)

请注意,虽然我们可以使用在 Cilk++ 实施中的方法(消除行表指针),但是我将选择一种在稍后非常有用的替代方法。首先,让我们将这种方法应用到简单的串行和简单的并行示例中,看看会发生什么。

为了帮助编译器对矩阵乘法进行向量化(并且最大限度地降低高速缓存移出),您可以转置矩阵 m2(得到 m2t),然后您可以颠倒内循环中索引的顺序。最初,这听起来很荒谬。为了转置矩阵 m2,您需要额外的开销:

执行分配(实际上分配 size+1 的空间,1 个空间用作行指针)
运行转置循环
读取并重写一个数组(将数组 m2 变为 m2t)

您可能会认为,这样做肯定会比简单的矩阵乘法慢。
那您就想错了。矩阵 m2 中的每个元素使用 N 次,矩阵 m1 中的每个元素也使用 N 次。我们可能会将 2N**2 次读 + N*N 次写出现的高速缓存未命中降低到 N 次读 + N/2 次写 + (2N**2)/2 次读 + N/2 次写。“/2”是因为现在内循环变为一个对两个顺序数组进行操作的点(DOT)函数,编译器能够非常高效地对其进行向量化。重写的串行循环(并且将内循环转换为内联函数)如下所示:

// compute DOT product of two vectors, return result as double
inline double DOT(double* v1, double* v2, size_t size)
{
    double temp = 0.0;
    for(size_t i = 0; i < size; i++)
    {
       temp += v1[i] * v2[i];
    }
    return temp;
}

现在,带有转置功能的矩阵相乘函数,如下所示:

void matrix_multiplyTranspose(
double** m1, double** m2, double** result, size_t size)
{
    // create a matrix to hold the transposition of m2
    double** m2t = create_matrix(size);
    // perform transposition m2 into m2t
    for (size_t i = 0; i < size; i++)
    {
        for (size_t j = 0; j < size; j++)
        {
            m2t[j][i] = m2[i][j];
	  }
    }
    // Now matrix multiply with the transposed matrix m2t
    for (size_t i = 0; i < size; i++)
    {
        for (size_t j = 0; j < size; j++)
        {
            result[i][j] = DOT(m1[i], m2t[j], size);
        }
    }
    // remove the matrix that holds the transposition of m2
    destroy_matrix(m2t, size);
}

请注意,在实际实施过程中,您可能希望保留在第一次调用分配的缓冲区。然后,在后续调用时,您可以检查之前分配能否满足当前调用。如果能够满足,则绕过分配步骤。如果不能满足,则删除缓冲区并分配新的缓冲区。这些进一步的优化细节留给读者去完善。

在我们的图表中将包括分配/释放开销。另外请注意,如果 m2 矩阵将使用多次(例如在过滤系统中),您可以执行一次转置,然后多次重复使用转置的矩阵。您应该记住这些实施策略,并且在您自己的代码中进行应用。

让我们看看采用转置技术(ST)的串行方法与未采用转置的串行(S)和并行(P)代码的性能对比情况:

图 5(上图)

图 6


改进的技巧,仍然使用行表,但是与串行和并行方法相比,在性能上取得了很大提升。然而,改进的技巧达到的性能提升远远低于 Cilk++ 演示程序消除行表方法得到的性能提升。我们可以从中学到什么?

改进的串行代码,使用一条线程,在 N 值较大的区间超过了两种不同处理器上使用 4 条或 8 条线程的并行代码。

虽然这种方法比之前的并行性能提升了 4 倍,但它仍然低于 Cilk++ 方法,这些信息是否有助于获得更好的算法?

为什么串行转置方法在 N 为 513 附近性能出现大幅提升?
虽然并行转置(PT)比并行(P)方法快 4 倍,但实际上前者的性能低于串行转置方法!为什么!?!?

这些问题的答案在于该图表线条表示的是与原始(非转置)串行性能的对比。让我们了解一下(非转置)串行方法实际的运行时间,看看能不能发现什么线索:

图 7


问题就在这儿。当 N = 513 时,我们发现串行时间以更陡的斜率变化。(Q6600 的图表显示早在 N 约为 350 时就已经出现了这种斜率变化)。

酷睿i7 920 上 N = 513 有何重要意义?我们主要访问两个维度为 513 x 513 的 double 型数组。使用的字节数为 2 x 513 x 513 x 8 = 4,210,704。嗯......

从网络上搜集的酷睿 i7 920 规格,我们知道这款处理器有

32 KB 一级指令高速缓存 + 32 KB 一级数据高速缓存(每个核心一个)
256 KB 二级高速缓存(每个核心一个)
8 MB 三级高速缓存(共享)
那么,为什么在大约 4 MB 而不是 8 MB 时出现性能下降?老实说,我没有对这一点进行透彻分析,我假设这是一个 TLB 问题,因此我可以直接说明在这个时候,似乎高速缓存系统已经达到了使用极限。重要的是,该曲线在这一点的急剧上升是由于与未优化的串行技巧相比,采用更出色的串行技巧获得了巨大的超标量性能提升。

我们从中得到的重要信息是,与串行转置方法相比,使用非转置技巧的并行方法获得的性能提升微不足道。(虽然使用 Cilk++ 技巧仍然可以快 3 倍多)

现在,我们将 Cilk++ 行表消除方法与串行转置方法进行比较,得到:

图 8

图 9


Cilk++ 技巧获得的性能提升大约是串行转置方法的 3.25 倍。
更值得一提的是,并行转置方法并不比串行转置方法快。为什么?而且更重要的是,性能出色的原因是什么?

一般来说,当您发现并行代码的性能不如串行代码的性能时,很可能表示您已经遇到了内存带宽问题。Cilk++ 代码的性能曲线清楚地表明,这种内存带宽问题是由于并行转置技巧的高速缓存利用率较差导致的。

Cilk++ 方法是我们能够使用的最佳方法吗?

答案是:不是。

我们将在下一篇文章讨论这个问题(第 3 部分)。

Jim Dempsey
jim@quickthreadprogramming.com
www.quickthreadprogramming.com

有关编译器优化的更完整信息,请参阅优化通知
类别: