循环体并行优化（二）——多维循环迭代空间的仿射变换及循环上下界不等式的矩阵表示法

本文链接：https://blog.csdn.net/u014038143/article/details/110622193

（重新整理于2020年冬。）

上回书说到（忘了我不是说书的了，习惯性口语，见谅！），我们可以通过一个简单的仿射变换将一维的“不连续”循环下标空间变换到一个“连续”的下标空间中。这次，我们则继续来看看如何将一个高维的“不连续”循环下标变换到“连续”的循环下标空间中。如果非要为这种变换加上一个理由的话，我认为那就是不要在我们的整数向量空间中留下太多的“洞”或者“空泡”。因为对于一个增量不为1的循环变量来说（包括绝对值大于1的负增量也就是递减量），跳过的坐标值形成的那些点就不能在这个迭代空间中，因而形成了这个迭代空间中的“洞”，这为我们后续的不等式分析带来了极大的不便，因为我们除了要写出循环下标的上下界之外，还要写一堆不等式来说明这个循环下标不能等于那些被跳过的值，对于分析来说这是极其不方便的。

大约在1800年前的丢番图他老人家就研究了整数方程，当然我想他也不喜欢这些“洞”。而我们后面其实大多数的内容都是源自于他老人家在1800年前发明的方法，现在被称之为“丢番图方程”，当然那时还没有程序员这种职业。介绍这个只是简单的让大家了解下我们这系列文章究竟是在说什么。

基于此，在进一步分析迭代空间之前，我们有必要将n层嵌套循环中的每一层循环的下标做仿射变换，使其增量均为1。要做到这一点，就需要我们将前一篇文章中的方法自然推广到n重循环上。

首先，我们来看一个2重循环的例子：

for( i = 3; i<10; i+=2)
	for( j = i; j<15; j+=2)
		a[i,j] = ...;

这里我们先对i做变换：

令 $\times ii + 3$ ; 变换外层循环为 for(ii=0; ii<3; ii++)

接着令 $\times jj + i$ ; 注意这里有外层循环的旧下标 $i$ ，所以继续带入 $i$ 的替换得到 $2\times jj + 2 \times ii + 3$ ; 此时要求 $\times jj + 2 \times ii + 3 < 15$ ; 得到 $j j + i i < 6$ ; 于是内层循环变成 for(jj=0;jj+ii<6;jj++)；最终这个双层循环就变为：

for( ii=0; ii<3; ii++ )
	for( jj=0; jj+ii<6; jj++ )
		a[2*ii+3,2*jj+2*ii+3] = ...;//看上去有点点复杂了

变换前后循环下标变化如下表所示：

i	ii	*2ii+3**	j（i=3）	j（i=5）	jj	*2jj+2ii+3（ii=0）*	*2jj+2ii+3（ii=1）*
3	0	3	3	5	0	3	5
5	1	5	5	7	1	5	7
7	2	7	7	9	2	7	9
9	3	9	9	11	3	9	11
			11	13	4	11	13
			13	15	5	13	15
			15		6	15

表中只列出了 $i = 3 、 5$ 时的的 $j$ 值，这个很容易从原循环下标中验证得到，同时也只列了 $i i = 0 、 1$ 时的 $\times ii + 3$ 及 $\times jj + 2 \times ii + 3$ 的对应值，同样这也可以很容易的从变换后的循环中推出。

(为了文章结构性果断的割一下)

到这里我们完整的介绍完了1、2维“非连续”循环迭代空间到“连续”循环迭代空间的变换方法，以此类推更多层循环的变换方法就不在啰嗦了，依葫芦画瓢即可得到，无非就是简单的等式代换即可完成，当然最终数组下标表达式将复杂到有点没法一下看懂。

接下来我们将进入更烧脑的不等式矩阵化的方法介绍。希望到这里你还继续保持清醒。同时也希望大家根据上篇文章的提示及时的复习了矩阵运算尤其是矩阵乘法的相关知识，最好也复习了线性代数的内容，这样后面的分析方法对你而言将没有什么难度。

为了便于理解，让我们首先来分析一个简单的for循环，以便理解for语句背后的数学含义，比如：

for( i=1; i<=10; i++)

很容易我们可以发现这个for循环其实规定了 $i$ 的取值范围是

$\leqslant i \leqslant 10$ 之间的整数。通常这个不等式被理解为两个独立的不等式： $\leqslant i 和 i \leqslant 10$ 。对于一个嵌套更深的多层循环，比如：

for(i=0; i<=10; i++)
	for(j=i; j<=10; j++)
//（…这可能是一个典型的冒泡排序法的双层循环写法…）

这样的嵌套循环，根据前面简单的例子，可以理解为4个不等式 $\leqslant i,i \leqslant 10,i \leqslant j,j \leqslant 10$ 的组合。当然这里 $i$ 和 $j$ 都是整数变量，同时注意其增量是 1（你应该已经明白怎么把增量不是 1 的循环仿射变换成 1 的方法了，所以这个地方为了便于后续分析的理解，我们就不再举增量不为 1 的循环的例子了），所以这些等式精确的描述了原来的双层循环的循环迭代空间的边界。

当然这样的独立的不等式对于整体化的分析来说是不方便的，虽然在数学上以及对原来程序的理解上这样都是正确的，但是对于进一步的分析来说，在描述迭代空间的时候，这四个独立的不等式都需要同时列出，才能最终明白 $i$ 和 $j$ 之间的关系。因此这时我们需要威力更强大的数学工具——矩阵出场来hold这种局面了。

首先我们教条式的先来看下在数学上如何用矩阵化的方式来描述一个迭代空间（改编自《编译原理》紫龙版）：

在一个深度为d的循环嵌套中的迭代可以用数学方式表示为
$\{ \vec x \ 在 \ \mathbf{Z}^{d} \ 中 \ | \ B \times \vec x + \vec b \geqslant \vec 0 \}$
其中

$\mathbf{Z}$ 为整数集合——包括正整数、负整数和零； $\mathbf{Z}^{d}$ 是 $\mathbf{Z}$ 上的 $d$ 维向量空间;
$B$ 是一个 $d$ 列的整数矩阵；
$\vec b$ 是一个有 $d$ 个分量的整数向量；
$\vec 0$ 是中 $\mathbf{Z}^{d}$ 的零向量；

这个定义我想学过线性代数的同学都能看的明白，无非将 “ $=$ ” 换成了 “ $\geqslant$ ” 。其中大部分的知识基本都是高中数学的内容，比如 $\mathbf{Z}$ 表示整数集合之类、 $d$ 维向量、 $d$ 维矩阵等。

以前面的循环为例：

for(i=0; i<=10; i++)
	for(j=i; j<=10; j++)

我们可以得到 $d = 2$ ，向量 $\vec x = (i,j)$ 接下来比较难的就是如何将之前等价的4个不等式 $\leqslant i,i \leqslant 10,i \leqslant j,j \leqslant 10$ 综合变换为矩阵 $\mathbf{B}$ 和向量 $\vec b$ 了。这时我们采取的策略就是逐个按照行来整理的方法，然后拼装出最终的矩阵和常数向量。

比如对于 $\leqslant i$ 我们整理得到 $\times i + 0 \times j + 0 \geqslant 0$ ; 以此类推有：

$\begin{aligned} & i \leqslant 10 \Longrightarrow -1 \times i + 0 \times j + 10 \geqslant 0 \\[1ex] & i \leqslant j \Longrightarrow -1 \times i + 1 \times j + 0 \geqslant 0 \\[1ex] & j \leqslant 10 \Longrightarrow 0 \times i + -1 \times j + 10 \geqslant 0 \end{aligned}$
需要注意的是不等式两边乘以负数时不等号要改变方向，这里用到的基本都还是初中时的数学知识，但愿你还没有还给老师。

接下来就是要将系数和变量分离了，线性代数成绩 $\geqslant 90$ 分的同学肯定已经觉得我太啰嗦了。其实我只是在这里想把这个重要的变换技巧交给那些“忘记”了的同学，以便更多的人都能看懂并记住这个方法。

首先复习下矩阵乘法的口诀“列等行、行乘列加、等行等列”（注：该口诀为本人高中时被逼之后的大胆“创造”，本人保留原创权。这个口诀我一直牢记至今，现在再次分享给大家，原版是“列等行、行乘列、等行等列”，因为考虑到此处的乘法是类似向量的点乘，所以多了一个加号，防止被简单的理解为乘积，而忘记了加）。具体含义就是说做乘法的两个矩阵，其中乘号左边的矩阵列数一定要等于乘号右边的矩阵行数（注意矩阵乘法不可交换，必须区分左右。同时这个“列等行”的要求是矩阵乘法的先决条件，不等就不能做乘法），然后用左边矩阵的行点乘右边矩阵的列（点乘就不解释了，你懂的），放在结果矩阵对应左边矩阵的行号和右边矩阵的列号上，最后的结果矩阵的行数就等于左边矩阵的行数，而结果矩阵的列数就等于右边矩阵的列数。简单举例如下：
$\begin{aligned} & \begin{bmatrix} 1 & 3 \\ 2 & 2 \\ 3 & 1 \end{bmatrix} \times \begin{bmatrix} 1 & 2 & 3 \\ 3 & 2 & 1 \end{bmatrix} \\ = & \begin{bmatrix} 1 \times 1 + 2 \times 2 + 3 \times 3 & 3 \times 1 + 2 \times 2 + 1 \times 3 \\ 1 \times 3 + 2 \times 2 + 3 \times 1 & 3 \times 3 + 2 \times 2 + 1 \times 1 \end{bmatrix} \\ = & \begin{bmatrix} 14 & 10 \\ 10 & 14 \end{bmatrix} \end{aligned}$
上面就是2行3列的矩阵乘以一个3行2列的矩阵，对应口诀的含义来说，左边矩阵的列数3等于右边矩阵的行数3，所以两个矩阵能够相乘，就是“列等行”；然后用左边矩阵的每一行点乘右边矩阵的每一列，就是“行乘列加”；最后结果矩阵为2行2列等于左边矩阵的行数，等于右边矩阵的列数，就是“等行等列”。

接着我们返回到那一堆推导出的不等式组上：
$\begin{aligned} 1 & \times i + 0 \times j + 0 \geqslant 0 \\[1ex] -1 & \times i + 0 \times j + 10 \geqslant 0 \\[1ex] -1 & \times i + 1 \times j + 0 \geqslant 0 \\[1ex] 0 & \times i + -1 \times j + 10 \geqslant 0 \end{aligned}$
这里先遮住所有的常数及之后的 $\geqslant 0$ ，然后变成下面这样：
$\begin{aligned} 1 & \times i + 0 \times j \\[1ex] -1 & \times i + 0 \times j \\[1ex] -1 & \times i + 1 \times j \\[1ex] 0 & \times i + -1 \times j \end{aligned}$
这时不要管每行表达式中复杂的公式，直观看上去这就是一个4行1列的矩阵，然后反过来用之前的口诀“等行等列”，得知这必定是一个4行n列矩阵和一个n行1列矩阵的乘积。

接着我们看到每行的公式都是简单的乘积和加法组成的多项式（线性！），此时反用口诀“行乘列加”，因为只有一个加号，所以这是个二项式，此时立刻可以知道 $n = 2$ ，即左边的矩阵是个4行2列的矩阵，右边的矩阵是一个2行1列的矩阵，看上去像下面这样：
$\begin{aligned} & \begin{bmatrix} \cdots & \cdots \\ \cdots & \cdots \\ \cdots & \cdots \\ \cdots & \cdots \end{bmatrix} \times \begin{bmatrix} \cdots \\ \cdots \end{bmatrix} \end{aligned}$
接着从加号前后的每个同类项中提取公因式作为右边矩阵的对应行，此处第一行的公因式是i，第二行的公因式是 j，剩余的每列因式就是左边矩阵的对应列，于是得到下面的矩阵乘法的公式：
$\begin{aligned} & \begin{bmatrix} 1 & 0 \\ -1 & 0 \\ -1 & 1 \\ 0 & -1 \end{bmatrix} \times \begin{bmatrix} i \\ j \end{bmatrix} \end{aligned}$
至此我们就知道了之前公式中的矩阵 $\mathbf{B}$ 和向量 $\vec b$ 。此处用的是列向量形式，因为为了还原为矩阵乘法形式的方便性我们将不等式竖着排列到了一起，大家可以试着将原来的不等式组横着排列，然后重复之前描述的步骤，就可以得到行向量的形式。不过那时公式就变成了 $\vec{x^*} \times \mathbf{B^*}$ 其中 $\mathbf{B^*}$ 表示矩阵 $\mathbf{B}$ 的转置。熟悉线性代数的读者应该明白这是矩阵乘法的特殊交换律，即 $\mathbf{A} \times \mathbf{B}=\mathbf{B^*}\times\mathbf{A^*}$ 其中 $\mathbf{A}、\mathbf{B}$ 表示矩阵。当然如果你记住了这个规律，将列向量表示法变成行向量的表示法就不需要重复之前的那个步骤了，直接转置矩阵，并颠倒乘法的左右顺序即可。还要注意的就是这里将向量也理解为特殊的 1 行或 1 列的矩阵。

接着将之前遮住的常数列变成单独的一个列向量，并把 $\geqslant 0$ 简单挪回来就得到：

$\begin{aligned} & \begin{bmatrix} 1 & 0 \\ -1 & 0 \\ -1 & 1 \\ 0 & -1 \end{bmatrix} \times \begin{bmatrix} i \\ j \end{bmatrix} + \begin{bmatrix} 0 \\ 10 \\ 0 \\ 10 \end{bmatrix} \geqslant 0 \end{aligned}$
当然要注意这里的0就是公式中说的向量 $\vec 0$ 。

至此我们就把循环迭代空间的矩阵化的主要方法讲解完了。进一步的还可以将这个表达式“齐次化”扩展成：
$\begin{aligned} & \begin{bmatrix} 1 & 0 & 0 \\ -1 & 0 & 10 \\ -1 & 1 & 0 \\ 0 & -1 & 10 \end{bmatrix} \times \begin{bmatrix} i \\ j \\ 1 \end{bmatrix} \geqslant 0 \end{aligned}$
这样我们看到左右的矩阵分别多了一列和一行，大家可以通过矩阵乘法验证这个矩阵乘积的结果和之前的矩阵乘积结果以及之前得到的不等式组是一样的。

当然“紫龙书”（指《编译原理》中文第二版对应英文原版第三版）上没有说到这种“齐次化”的扩展，我只是根据多年学习3D数学处理的经验做了这个扩展，不知道会不会对后续的处理有帮助，待我继续学习“紫龙书”后再回头来看这个问题。当然无论如何，掌握这种“齐次化”的变换方式，总是有好处的。实际上这也就是将一个仿射变换变成了高一维空间中“超平面”上的线性变换。

当然熟练掌握这个方法之后，对于这些与循环体等价的不等式，我们可以直接通过“观察”的方法，大胆的写出矩阵变换后的形式。例如：

for(i=0; i<=6; i++)
	for(j=i; j>=0; j--)

来说，虽然我们要先进行一下仿射变换才能变成我们要求的“连续”空间，但是此处我们发现j的增量绝对值为1，所以不变换也行。因此“观察”后得到：
$\begin{aligned} & \begin{bmatrix} 1 & 0 \\ -1 & 0 \\ -1 & 1 \\ 0 & 1 \end{bmatrix} \times \begin{bmatrix} i \\ j \end{bmatrix} + \begin{bmatrix} 0 \\ 6 \\ 0 \\ 0 \end{bmatrix} \geqslant 0 \end{aligned}$
再观察一个例子：

for(i=0; i<=9; i++)
	for(j=0; j<=8; j++)

变换后得到：

$\begin{aligned} & \begin{bmatrix} 1 & 0 \\ -1 & 0 \\ 0 & 1 \\ 0 & -1 \end{bmatrix} \times \begin{bmatrix} i \\ j \end{bmatrix} + \begin{bmatrix} 0 \\ 9 \\ 0 \\ 8 \end{bmatrix} \geqslant 0 \end{aligned}$
其中规律就是，对于 $d$ 层循环，就先把 $d$ 层循环的每个下标变量列为列向量，矩阵 $\mathbf{B}$ 的大小就是有 $2 d$ 行、 $d$ 列、因为每重循环有一个上界和一个下界，所以有 $2 d$ 行。对于矩阵的每一行来说，每个列的位置就对应该重循环下标变量的系数，如果循环下标变量出现在 $\geqslant$ 或者 $\leqslant$ 符号尖端所指的一侧那么系数就是 $- 1$ （更正确的说法应该是变量系数 $\times (-1)$ ），比如 $i\leqslant10$ 尖端指向 $i$ ，对应系数就是 $- 1$ 。否则开口一侧的变量系数就取 $1$ （或原系数）。当然未出现的变量系数就是 $0$ 。对于 $\leqslant 0$ ，则对应的 $i$ 和 $j$ 的系数都是 $- 1$ 。若是 $\leqslant 0$ 则 $i$ 的系数是 $- 1$ ，而 $j$ 的系数就是 $\times -1 = 1$ 。

常数列向量就将所有循环下标变量不等式中的常数列在一列。for 循环中初值表达式和判断表达式中都有可能出现常数，复杂情况下可能要对上下界中的常数数进行带符号求和运算。常数的符号按照出现在 $\geqslant$ 或者 $\leqslant$ 符号尖端一侧乘以 $- 1$ ，或者出现在开口一侧不变号的方式填写即可。 $0$ 常数或没有常数的就填 $0$ 。最终整个矩阵表达式的结果都是 $\geqslant \vec 0$ 。（希望你没有被绕晕。）

至此我们可以先小结一下，将循环迭代空间的分离不等式“综合”为一个矩阵表达式的好处如下：
1. 首先对于数学上的向量空间处理来说，矩阵化带来的好处就是能够整体的去看待和分析一个空间，更进一步可以通过连续的矩阵乘法及加法将一个向量空间变换到另一个等价的向量空间，也就是常说的线性变换（因为有常数列，其实应该是仿射变换更准确）。
2. 其次对于此系列文章的主题——循环体并行优化来说，将循环迭代空间进行矩阵表示之后，就可以与之后的数据空间、处理器空间（均做类似的矩阵化表示）等，通过解线性方程组（丢番图方程）的方式来分析得到数据访问冲突，也就是之前说的读/写冲突，通常就是丢番图方程组的解集，或者处理器分配、数据连续放置高速缓存等问题的答案，从而最终可以高效高速的将循环体并行化执行，以最大化利用多处理器系统的优势。这也就是为什么这系列文章最终绕到矩阵表达方式的根本原因。

(为了文章结构性果断的再割一下)

如果看到这里你还没有完全晕厥的话，就让我们进入谋杀脑细胞的阶段。现在我们要开始打通任督二脉了，请系好安全带！Let’s go!

首先我们再折回头来看最开始我们用上一篇文章中的直接仿射变换迭代空间的循环例子：

for(i = 3; i<10; i+=2)
	for(j = i; j<15; j+=2)

其实按照我们对 for 循环的理解，这个循环体实际是在说这样的数学事实：

$\begin{aligned} & i \geqslant 3; & 2 \times i < 10; \\[1ex] & j \geqslant i; & 2 \times j < 15; \end{aligned}$
其中每个循环上界表达式中的系数 $2$ ，即为各自循环下标的增量。

再结合线性代数的知识，我们知道仿射变换都可以用矩阵来表达，那么既然已经知道了该循环体对应的这些不等式，那么是不是就不需要用我们最开始说的仿射变换方法变到“连续”的空间，而直接矩阵化表示呢？答案显然是yes！

那就让我们更直接的用“观察”法来做一下吧，如果你还不适应观察法，那么就按照我啰嗦了半天的那个“倒推法”一步步自己试试了。

首先，不要太激动，细节决定成败，发现那个该死的 $<$ 没有？对头，我在讲矩阵化的时候可都是用的 $\geqslant$ 或 $\leqslant$ 哦，二者显然是有区别的。聪明的你肯定也想到该怎么办了，对头，那就是改成： $\times i \leqslant 9 , 2 \times j \leqslant 14$ 因为我们一直玩的是整数循环！