高性能计算应用优化之代码实现调优(二)
内存优化
程序编写过程中对内存区域的分配、使用能够影响程序的性能，主要有以下几个方面：内存别名；内存引用；内存分配；内存访问；内存对齐；内存预取。
内存别名
当 xp 和 yp 指向同样的内存（memory aliasing）时，func1 和 func2 是两个完全不同的函数，所以编译器不会假定指针相同或不同，更不会尝试将 func1 优化为 func2。
void func1(long *xp, long *yp) {*xp += *yp;*xp += *yp;
}可以显式使用__restrict__修饰指针，表明不存在和被修饰的指针指向同一块内存的指针，此时编译器会将func1优化为和 func2 等效。
void func1(long *__restrict xp,long *__restrict yp) {*xp += *yp;*xp += *yp;
}内存引用
在程序计算中直接使用指针指向的值被称为“内存引用”，其过程为从内存中取出乘数 movss (%rdx), %xmm0，计算完之后，有需要更新，将其写回内存movss %xmm0, (%rdx)。这种操作涉及到了CPU从数据总线中向内存中取值，通常速度远远慢于CPU本身的计算操作，也慢于CPU取出内部寄存器值的操作，很多时候，一个程序的计算瓶颈就在这些去内存的操作中，因此要尽量避免不必要的内存引用。
内存引用会降低程序性能，所以可以按照函数的实际语义进行特定的优化。在函数multiply1中由于存在内存引用，编译器不能将数据存储在寄存器中，而必须将其存储在内存中。不开优化时生成的汇编中，每次循环都有三个操作：读内存一次（*dest）、执行乘法、写内存一次（*dest）。
void multiply1(double arr[], long length, double* dest) {
long i;
*dest = 1;
for (i = 0; i < length; i++) {
*dest = *dest * arr[i];
}
}函数multiply2是编译器开了-O2优化后对应的 C 语言代码，我们可以看到由于使用了一个临时变量，所以乘法操作不需要再读取内存，每次循环只有两个操作：执行乘法、写内存一次。
void multiply2(double arr[], long length, double* dest) {
long i;
double res = 0;
*dest = 1;
for (i = 0; i < length; i++) {
res = res * arr[i];
*dest = res;
}
}在函数multiply3中我们根据函数multiply1的语义进行了人工优化，代码不需要额外读内存，也不需要写内存，所以生成的汇编代码中，每次循环只有一个操作：执行乘法。
void multiply3(double arr[], long length, double* dest) {
long i;
double res = 0;
for (i = 0; i < length; i++) {
res = res * arr[i];
}
*dest = res;
}但这种优化存在内存别名问题，如果arr = [2, 3, 5]，那么multiply1(arr, 3, arr + 2)最终会得到arr[2] = 36，而multiply3(arr, 3, arr + 2)最终会得到arr[2] = 30，两者结果不一致。
内存分配
内存分配需要由一系列系统调用完成，合理的分配和回收策略能够提高程序的性能。在应用程序中有较多malloc（或allocate）时，可以尝试使用其他的内存分配库代替默认的内存分配库以提高性能。在线程数较少，或释放强度较低的情况下，较为简洁的tcmalloc性能稍胜过jemalloc。而在核心数较多、申请释放强度较高的情况下，jemalloc因为锁竞争强度远小于tcmalloc，会表现出较强的性能优势。若使用Intel编译器，也可尝试Intel的qkmalloc。
尽量不要在循环中执行动态数组分配、内存拷贝、数组复制等耗时的操作，减少不必要的系统调用。某些应用中每次迭代均需进行动态数组分配，可以将这些数组提取到module中，在初始化时分配，这样虽然增加了一些内存开销，但能极大降低运行时间开销。
内存拷贝、数组赋值等操作可能会调用memset函数，也需要尽量避免。例如下面的语句会首先赋初值为0，再进行计算，这时对a的初始化是不必要的，可以省略（需要谨慎地使用该方法）。
a = 0
……
do i = 1, 10a(i) = ……
enddo内存访问
系统对内存的访问一般是连续的。内存访问的一个经典问题是 AOS v.s. SOA。所谓 Array of structures (AOS)，即同一个结构体（以N-body问题中的粒子举例）中的属性在内存中连续。在 C++ 中，可以这么表示：
// Array of structures
struct Particle {float x, y, z, w};
Particle particles[1000];对于 Structure of array（SOA），可以这么表示：
// Structure of arrays
struct Particles {float x[1000];float y[1000];float z[1000];float w[1000];
};这样的好处是一个属性（如 x）的数据在内存中是连续的，如果有一些操作只用到某几个属性，那么 SOA 的内存布局可以完全不浪费 w 造成的内存带宽。在 AOS 的 数据布局下，由于一个结构体的 x, y, z, w 总是在同一个 cacheline 内，没有办法只取一部分元素。导致内存带宽的浪费。在 GPU 上，SOA 也能发挥硬件的合并机制，硬件会把一个 warp （32 个 thread）的数据访问打包以提高效率；在 CPU 上，SOA 对于 SIMD 模式的编程也更加友好，如 x86 架构下可以使用_mm_load_ps等 SIMD 指令对数据进行向量化加载。
内存访问的另一个问题是访问的连续与不连续。通常我们对内存地址的访问都是连续的，但在某些应用中会有如下的代码，在很多时候这种代码是很低效的，因为不连续内存地址的访问使编译器难以进行有效的优化，在实际编写代码时要注意避免此类实现。
a(:, jms+1:jme-1, 1:layer) = b下面这种传参方式会成为代码的瓶颈之一，编译器会首先将a数组的对应元素拷贝到临时数组，再进行传参。可以首先将a(:, jms+1:jme-1, 1:layer)拷贝到临时数组再进行传参。
call somefunc(a(:, jms+1:jme-1, 1:layer))内存对齐
内存对齐可以方便计算机去读写数据，提高对内存的读写速度。对齐的地址一般都是 n（n = 2、4、8）的倍数。例如4 个字节的变量，例如 float、int 类型的变量，放在 4 的整数倍地址上，8 个字节的变量，例如 long long、double 类型的变量，放在 8 的整数倍地址上。
假如没有内存对齐机制，数据可以任意存放，现在一个int变量存放在从地址1开始的联系四个字节地址中，该处理器去取数据时，需要读取多次才能获取到变量。内存对齐后，int类型数据只能存放在按照对齐规则的内存中，在取数据时一次性就能将数据读出来。
对C++而言，我们可以通过预编译命令#pragma pack(n)，其中n= 1,2,4,8,16来改变对齐系数，使用#pragma pack ()取消自定义字节对齐方式。
对Fortran而言，可以通过编译选项指定对齐，也可以通过在代码中加入制导语句指导编译器对齐。可以使用!DIR$ ASSUME_ALIGNED address1:n1告诉编译器，某个内存中的数组变量address是已经以n1对齐的，该语句必须放置在可执行代码中，例如
TYPE NODE
REAL(KIND=8), POINTER :: A(:,:)
END TYPE NODE
TYPE(NODE), POINTER :: NODES ALLOCATE(NODES)
ALLOCATE(NODES%A(1000,1000))
!DIR$ ASSUME_ALIGNED NODES%A(1,1) : 16
DO I=1,N
NODES%A(1,I) = NODES%A(1,I)+1
ENDDOFortran中也可使用ATTRIBUTES指令选项!DIR$ ATTRIBUTES ALIGN: n:: object为变量以及派生类型的可分配或指针变量object指定以n字节对齐。例如下边的代码中内存是以64字节对齐的。
TYPE EXAMPLE
!DIR$ ATTRIBUTES ALIGN : 64 :: R_alloc
REAL, ALLOCATABLE :: R_alloc ( : )
END TYPE EXAMPLE
TYPE (EXAMPLE) :: MyVar
ALLOCATE (MyVar%R_alloc(1000))GCC中提供了一些buildin函数，以实现对齐内存分配。__builtin_alloca函数在调用函数的堆栈上分配一个大小为size字节的对象。该对象在由__BIGGEST_ALIGNMENT__宏确定的目标默认堆栈对齐边界上对齐。__builtin_alloca函数返回分配对象的第一个字节的指针。分配对象的生命周期在调用函数返回给其调用者之前结束，即使__builtin_alloca在嵌套块内调用也是如此。例如，下面的函数在堆栈上分配了n字节大小的八个对象，并将每个对象的指针存储在数组a的连续元素中。然后，它将该数组传递给函数g，该函数可以安全地使用数组元素指向的存储空间。由于__builtin_alloca函数不验证其参数，因此调用者有责任确保参数不会导致超出堆栈大小限制。
void f (unsigned n)
{void *a [8];for (int i = 0; i != 8; ++i)a [i] = __builtin_alloca (n);
}__builtin_alloca_with_align函数在调用函数的堆栈上分配一个大小为size字节的对象。分配的对象在由参数alignment指定的边界上对齐，该参数的单位为位（而不是字节）。size参数必须为正，并且不能超过堆栈大小限制。该函数返回分配对象的第一个字节的指针。分配对象的生命周期在调用该函数的块的末尾结束。例如：在下面的函数中，对函数g的调用是不安全的，因为当overalign非零时，由__builtin_alloca_with_align分配的空间可能已经在调用它的if语句结束时释放了。
void f (unsigned n, bool overalign)
{void *p;if (overalign)p = __builtin_alloca_with_align (n, 64 /* bits */);elsep = __builtin_alloc (n);g (p, n); // unsafe
}内存预取
预取可以进一步帮助隐藏读取内存带来的延迟，分为硬件预取和软件预取。硬件预取由硬件预取器完成，根据以往数据的需求，来猜测以后需要的数据，并且提前加载进来，一般来说只有线性的地址访问规律（包括顺序、逆序；连续、跨步）能被识别出来；软件预取由程序通过宏或函数明确请求。
软件预取可以通过mm_prefetch实现，对应的Fortran程序如下所示，在循环中首先计算a(i,j)之后发出预取指令，将后续的元素提前取入。
subroutine spread_lf (a, b)
PARAMETER (n = 1025)
real*8 a(n,n), b(n,n), c(n)
do j = 1,n
do i = 1,100
a(i, j) = b(i-1, j) + b(i+1, j)
call mm_prefetch (a(i+20, j), 1)
call mm_prefetch (b(i+21, j), 1)
enddo
enddo
end也可以在do循环前使用制导语句指定预取。例如，以下的制导语句指定了不对A1进行预取，为B从L2及更高缓存发出向量化预取，以16为向量化的长度，
!DIR$ NOPREFETCH A1
!DIR$ PREFETCH B:1:16
DO J = 1,N
A1(J) = B(J-1) + B(J+1)
END DO对于C++来说，预取可以以如下方式实现，这里的提前量是64，不能太多，否则会污染缓存，也不能太少，无法隐藏计算延迟。
#pragma omp parallel for
for (size_t i = 0; i < n / 16; i++) {size_t next_r = randomize(i + 64) % (n / 16);_mm_prefetch(&a[next_r * 16], _MM_HINT_T0);size_t r = randomize(i) % (n / 16);for (size_t j = 0; j < 16; j++) {benchmark::DoNotOptimize(a[r * 16 + j]);}
}计算优化
消除伪共享
由于CPU之间通信的最小单位是缓存行（64 字节），在多线程环境下，两个核心可能会访问到同一缓存行（如下C++代码段）。假设一个核心修改了该缓存行的前 32 字节，另一个修改了后 32 字节，同时写回的话，不能两个都正确写入。CPU为了安全起见，同时只能允许一个核心写入同一地址的缓存行。从而导致读写这个变量的速度受限于三级缓存的速度，而不是一级缓存的速度。
int n = 1<<23;
std::vectora(n);
std::vectortmp(omp_get_max_threads());
#pragma omp parallel for
for (int i = 0; i < n; i++) {tmp[omp_get_thread_num()] += a[i];
}如果需要消除错误共享，需要把每个核心写入的地址尽可能分散。比如这里，我们把每个核心访问的地方跨越 16 KB，每个核心之间不会发生冲突，从而可以放心地放在一级缓存里进行操作。（该部分代码仅为示例，用户实际编写代码时需要自行设计避免冲突。）
std::vectortmp(omp_get_max_threads() * 4096);
#pragma omp parallel for
for (int i = 0; i < n; i++) {tmp[omp_get_thread_num() * 4096] += a[i];
}消除分支
消除分支能够提高性能，因为它减少了误预测的可能性，也减少了所要求的分支目标缓冲(Branch Target Buffer, BTB)项。在向量化中，分支也增加了编译器对代码进行向量化的难度，降低了向量化程度。
在实际编写代码时，有如下几种方法可以消除分支。例如使用mask数组代替if，进行部分冗余计算，安排代码使得基本块连续等。
消除计算
下面代码中的除法可以使用更廉价的乘法代替，改为 a*0.5f。
float func(float a) {return a / 2;
}进一步地，可以在循环中分离公共除数，提前计算好 b 的倒数避免重复求除法。
void func(float *a, float b) {for (int i = 0; i < 1024; i++) {a[i] /= b;}
}
void func(float *a, float b) {float inv_b = 1 / b;for (int i = 0; i < 1024; i++) {a[i] *= inv_b;}
}内联汇编
汇编语言的种类有两种风格：Intel Style和AT&T。GUN C compiler如GCC使用AT&T语法。使用内联汇编主要目的是为了提高效率，同时还是为了实现 C 语言无法实现的部分。在一些对性能要求较高的场景，使用内联汇编可以直接插入优化后的汇编代码，减少函数调用开销和其他额外开销。
内联汇编的基本格式如下，共四个部分：汇编语句，输出部分，输入部分，会被修改的部分。
asm("汇编语句"
: 输出部分
: 输入部分
: 会被修改的部分);如果汇编语句的声明和执行必须在我们指定的位置，（即不得从循环移出作为优化），需要使用关键字“volatile”或“__volatile__”声明。用法如下：
asm volatile ( "...;""...;": ... );
或
__asm__ __volatile__ ( "...;""...;": ... );