C# 使用SIMD向量类型加速浮点数组求和运算（3）：循环展开 ...

作者： zyl910

目录

• 一、背景
  
  
  
  • 1.1 循环展开简介
    
  • 1.2 测试准备
    
  
  
• 二、在C#中使用
  
  
  
  • 2.1 对基础算法做循环展开
    
    
    
    • 2.1.1 测试结果:
      
    
    
  • 2.2 对 Vector4 版算法做循环展开
    
    
    
    • 2.2.1 测试结果:
      
    
    
  • 2.3 对 Vector 版算法做循环展开
    
    
    
    • 2.3.1 测试结果:
      
    
    
  • 2.4 对 Avx版算法做循环展开
    
    
    
    • 2.4.1 测试结果:
      
    
    
  • 2.5 对 Avx版算法做循环展开16次
    
    
    
    • 2.5.1 测试结果:
      
    
    
  • 2.6 尝试用数组来存储循环展开的临时变量
    
    
    
    • 2.6.1 测试结果:
      
    
    
  • 2.7 尝试用栈数组来减少相关性
    
    
    
    • 2.7.1 测试结果:
      
    
    
  • 2.8 测试结果汇总
    
  
  
• 三、在C++中使用
  
  
  
  • 3.1 修改代码
    
  • 3.2 测试结果
    
  
  
• 四、小结
  
• 参考文献
  

一、背景

先前的2篇文章，说了向量类型的类型选择问题。本文讨论一个使用方面的问题——循环展开。
现在的CPU采用了流水线、超标量等机制来提高运算性能。如果完全是顺序代码，那么流水线的效果会非常好。
但是程序中不可避免的需要 分支 与 循环来处理各种复杂的逻辑。分支与循环会被编译为跳转指令，而跳转指令会导致CPU流水线失效，对性能的影响很大。虽然现代处理器增加了分支预测技术，但总会有预测失败的概率。
尤其是在使用向量类型进行SIMD运算时，因向量类型仅尽可能榨干CPU内部的ALU（算术逻辑单元），于是在跳转时的性能损失更大。
故在使用向量类型处理大规模数学计算时，应尽可能的避免分支与循环。
对于分支，最好尽量将分支挪到内循环外。若是内循环中必须的分支，可尽量用位掩码等办法来写无分支代码。
对于循环，一般可使用循环展开技术，来避免短的循环。
1.1 循环展开简介

摘录——

1. 循环展开（Loop unrolling）技术是一种提升程序执行速度的非常有效的优化方法，它可以由程序员手工编写，也可由编译器自动优化。循环展开的本质是，利用CPU指令级并行，来降低循环的开销，当然，同时也有利于指令流水线的高效调度。
2. ……
4. 循环展开的优点：
5. 第一，减少了分支预测失败的可能性。
6. 第二，增加了循环体内语句并发执行的可能性，当然，这需要循环体内各语句不存在数据相关性。
8. 循环展开的缺点：
9. 第一，造成代码膨胀，导致ELF文件（或Windows PE文件）尺寸增大。
10. 第二，代码可读性显著降低，前一个人写的循环展开代码，很可能被不熟悉的后续维护人员改回去。

复制代码

1.2 测试准备

注意，循环展开提高的是流水线性能，对小循环效果明显。此时分支造成的延时，大多与内循环的运算耗时差不多。
对于有些复杂的大循环，内循环的运算耗时已经很大了，而分支造成的延时仍是常数值，比例下降了很多。此时循环展开的收益就少了。
由于循环展开是程序员手工编写的，故必须在编码前就确定好展开次数。
本文就来探讨一下大多数时候的展开次数选择。
展开2倍的话，性能最多为原来2倍，即大多数情况下只有1倍多的性能提升，提升不大。
展开2倍的话，性能最多为原来4倍。区间大了，很多时候能达到2倍以上的提升。
故一开始可以用展开4倍来测试。下面将进行测试。
测试电脑的配置信息为：lntel(R) Core(TM) i5-8250U CPU @ 1.60GHz、Windows 10。
二、在C#中使用

为了对比测试 Avx指令的效果，故可在 BenchmarkVectorCore30 工程里进行测试。因是64位操作系统，故选取 x64、Release版的测试结果.
2.1 对基础算法做循环展开

回顾一下基础算法:

1. private static float SumBase(float[] src, int count, int loops) {
2.     float rt = 0; // Result.
3.     for (int j=0; j< loops; ++j) {
4.         for(int i=0; i< count; ++i) {
5.             rt += src[i];
6.         }
7.     }
8.     return rt;
9. }

复制代码

改为循环展开4倍后，代码为:

1. private static float SumBaseU4(float[] src, int count, int loops) {
2.     float rt = 0; // Result.
3.     float rt1 = 0;
4.     float rt2 = 0;
5.     float rt3 = 0;
6.     int nBlockWidth = 4; // Block width.
7.     int cntBlock = count / nBlockWidth; // Block count.
8.     int cntRem = count % nBlockWidth; // Remainder count.
9.     int p; // Index for src data.
10.     int i;
11.     for (int j = 0; j < loops; ++j) {
12.         p = 0;
13.         // Block processs.
14.         for (i = 0; i < cntBlock; ++i) {
15.             rt += src[p];
16.             rt1 += src[p + 1];
17.             rt2 += src[p + 2];
18.             rt3 += src[p + 3];
19.             p += nBlockWidth;
20.         }
21.         // Remainder processs.
22.         //p = cntBlock * nBlockWidth;
23.         for (i = 0; i < cntRem; ++i) {
24.             rt += src[p + i];
25.         }
26.     }
27.     // Reduce.
28.     rt = rt + rt1 + rt2 + rt3;
29.     return rt;
30. }

复制代码

之前内循环只处理1个数据，现在内循环处理了4个数据。
注意内循环在处理者4个数据时，并不是直接将结果全部累加到 rt 变量，而是使用新增的 rt1、rt2、rt3 变量来临时存储累加值。这是为了消除变量之间的相关性，因为变量之间的相关性会影响流水线性能，故分别使用独立的变量就好了。
最后在 Reduce 阶段，将 rt1、rt2、rt3 的值累加到 rt。
2.1.1 测试结果:

测试结果摘录如下：

1. SumBase:        6.871948E+10    # msUsed=4938, MFLOPS/s=829.485621709194
2. SumBaseU4:      2.748779E+11    # msUsed=1875, MFLOPS/s=2184.5333333333333, scale=2.6336

复制代码

可以发现，基础算法使用4倍循环展开后，性能是原先的 2.6336 倍。
2.2 对 Vector4 版算法做循环展开

回顾一下Vector4 版算法:

1. private static float SumVector4(float[] src, int count, int loops) {
2.     float rt = 0; // Result.
3.     const int VectorWidth = 4;
4.     int nBlockWidth = VectorWidth; // Block width.
5.     int cntBlock = count / nBlockWidth; // Block count.
6.     int cntRem = count % nBlockWidth; // Remainder count.
7.     Vector4 vrt = Vector4.Zero; // Vector result.
8.     int p; // Index for src data.
9.     int i;
10.     // Load.
11.     Vector4[] vsrc = new Vector4[cntBlock]; // Vector src.
12.     p = 0;
13.     for (i = 0; i < vsrc.Length; ++i) {
14.         vsrc[i] = new Vector4(src[p], src[p + 1], src[p + 2], src[p + 3]);
15.         p += VectorWidth;
16.     }
17.     // Body.
18.     for (int j = 0; j < loops; ++j) {
19.         // Vector processs.
20.         for (i = 0; i < cntBlock; ++i) {
21.             // Equivalent to scalar model: rt += src[i];
22.             vrt += vsrc[i]; // Add.
23.         }
24.         // Remainder processs.
25.         p = cntBlock * nBlockWidth;
26.         for (i = 0; i < cntRem; ++i) {
27.             rt += src[p + i];
28.         }
29.     }
30.     // Reduce.
31.     rt += vrt.X + vrt.Y + vrt.Z + vrt.W;
32.     return rt;
33. }

复制代码

改为循环展开4倍后，代码为:

1. private static float SumVector4U4(float[] src, int count, int loops) {
2.     float rt = 0; // Result.
3.     const int LoopUnrolling = 4;
4.     const int VectorWidth = 4;
5.     int nBlockWidth = VectorWidth * LoopUnrolling; // Block width.
6.     int cntBlock = count / nBlockWidth; // Block count.
7.     int cntRem = count % nBlockWidth; // Remainder count.
8.     Vector4 vrt = Vector4.Zero; // Vector result.
9.     Vector4 vrt1 = Vector4.Zero;
10.     Vector4 vrt2 = Vector4.Zero;
11.     Vector4 vrt3 = Vector4.Zero;
12.     int p; // Index for src data.
13.     int i;
14.     // Load.
15.     Vector4[] vsrc = new Vector4[count / VectorWidth]; // Vector src.
16.     p = 0;
17.     for (i = 0; i < vsrc.Length; ++i) {
18.         vsrc[i] = new Vector4(src[p], src[p + 1], src[p + 2], src[p + 3]);
19.         p += VectorWidth;
20.     }
21.     // Body.
22.     for (int j = 0; j < loops; ++j) {
23.         p = 0;
24.         // Vector processs.
25.         for (i = 0; i < cntBlock; ++i) {
26.             vrt += vsrc[p]; // Add.
27.             vrt1 += vsrc[p + 1];
28.             vrt2 += vsrc[p + 2];
29.             vrt3 += vsrc[p + 3];
30.             p += LoopUnrolling;
31.         }
32.         // Remainder processs.
33.         p = cntBlock * nBlockWidth;
34.         for (i = 0; i < cntRem; ++i) {
35.             rt += src[p + i];
36.         }
37.     }
38.     // Reduce.
39.     vrt = vrt + vrt1 + vrt2 + vrt3;
40.     rt += vrt.X + vrt.Y + vrt.Z + vrt.W;
41.     return rt;
42. }

复制代码

跟刚才的办法一样，使用新增的 rt1、rt2、rt3 变量来临时存储累加值，消除变量之间的相关性。
最后在 Reduce 阶段，将 vrt1、vrt2、vrt3 的值累加到 vrt。
2.2.1 测试结果:

测试结果摘录如下：

1. SumBase:        6.871948E+10    # msUsed=4938, MFLOPS/s=829.485621709194
2. SumBaseU4:      2.748779E+11    # msUsed=1875, MFLOPS/s=2184.5333333333333, scale=2.6336SumVector4:     2.748779E+11    # msUsed=1218, MFLOPS/s=3362.8899835796387, scale=4.054187192118227SumVector4U4:   1.0995116E+12   # msUsed=532, MFLOPS/s=7699.248120300752, scale=9.281954887218046

复制代码

SumVector4U4对比基础算法（SumBase），性能倍数是 9.281954887218046。
SumVector4U4对比未循环展开的算法（SumVector4），倍数是 9.281954887218046/4.054187192118227=2.2894736842105263092984587836542
2.3 对 Vector 版算法做循环展开

回顾一下 Vector 版算法:

1. private static float SumVectorT(float[] src, int count, int loops) {
2.     float rt = 0; // Result.
3.     int VectorWidth = Vector<float>.Count; // Block width.
4.     int nBlockWidth = VectorWidth; // Block width.
5.     int cntBlock = count / nBlockWidth; // Block count.
6.     int cntRem = count % nBlockWidth; // Remainder count.
7.     Vector<float> vrt = Vector<float>.Zero; // Vector result.
8.     int p; // Index for src data.
9.     int i;
10.     // Load.
11.     Vector<float>[] vsrc = new Vector<float>[cntBlock]; // Vector src.
12.     p = 0;
13.     for (i = 0; i < vsrc.Length; ++i) {
14.         vsrc[i] = new Vector<float>(src, p);
15.         p += VectorWidth;
16.     }
17.     // Body.
18.     for (int j = 0; j < loops; ++j) {
19.         // Vector processs.
20.         for (i = 0; i < cntBlock; ++i) {
21.             vrt += vsrc[i]; // Add.
22.         }
23.         // Remainder processs.
24.         p = cntBlock * nBlockWidth;
25.         for (i = 0; i < cntRem; ++i) {
26.             rt += src[p + i];
27.         }
28.     }
29.     // Reduce.
30.     for (i = 0; i < VectorWidth; ++i) {
31.         rt += vrt[i];
32.     }
33.     return rt;
34. }

复制代码

改为循环展开4倍后，代码为:

1. private static float SumVectorTU4(float[] src, int count, int loops) {
2.     float rt = 0; // Result.
3.     const int LoopUnrolling = 4;
4.     int VectorWidth = Vector<float>.Count; // Block width.
5.     int nBlockWidth = VectorWidth * LoopUnrolling; // Block width.
6.     int cntBlock = count / nBlockWidth; // Block count.
7.     int cntRem = count % nBlockWidth; // Remainder count.
8.     Vector<float> vrt = Vector<float>.Zero; // Vector result.
9.     Vector<float> vrt1 = Vector<float>.Zero;
10.     Vector<float> vrt2 = Vector<float>.Zero;
11.     Vector<float> vrt3 = Vector<float>.Zero;
12.     int p; // Index for src data.
13.     int i;
14.     // Load.
15.     Vector<float>[] vsrc = new Vector<float>[count / VectorWidth]; // Vector src.
16.     p = 0;
17.     for (i = 0; i < vsrc.Length; ++i) {
18.         vsrc[i] = new Vector<float>(src, p);
19.         p += VectorWidth;
20.     }
21.     // Body.
22.     for (int j = 0; j < loops; ++j) {
23.         p = 0;
24.         // Vector processs.
25.         for (i = 0; i < cntBlock; ++i) {
26.             vrt += vsrc[p]; // Add.
27.             vrt1 += vsrc[p + 1];
28.             vrt2 += vsrc[p + 1];
29.             vrt3 += vsrc[p + 1];
30.             p += LoopUnrolling;
31.         }
32.         // Remainder processs.
33.         p = cntBlock * nBlockWidth;
34.         for (i = 0; i < cntRem; ++i) {
35.             rt += src[p + i];
36.         }
37.     }
38.     // Reduce.
39.     vrt = vrt + vrt1 + vrt2 + vrt3;
40.     for (i = 0; i < VectorWidth; ++i) {
41.         rt += vrt[i];
42.     }
43.     return rt;
44. }

复制代码

跟刚才的办法一样，使用新增的 rt1、rt2、rt3 变量来临时存储累加值，消除变量之间的相关性。
最后在 Reduce 阶段，将 vrt1、vrt2、vrt3 的值累加到 vrt。
2.3.1 测试结果:

测试结果摘录如下：

1. SumBase:        6.871948E+10    # msUsed=4938, MFLOPS/s=829.485621709194
2. SumBaseU4:      2.748779E+11    # msUsed=1875, MFLOPS/s=2184.5333333333333, scale=2.6336SumVectorT:     5.497558E+11    # msUsed=609, MFLOPS/s=6725.7799671592775, scale=8.108374384236454SumVectorTU4:   2.1990233E+12   # msUsed=203, MFLOPS/s=20177.339901477833, scale=24.32512315270936

复制代码

SumVectorTU4对比基础算法（SumBase），性能倍数是 24.32512315270936。
SumVectorTU4对比未循环展开的算法（SumVectorT），倍数是 24.32512315270936/8.108374384236454=2.9999999999999997533414337788579
初步发现 Vector循环展开（2.9999）带来的性能提升, 比VectorT循环展开（2.2894）更高一些。
2.4 对 Avx版算法做循环展开

先前分别尝试用 数组、Span、指针 的办法来操纵数据、使用Avx指令集。现在对这3种办法，均写一套循环展开4次的代码:

1. /// <summary>
2. /// Sum - Vector AVX.
3. /// </summary>
4. /// <param name="src">Soure array.</param>
5. /// <param name="count">Soure array count.</param>
6. /// <param name="loops">Benchmark loops.</param>
7. /// <returns>Return the sum value.</returns>
8. private static float SumVectorAvx(float[] src, int count, int loops) {
9. #if Allow_Intrinsics
10.     float rt = 0; // Result.
11.     //int VectorWidth = 32 / 4; // sizeof(__m256) / sizeof(float);
12.     int VectorWidth = Vector256<float>.Count; // Block width.
13.     int nBlockWidth = VectorWidth; // Block width.
14.     int cntBlock = count / nBlockWidth; // Block count.
15.     int cntRem = count % nBlockWidth; // Remainder count.
16.     Vector256<float> vrt = Vector256<float>.Zero; // Vector result.
17.     int p; // Index for src data.
18.     int i;
19.     // Load.
20.     Vector256<float>[] vsrc = new Vector256<float>[cntBlock]; // Vector src.
21.     p = 0;
22.     for (i = 0; i < cntBlock; ++i) {
23.         vsrc[i] = Vector256.Create(src[p], src[p + 1], src[p + 2], src[p + 3], src[p + 4], src[p + 5], src[p + 6], src[p + 7]); // Load.
24.         p += VectorWidth;
25.     }
26.     // Body.
27.     for (int j = 0; j < loops; ++j) {
28.         // Vector processs.
29.         for (i = 0; i < cntBlock; ++i) {
30.             vrt = Avx.Add(vrt, vsrc[i]);    // Add. vrt += vsrc[i];
31.         }
32.         // Remainder processs.
33.         p = cntBlock * nBlockWidth;
34.         for (i = 0; i < cntRem; ++i) {
35.             rt += src[p + i];
36.         }
37.     }
38.     // Reduce.
39.     for (i = 0; i < VectorWidth; ++i) {
40.         rt += vrt.GetElement(i);
41.     }
42.     return rt;
43. #else
44.     throw new NotSupportedException();
45. #endif
46. }
48. /// <summary>
49. /// Sum - Vector AVX - Loop unrolling *4.
50. /// </summary>
51. /// <param name="src">Soure array.</param>
52. /// <param name="count">Soure array count.</param>
53. /// <param name="loops">Benchmark loops.</param>
54. /// <returns>Return the sum value.</returns>
55. private static float SumVectorAvxU4(float[] src, int count, int loops) {
56. #if Allow_Intrinsics
57.     float rt = 0; // Result.
58.     const int LoopUnrolling = 4;
59.     int VectorWidth = Vector256<float>.Count; // Block width.
60.     int nBlockWidth = VectorWidth * LoopUnrolling; // Block width.
61.     int cntBlock = count / nBlockWidth; // Block count.
62.     int cntRem = count % nBlockWidth; // Remainder count.
63.     Vector256<float> vrt = Vector256<float>.Zero; // Vector result.
64.     Vector256<float> vrt1 = Vector256<float>.Zero;
65.     Vector256<float> vrt2 = Vector256<float>.Zero;
66.     Vector256<float> vrt3 = Vector256<float>.Zero;
67.     int p; // Index for src data.
68.     int i;
69.     // Load.
70.     Vector256<float>[] vsrc = new Vector256<float>[count / VectorWidth]; // Vector src.
71.     p = 0;
72.     for (i = 0; i < vsrc.Length; ++i) {
73.         vsrc[i] = Vector256.Create(src[p], src[p + 1], src[p + 2], src[p + 3], src[p + 4], src[p + 5], src[p + 6], src[p + 7]); // Load.
74.         p += VectorWidth;
75.     }
76.     // Body.
77.     for (int j = 0; j < loops; ++j) {
78.         p = 0;
79.         // Vector processs.
80.         for (i = 0; i < cntBlock; ++i) {
81.             vrt = Avx.Add(vrt, vsrc[p]);    // Add. vrt += vsrc[p];
82.             vrt1 = Avx.Add(vrt1, vsrc[p + 1]);
83.             vrt2 = Avx.Add(vrt2, vsrc[p + 2]);
84.             vrt3 = Avx.Add(vrt3, vsrc[p + 3]);
85.             p += LoopUnrolling;
86.         }
87.         // Remainder processs.
88.         p = cntBlock * nBlockWidth;
89.         for (i = 0; i < cntRem; ++i) {
90.             rt += src[p + i];
91.         }
92.     }
93.     // Reduce.
94.     vrt = Avx.Add(Avx.Add(vrt, vrt1), Avx.Add(vrt2, vrt3)); // vrt = vrt + vrt1 + vrt2 + vrt3;
95.     for (i = 0; i < VectorWidth; ++i) {
96.         rt += vrt.GetElement(i);
97.     }
98.     return rt;
99. #else
100.     throw new NotSupportedException();
101. #endif
102. }
104. /// <summary>
105. /// Sum - Vector AVX - Span.
106. /// </summary>
107. /// <param name="src">Soure array.</param>
108. /// <param name="count">Soure array count.</param>
109. /// <param name="loops">Benchmark loops.</param>
110. /// <returns>Return the sum value.</returns>
111. private static float SumVectorAvxSpan(float[] src, int count, int loops) {
112. #if Allow_Intrinsics
113.     float rt = 0; // Result.
114.     int VectorWidth = Vector256<float>.Count; // Block width.
115.     int nBlockWidth = VectorWidth; // Block width.
116.     int cntBlock = count / nBlockWidth; // Block count.
117.     int cntRem = count % nBlockWidth; // Remainder count.
118.     Vector256<float> vrt = Vector256<float>.Zero; // Vector result.
119.     int p; // Index for src data.
120.     ReadOnlySpan<Vector256<float>> vsrc; // Vector src.
121.     int i;
122.     // Body.
123.     for (int j = 0; j < loops; ++j) {
124.         // Vector processs.
125.         vsrc = System.Runtime.InteropServices.MemoryMarshal.Cast<float, Vector256<float> >(new Span<float>(src)); // Reinterpret cast. `float*` to `Vector256<float>*`.
126.         for (i = 0; i < cntBlock; ++i) {
127.             vrt = Avx.Add(vrt, vsrc[i]);    // Add. vrt += vsrc[i];
128.         }
129.         // Remainder processs.
130.         p = cntBlock * nBlockWidth;
131.         for (i = 0; i < cntRem; ++i) {
132.             rt += src[p + i];
133.         }
134.     }
135.     // Reduce.
136.     for (i = 0; i < VectorWidth; ++i) {
137.         rt += vrt.GetElement(i);
138.     }
139.     return rt;
140. #else
141.     throw new NotSupportedException();
142. #endif
143. }
145. /// <summary>
146. /// Sum - Vector AVX - Span - Loop unrolling *4.
147. /// </summary>
148. /// <param name="src">Soure array.</param>
149. /// <param name="count">Soure array count.</param>
150. /// <param name="loops">Benchmark loops.</param>
151. /// <returns>Return the sum value.</returns>
152. private static float SumVectorAvxSpanU4(float[] src, int count, int loops) {
153. #if Allow_Intrinsics
154.     float rt = 0; // Result.
155.     const int LoopUnrolling = 4;
156.     int VectorWidth = Vector256<float>.Count; // Block width.
157.     int nBlockWidth = VectorWidth * LoopUnrolling; // Block width.
158.     int cntBlock = count / nBlockWidth; // Block count.
159.     int cntRem = count % nBlockWidth; // Remainder count.
160.     Vector256<float> vrt = Vector256<float>.Zero; // Vector result.
161.     Vector256<float> vrt1 = Vector256<float>.Zero;
162.     Vector256<float> vrt2 = Vector256<float>.Zero;
163.     Vector256<float> vrt3 = Vector256<float>.Zero;
164.     int p; // Index for src data.
165.     ReadOnlySpan<Vector256<float>> vsrc; // Vector src.
166.     int i;
167.     // Body.
168.     for (int j = 0; j < loops; ++j) {
169.         p = 0;
170.         // Vector processs.
171.         vsrc = System.Runtime.InteropServices.MemoryMarshal.Cast<float, Vector256<float>>(new Span<float>(src)); // Reinterpret cast. `float*` to `Vector256<float>*`.
172.         for (i = 0; i < cntBlock; ++i) {
173.             vrt = Avx.Add(vrt, vsrc[p]);    // Add. vrt += vsrc[p];
174.             vrt1 = Avx.Add(vrt1, vsrc[p + 1]);
175.             vrt2 = Avx.Add(vrt2, vsrc[p + 2]);
176.             vrt3 = Avx.Add(vrt3, vsrc[p + 3]);
177.             p += LoopUnrolling;
178.         }
179.         // Remainder processs.
180.         p = cntBlock * nBlockWidth;
181.         for (i = 0; i < cntRem; ++i) {
182.             rt += src[p + i];
183.         }
184.     }
185.     // Reduce.
186.     vrt = Avx.Add(Avx.Add(vrt, vrt1), Avx.Add(vrt2, vrt3)); // vrt = vrt + vrt1 + vrt2 + vrt3;
187.     for (i = 0; i < VectorWidth; ++i) {
188.         rt += vrt.GetElement(i);
189.     }
190.     return rt;
191. #else
192.     throw new NotSupportedException();
193. #endif
194. }
196. /// <summary>
197. /// Sum - Vector AVX - Ptr.
198. /// </summary>
199. /// <param name="src">Soure array.</param>
200. /// <param name="count">Soure array count.</param>
201. /// <param name="loops">Benchmark loops.</param>
202. /// <returns>Return the sum value.</returns>
203. private static float SumVectorAvxPtr(float[] src, int count, int loops) {
204. #if Allow_Intrinsics && UNSAFE
205.     unsafe {
206.         float rt = 0; // Result.
207.         int VectorWidth = Vector256<float>.Count; // Block width.
208.         int nBlockWidth = VectorWidth; // Block width.
209.         int cntBlock = count / nBlockWidth; // Block count.
210.         int cntRem = count % nBlockWidth; // Remainder count.
211.         Vector256<float> vrt = Vector256<float>.Zero; // Vector result.
212.         Vector256<float> vload;
213.         float* p; // Pointer for src data.
214.         int i;
215.         // Body.
216.         fixed(float* p0 = &src[0]) {
217.             for (int j = 0; j < loops; ++j) {
218.                 p = p0;
219.                 // Vector processs.
220.                 for (i = 0; i < cntBlock; ++i) {
221.                     vload = Avx.LoadVector256(p);    // Load. vload = *(*__m256)p;
222.                     vrt = Avx.Add(vrt, vload);    // Add. vrt += vsrc[i];
223.                     p += nBlockWidth;
224.                 }
225.                 // Remainder processs.
226.                 for (i = 0; i < cntRem; ++i) {
227.                     rt += p[i];
228.                 }
229.             }
230.         }
231.         // Reduce.
232.         for (i = 0; i < VectorWidth; ++i) {
233.             rt += vrt.GetElement(i);
234.         }
235.         return rt;
236.     }
237. #else
238.     throw new NotSupportedException();
239. #endif
240. }
242. /// <summary>
243. /// Sum - Vector AVX - Ptr - Loop unrolling *4.
244. /// </summary>
245. /// <param name="src">Soure array.</param>
246. /// <param name="count">Soure array count.</param>
247. /// <param name="loops">Benchmark loops.</param>
248. /// <returns>Return the sum value.</returns>
249. private static float SumVectorAvxPtrU4(float[] src, int count, int loops) {
250. #if Allow_Intrinsics && UNSAFE
251.     unsafe {
252.         float rt = 0; // Result.
253.         const int LoopUnrolling = 4;
254.         int VectorWidth = Vector256<float>.Count; // Block width.
255.         int nBlockWidth = VectorWidth * LoopUnrolling; // Block width.
256.         int cntBlock = count / nBlockWidth; // Block count.
257.         int cntRem = count % nBlockWidth; // Remainder count.
258.         Vector256<float> vrt = Vector256<float>.Zero; // Vector result.
259.         Vector256<float> vrt1 = Vector256<float>.Zero;
260.         Vector256<float> vrt2 = Vector256<float>.Zero;
261.         Vector256<float> vrt3 = Vector256<float>.Zero;
262.         Vector256<float> vload;
263.         Vector256<float> vload1;
264.         Vector256<float> vload2;
265.         Vector256<float> vload3;
266.         float* p; // Pointer for src data.
267.         int i;
268.         // Body.
269.         fixed (float* p0 = &src[0]) {
270.             for (int j = 0; j < loops; ++j) {
271.                 p = p0;
272.                 // Vector processs.
273.                 for (i = 0; i < cntBlock; ++i) {
274.                     vload = Avx.LoadVector256(p);    // Load. vload = *(*__m256)p;
275.                     vload1 = Avx.LoadVector256(p + VectorWidth * 1);
276.                     vload2 = Avx.LoadVector256(p + VectorWidth * 2);
277.                     vload3 = Avx.LoadVector256(p + VectorWidth * 3);
278.                     vrt = Avx.Add(vrt, vload);    // Add. vrt += vsrc[i];
279.                     vrt1 = Avx.Add(vrt1, vload1);
280.                     vrt2 = Avx.Add(vrt2, vload2);
281.                     vrt3 = Avx.Add(vrt3, vload3);
282.                     p += nBlockWidth;
283.                 }
284.                 // Remainder processs.
285.                 for (i = 0; i < cntRem; ++i) {
286.                     rt += p[i];
287.                 }
288.             }
289.         }
290.         // Reduce.
291.         vrt = Avx.Add(Avx.Add(vrt, vrt1), Avx.Add(vrt2, vrt3)); // vrt = vrt + vrt1 + vrt2 + vrt3;
292.         for (i = 0; i < VectorWidth; ++i) {
293.             rt += vrt.GetElement(i);
294.         }
295.         return rt;
296.     }
297. #else
298.     throw new NotSupportedException();
299. #endif
300. }

复制代码

2.4.1 测试结果:

测试结果摘录如下：

1. SumBase:        6.871948E+10    # msUsed=4938, MFLOPS/s=829.485621709194
2. SumBaseU4:      2.748779E+11    # msUsed=1875, MFLOPS/s=2184.5333333333333, scale=2.6336SumVectorAvx:   5.497558E+11    # msUsed=609, MFLOPS/s=6725.7799671592775, scale=8.108374384236454SumVectorAvxSpan:       5.497558E+11    # msUsed=625, MFLOPS/s=6553.6, scale=7.9008SumVectorAvxPtr:        5.497558E+11    # msUsed=610, MFLOPS/s=6714.754098360656, scale=8.095081967213115SumVectorAvxU4: 2.1990233E+12   # msUsed=328, MFLOPS/s=12487.80487804878, scale=15.054878048780488SumVectorAvxSpanU4:     2.1990233E+12   # msUsed=312, MFLOPS/s=13128.205128205129, scale=15.826923076923078SumVectorAvxPtrU4:      2.1990233E+12   # msUsed=157, MFLOPS/s=26089.171974522294, scale=31.452229299363058

复制代码

未做循环展开时，这3钟办法的性能拉不开差距，都是8倍左右。
而现在用了循环展开后，数组版（SumVectorAvxU4）、Span版（SumVectorAvxSpanU4）只有15倍左右的性能提升。而指针版有 31倍性能提升，是 数组版、Span版 的2倍。
可能是因为指针更贴近底层硬件、更易于编译器优化。故当使用内在函数时，推荐优先使用指针。
SumVectorAvxPtrU4 对比基础算法（SumBase），性能倍数是 31.452229299363058。
SumVectorAvxPtrU4 对比未循环展开的算法（SumVectorAvxPtr），倍数是 31.452229299363058/8.095081967213115=3.8853503184713375449974589366746。
2.5 对 Avx版算法做循环展开16次

刚才尝试了4倍循环展开，故理论上限是4倍。而SumVectorAvxPtrU4版有约 3.8853 倍性能提升，故可考虑进一步加大，于是可测试一下 4*4=16 次的循环展开。
将 SumVectorAvxPtr 改造为循环展开16次的，代码如下:

1. private static float SumVectorAvxPtrU16(float[] src, int count, int loops) {
2. #if Allow_Intrinsics && UNSAFE
3.     unsafe {
4.         float rt = 0; // Result.
5.         const int LoopUnrolling = 16;
6.         int VectorWidth = Vector256<float>.Count; // Block width.
7.         int nBlockWidth = VectorWidth * LoopUnrolling; // Block width.
8.         int cntBlock = count / nBlockWidth; // Block count.
9.         int cntRem = count % nBlockWidth; // Remainder count.
10.         Vector256<float> vrt = Vector256<float>.Zero; // Vector result.
11.         Vector256<float> vrt1 = Vector256<float>.Zero;
12.         Vector256<float> vrt2 = Vector256<float>.Zero;
13.         Vector256<float> vrt3 = Vector256<float>.Zero;
14.         Vector256<float> vrt4 = Vector256<float>.Zero;
15.         Vector256<float> vrt5 = Vector256<float>.Zero;
16.         Vector256<float> vrt6 = Vector256<float>.Zero;
17.         Vector256<float> vrt7 = Vector256<float>.Zero;
18.         Vector256<float> vrt8 = Vector256<float>.Zero;
19.         Vector256<float> vrt9 = Vector256<float>.Zero;
20.         Vector256<float> vrt10 = Vector256<float>.Zero;
21.         Vector256<float> vrt11 = Vector256<float>.Zero;
22.         Vector256<float> vrt12 = Vector256<float>.Zero;
23.         Vector256<float> vrt13 = Vector256<float>.Zero;
24.         Vector256<float> vrt14 = Vector256<float>.Zero;
25.         Vector256<float> vrt15 = Vector256<float>.Zero;
26.         float* p; // Pointer for src data.
27.         int i;
28.         // Body.
29.         fixed (float* p0 = &src[0]) {
30.             for (int j = 0; j < loops; ++j) {
31.                 p = p0;
32.                 // Vector processs.
33.                 for (i = 0; i < cntBlock; ++i) {
34.                     //vload = Avx.LoadVector256(p);    // Load. vload = *(*__m256)p;
35.                     vrt = Avx.Add(vrt, Avx.LoadVector256(p)); // Add. vrt[k] += *((*__m256)(p)+k);
36.                     vrt1 = Avx.Add(vrt1, Avx.LoadVector256(p + VectorWidth * 1));
37.                     vrt2 = Avx.Add(vrt2, Avx.LoadVector256(p + VectorWidth * 2));
38.                     vrt3 = Avx.Add(vrt3, Avx.LoadVector256(p + VectorWidth * 3));
39.                     vrt4 = Avx.Add(vrt4, Avx.LoadVector256(p + VectorWidth * 4));
40.                     vrt5 = Avx.Add(vrt5, Avx.LoadVector256(p + VectorWidth * 5));
41.                     vrt6 = Avx.Add(vrt6, Avx.LoadVector256(p + VectorWidth * 6));
42.                     vrt7 = Avx.Add(vrt7, Avx.LoadVector256(p + VectorWidth * 7));
43.                     vrt8 = Avx.Add(vrt8, Avx.LoadVector256(p + VectorWidth * 8));
44.                     vrt9 = Avx.Add(vrt9, Avx.LoadVector256(p + VectorWidth * 9));
45.                     vrt10 = Avx.Add(vrt10, Avx.LoadVector256(p + VectorWidth * 10));
46.                     vrt11 = Avx.Add(vrt11, Avx.LoadVector256(p + VectorWidth * 11));
47.                     vrt12 = Avx.Add(vrt12, Avx.LoadVector256(p + VectorWidth * 12));
48.                     vrt13 = Avx.Add(vrt13, Avx.LoadVector256(p + VectorWidth * 13));
49.                     vrt14 = Avx.Add(vrt14, Avx.LoadVector256(p + VectorWidth * 14));
50.                     vrt15 = Avx.Add(vrt15, Avx.LoadVector256(p + VectorWidth * 15));
51.                     p += nBlockWidth;
52.                 }
53.                 // Remainder processs.
54.                 for (i = 0; i < cntRem; ++i) {
55.                     rt += p[i];
56.                 }
57.             }
58.         }
59.         // Reduce.
60.         vrt = Avx.Add( Avx.Add( Avx.Add(Avx.Add(vrt, vrt1), Avx.Add(vrt2, vrt3))
61.             , Avx.Add(Avx.Add(vrt4, vrt5), Avx.Add(vrt6, vrt7)) )
62.             , Avx.Add( Avx.Add(Avx.Add(vrt8, vrt9), Avx.Add(vrt10, vrt11))
63.             , Avx.Add(Avx.Add(vrt12, vrt13), Avx.Add(vrt14, vrt15)) ) )
64.         ; // vrt = vrt + vrt1 + vrt2 + vrt3 + ... vrt15;
65.         for (i = 0; i < VectorWidth; ++i) {
66.             rt += vrt.GetElement(i);
67.         }
68.         return rt;
69.     }
70. #else
71.     throw new NotSupportedException();
72. #endif
73. }

复制代码

2.5.1 测试结果:

测试结果摘录如下：

1. SumBase:        6.871948E+10    # msUsed=4938, MFLOPS/s=829.485621709194
2. SumBaseU4:      2.748779E+11    # msUsed=1875, MFLOPS/s=2184.5333333333333, scale=2.6336SumVectorAvxPtr:        5.497558E+11    # msUsed=610, MFLOPS/s=6714.754098360656, scale=8.095081967213115SumVectorAvxPtrU4:      2.1990233E+12   # msUsed=157, MFLOPS/s=26089.171974522294, scale=31.452229299363058SumVectorAvxPtrU16:     8.386202E+12    # msUsed=125, MFLOPS/s=32768, scale=39.504

复制代码

SumVectorAvxPtrU16 对比基础算法（SumBase），性能倍数是 39.504。
SumVectorAvxPtrU16 对比未循环展开的算法（SumVectorAvxPtr），倍数是 39.504/8.095081967213115=4.8799999999999998517618469015796。
SumVectorAvxPtrU16 对比循环展开4次的算法（SumVectorAvxPtrU4），倍数是 39.504/31.452229299363058=1.2559999999999999730384771162414。
从循环展开4次，改为循环展开16次，性能倍数只是从 31倍多，提升到 39 倍多，仅提升 25% 左右。
性能提升的少，但编码麻烦多了。看来循环展开16次的性价比很低，故一般情况下用循环展开4次就行了。
2.6 尝试用数组来存储循环展开的临时变量

使用循环展开N次时，将会导致临时变量数量是非循环展开版的N倍。例如刚才的 SumVectorAvxPtrU16 函数，因循环展开16次，导致临时变量是非循环展开版的16倍，写起了很啰嗦。
这些变量的类型是一样的，放到数组中的话，代码会清晰不少，但会不会影响性能呢？
于是做了一个测试，代码如下：

1. private static float SumVectorAvxPtrU16A(float[] src, int count, int loops) {
2. #if Allow_Intrinsics && UNSAFE
3.     unsafe {
4.         float rt = 0; // Result.
5.         const int LoopUnrolling = 16;
6.         int VectorWidth = Vector256<float>.Count; // Block width.
7.         int nBlockWidth = VectorWidth * LoopUnrolling; // Block width.
8.         int cntBlock = count / nBlockWidth; // Block count.
9.         int cntRem = count % nBlockWidth; // Remainder count.
10.         int i;
11.         //Vector256<float>[] vrt = new Vector256<float>[LoopUnrolling]; // Vector result.
12.         Vector256<float>* vrt = stackalloc Vector256<float>[LoopUnrolling]; ; // Vector result.
13.         for (i = 0; i< LoopUnrolling; ++i) {
14.             vrt[i] = Vector256<float>.Zero;
15.         }
16.         float* p; // Pointer for src data.
17.         // Body.
18.         fixed (float* p0 = &src[0]) {
19.             for (int j = 0; j < loops; ++j) {
20.                 p = p0;
21.                 // Vector processs.
22.                 for (i = 0; i < cntBlock; ++i) {
23.                     //vload = Avx.LoadVector256(p);    // Load. vload = *(*__m256)p;
24.                     vrt[0] = Avx.Add(vrt[0], Avx.LoadVector256(p)); // Add. vrt[k] += *((*__m256)(p)+k);
25.                     vrt[1] = Avx.Add(vrt[1], Avx.LoadVector256(p + VectorWidth * 1));
26.                     vrt[2] = Avx.Add(vrt[2], Avx.LoadVector256(p + VectorWidth * 2));
27.                     vrt[3] = Avx.Add(vrt[3], Avx.LoadVector256(p + VectorWidth * 3));
28.                     vrt[4] = Avx.Add(vrt[4], Avx.LoadVector256(p + VectorWidth * 4));
29.                     vrt[5] = Avx.Add(vrt[5], Avx.LoadVector256(p + VectorWidth * 5));
30.                     vrt[6] = Avx.Add(vrt[6], Avx.LoadVector256(p + VectorWidth * 6));
31.                     vrt[7] = Avx.Add(vrt[7], Avx.LoadVector256(p + VectorWidth * 7));
32.                     vrt[8] = Avx.Add(vrt[8], Avx.LoadVector256(p + VectorWidth * 8));
33.                     vrt[9] = Avx.Add(vrt[9], Avx.LoadVector256(p + VectorWidth * 9));
34.                     vrt[10] = Avx.Add(vrt[10], Avx.LoadVector256(p + VectorWidth * 10));
35.                     vrt[11] = Avx.Add(vrt[11], Avx.LoadVector256(p + VectorWidth * 11));
36.                     vrt[12] = Avx.Add(vrt[12], Avx.LoadVector256(p + VectorWidth * 12));
37.                     vrt[13] = Avx.Add(vrt[13], Avx.LoadVector256(p + VectorWidth * 13));
38.                     vrt[14] = Avx.Add(vrt[14], Avx.LoadVector256(p + VectorWidth * 14));
39.                     vrt[15] = Avx.Add(vrt[15], Avx.LoadVector256(p + VectorWidth * 15));
40.                     p += nBlockWidth;
41.                 }
42.                 // Remainder processs.
43.                 for (i = 0; i < cntRem; ++i) {
44.                     rt += p[i];
45.                 }
46.             }
47.         }
48.         // Reduce.
49.         for (i = 1; i < LoopUnrolling; ++i) {
50.             vrt[0] = Avx.Add(vrt[0], vrt[i]); // vrt[0] += vrt[i]
51.         }
52.         for (i = 0; i < VectorWidth; ++i) {
53.             rt += vrt[0].GetElement(i);
54.         }
55.         return rt;
56.     }
57. #else
58.     throw new NotSupportedException();
59. #endif
60. }

复制代码

2.6.1 测试结果:

测试结果摘录如下：

1. SumBase:        6.871948E+10    # msUsed=4938, MFLOPS/s=829.485621709194
2. SumBaseU4:      2.748779E+11    # msUsed=1875, MFLOPS/s=2184.5333333333333, scale=2.6336SumVectorAvxPtr:        5.497558E+11    # msUsed=610, MFLOPS/s=6714.754098360656, scale=8.095081967213115SumVectorAvxPtrU4:      2.1990233E+12   # msUsed=157, MFLOPS/s=26089.171974522294, scale=31.452229299363058SumVectorAvxPtrU16:     8.386202E+12    # msUsed=125, MFLOPS/s=32768, scale=39.504SumVectorAvxPtrU16A:    8.3862026E+12   # msUsed=187, MFLOPS/s=21903.74331550802, scale=26.406417112299465

复制代码

可以发现 SumVectorAvxPtrU16A 的性能比 SumVectorAvxPtrU16 差。
曾经以为是因为数组是在堆中分配的（new Vector256）引起的，有堆内存分配的开销，且需要多次寻址才能定位变量。
随后改为栈中分配的数组（stackalloc Vector256），且用最贴近硬件的指针来操作，可性能几乎一致。故猜测可能是编译优化时难以将它们优化为寄存器变量。
故在使用循环展开时，临时变量不要用数组来存，还是逐个定义局部变量比较好。
2.7 尝试用栈数组来减少相关性

还尝试了用栈数组来减少相关性，代码如下：
[code]private static float SumVectorAvxPtrUX(float[] src, int count, int loops, int LoopUnrolling) {#if Allow_Intrinsics && UNSAFE    unsafe {        float rt = 0; // Result.        //const int LoopUnrolling = 16;        if (LoopUnrolling