4k aliasing對分塊算法的影響和lapacke中行列主序的問題
由於程序運行時間測量的不準確性，雖然測量時間時已經採取了運行多次取中位值的方法，我不保證能夠重現結果。
下圖是使用dgemm()矩陣乘法用時計算所得CPU浮點計算能力（Gflops，y軸，越大越好）隨運行時間（x軸）變化的曲線。可見1950x(3.4G,3.75G最上面紅藍黑三條)的波動是非常劇烈的，偶爾還會來個尖峯或者波谷，最上面的紅線要穩定些應該是因爲排除了4K alisaing的影響。倒數第二的紅線是筆記本低壓U i7-8550U（外接電源），短暫睿頻結束後很快就穩定到110Gflops的位置。
下圖是AMD Thread Ripper 1950x（後面簡寫爲TR 1950x）在不同配置下的Gflops結果，可見正確的配置是很重要的。後面的測試中由於失誤OpenBLAS在TR 1950x下運行的是HASWELL代碼，這可能會導致運行時間變長。
所有文章都告訴我們要對齊，struct要對齊 padding，棧要對齊，數據要對齊。但沒人告訴你過度的對齊會帶來性能損失，嚴格的說不是對齊，而是2整數次冪強迫症。爲充分利用cache大量數值算法採用的分塊的方法，分塊內部的kernel經常會同時對相差若干lda的數據進行操作（SSE,AVX指令），如果lda恰好爲某個2的整數次冪左右（這與cache有關，一般爲4K字節），那麼會映射到同一個cache地址，當分塊內映射到同一cache地址的數據超過了ways大小就會起衝突。不過一般來說由於現代x86 cache的ways還是比較多的，分塊kernel也不容易超過ways大小。嚴格的驗證需要使用vtune,perf等工具檢測cache miss。
對了，如果不知道lda，可以去了解下blas，lapack的函數接口。lda=leading dimension of A
下圖是若干矩陣乘法函數的性能比較。primitive那個是我手寫的，按A*B爲A的列向量組合的思路+FMA指令，並不真的是單純for循環的那種primitive。在i5-7500上除了primitive不受lda，ldc設置影響外其他函數包括openblas和mkl都受ldc影響，nopacking（做了分塊但沒做packing）受ldc和lda影響特別大，或者說，packing後內層for循環的運算都在固定的一個小工作內存區域內，減弱了lda的影響，但ldc依然對dgemm的性能，哪怕是mkl版本有可觀的影響。
很多人簡單創建兩個1024*1024/2048*2048/4096*4096的矩陣然後計算矩陣乘法估算gflops，卻沒想過這樣得到的結果是不準確的。
intel的文檔上也提到了這個問題，不過重點說的是地址相差4K整數倍下的寫後讀問題。
Intel 64 and IA-32 Architectures Optimization Reference Manual, June 2016.11.8 4K ALIASING
4-KByte memory aliasing occurs when the code stores to one memory location and shortly after that it
loads from a different memory location with a 4-KByte offset between them. For example, a load to linear
address 0x400020 follows a store to linear address 0x401020.
The load and store have the same value for bits 5 - 11 of their addresses and the accessed byte offsets
should have partial or complete overlap.
4K aliasing may have a five-cycle penalty on the load latency. This penalty may be significant when 4K
aliasing happens repeatedly and the loads are on the critical path. If the load spans two cache lines it
might be delayed until the conflicting store is committed to the cache. Therefore 4K aliasing that happens
on repeated unaligned Intel AVX loads incurs a higher performance penalty.
To detect 4K aliasing, use the LD_BLOCKS_PARTIAL.ADDRESS_ALIAS event that counts the number of
times Intel AVX loads were blocked due to 4K aliasing.
To resolve 4K aliasing, try the following methods in the following order:
• Align data to 32 Bytes.
• Change offsets between input and output buffers if possible.
• Use 16-Byte memory accesses on memory which is not 32-Byte aligned.
另外文檔上還強調如果load,store的對齊不可兼得則優先考慮store操作的對齊。在故意設置lda,ldb,ldc使之不對齊後我發現真正對dgemm速度（i5-7500, n=2048,4096）有顯著影響的只有ldc，這也與文檔說法相符，但不確定是否有packing操作的影響。
ZEN等AMD架構上的4K間隔寫後讀也有同樣的問題。
談到amd不得不說intel家產品對其不友好的態度，雖然mkl文檔裏到處是提示。但是在TR1950x上連mkl_cbwr_set都不讓用還是挺過分的。嗯，我並沒有安裝mkl，是用matlab自帶的那個mkl.dll，手動導出函數列表發現了mkl_cbwr_set這個玩意，使用時也被坑過，看的c版本說好是int mkl_cbwr_set(int)，結果程序崩了，容我三思後意識到這玩意是Fortran的，正確姿勢是int mkl_cbwr(int*). 還有個mkl_enable_instructions ，但沒能在mkl.dll的導出函數表裏發現它。mkl在TR 1950x上的表現也沒讓我失望，一如既往的爛，後面有幾個測試爲節省時間就直接跳過mkl了。所以，amd上就忘了mkl吧。
下面這個是dposv的測試結果。我知道你們不會看圖的，只說結果。自己手寫的my code針對column/row major,upper/lower使用不同的處理方式，所以沒有明顯規律。而openblas和mkl對column major/row major非常敏感，並且受lda影響（因爲內部也是分塊的），這點並沒有寫到手冊上！雖然這是實現相關的內容，但庫開發者應提醒用戶lapacke函數的性能與輸入數據的layout有關。openblas的lapacke直接用的netlib實現，看過的就知道lapacke只不過是給lapack套了個C語言的殼子，內部實現完全依賴lapack，就是函數輸入參數加了個major，當其爲row major時我看過的那些函數都是將輸入矩陣做個轉置，調用lapack函數，算完後再轉置，釋放內存。因爲大稠密矩陣轉置必然帶來不連續的內存訪問，實際效率遠低於普通的memcpy，而且這樣內存使用量直接翻倍。