矩阵并行计算

一、矩阵卷帘存储

串行计算中无需考虑数据在处理中怎样存储，但是在并行计算中需要进行关注，以便进行负载均衡等操作。

卷帘存储：矩阵A，矩阵中的元素为$A_{ij}$，处理机阵列$p\times q$，小块$A_{ij}$存放在处理器$P_{kl}$处（k = i mod p, l = j mod q)。
如下图所示：

二、矩阵向量并行乘法

$y=Ax,A\in {\mathbb{R}}^{mn},x\in {\mathbb{R}}^n$

2.1 串行算法

2.1.1 i-j实现

矩阵A的一行乘以向量x，外层循环i循环一次确定结果y中的一个值，内层循环j负责控制遍历矩阵A的一行，x的每一行（即为每一个元素）。

void matvec_i_j(float *y, float **A, float *x, int m, int n)
{
	float *arr = (float*)calloc(m, sizeof(float));
	memset(arr, 0, m);
	for (int i = 0;i < m;i++)
		for (int j = 0;j < n;j++)
			arr[i] += A[i][j] * x[j];

	printf("Multiply the matrix by the vector: \n");
	for (int i = 0;i < m;i++)
		printf("%f  \n", arr[i]);
}

2.1.2 j-i实现

矩阵A的一列乘向量x的一个值，即为外层循环j循环一次计算出结果y中每个值的第j个分量，再对每次算出的值相加即可。

void matvec_j_i(float *y, float **A, float *x, int m, int n)
{
    float *arr = (float *)calloc(m, sizeof(float));
    memset(arr, 0, m);
    for (int j = 0; j<n; j++){
        for (int i = 0; i<m; i++){
            arr[i] += A[i][j] * x[j];
        }
    }
    printf("Multiply the matrix by the vector j-i: \n");
    for (int i = 0; i<m; i++)
        printf("%f  \n", arr[i]);
}

二者运行结果相同，如下图所示：

2.2 并行算法

2.2.1 行分块划分矩阵

将矩阵 A 按行划分成如下的行块子矩阵
$A={\left[A_0^T,A_1^T,\dots ,A_{p-1}^T\right]}^T$
则：
$\begin{array}{rl}Ax& ={\left[A_0^T,A_1^T,\dots ,A_{p-1}^T\right]}^{T}x\\ & =\left[{\left(A_{0}x\right)}^T,{\left(A_{1}x\right)}^T,\dots ,{\left(A_{p-1}x\right)}^T\right]\end{array}$

将$A_i$存放在结点$P_i$中，每个结点计算$A_{i}x$，最后调用 MPI_GATHER 或 MPI_GATHERV 即可。

2.2.2 列分块划分矩阵

将矩阵 A 按列划分，并对 x 也做相应的划分
$A=\left[A_0,A_1,\dots ,A_{p-1}\right],x={\left(x_0^T,x_1^T,\dots ,x_{p-1}^T\right)}^T$
则有结果为：
$Ax={\sum }_{i=0}^{p-1}{A}_i{x}_i^T$
将 $A_i$ 和 $x_i$ 存放在结点$P_i$ 中，每个结点计算 $A_i{x}_i^T$，最后调用 MPI_REDUCE 或 MPI_ALLREDUCE 即可。

三、矩阵乘并行计算

本节要考虑的计算问题为：
$$C=A\times B$$

3.1 串行矩阵乘法

3.1.1i-j-k实现

矩阵A的行和矩阵B的列相乘相加即可。

for (int i = 0;i < m;i++){
	for (int j = 0;j < n;j++){
		for (int r=0;r<k;r++)
			arr[i][j] += A[i][r] * B[r][j];
	}
}

3.2 矩阵并行分块乘法

3.2.1 行列分块

矩阵A按行进行划分，矩阵B按列进行划分。
$$A={\left[\begin{array}{llll}{A}_0^T& A_1^T& \cdots & A_{p-1}^T\end{array}\right]}^T,B=\left[\begin{array}{llll}{B}_0& B_1& \cdots & B_{p-1}\end{array}\right]$$
第一轮计算后得到
$$C_{ij}=A_i\times B_j$$
然后第i号处理器将其所有的矩阵B的列块传递至第i-1号处理器中，再次进行计算。B的列块循环传送p-1次后得到C的最终结果。此时编号为i的处理器中存储有子矩阵$C_{i0},C_{i1},...,C_{i(p-1)}$。
循环过程如下图所示，第一次计算时为图(a)，传递后为图(b)。

其实整个计算的过程是从对角线开始斜向上进行计算的，计算过程如下图所示：

for i=0 to p-1
	j=(i+myid) mod p
	Cj=A*B
	src = (myid+1) mod p
	dest = (myid-1+p) mod p
	if (i!=p-1) 
		send(B,dest)
		recv(B,src)
	end if
end for

本算法中，Cj = Cmyid, j ,A = Amyid ,B 在处理器中每次循环向前移动一个处理器，即每次交换一个子矩阵数据块，共交换 p-1 次。

3.2.2 行行分块算法

将矩阵A、B均按行进行划分，如下：
$$A={\left[A_0^T,A_1^T,\dots ,A_{p-1}^T\right]}^T,B=\left[B_0^T,B_1^T,\dots ,B_{p-1}^T\right]$$
再对矩阵A的每一行再次划分，划分为单个元素：
$$A_i=\left[A_{i0},A_{i1},\dots ,A_{i,p-1}\right],i=0,1,\dots ,p-1$$
要计算出矩阵C的第一行，需要用矩阵A中的第一行去乘矩阵B的每一列，也就是
$$C_i=A_i\times B$$
此时真正需要相乘的是矩阵A的第一行的第j个元素和矩阵B的第j行，每次相乘会得到矩阵$C_{ij}$的一部分，j++依次计算即可得到矩阵C的第i行的值，即为 ：
$$C_i=\sum _{j=0}^{p-1}{A}_{ij}\times B_j$$

mp1 ≡ (myid+1) mod p, mm1 ≡ (myid−1 + p) mod p
for i = 0 to p − 1 do
	l ≡ (i+myid) mod p
	C = C + Al × B
	if i ≠ p − 1, send(B, mm1), recv(B, mp1)
end{for}

数据交换方式与3.2.1相同，不同的只是计算C的方式。

3.2.3 列列分块

矩阵A和矩阵B均按列进行分块，C也按列进行划分。
与3.2.2类似此处对B的一列再次进行划分，划分为单个元素与A相乘，计算方法如下：
$$C_j=A\ast B_j=\sum _{i=0}^{p-1}{A}_i{B}_{ij}$$
此时唯一不同在于要进行交换的是A矩阵，不是B矩阵。A矩阵在每次循环后，发送$A_j$至上一个处理器。

mp1 ≡ (myid+1) mod p, mm1 ≡ (myid−1 + p) mod p
for i = 0 to p − 1 do
	l ≡ (i+myid) mod p
	C = C + A × Bl
	if i ≠ p − 1, send(A, mm1), recv(A, mp1)
end{for}

3.2.4 列行分块

矩阵A按列分块，矩阵B按行进行分块.
$$A=\left[\begin{array}{llll}{A}_0& A_1& \cdots & A_{p-1}\end{array}\right],B_i=\left[\begin{array}{llll}{B}_{i0}& B_{i1}& \dots & B_{i,p-1}\end{array}\right]$$
此时计算C时按照下式进行：
$$C=A\ast B=\sum _{i=0}^{p-1}{A}_i{B}_i$$
每次计算都可以算出C的一个子矩阵，一共算出P个与C相同大小的子矩阵，再进行相加即可得到C。
假设C也是按列存放在处理器中的则有将B的一行再次进行了划分（按行存放与3.2.2类似）：
$$C_j=\sum _{i=0}^{p-1}{A}_i{B}_{ij}$$

C = A × Bmyid
for i = 1 to p − 1 do
	l ≡ (i+myid) mod p, k ≡ (p − i+myid) mod p
	T = A × Bl
	send(T, l), recv(T, k)
	C = C + T
end{for}

3.3 Cannon乘法

Cannon乘法和并行分块乘法不同，不是仅仅让 B 矩阵的各列块循环移动，而是有目的地让 A 的各行块以及 B 的各列块皆施行循环移位，从而实现对 C 的子块的计算。
将矩阵 A 和 B 分成 p 个方块$A_{ij}$和$B_{ij}$$(0\le i,j\le \sqrt{p}-1)$，每块大小为$\lceil n/\sqrt{p}\rceil \times \lceil n/\sqrt{p}\rceil$，并将它们分配给$\sqrt{p}\times \sqrt{p}$ 个处理器$\left(P_{00},P_{01},\dots ,P\sqrt{p}-1\sqrt{p}-1\right)$。

算法执行分三步：

整个计算过程如下：