为了完成羞涩教授猛哥的将Multigrid算法的速度发挥到极致的目标， 顺着这些专门研究利用cache数据的德国人们（
DiME）的脚步， 作了一些笔记， 虽然自己不姓后藤…… 中文毕竟比英文好看一些， 要我看我也愿意看中文的。
什么是内存壁垒（Memory Wall） ？
用70年代后期的Multigrid方法解决泊松问题， 一个未知数需要三十个操作， 对于更一般的椭圆方程问题， 一个未知需要一百个左右的操作， 对于一个现在的CPU，一秒钟可以有1Gflop， 如果一切完美， 只考虑CPU频率， 我们可以在一秒钟内解决1000万个未知， 需要O（100）Mbyte的内存。 但是实际上只要多于一万到十万个未知， 你的代码就会崩溃， 需要一个小时左右， 需要不可计数的内存。
问题不在于处理数据， 而在于转移数据。 访问数据有两个指标， 延迟和带宽。 延迟--读写内存的时间太长， CPU内的数据访问需要0.5ns， 而内存需要50ns。 带宽--每秒可以读写的byte数， 对于1GFLOP的CPU需要24Gbyte/秒的带宽。 而现在CPU所拥有的峰值带宽是5Gbyte/s， 而实际应用上远远达不到这个数字。 解决的办法除了可以提高带宽无法提高延迟的Interleaving以外， 就是两方面都可以提高的cache， 通常现在会有两到三层的cache， 但是只有在数据在本地的cache的时候才能用。 载入cache的原则有根据时间的还有根据空间的， 每次载入的是一个cache line， 通常是32-128byte， 现在最长时1kbyte。 这里需要了解的是替换的策略， 关联， cache line的大小， hit ratio和miss ratio的概念。 这是一些现在电脑的cache数据。
•
IBM Power 3: • L1 = 64 KB, 128-way set associative
L2 = 4 MB, direct mapped, line size = 128, write back
•
Compaq EV6 (Alpha 21264):
• L1 = 64 KB, 2-way associative, line size= 32
• L2 = 4 MB (or larger), direct mapped, line size = 64
•
HP PA: no L2
• PA8500, PA8600: L1 = 1.5 MB, PA8700: L1 = 2.25 MB
•
AMD Athlon: L1 = 64 KB, L2 = 256 KB
•
Intel Pentium 4: L1 = 8 KB, L2 = 256 KB
•
Intel Itanium:
• L1 = 16 KB, 4-way associative
• L2 = 96 KB
6-way associative, L3 = off chip, size varies
让程序更快的原则
- 选择最好的算法。这是废话， 不过意思是告诉你， 破算法没有什么优化的价值。要优化就要优化MG这样优秀的算法。
-on
-fast
-unroll
-arch
阅读编译器的man，w ww.specbench.org特定平台的选项， 试验和比较
- 使用合适的数据布局。 按顺序访问数据。例子：访问三个向量 double a[n],b[n],c[n]。 这样的效率通常更高 doubel abc[n][3]
评估- 子程序水平的评估： 编译器在开始和结束的时候插入timing calls， 粗略的评估， 评估的代价高。
- 时钟水平的评估： 定时的执行OS的中断指令， 监视代码的存放位置， 更详细的评估， 代价仍然很高。
- Data cache misses (for different levels)
- Instruction cache misses
- TLB misses
- Branch mispredictions
- Floating-point and/or integer operations
- Load/store instructions
- Etc.
工具
PCL =Performance Counter Library, 跨平台，在代码内或代码外使用。
PAPI =Peformance API，跨平台， 高级和低级的接口。
DCPI ：基于Alpha的Compaq Tru64 Unix
例子： 代码外
PCL Hardware performance monitor: hpm
% hpm –-events PCL_CYCLES, PCL_MFLOPS ./mg
hpm: elapsed time: 5.172 s
hpm: counter 0 : 2564941490 PCL_CYCLES
hpm: counter 1 : 19.635955 PCL_MFLOPS
或者是代码内 ……
可见要查看cache的miss ratio， 就必须要用第三种方法了。
浮点数计算的优化- Loop unrolling
- Fused Multiply-Add (FMA) instructions
- Exposing instruction-level parallelism (ILP)
- Software pipelining (again: exploit ILP)
- Aliasing
- Special functions
- Eliminating overheads
• if statements
• Loop overhead
• Subroutine calling overhead
Loop unrolling
更少的判断/跳转指令， （loop更大， 代价就越小）
更少的load指令
--让loop更大，比如i：n步进为1的循环，改成i：n步进为4的循环。
--flop-to-load-ration
--重命名， Fortran不允许， c可以。 F更快。c可以使用编译选项 -noalias或c的关键字restrict
--特别的函数，比如除， 开根号， 指数，对数， 三角运算代价更高， 会用掉几个周期。
--避免在最内层的循环使用if， 但是不存在这样的技术。 不要在判断时用到上面说的那些代价高的函数。
--循环的代价：cpu需要释放一定的寄存器， 比如说循环寄存器，地址寄存器来开始循环， 所以短的循环相对代价更高。 如果n>m, 这样的循环
do i=1，n
do j=1，m
就没有下面这个循环有效率
do i=1，m
do j=1，n
--
call子函数的代价。子函数对程序的良好结构很重要， 但是呼叫的代价很高， 大概是100个cycles。 传递值参数需要复制数据，
如果使用不当代价会非常高。 但是使用指针可能会不安全， 可以在使用时用const的声明。 通常， 在很小的循环中， 不应该呼叫子函数。
在c中使用incline声明， 使用宏。
针对内存的数据布局优化
stride-1内存访问。
用合并数组的方式创造基于cache优化的数据结构。
数组填充（padding）
4*3数组A[i][j]
内存地址
8 A[][] A[][] A[][] A[][]
4 A[][] A[][] A[][] A[][]
0 A[][] A[][] A[][] A[][]
stride-1（一步数组访问）：数组变换， 和上面的循环优化一样， 一般由编译器完成。 重用cache line的内容。
cache-aware的数据结构：
将需要的数据放在一起， 提高数据在内存空间的集中性， 将其写到一个struct中
数组填充：
建立比需要更大的数组， 改变相对的内存距离， 避免严重的cache匹配失败问题。 例如将(FORTRAN中以列为顺序)double precision u(1024, 1024)改为double precision u(1024+pad, 1024)
如何padding
C-W Tseng（UMD） 数组内的intra-variable pad和数组间的inter-variable pad.
cache
的映射规则有直接映射， 部分关联映射和全关联映射， 对于前两者， 两块不同的内存位置是可能会竞争同一块cache行的。
从而会引起匹配cache失败。 通过在数组之间或者数组里面增加pad， 可以让竞争的内存分配到cache中不同的地方。 但是，
这个padding的计算是非常困难的。
数据访问的优化。-loop unrolling
-loop interchange
-loop fusion
-loop split
-loop skewing
-loop blocking
循环合并：将几个循环合并为一个，减少cache miss， 增加了时间上的集中性。 当数据集合被重复使用的时候可以使用， 比如说叠代运算。
do i= 1,N
a(i)= a(i)+b(i)
enddo
do i= 1,N
a(i)= a(i)*c(i)
enddo
do i= 1,N
a(i)= (a(i)+b(i))*c(i)
enddo
改为
a(i)= a(i)+b(i)
do i= 1,N
a(i)= (a(i)+b(i))*c(i)
a就从load到cache两次变为load到cache一次。
例子：红黑GS叠代
传统的方法--先松弛全部的红点， 再松弛全部的黑点。
但
是实际上只要黑点上面一行的红点更新后， 更新这个黑点所需要的四个红点就都已经都更新了， 所以可以在更新上密集按那一行红点的同时 ， 更新黑点。
这样比传统的方法将矩阵两次载入cache变为之需要一次！ 当然， 第一行的红点和倒数第二行的黑点需要特殊的处理。
可以看到这个gif图片上的动态显示， 点外围的圈数代表松弛的次数。
循环分割： 循环合并的逆操作。 提升编译器的优化， 提升指令cache的利用率。 （尽量的利用一次移入cache的数据集合）
分块循环（loop blocking, loop tiling）： 将数据集合分割成适合cache尺寸大小的小块子集。 在转移到下一块数据子集的时候尽可能的进行更多的运算。 因为数据的相互依赖关系的存在， 这有时并不是那么容易实现。
应用一维分块技术后的代码
B 是GS叠代分的块数
for it= 1 to numIter/B do
// Special handling: rows 1, …, 2B-1
// Not shown here …
// Inner part of the 2D grid
for k= 2*B to n- 1 do
for i= k to k- 2*B+1 by ?2 do
for j= 1+(k+1)%2 to n- 1 by 2 do
relax(u(j,i))
relax(u(j,i- 1))
end for
end for
end for
// Special handling: rows n-2B+1, …, n-1
// Not shown here …
end for
数据被B次载入cache， 如果2*B+2正好合适cache的大小。
数据如果过大， 可以采用2D的分块。 或者更复杂的分块。
介绍了根据cache优化的Multigrid方法的一个包， DiMEPACK。 用的是C++的接口， F77的子程序， 接下来看看这个包。 不过里面子程序用FORTRAN写的， 用M4宏处理器， 有大量的宏， 跟以前看的程序还不一样。
有些事情其实本应该是编译器去做的， 但是编译器做得还没有那么好， 自己写的时候注意一下， 考虑到内部的结构， 就可以很方便的提高速度。