Optimizing Parallel Multiplication Operation for Rectangular and Transposed Matrices
===
**[論文連結](http://hpc.pnl.gov/srumma/srumma-icpads.pdf)**
[youtube連結](https://www.youtube.com/watch?v=HGWQbTPBA54)
## Index
* paper study
* OpenMPI
## paper study
### 平行矩陣乘法
- **1D-systolic, 2D-systolic**
- **Cannon's algorithm**
- **Fox's algorithm**
- **Berntsen’s algorithm**
- **DNS algorithm**
- **SUMMA**: Scalable Universal Matrix Multiplication Algorithm
- **SRUMMA**: Shared and Remote-memory based Universal Matrix Multiplication Algorithm
- ...還有很多延伸的演算法，適合在伺服器或分散式架構上運算。
:::success
SRUMMA就是一開始<`0140454`>, <`linachiu`>寫的[submatrix_multiply](https://github.com/0140454/matrix-multiplication/blob/master/impl.c#L25)的分散式版本加上`communication`，`communication`可以想成是在計算的同時進行下一次資料的傳輸，與`prefetch`想法很像只是記憶體層級不一樣，且`communication`通常是透過 MPI 介面的點對點傳輸。
:::
### shared-memory and distributed-memory
最常見的 shared-memory system 就是`OpenMP`了，只要幾行指令 compiler 就會幫我們將程式平行化，但這篇論文想探討的是 remote memory access (`RMA`)，它基本上就是`DMA`的翻版，它會有額外的硬體去處理 cluster 對 cluster 記憶體之間資料傳送。
而 distributed-memory 就是透過`message passing interface(MPI)`傳遞資料，雖然必須透過網路傳遞資料，但其擴充性較高。
### SRUMMA Baseline Efficiency
考慮一個 $C=AB$ 的矩陣乘法，其中A, B, C分別是`mxk, kxn, mxn`的矩陣。
$t_w =每個\ element\ 傳輸時間$
$t_s =latency$
$p*q,\ process\ grid\ in\ two\ dimension$
$P,\ processors\ number$
上圖是一個$p=q=3$的例子，經過分塊處理後A, B, C的區塊大小分別為$\left(\frac{m}{p}\right)*\left(\frac{k}{q}\right),\left(\frac{k}{p}\right)*\left(\frac{n}{q}\right),\left(\frac{m}{p}\right)*\left(\frac{n}{q}\right)$
假設 m=n=k=N，P=p^2^=q^2^，而每一次運算，`RMA`會取得`q`個矩陣A的 blocks，矩陣B要取`p`個。
$T_{par\_rma} = \text{Compute Time}(T_{comp})+ \text{Communication Time}(T_{comm})$
$T_{comm} = T_{A\_row\_comm} + T_{B\_col\_comm}$
$\begin{align}T_{A\_row\_comm} & = \{\text{data transfer time of size} \frac{mk}{pq}\}+\{\text{latancy}\} \\
& = \left(\left(\frac{mk}{pq}\right)t_w+t_s\right)q
\end{align}$
$T_{B\_col\_comm} = \left(\left(\frac{nk}{pq}\right)t_w+t_s\right)p$
總後可以得出結論，平行矩陣乘法的複雜度為。
$\begin{align}T_{par\_rma} &= \frac{mnk}{P} + \left(\left(\frac{mk}{pq}\right)t_w+t_s\right)q + \left(\left(\frac{nk}{pq}\right)t_w+t_s\right)p \\
&= O\left(\frac{N^3}{P}\right)+O\left(\frac{N^2}{\sqrt{P}}\right)+O\left(\sqrt{P}\right)\end{align}$
若在執行運算的同時進行傳輸，則第二項的傳輸時間$t_w$可以視為0，因此新的複雜度為。
$T_{par\_rma} =O\left(\frac{N^3}{P}\right)+O\left(\sqrt{P}\right)$
### diagonal shift
論文整理了幾點平行矩陣乘法要注意的細節，其中一點是要注意`computation`及`communication`重疊的問題。
上圖(a)是一個矩陣被分配給4個 cluster，每個 cluster 有4個CPU，因此矩陣被分成16份，每份由一個CPU負責。
如果沒有 diagonal shift 如圖(b)，node 1有`P~00~、P~10~、P~20~、P~30~`，node 2有`P~01~、P~11~、P~21~、P~31~`，我們知道P~00~要取得`A 的P~00~、P~10~、P~20~、P~30~`以及`B 的P~00~、P~01~、P~02~、P~03~`的資料，但如果程式一開始，node 1的CPU 會同時存取node 2的CPU，會吃掉許多頻寬。
diagonal shift可以減輕這個問題，如圖(3)所示，如果我把一開始的位置分均分配在4個nodes上，P~00~一開始會分別取得node 2, node 3, node 4的資料，而P~11~會去取node 1, node 3, node 4的資料，如此可以錯開存取的順序。
### Rectangular Matrices
SRUMMA有3種分割矩陣的方式，依照 row-base, column-base 或 square-base的不同可分成以下的形式。若 matrix A 與 matrix B 分別分割成 s 個區塊，first-order 複雜度是 O(s)，second-order是 O(s^2^)，third-order 是 O(s^3^)。
> 但這裡我覺得很奇怪，first-order及second-order的複雜度應該調過來才對
### 矩形矩陣乘法
### Impact of communication protocols
論文最後做三個實驗，前面兩個是比較不同形狀的矩陣在不同 order 分割下的效能，以及轉置矩陣的效能。那這邊直接節錄第三個實驗，想要知道 MPI 與 RMA 對效能有何影響。
## OpenMPI
看paper的時候看到一些關鍵字`numa` `mpi`，在自己的筆電上下`lscpu`也有出現numa nodes，於是開始找`numa`相關資料。
```
Architecture: x86_64
CPU 作業模式： 32-bit, 64-bit
Byte Order: Little Endian
CPU(s): 4
On-line CPU(s) list: 0-3
Socket(s): 1
NUMA 節點： 1
Model name: Intel(R) Core(TM) i5-5200U CPU @ 2.20GHz
CPU MHz： 1436.789
CPU max MHz: 2700.0000
CPU min MHz: 500.0000
BogoMIPS: 4389.54
虛擬： VT-x
L1d 快取： 32K
L1i 快取： 32K
L2 快取： 256K
L3 快取： 3072K
NUMA node0 CPU(s): 0-3
```
### numactl, numastat
`numactl`可以讓 processes 跑在 `NUMA`架構上或記憶體位置 (使用`--membind=nodes`, `--pyhscpubind=cpus`)，可以用`numactl --show`查看目前機器上的 numa nodes 及 CPU數，像我的筆電就有4顆，但還沒有連接其它機器所以是 node 0。
而`numastat`可以查看 numa nodes
```
kai@kai-TravelMate:~$ numactl --show
policy: default
preferred node: current
physcpubind: 0 1 2 3
cpubind: 0
nodebind: 0
membind: 0
```
`numactl --physcpubind=+0-4,8-12 ./myapplic`是一個簡單的例子，讓特定的CPU執行程式。
`numactl --cpubind=0 --membind=0,1`在node 0上面執行程式，而記憶體綁定在node 0及node 1。
還有很多例子可上man page查看。
### [hwloc](https://www.open-mpi.org/projects/hwloc/)
[hwloc install](https://www.open-mpi.org/projects/hwloc/tutorials/20120702-POA-hwloc-tutorial.html)
後來才發現hwloc只是幫使用者取得整個分散式架構的資料，不過有幾個有趣的指令如`lstopo`，另外，可以跨平台地讀取機器的資訊，如CPU數量。
### OpenMPI
參考 [官網](https://www.open-mpi.org/) 及 [MPI Hello World](http://mpitutorial.com/tutorials/mpi-hello-world/)。
另外，google搜尋`mpi matrix multiplication`可以找到很多範例，只是目前測試的結果還是很怪，需要再做修改。
#### [MPI Hello World](http://mpitutorial.com/tutorials/mpi-hello-world/)
用這個簡單的例子來說明`mpi`
```clike=
#include <mpi.h>
#include <stdio.h>
int main(int argc, char** argv) {
// Initialize the MPI environment
MPI_Init(&argc, &argv);
// Get the number of processes
int world_size;
MPI_Comm_size(MPI_COMM_WORLD, &world_size);
// Get the rank of the process
int world_rank;
MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
// Get the name of the processor
char processor_name[MPI_MAX_PROCESSOR_NAME];
int name_len;
MPI_Get_processor_name(processor_name, &name_len);
if (world_rank == 0) {
printf("This is master task\n");
}
if (world_rand > 0) {
// Print off a hello world message
printf("Hello world from processor %s, rank %d out of %d processors\n",
processor_name, world_rank, world_size);
}
// Finalize the MPI environment.
MPI_Finalize();
}
```
MPI_Init(&argc, &argv)
: 初始化`MPI`環境，可以傳入參數NULL
MPI_Comm_size(MPI_COMM_WORLD, &world_size)
: 取得所有的CPU數目，一般情況`MPI_COMM_WORLD`不用改
MPI_Comm_rank(MPI_COMM_WORLD, &world_rank)
: 取得目前CPU的順序，不重覆的，假如總共有8個CPU，它會回傳 0~7 其中一個數。
MPI_Finalize()
: 它會結束`MPI`環境
`MPI`一開啟，所有CPU都會跑這段程式，就要用`MPI_Comm_rank()`來區別task 0及其它的task，來達到平行化的環境。
另外還有
`MPI_Send (buf, count, datatype, dest, tag, communicator, ierr)`
`MPI_Recv (buf, count, datatype, source, tag, communicator, status, ierr)`
用來傳送及接收資料，以上只是`OpenMPI`簡單的用法，我在寫矩陣乘法時也只用到上面的api，只是我的程式還有問題，還要再研究
###### tags: `paper_study`
Sign in with Wallet
Connect another wallet