矩阵转置向量优化

问题描述

使用simd向量优化转置算子

如果我们要对$n\ast n$的矩阵$m$实现转置，可以用以下代码实现

for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) {
        m_out[i][j] = m[j][i];
    }
}

不考虑cache的局部性原理的话，需要$n^2$次访存去读数据以及$n^2$次访存去写数据，而在使用simd向量对计算进行优化时，每次load的数据块是连续的，无法将不连续的数据块一次访存load为一个向量，因此我们可以使用向量的平移和翻转操作间接的实现矩阵的转置操作。

向量平移

simd_concat

也可以视为是向量的拼接操作，比如1 2 3 4右移一位得到2 3 4 1

向量翻转

simd_seleq

也就是向量的选择操作，可以实现对两个向量指定位置的翻转，比如a={1,2,3,4}，b={5,6,7,8}，翻转第一个和第三个位置可以得到a={5,2,7,4}，b={1,6,3,8}。

注：不同编译器使用的平移和翻转指令名称可能不同，往往需要两条simd_seleq指令才能实现两个向量的交换操作，下文为了说明方便统一将两次simd_seleq指令记为一次翻转操作。

以4*4的矩阵转置操作为例

1 5 9  13          1  2  3  4
2 6 10 14    -->   5  6  7  8
3 7 11 15          9  10 11 12
4 8 12 16          13 14 15 16

平移操作可以使元素在行内位置发生改变，翻转操可以试元素所在行号发生改变。

如果我们将第2、3、4行分别右移1、2、3个位置，可以得到

1  5  9  13
14 2  6  10
11 15 3  7
8  12 16 4

然后将1，2行第二个元素，1，3行第三个元素，1，4行第四个元素进行翻转，第一行元素就变成了想要转置的最终元素。

原理也很简单，原矩阵元素2在第二行第一列，要想到达最终位置第一行第二列，只能选择右移一位然后向上翻转。如果先翻转再右移会将1翻转下来，所以我们选择将2右移后再翻转，其它元素也是一样。

基于上面平移后的矩阵，我们想要变换为最终的转置矩阵有两种思路：减而治之和分而治之

减而治之

减而治之是我们最容易想到的思路，也就是先让第一行元素到达最终状态，然后是第2，3，4行。

翻转的原则是保留应该在本行的元素，换出去不属于本行的元素。

// 平移后的矩阵
1  5  9  13
14 2  6  10
11 15 3  7
8  12 16 4
// 第1行元素与2，3，4行元素进行翻转
1  2  3  4
14 5  6  10
11 15 9  7
8  12 16 13
// 第2行元素与3，4行元素进行翻转
1  2  3  4
8  5  6  7
11 15 9  10
14 12 16 13
// 第3行元素与第4行元素进行翻转
1  2  3  4
8  5  6  7
11 12 9  10
14 15 16 13
// 第2，3，4行元素平移
1  2  3  4
5  6  7  8
9  10 11 12
13 14 15 16

开头和末尾都进行了3次平移，一共6次平移操作；第1行与下面3行进行了翻转，第2行与下面2行进行了翻转，第3行与下面1行进行了翻转，翻转操作一共：3 + 2 + 1 = 6次，可以看出采用减而治之的方式对$n\ast n$的矩阵进行转置需要2(n-1)次平移和(n-1) + (n-2) +…+1 = n(n-1)/2次翻转。

分而治之

减而治之的方案在$4\ast 4$矩阵的转置上似乎效率还可以，但是矩阵规模一旦增加，平方级别的复杂度会让指令条数飞速增长，比如$16\ast 16$的矩阵，就需要30次平移和120次翻转操作才能够实现转置。

分而治之和减而治之的思想都是通过翻转操作，让不属于本行的元素翻转出去，也就是初始平移后，只考虑将元素翻转到他最终的行里，最后通过平移操作一起归位。分而治之的思想是让矩阵前一半的元素和后一半的元素分开，这样问题的规模就会减半，以$2\ast 2$的矩阵为例：

// 平移后的矩阵
1  5  9  13
14 2  6  10
11 15 3  7
8  12 16 4
// 行跨步为2进行翻转，对1，3和2，4行元素进行翻转
1  5  3  7
8  2  6  4
11 15 9  13
14 12 16 10
// 可以发现前两行的元素都属于前两行，后两行的元素都属于后两行
//对上下两个子问题分别进行行跨步为1的翻转，即翻转1，2和3，4
1  2  3  4
8  5  6  7
11 12 9  10
14 15 16 13
// 第2，3，4行元素平移
1  2  3  4
5  6  7  8
9  10 11 12
13 14 15 16

分而治之的思路平移操作也是6次，翻转操作行跨步为2和行跨步为1的翻转各2次，一共4次，在$4\ast 4$的矩阵转置中翻转虽然只少了两次，但是随着矩阵规模的增大，操作次数就不是一个量级的了。

对于$n\ast n$的矩阵，前$n/2$行元素每行恰好有一半元素属于后$n/2$行，所以我们可以通过行跨步为$n/2$的翻转操作将不输于上面矩阵的元素翻转下去，一共做$n/2$次翻转可以将规模为$n$的问题划分为两个规模为$n/2$的子问题；对于$n/2$的子问题，也需要做$n/4$次翻转进一步划分子问题，直到问题规模为1。翻转操作次数为：

$$n/2+2\ast (n/4)+4\ast (n/8)+...+n/2\ast 1=n/2\ast logn$$

总结

总结下两种算法对$n\ast n$矩阵转置的操作次数：

算法	平移次数	翻转次数
减而治之	2(n-1)	n(n-1)/2
分而治之	2(n-1)	n/2*logn

以$16\ast 16$矩阵为例，减而治之的思路需要120次翻转操作，分而治之的思路则只需要32次翻转操作。

下面代码模拟了$16\ast 16$矩阵转置分而治之的翻转过程：

#include <iostream>
#include <algorithm>
using namespace std;
const int N = 16;
int m[N][N];
int tmp[N];
int w[N] = {1, 2, 4, 8, 16, 32, 64, 128, 256, 512,
			1024, 2048, 4096, 8192, 16384, 32768};
int pos1[15][N] = {
	{0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1},
	{1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1},
	{1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1},
	{1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1},
	{1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1},
	{1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1},
	{1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1},
	{1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1},
	{0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1},
	{1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1},
	{1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1},
	{1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1},
	{0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1},
	{1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1},
	{0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1}
};
int pos2[N];
void print() {
	for (int i = 0; i < N; i++) {
		cout<<i<<":";
		if (i < 10) cout<<"  ";
		else cout<<" ";
		for(int j = 0; j < N; j++) {
			cout<<m[i][j];
			if (m[i][j] < 10) cout<<"    ";
			else if(m[i][j] < 100) cout<<"   ";
			else cout<<"  ";
		}
		cout<<endl;
	}
	cout<<"-----------------------------------------";
	cout<<"-----------------------------------------"<<endl;
}
void move(int row, int step) {
	if (!step) return;
	for (int i = 0; i < N;i++) {
		int idx = (i + step) % N;
		tmp[idx] = m[row][i];
	}
	memcpy(m[row], tmp, sizeof(tmp));
}
void swap_row(int a, int b, int st) {
	for (int i = 0; i < N;i++) {
		if (st >> i & 1) {
			swap(m[a][i], m[b][i]);
		}
	}
	print();
}
int main() {
	for (int i = 0; i < N; i++) {
		for (int j = 0; j < N; j++) {
			int t = i * N + j;
			m[j][i] = t;
		}
	}
	for (int i = 0; i < N; i++) move(i, i); // 先逐行平移
	for (int i = 0; i < 15;i++) {
		for (int j = 0; j < N;j++) {
			pos2[i] += w[j] * pos1[i][j];
		}
	}
	print();
	cout<<"进行行跨步为8的翻转操作"<<endl;
	swap_row(0, 8, pos2[0]);
	swap_row(1, 9, pos2[1]);
	swap_row(2, 10, pos2[2]);
	swap_row(3, 11, pos2[3]);
	swap_row(4, 12, pos2[4]);
	swap_row(5, 13, pos2[5]);
	swap_row(6, 14, pos2[6]);
	swap_row(7, 15, pos2[7]);
	cout<<"进行行跨步为4的翻转操作"<<endl;
	swap_row(0, 4, pos2[8]);
	swap_row(8, 12, pos2[8]);
	swap_row(1, 5, pos2[9]);
	swap_row(9, 13, pos2[9]);
	swap_row(2, 6, pos2[10]);
	swap_row(10, 14, pos2[10]);
	swap_row(3, 7, pos2[11]);
	swap_row(11, 15, pos2[11]);
	cout<<"进行行跨步为2的翻转操作"<<endl;
	swap_row(0, 2, pos2[12]);
	swap_row(4, 6, pos2[12]);
	swap_row(8, 10, pos2[12]);
	swap_row(12, 14, pos2[12]);
	swap_row(1, 3, pos2[13]);
	swap_row(5, 7, pos2[13]);
	swap_row(9, 11, pos2[13]);
	swap_row(13, 15, pos2[13]);
	cout<<"进行行跨步为1的翻转操作"<<endl;
	swap_row(0, 1, pos2[14]);
	swap_row(2, 3, pos2[14]);
	swap_row(4, 5, pos2[14]);
	swap_row(6, 7, pos2[14]);
	swap_row(8, 9, pos2[14]);
	swap_row(10, 11, pos2[14]);
	swap_row(12, 13, pos2[14]);
	swap_row(14, 15, pos2[14]);
	return 0;
}