优化任务：

4.1 原始程序及说明：

运行环境：Linux
程序框架部分如下：
说明：生成数组时，令其依次等于 1, 2, ...，这样进行平滑处理后的数组与原数组相同，方便判断程序是否正确运行；

const int N = 1920;
const int M = 1080;
int img[N][M], tmp[M]; // [0, N), [0, M)

void init(int A[N][M])
{
    int c = 0;
    for (int i=0; i< N; i++)
        for (int j=0; j< M; j++) A[i][j] = ++c;
}
bool judge(int A[N][M])
{
    int c = 0;
    for (int i=0; i< N; i++)
        for (int j=0; j< M; j++) if (A[i][j]!=(++c)) return 0;
    return 1;
}
void work0(int A[N][M])
{
    for (int j=1; j< M-1; j++) tmp[j] = A[0][j];
    int t;
    for (int i=1; i< N-1; i++) {
        for (int j=1; j< M-1; j++) {
            t = (A[i-1][j] + A[i+1][j] + A[i][j-1] + A[i][j+1]) >> 2;
            A[i-1][j] = tmp[j], tmp[j] = t;
        }
    }
    for (int j=1; j< M-1; j++) A[N-2][j] = tmp[j];
}
void call(void(*f)(int(*)[M]))
{
    clock_t s1, e1;
    double tim = 0.0;
    int T = 1000;
    for (; T; T--) {
        init(img);
        s1 = clock();
        f(img);
        e1 = clock();
        tim += (double)(e1 - s1)/ CLOCKS_PER_SEC;
        if (!judge(img)) puts("Error"), getchar();
    }
    printf(" %lf ms\n", tim);
}
int main()
{
    int cas = 1;
    for (; cas; cas--) {
        printf("func0 : "), call(work0);
        printf("func1 : "), call(work1);
        printf("func2 : "), call(work2);
        printf("func3 : "), call(work3);
        printf("func4 : "), call(work4);
        printf("func5 : "), call(work5);
        printf("func6 : "), call(work6);
        printf("func7 : "), call(work7);
        printf("-------------------\n");
    }
}

初始函数 work0 ，说明：
采用逐行遍历的方式。tmp 数组存储上一行的运算结果，在当前行计算完成后更新上一行。
可能的瓶颈：引用数组较为频繁；数组的每个位置被多次引用；数组引用地址的计算存在重复；cache 命中率可能存在提升空间，等等。

4.2 优化后的程序及说明

面向 CPU 的优化：work0, work1, work2, work3, work4, work5, work6, work7
面向 cache 的优化：work0, work8

work0

需要说明的是，实际上 work0 已经是相对某些实现方式的优化，包括但不限于：采用 1x1080 的额外空间而不是 2x1080 或更多；采用地址上连续的遍历方式能够提升空间局部性和时间局部性。

work1（2x1循环展开，倒序枚举）

void work1(int A[N][M])
{
    for (int j=M-2; j> 0; j--) tmp[j] = A[N-1][j];
    int t0, t1; //bool flg = (M&1)>0;
    for (int i=N-2; i> 0; i--) {
        for (int j=M-2; j> 0; j-=2) {
            t0 = (A[i-1][j] + A[i+1][j] + A[i][j-1] + A[i][j+1]) >> 2;
            t1 = (A[i-1][j-1] + A[i+1][j-1] + A[i][j-2] + A[i][j]) >> 2;
            A[i+1][j] = tmp[j], tmp[j] = t0;
            A[i+1][j-1] = tmp[j-1], tmp[j-1] = t1;
        }
        /*if (flg)
            t0 = (A[i-1][1] + A[i+1][1] + A[i][0] + A[i][2]) >> 2,
            A[i-1][1] = tmp[1], tmp[1] = t0;*/
    }
    for (int j=M-2; j> 0; j--) A[1][j] = tmp[j];
}

2x1 循环展开：但由于不存在乘法这类时钟周期大的运算，估计该优化无效果，仅做对比；
倒序枚举：去除了正序枚举中每次判断循环结束时需要额外计算 N-1 和 M-1 ，但由于 N 和 M 已申明为常量，估计编译后已经不需要额外计算，仅做对比。

work2（3x1 循环展开+指针调用，临时变量）

void work2(int A[N][M])
{
    int t, a0, a1, a2, *u, *d, *l, *r, *p, i, j;
	for (j=0; j< M; j++) tmp[j] = A[N-1][j];
    u = &A[N-3][M-2], d = &A[N-1][M-2];
    l = &A[N-2][M-3], r = &A[N-2][M-1];
    p = &tmp[M-2];
    for (i=N-2; i> 0; i--) {
        a0 = *r, a1 = *(r-1), a2 = *l;
        for (j=M-2; j> 1; j-=3) {
            t = (a0 + a2 + (*u) + (*d)) >> 2;
            *d = *p, *p = t;
            u--, d--, l--, r--, p--; a0 = *l;
            t = (a0 + a1 + (*u) + (*d)) >> 2;
            *d = *p, *p = t;
            u--, d--, l--, r--, p--; a1 = *l;
            t = (a1 + a2 + (*u) + (*d)) >> 2;
            *d = *p, *p = t;
            u--, d--, l--, r--, p--; a2 = *l;
        }
        // 1078 % 3 = 1
        t = (a0 + a2 + (*u) + (*d)) >> 2;
        *d = *p, *p = t;
        u -= 3, d -= 3, l -= 3, r -= 3, p = &tmp[M-2];
    }
    for (j=M-2; j> 0; j--) A[1][j] = tmp[j];
}

3x1 循环展开+指针调用：注意到引用数组时需要通过行乘列大小加列的方式计算地址，可以采用指针 *u, *d, *l, *r, *p 优化掉重复计算，分别表示上、下、左、右和 tmp 对应位置；注意到对于同一行遍历时，每个位置会被其左、右各引用一次，我们可以用临时变量 a0, a1, a2 存储它；t 用于临时存储计算结果，i, j 为循环变量；共计 11 个临时变量；
若这些临时变量有部分存于寄存器时，预计将有非常优秀的优化效果！

work3（添加 register）

在 work2 的基础上，在临时变量声明前加 register，建议编译器将其存于寄存器中。

work4（减少 2 个临时变量）

void work4(int A[N][M])
{
    register int a0, a1, a2, *u, *d, *l, *p, i, j;
	for (j=0; j< M; j++) tmp[j] = A[N-1][j];
    u = &A[N-3][M-2], d = &A[N-1][M-2];
    l = &A[N-2][M-3];
    p = &tmp[M-2];
    for (i=N-2; i> 0; i--) {
        a0 = *(l+2), a1 = *(l+1), a2 = *l;
        for (j=M-2; j> 1; j-=3) {
            a0 = (a0 + a2 + (*u) + (*d)) >> 2;
            *d = *p, *p = a0;
            u--, d--, l--, p--; a0 = *l;
            a1 = (a0 + a1 + (*u) + (*d)) >> 2;
            *d = *p, *p = a1;
            u--, d--, l--, p--; a1 = *l;
            a2 = (a1 + a2 + (*u) + (*d)) >> 2;
            *d = *p, *p = a2;
            u--, d--, l--, p--; a2 = *l;
        }
        // 1078 % 3 = 1
        a0 = (a0 + a2 + (*u) + (*d)) >> 2;
        *d = *p, *p = a0;
        u -= 3, d -= 3, l -= 3, p = &tmp[M-2];
    }
    for (j=M-2; j> 0; j--) A[1][j] = tmp[j];
}

仔细分析代码，可以注意到我们能完美地将 t, *r 省去，用其他变量代替其在计算式中的位置；但是有可能降低了并行性。

work5（减少 2 个临时变量）

去除 *u, *d，而用 *(l-1079), *(l+1081) 代替之；
由于之前每次枚举时需要执行 u--, d-- ，故修改后运算指令数大致不变。

work6（修改运算顺序，减少 1 个临时变量）

void work6(int A[N][M])
{
    register int a0, a2, *l, *p, i, j;
	for (j=0; j< M; j++) tmp[j] = A[N-1][j];
    l = &A[N-2][M-3], p = &tmp[M-2];
    for (i=N-2; i> 0; i--) {
        a0 = *(l+2), a2 = *l;
        for (j=M-2; j> 1; j-=3) {
            a0 = (a0 + a2 + (*(l-1079)) + (*(l+1081))) >> 2;
            *(l+1081) = *p, *p = a0;
			a0 = (*(l-2) + a2 + (*(l-1081)) + (*(l+1079))) >> 2;
			p -= 2, *(l+1079) = *p, *p = a0;
			a0 = *(l-1);
			a2 = (a0 + (*(l+1)) + (*(l-1080)) + (*(l+1080))) >> 2;
			p += 1, *(l+1080) = *p, *p = a2;
			l -= 3, p -= 2, a2 = *l;
        }
        // 1078 % 3 = 1
        a0 = (a0 + a2 + (*(l-1079)) + (*(l+1081))) >> 2;
        *(l+1081) = *p, *p = a0;
        l -= 3, p = &tmp[M-2];
    }
    for (j=M-2; j> 0; j--) A[1][j] = tmp[j];
}

去除 a1 ，为了仍然能够减少同一行每个地址的重复引用，对于每个循环展开的 3x1，原本采用的是逐个计算 1、2、3 的顺序，现在采用按 1、3、2 的顺序计算，具体细节见代码。

work7（减少 1 个临时变量）

void work7(int A[N][M])
{
    register int a0, *l, *p, i, j;
	for (j=0; j< M; j++) tmp[j] = A[N-1][j];
    l = &A[N-2][M-3], p = &tmp[M-2];
    for (i=N-2; i> 0; i--) {
        a0 = *(l+2);
        for (j=M-2; j> 1; j-=3) {
            a0 = ((*l) + a0 + (*(l-1079)) + (*(l+1081))) >> 2;
            *(l+1081) = *p, *p = a0;
			a0 = ((*l) + (*(l-2)) + (*(l-1081)) + (*(l+1079))) >> 2;
			p -= 2, *(l+1079) = *p, *p = a0;
			a0 = ((*(l-1)) + (*(l+1)) + (*(l-1080)) + (*(l+1080))) >> 2;
			p += 1, *(l+1080) = *p, *p = a0;
			l -= 3, p -= 2, a0 = *(l+2);
        }
        // 1078 % 3 = 1
        a0 = ((*l) + a0 + (*(l-1079)) + (*(l+1081))) >> 2;
        *(l+1081) = *p, *p = a0;
        l -= 3, p = &tmp[M-2];
    }
    for (j=M-2; j> 0; j--) A[1][j] = tmp[j];
}

去除 a2，目的主要是使这 5 个临时变量均存储在寄存器中。

work8（关于 cache）

结论：按当前的二重循环方式遍历，已经几乎达到 cache 优化的上限。以下试说明之：
通过 CPU-Z 查询一级数据缓存（4个）的参数：
data cache: 48-KB, 12-way set associative, 64-byte line size
计算得知：cache set : 2^6 = 64 , cache line : 12 , block : 2^6 = 64
故 b = 6, s = 6, t = m-b-s
此实验中 img[1920][1080] 数组的数据类型为 int，每个元素占 4 个字节，故每个块可以存储 16 个元素；每遍历过 16 个元素时产生一次 cache miss，且当前元素所在的新放置的块所映射的组恰好是上一个元素所映射的组的下一个（模64意义下）。当直接按地址升序访问数组（逐行遍历）时，由于 1080 / 16 = 67，每次遍历完数组的一行，对于 cache 的放置情况如图：

以 work4 函数为例，对于当前的位置 [i, j]，访问了数组 [i-1, j], [i+1, j] 位置（[i, j-1], [i, j+1] 位于临时变量中，其理想状态下由寄存器存储，或者位于栈内，依然具有良好的空间局部性）：我们可以认为 [i-1, j] 通过之前的访问和放置，目前仍在 cache 中（因为相对同一组内的其他 11 行，它是最近加入的，假定它更不容易被替换，况且有 4 个 L1）；只有对于 [i+1, j]，每遍历 16 个元素发生一次 miss，故总共不命中次数是 \(\frac{NM}{16}\)，很显然这就是不命中次数的下界。故逐行遍历数组已经达到了 cache 优化的上限。

4.3 优化前后的性能测试

见 4.1，利用 clock 函数得到程序运行前后的时间，相减得之。重复执行 T=1000 次并取平均，尽量消除运行程序时的效率波动因素。
整理得下表如下：

*v : Number of declared temporary variable

4.4 面向泰山服务器优化后的程序与测试结果

需要考虑的服务器cpu的cache大小：

stu_120L021220@node210:~$ lscpu
Architecture:                    aarch64
CPU op-mode(s):                  64-bit
...
L1d cache:                       6 MiB
L1i cache:                       6 MiB
L2 cache:                        48 MiB
L3 cache:                        192 MiB

这比本地的 L1 大多了... 因此我们仍按当前的二重循环方式遍历。
无 O2 编译优化：

stu_120L021220@node210:~/lab3$ gcc lab3.cpp
stu_120L021220@node210:~/lab3$ ./a.out
func0 :  23.860822 ms
func1 :  23.003251 ms
func2 :  14.036629 ms
func3 :  3.632750 ms
func4 :  3.323216 ms
func5 :  4.140020 ms
func6 :  4.201077 ms
func7 :  4.353720 ms
-------------------

加 O2 编译优化：

stu_120L021220@node210:~/lab3$ gcc lab3-o2.cpp
stu_120L021220@node210:~/lab3$ ./a.out
func0 :  3.219893 ms
func1 :  4.054577 ms
func2 :  2.793952 ms
func3 :  2.722771 ms
func4 :  2.681660 ms
func5 :  2.693084 ms
func6 :  2.808298 ms
func7 :  3.206895 ms
-------------------

可见无 O2 下 work4 的运行效果接近 O2 下的初始程序，而 O2 下对程序的优化其效果也比较显而易见。

4.5 还可以采取的进一步的优化方案

SSE/AVX：通过使用寄存器来进行并行加速；
多进程优化：fork，每个进程负责各自的工作任务，通过mmap共享内存或磁盘等进行交互。