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

作者:Jim Dempsey

前面的文章(第 2 部分)中,我们了解了通过重组循环以及使用临时数组可以通过 SSE 小型向量优化(由编译器完成)获得的性能提升,然而通过更改布局以及数组访问顺序提高缓存利用率,我们可以获得更大性能提升。性能提升促使我们研究内存带宽的限制,这些限制使串行方法在性能上超过了(串行方法的)并行版本。

内存带宽限制是一项重要的考虑因素,并且会引起编程策略的变化:为了获得更高的性能提升,我们将从使数据访问模式尽可能满足最接近的高速缓存级别的角度解决这一问题。

然而我们应该如何做?

在支持超线程技术的系统上(例如酷睿 i7 920),一个核心内的超线程技术可共享 256 KB 二级高速缓存,视内部架构而定,可共享 32 KB 一级数据高速缓存。而插槽内的所有核心将共享 8 MB 三级高速缓存(在其它处理器型号上可能为 6 MB 或 12 MB)。

Cilk++ 方法使用常见的“分治方法”(或区块方法)将矩阵分割成较小的工作集,恰好适合线程可用的高速缓存。这是一个很好的起点。接下来我们需要做的是,如何在系统处理器内的高速缓存之间协调相关活动。

虽然您可以使用嵌套循环进行分割,但是请考虑下列因素:当您减小区块的大小时,您也增加了与线程调度程序交互的次数(进出最内层循环)。线程调度耗费的资源并不少。在第一个图表中,开销在 50 x 50 区块大小达到平衡。能否使用其它编程方法,使您能够利用 2x2、2x1 或 1x2 区块,获得卓越的性能?

答案是:可以

使用大多数线程工具包(例如 MS Parallel Collections、OpenMP 或 Cilk++),很难有效地将任务分配给特定线程或者一组线程。高速缓存级别任务分组是 QuickThread 的一项内置设计特性。例如:

 

// divide the iteration space across each multi-core processor
// L3$ specifies by L3 cache
//		.OR.
// within processor (Socket) for processors without L3 cache
parallel_for( OneEach_L3$, intptr_t(0), size,
    [&](intptr_t iBeginL3, intptr_t iEndL3)
    {
        // divide our L3 iteration space by L2 within this threads L3
        parallel_for( OneEach_Within_L3$ + L2$, iBeginL3, iEndL3,
            [&](intptr_t iBeginL2, intptr_t iEndL2)
            {
                // Here we are running as the Master thread of a
                // 2 team member team (or 1 in the event of older
                // processor)
                // 
                // Now bring in our other team member(s)
                // form team of threads sharing this threads L2 cache
                parallel_distribute( L2$,
                    [&](intptr_t iTMinL2, intptr_t nTMinL2)
                    {
                        // ... (Do Work)
                    } // [&](intptr_t iTMinL2, intptr_t nTMinL2)
                ); // parallel_distribute( L2$,
            } // [&](intptr_t iBeginL2, intptr_t iEndL2)
        ); // parallel_for( OneEach_Within_L3$ + L2$, iBeginL3, iEndL3,
    } // [&](intptr_t iBeginL3, intptr_t iEndL3)
); // parallel_for( OneEach_L3$, intptr_t(0), size,

 

外层循环对应每个插槽部分,更深一层循环对应每个二级高速缓存,内层循环对应按共享二级高速缓存的线程数量分割的 n 路(n-Way)(酷睿 i7 920 有 2 条线程,Q6600 也有 2 条线程)

这里介绍的技巧是使用 parallel_distribute 作为最内层部分,然后使用循环分割 parallel_distribute 内层状态机循环中的外层循环。什么意思?

QuickThread 中的 parallel_for 执行线程集选择和分割迭代空间的任务,但是没有具体执行迭代。由于 QuickThread 的这一特性,我们可以将迭代域传递到更深层的循环嵌套中。第二个循环也完成相同的任务(将迭代空间传递到最内层的 parallel_distribute)。现在,让我们了解一下为什么这样做非常重要。

这种替代技巧创造了一个让智能线程群合作解决问题的环境。线程将能够知道,哪些线程共享速度最高的缓存,哪些线程共享容量最大的缓存,哪些线程与当前下称不共享缓存。外层两个循环的子范围以及 parallel_distribute 中的团队成员编号(iTMinL2)都隐含了这种对线程的确认。

您可以很容易地通过下列代码,添加每个 NUMA 节点的其它外层循环:

parallel_for( OneEach_M0$, intptr_t(0), size,
...

然而,对于矩阵乘法这样做不会带来任何好处,除非您的系统上配备了在处理器外的四级高速缓存,并且处理器与同一 NUMA 节点共享该缓存。NUMA 感知问题可以在内存分配例程中得到很好的解决,可以将其放在循环结构上方的三级高速缓存分配层内。

内部状态机循环将分割输出数组:

插槽(三级高速缓存)
内核对/超线程逻辑处理器
之间的内核对(Amongst Core Pair)/超线程逻辑处理器


为了实现这一点,我们将矩阵分为 2x2 区块,每个 2x2 区块由一个双线程组提供服务。我将这个双线程组称为“双打(Tag Team)”(正像双打摔跤手一样)。在酷睿 i7 上,我们有 4 组双打,每个核心一组,如果算上同一核心上的超线程逻辑处理器,每个组有两条线程。在 Q6600 上双打变为共享相同二级高速缓存的两条线程(Q6600 有 2 块二级高速缓存,2 个核心共享一块)。

如果您曾将看过电视摔跤节目,您会知道这是一种摔跤形式,每组有两名成员。任何时刻,每组都应该让一名队员登上摔跤场(方台),当他们希望让队友上场与对手比赛时,必须和自己的队友拍手换人。我用斜体写应该是因为,拍手换人通常使换人一方在某一时刻有两名摔跤手站在摔跤场上。在这种情况下,他们会完全压倒对手。我将向您展示,如何利用 QuickThread 中提供的高速缓存指挥的线程调度实现同样的效果。


在同一核心内的超线程逻辑处理器有两条整数执行路径,但只有一条浮点执行路径。我们可以通过编码利用这一点,让双打的一条线程执行转置,同时另一条线程对 2x2 区块较低的 2x1 部分执行点积,并且在转置时对右上角的 1x1 执行点积。执行转置的线程,稍后将对区块 2x1 部分的左侧单元执行点积。这样双打的两条线程会执行不同的工作。由于我们使用 parallel_distribute,首先完成的线程(通常是不执行转置的线程),将承担下一个区块的转置工作。也就是说,您通过角色转换实现负载平衡。

输出矩阵被分为 2x2 区块。线程组从线程池派生为双线程组,每条线程共享级别最接近的高速缓存。在酷睿 i7 920 上,它们是每个核心内的超线程逻辑处理器;在英特尔 Q6600 上,它们是共享相同二级高速缓存的两个核心;在 AMD Opteron 270 上,它们是在同一 NUMA 节点上相同处理中的两个核心……


图 10

图 10 展示每个插槽区块初始的工作分配(彩色背景,按处理器线程组迭代(较重的轮廓),以及 2x2 单元区块最轻的轮廓。

位于对角线上的 2x2 区块需要进行额外的工作,以便将 m2 矩阵转职为 m2t 矩阵。主线程(双线程组中的偶数线程)使用超线程核心的整数执行单元执行转置(需要大量内存访问),另一条线程使用超线程核心的浮点执行单元对 2x2 区块的 3 个单元执行点积。通过探听双打中主线程的转置进度,并发执行点积和转置。然后双线程组中的主线程执行最终的点积。还需要进行一些实验,才能了解双线程组中的从线程能否执行所有 4 个点积。参见下文:



第一个双线程组对指定区域(如上所示)最左上角的 2x2 区块进行操作(以及转置)。在并发进行此操作的同时,第二组对指定区域最左上角区块进行操作(以及转置):



以此类推到达最后一个插槽,最后一组对指定区域最左上角的 2x2 区块进行操作(以及转置):



完成并行结构并不意味着完成了对这些对角区块的操作。相反,向邮箱写入完成状态意味着完成。这种完成将与其它线程的活动异步。而且由于对共享的内存邮箱进行非互锁写入会产生一定“开销”。

完成对角线最左上角的操作后,双线程组会查看一个标志,指示是否有处于工作饥饿状态的线程(初期该标识将指示没有处于工作饥饿状态的线程)。如果没有处于工作饥饿状态的线程,双线程组会继续刚完成的对角线区块所在列中的其它 2x2 区块进行处理:


由上表中绿色的列表示。然而,这时每条线程会计算 2x2 区块中 1x2 的双点积。完成对 2x2 区块的操作后,双线程组会查看一个标志,指示是否有处于工作饥饿状态的线程,如果没有处于工作饥饿状态的线程,双线程组会继续执行该列的其它区块,如果有处于工作饥饿状态的线程,它会沿着对角线进行处理。

按列依次处理的目的在于,在转置后刚刚转置的两列最有可能存储于高速缓存中(一级高速缓存、二级高速缓存或三级高速缓存)。

当双线程组完成对角线以及在对角线上/下方的列单元(位于其双线程组区域内)操作后,该线程组将查看在其三级高速缓存中其它线程组完成转置工作的状态(通过邮箱)



当第一组(插槽 0 中的线程 0/1)完成其工作,而且第二组(2/3)完成了对角线最左上角的操作后,第一组将处理浅绿色列。从统计角度来说,这些工作将在查看邮箱后完成(在利用点积进行 4 次转置,并且加上 12 次双点积时间延迟)。

每个线程组都相同。在线程组处理与其共享三级高速缓存的线程组对角线上/下方的列时,如果它发现三级高速缓存中指定的列不完整,它会发出“工作饥饿”线程通知,以便其三级高速缓存共享线程组能够中断列处理,在完成当前列之前先对下一个对角线进行处理。

除了处理处于饥饿状态的线程标志以外,每个线程组都与其插槽中的其它线程组保持独立。

这一迭代过程将一直继续,直到某个线程组完成在其插槽区块中指定的所有工作:



这时,它将查看其它插槽对角线的对角线完成状态邮箱。因此这是一个状态机而不是嵌套的内部循环,它会以任何顺序完成其它插槽的对角线,但是很可能会以每个其它插槽的每个双线程组中对角线 2x2 区块的升序完成。如单独的绿色条所示,在插槽 0 组 0/1 区域右侧的插槽 1 区域,位于下面第二组红色区域(插槽 1)已完成对角线的上方。




下面,我们来制作图表,核实这些结果数据:




当 N = 128 : 1024 时,平均扩展系数(与串行转置相比)为:

采用并行转置双打方法时,Q6600 的扩展系数为 4.96 倍,采用 Cilk++ 时为 3.06 倍
采用并行转置双打方法时,酷睿 i7 920 的扩展系数为 4.69 倍,采用 Cilk++ 时为 3.20 倍

指定部分的峰值,大约 880 个 Q6600 显示采用并行双打方法时大约为 7.5 倍,采用 Cilk++ 方法是为 3 倍。

并行双打方法绝对优于 Cilk++ 方法(至少在该 N 范围内)。

为丰富该图,我们将其与不带转置的原始串行方法进行比较



现在,我们可以看到,在图表的指定点,我们在 Q6600 上实现了高出串行方法 38 倍的性能提升,在酷睿 i7 920 上实现了高出串行方法近 80 倍的性能提升。

并行转置双打的完整描述可以在 QuickThread 演示套件中包含的 MatrixMultiply.cpp 样本代码中找到。www.quickthreadprogramming.com

这一点能否改进……

我的答案是“可以改进”。

我在这篇博文的前面提到过:“永久保存内置小型矢量函数”。

在矩阵相乘的内部,最内层循环执行点乘积。所有编程变量都使用内联函数执行点乘积。QuickThread 并行双打(PTT)使用该点函数上的额外变量同时生成两个点乘积。2 线程组的各条线程都在 2x2 区块上运行,为 1/2 区块生成结果(1x2)。同时在一条线程内执行 2 个点乘积可以提高大型矩阵的高速缓存命中率。

以下是两个点函数,它们经过修改可以使用 xmm 内置函数(我不擅长使用 xmm 内置函数)和两个选择函数(可以调用测试程序中相应的点)。

// compute DOT product of two vectors, returns result as double
// used in QuickThread variations
double DOT(double v1[], double v2[], intptr_t size)
{
    double temp = 0.0;
    for(intptr_t i = 0; i < size; i++)
    {
       temp += v1[i] * v2[i];
    }
	return temp;
}

double xmmDOT(double v1[], double v2[], intptr_t size)
{
	// __declspec(align(16)) not working reliably for me
	double temp[4];
	intptr_t alignedTemp = (((intptr_t)&temp[0]) & 8) >> 3;

	__m128d	_temp = _mm_set_pd(0.0, 0.0);
	__m128d *_v1 = (__m128d *)v1;
	__m128d *_v2 = (__m128d *)v2;

	intptr_t	halfSize = size / 2;

	for(intptr_t i = 0; i < halfSize; i++)
	{
	   _temp = _mm_add_pd(_temp, _mm_mul_pd(_v1[i], _v2[i]));
	}
	// fix code to remove temp[4] array
	_mm_store_pd( &temp[alignedTemp], _temp);
	if(size & 1)
		temp[alignedTemp] += v1[size-1] * v2[size-1];

	return temp[alignedTemp] + temp[alignedTemp+1];
}

// compute two DOT products at once
// effectively
// r[0] = DOT(v1, v2, size);
// r[1] = DOT(v1, v3, size);
// except running both results at the same time
void DOTDOT(double v1[], double v2[], double v3[], double r[2],  intptr_t size)
{
	double	temp[2];
	temp[0] = 0.0;
	temp[1] = 0.0;
	for(int i=0; i < size; ++i)
	{
		temp[0] += v1[i] * v2[i];
		temp[1] += v1[i] * v3[i];
	} // for(int i=0; i < size; ++i)
	r[0] = temp[0];
	r[1] = temp[1];
}

// compute two DOT products at once
// effectively
// r[0] = DOT(v1, v2, size);
// r[1] = DOT(v1, v3, size);
// except running both results at the same time
void xmmDOTDOT(double v1[], double v2[], double v3[], double r[2],  intptr_t size)
{
	// __declspec(align(16)) not working reliably for me
	double	temp[6];
	intptr_t alignedTemp = (((intptr_t)&temp[0]) & 8) >> 3;
	__m128d	_temp0 = _mm_set_pd(0.0, 0.0);
	__m128d	_temp1 = _mm_set_pd(0.0, 0.0);
	__m128d *_v1 = (__m128d *)v1;
	__m128d *_v2 = (__m128d *)v2;
	__m128d *_v3 = (__m128d *)v3;

	intptr_t	halfSize = size / 2;

	for(intptr_t i = 0; i < halfSize; i++)
	{
	   _temp0 = _mm_add_pd(_temp0, _mm_mul_pd(_v1[i], _v2[i]));
	   _temp1 = _mm_add_pd(_temp1, _mm_mul_pd(_v1[i], _v3[i]));
	}
	_mm_store_pd( &temp[alignedTemp], _temp0);
	_mm_store_pd( &temp[alignedTemp+2], _temp1);
	if(size & 1)
	{
		temp[alignedTemp] += v1[size-1] * v2[size-1];
		temp[alignedTemp+2] += v1[size-1] * v3[size-1];
	}

	r[0] =  temp[alignedTemp] + temp[alignedTemp+1];
	r[1] =  temp[alignedTemp+2] + temp[alignedTemp+3];
}

bool UseXMM = false;

double doDOT(double v1[], double v2[], intptr_t size)
{
	if(UseXMM)
		return xmmDOT(v1, v2, size);
	return DOT(v1, v2, size);
}

void doDOTDOT(double v1[], double v2[], double v3[], double r[2],  intptr_t size)
{
	if(UseXMM)
		xmmDOTDOT(v1, v2, v3, r, size);
	else
		DOTDOT(v1, v2, v3, r, size);
}

如评论部分所述,我遇到过编译器指令对齐问题,而且必须稍稍调整代码以解决这一问题。将 xmm 内置函数融合到这两个例程中相对简单(对于像我这样没有经验的 xmm 程序员)。我们来看看结果:



在 Q6600(4 核心,无超线程技术)上,采用 xmm 内置函数的 S/PTTx 明显有优势,大约 30% 的改进。QuickThread 并行双打“x”方法主要使用 xmmDOTDOT(双点函数)。然而,使用 xmmDOT(单点)函数的 S/STx 稍逊于不带内置函数的 C++ 代码。这说明,编译器优化好于经验不足的内置函数程序员实施的手动优化(很常见)。

分析酷睿 i7 920 时,我们看到了不同的画面:



该处理器(采用超线程技术)使用手写 xmm 帮助函数时,性能下降了 12%。请注意,两个系统的可执行内容相同。

对于那些提供可在各种系统上使用的程序的厂商,这条信息非常重要。了解平台特定行为,意味着你可以在程序启动时查询系统,然后更改可以使用的代码变量。通常,你需要在代码中选择一个函子(函数指针),而不是使用选择函数。

对于未知行为,你可以使用启发式方法。在前几个调用中,应用可以选择代码的各个变量,测算各个方法的执行时间。获取足够的样本后,选择性能最好的函数变量,将相应的函子插入调度指针。

这一点能否改进……

我的答案是“可以改进”。

请注意,图表线条太过紧凑。额外调整可以让线条更加协调,因此请向上移动趋势线条。这相当于比当前并行转置双打方法预计改进了 15%。(无超线程技术的系统采用 xmmDOTDOT,有超线程技术的系统采用 DOTDOT)

15% 通常是微不足道的,然而,请注意并行双打方法和 Cilk++ 方法似乎在 1024 点出现了下降。对于大型矩阵,尤其是在多插槽系统上,应当运行更多测试。当我有机会收集此类数据时,我定会公布关于该性能测试的最新信息。

有改进这一点的高超技巧吗?
经验证明,的确有。

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

 

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