优化C++软件(12)

9.9. 顺序访问数据

在数据被顺序访问时，缓存的工作最高效。在前向访问数据时，效率稍低，而在以随机的方式访问数据时，效率要低得多。这适用于读与写数据。

应该在最里层循环中，通过改变最后一位索引来访问多维数组。这反映了元素在内存中保存的顺序。例子：

// Example 9.4

const int NUMROWS = 100, NUMCOLUMNS = 100;

int matrix[NUMROWS][NUMCOLUMNS];

int row, column;

for (row = 0; row < NUMROWS; row++)

     for (column = 0; column < NUMCOLUMNS; column++)

          matrix[row][column] = row + column;

不要交换这两个循环的次序（除了在Fortran里，那里的储存顺序是相反的）。

9.10. 在大数据结构中的缓存竞争

不总是可能顺序访问一个多维数组。某些应用（如线性代数）要求其他访问模式。如果在一个大矩阵中，行间距离恰好等于关键步长，这会导致严重的时延，如第71页所述。如果矩阵行的大小是2的一个大指数，这将出现。

下面的例子展示了这。我的例子是转置一个二次矩阵的函数，即每个matric[r][c]元素与matric[c][r]元素交换。

// Example 9.5a

const int SIZE = 64;                                        // number of rows/columns in matrix

void transpose(double a[SIZE][SIZE]) {      // function to transpose matrix

      // define a macro to swap two array elements:

      #define swapd(x,y) {temp=x; x=y; y=temp;}

      int r, c; double temp;

      for (r = 1; r < SIZE; r++) {                         // loop through rows

            for (c = 0; c < r; c++) {                       // loop columns below diagonal

                  swapd(a[r][c], a[c][r]);               // swap elements

            }

      }

}

void test () {

      __declspec(__align(64))                       // align by cache line size

      double matrix[SIZE][SIZE];                   // define matrix

      transpose(matrix);                                // call transpose function

}

转置一个矩阵与在对角线反射它一样。对角线下的每个元素matric[r][c]与对角线上镜像位置的元素matric[c][r]交换。例子9.5a中c循环从最左边的列到对角线。对角线上的元素不变。

这个代码的问题是，如果对角线下元素matric[r][c]是按行访问的，那么对角线上的镜像元素matric[c][r]就是按列访问的。

现在假定我们在Pentium 4计算机上使用一个64 x 64矩阵运行这个代码，其中1级数据缓存是8 kb = 8192字节，4路，行大小是64。每缓存行可以保存8个8字节的double。关键步长是8192/ 4 = 2048字节 = 4行。

让我们看一下在循环里发生了什么，例如在r = 28时。我们从对角线下的行28获取元素，并与对角线上列28交换。行28的前8个元素共享同一个缓存行，但这8个元素将进入列28的8个不同的缓存行，因为缓存行遵循行，但不遵循列。这些缓存行中的每4个属于同一个缓存集。在我们到达列28的元素16时，缓存将逐出由这个列元素0使用的缓存行。元素17将逐出元素1。元素18将逐出元素2，以此类推。这意味着，在我们交换列29与行29时，我们在对角线上使用的所有缓存行已经丢失了。每个缓存行必须被重载8次，因为在我们需要下一个元素前它被逐出了。通过测量在Pentium 4上，使用例子9.5a转置不同大小矩阵所需的时间，我确认。我实验的结果在下面给出。时间单位是时钟周期每数组元素。

矩阵大小	总k字节	每元素时间
63 x 63	31	11.6
64 x 64	32	16.4
65 x 65	33	11.8
127 x 127	126	12.2
128 x 128	128	17.4
129 x 129	130	14.4
511 x 511	2040	38.7
512 x 512	2048	230.7
513 x 513	2056	38.1

表9.1. 转置不同大小矩阵的时间，时钟周期每元素

这个表显示在矩阵大小是1级缓存大小的倍数时，转置该矩阵需要多40%的时间。这是因为关键步长是矩阵行的倍数。这个时延比从2级缓存重新载入1级缓存所需的时间短，因为乱序执行机制可以预取数据。

当在2级缓存中出现竞争时，影响要大得多。2级缓存是512kb，8路。2级缓存的关键步长是512kb/ 8 = 64kb。这对应一个512 x 512矩阵的16行。表9.1中我的实验结果显示，当在2级缓存中出现竞争时，比不出现竞争，需要6倍的时间转置一个矩阵。为什么2级缓存竞争的影响比1级缓存竞争大得多的原因是，2级缓存一次不能预取超过1行。

解决这个问题的一个简单方法是使矩阵中的行超过所需，以避免关键步长是矩阵行大小的倍数。我尝试使矩阵512 x 520，留下8列不使用。这消除了竞争，时间消耗降到36。

存在不可能向矩阵添加不使用列的情形。例如，一个数学函数库应该在所有矩阵大小上都有高效工作。在这个情形里，一个高效的解决方案是把矩阵分为较小的方块，一次处理一个方块。这称为方形分块或tiling。这技术在例子9.5b中展示。

// Example 9.5b

void transpose(double a[SIZE][SIZE]) {

      // Define macro to swap two elements:

      #define swapd(x,y) {temp=x; x=y; y=temp;}

      // Check if level-2 cache contentions will occur:

      if (SIZE > 256 && SIZE % 128 == 0) {

            // Cache contentions expected. Use square blocking:

            int r1, r2, c1, c2; double temp;

            // Define size of squares:

            const int TILESIZE = 8; // SIZE must be divisible by TILESIZE

            // Loop r1 and c1 for all squares:

            for (r1 = 0; r1 < SIZE; r1 += TILESIZE) {

                  for (c1 = 0; c1 < r1; c1 += TILESIZE) {

                        // Loop r2 and c2 for elements inside sqaure:

                        for (r2 = r1; r2 < r1+TILESIZE; r2++) {

                              for (c2 = c1; c2 < c1+TILESIZE; c2++) {

                                    swapd(a[r2][c2],a[c2][r2]);

                              }

                        }

                  }

                  // At the diagonal there is only half a square.

                  // This triangle is handled separately:

                  for (r2 = r1+1; r2 < r1+TILESIZE; r2++) {

                        for (c2 = r1; c2 < r2; c2++) {

                              swapd(a[r2][c2],a[c2][r2]);

                        }

                  }

            }

      }

      else {

            // No cache contentions. Use simple method.

            // This is the code from example 9.5a:

            int r, c; double temp;

            for (r = 1; r < SIZE; r++) {               // loop through rows

                  for (c = 0; c < r; c++) {              // loop columns below diagonal

                        swapd(a[r][c], a[c][r]);       // swap elements

                  }

            }

      }

}

在我的实验里，对于512x512矩阵，这个代码每元素需要50时钟周期。

2级缓存中的竞争代价是如此高，对它们做些什么十分重要。因此，你应该留意矩阵列数是2大指数的情形。1级缓存中的竞争代价不那么高。对1级缓存使用复杂的技术，像方形分块，也许不值得。

方形分块与类似的方法在. Goedecker与A. Hoisie的书《Performance Optimization of Numerically Intensive Codes》，SIAM 2001中有进一步描述。

9.11. 显式缓存控制

带有SSE与SSE2指令集的微处理器有允许你操作数据缓存的特定指令。这些指令可从支持固有函数的编译器（即Microsoft、Intel与Gnu）访问。其他编译器需要汇编代码使用这些指令。

函数	汇编名	固有函数名	指令集
Prefetch	PREFETCH	_mm_prefetch	SSE
Store 4 bytes without cache	MOVNTI	_mm_stream_si32	SSE2
Store 8 bytes without cache	MOVNTQ	_mm_stream_pi	SSE
Store 16 bytes without cache	MOVNTPS	_mm_stream_ps	SSE
Store 16 bytes without cache	MOVNTPD	_mm_stream_pd	SSE2
Store 16 bytes without cache	MOVNTDQ	_mm_stream_si128	SSE2

表9.2. 缓存控制指令

除了表9.2提到的，还有其他缓存控制指令，比如冲刷与屏障指令，但这些几乎与优化无关。

数据预取

预取指令可用于预取我们预期在程序流中稍后会使用的缓存行。不过，在任何我测试过的例子中，这没有提高执行速度。原因是归因于乱序执行以及先进的预测机制，现代处理器自动预取数据。对包含具有不同步长的多个流的规则访问模式，现代微处理器能够自动预取数据。因此，如果数据访问以具有固定步长的规则模式进行，你无需显式预取数据。

非缓存内存储存

非缓存写比非缓存读代价高，因为这个写会导致整个缓存行被读入，然后回写。

所谓的非临时写指令（MOVNT）设计来解决这个问题。这些指令直接写入内存，不载入缓存行。在写非缓存内存，以及不期望在该缓存行被逐出前从相同或附近地址读的情形里，这是有优势的。对同一个内存区域，不要混用非临时写与普通写或读。

非临时写指令不适用例子9.5，因为我们从相同的地址读写，因此缓存行无论如何都将被载入。如果我们修改例子9.5，使它只写，那么非临时写的影响变得显著。下面的例子转置一个矩阵，并将结果保存在另一个数组中。

// Example 9.6a

const int SIZE = 512; // number of rows and columns in matrix

// function to transpose and copy matrix

void TransposeCopy(double a[SIZE][SIZE], double b[SIZE][SIZE]) {

      int r, c;

      for (r = 0; r < SIZE; r++) {

            for (c = 0; c < SIZE; c++) {

                  a[c][r] = b[r][c];

            }

      }

}

这个函数按列模式写矩阵a，其中关键步长将导致所有的写在1级与2级缓存中载入新的缓存行。使用非临时写指令阻止了2级缓存为矩阵a载入任何缓存行：

// Example 9.6b.

#include "xmmintrin.h" // header for intrinsic functions

// This function stores a double without loading a cache line:

static inline void StoreNTD(double * dest, double const & source) {

      _mm_stream_pi((__m64*)dest, *(__m64*)&source); // MOVNTQ

      _mm_empty(); // EMMS

}

const int SIZE = 512; // number of rows and columns in matrix

// function to transpose and copy matrix

void TransposeCopy(double a[SIZE][SIZE], double b[SIZE][SIZE]) {

      int r, c;

      for (r = 0; r < SIZE; r++) {

            for (c = 0; c < SIZE; c++) {

                  StoreNTD(&a[c][r], b[r][c]);

            }

      }

}

在Pentium 4计算机上测量不同矩阵大小的每矩阵单元的执行时间。测量结果如下：

矩阵大小	例子9.6a每元素时间	例子9.6b每元素时间
64x64	14.0	80.8
65x65	13.6	80.9
512x512	378.7	168.5
513x513	58.7	168.3

表9.3. 转置与拷贝不同大小矩阵的时间，时钟周期每元素

如表9.3所示，仅在预期一个2级缓存不命中时，不使用缓存保存数据的方法是有优势的。64x64矩阵在1级缓存中导致不命中。这几乎对整体执行时间没有影响，因为在一个写操作上的缓存不命中不会推迟后续指令。512x512矩阵导致2级缓存不命中。这对执行时间有巨大的影响，因为内存总线饱和。这可以通过使用非临时写缓解。如果可以其他方式防止缓存竞争，如9.10节所述，非临时写指令不是最优的。

使用表9.2中列出的指令有某些限制。所有这些指令要求微处理器有SSE或SSE2指令集，如表中列出。16字节指令指令MOVNTPS、MOVNTPD与MOVNTDQ要求操作系统支持XMM寄存器；参考第125页。

在使用#pragma vector nontemporal时，Intel编译器可以在向量化代码中自动插入非临时写。不过，在例子9.6b中这不能工作。

在任何浮点指令之前，MOVNTQ指令必须后跟一条EMMS指令。如例子9.6b所示，这被编码为_mm_empty()。MOVNTQ指令不能用在Windows的64位设备驱动中。