登录  
 加关注
查看详情
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

Gamebaby Rock Sun的博客

我只知道一件事情,就是我一无所知。

 
 
 

日志

 
 
关于我

曾经拥有的,不要忘记, 已经得到的,更要珍惜, 属于自己的,不要放弃, 已经失去的,留着回忆, 想要得到的,必须努力, 但最重要的,是好好爱惜自己!

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

2016-11-20 14:23:17|  分类: 并行计算 |  标签: |举报 |字号 订阅

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

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

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

首先,我们来看一个2重循环的例子(顺便诅咒一下网易Blog糟糕的代码和公式页面排版!):

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

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

     a[i,j] = ...;

这里我们先对i做变换:

   i = 2*ii+3; 变换外层循环为for(ii=0;ii<3;ii++)

接着 令 j = 2*jj + i; 注意这里有外层循环的旧下标i,所以继续带入i的替换得到j=2*jj+2*ii+3; 此时要求 2*jj+2*ii+3<15; 得到jj+ii<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

2*ii+3

ji=3

ji=5

jj

2*jj+2*ii+3ii=0

2*jj+2*ii+3ii=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=35时的的j值,这个很容易从原循环下标中验证得到,同时也只列了ii=01时的2*ii+32*jj+2*ii+3的对应值,同样这也可以很容易的从变换后的循环中推出。

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

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

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

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

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

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

1i10之间的整数。通常这个不等式被理解为两个独立的不等式:1ii10。对于一个嵌套更深的多层循环,比如:

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

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

这可能是一个典型的冒泡排序法的双层循环写法

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

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

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

循环体并行优化(二) ——多维循环迭代空间的仿射变换及循环上下界不等式的矩阵表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的博客


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

以前面的循环为例:

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

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

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

比如对于0i我们整理得到1*i+0*j+00;以此类推有

i10 -1*i+0*j+100

ij  -1*i+1*j+00

j10 0*i+-1*j+100

需要注意的是不等式两边乘以负数时不等号要改变方向,这里用到的基本都还是初中时的数学知识,但愿你还没有还给老师。

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

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

循环体并行优化(二) ——多维循环迭代空间的仿射变换及循环上下界不等式的矩阵表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的博客 

上面就是23列的矩阵乘以一个32列的矩阵,对应口诀的含义来说,左边矩阵的列数3等于右边矩阵的行数3,所以两个矩阵能够相乘,就是“列等行”;然后用左边矩阵的每一行点乘右边矩阵的每一列,就是“行乘列加”;最后结果矩阵为22列等于左边矩阵的行数,等于右边矩阵的列数,就是“等行等列”。

接着我们返回到那一堆推导出的不等式组上:

1*i+0*j+00

-1*i+0*j+100

-1*i+1*j+00

0*i+-1*j+100

这里先遮住所有的常数及之后的≥0,然后变成下面这样:

1*i+0*j

-1*i+0*j

-1*i+1*j

0*i+-1*j

这时不要管每行表达式中复杂的公式,直观看上去这就是一个41列的矩阵,然后反过来用之前的口诀“等行等列”,得知这必定是一个4n列矩阵和一个n1列矩阵的乘积。

接着我们看到每行的公式都是简单的乘积和加法组成的多项式(线性!),此时反用口诀“行乘列加”,因为只有一个加号,所以这是个二项式,此时立刻可以知道n=2,即左边的矩阵是个42列的矩阵,右边的矩阵是一个21列的矩阵,看上去像下面这样:

循环体并行优化(二) ——多维循环迭代空间的仿射变换及循环上下界不等式的矩阵表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的博客

 

接着从加号前后的每个同类项中提取公因式作为右边矩阵的对应行,此处第一行的公因式i,第二行的公因式是j,剩余的每列因式就是左边矩阵的对应列,于是得到下面的矩阵乘法的公式:

循环体并行优化(二) ——多维循环迭代空间的仿射变换及循环上下界不等式的矩阵表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的博客


至此我们就知道了之前公式中的矩阵B和向量

。此处用的是列向量形式,因为为了还原为矩阵乘法形式的方便性我们将不等式竖着排列到了一起,大家可以试着将原来的不等式组横着排列,然后重复之前描述的步骤,就可以得到行向量的形式。不过那时公式就变成了

循环体并行优化(二) ——多维循环迭代空间的仿射变换及循环上下界不等式的矩阵表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的博客

其中B*

表示矩阵B的转置。熟悉线性代数的读者应该明白这是矩阵乘法的特殊交换律,即

循环体并行优化(二) ——多维循环迭代空间的仿射变换及循环上下界不等式的矩阵表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的博客

 

其中AB表示矩阵。当然如果你记住了这个规律,将列向量表示法变成行向量的表示法就不需要重复之前的那个步骤了,直接转置矩阵,并颠倒乘法的左右顺序即可。还要注意的就是这里将向量也理解为特殊的1行或1列的矩阵。

接着将之前遮住的常数列变成单独的一个列向量,并把≥0简单挪回来就得到:

循环体并行优化(二) ——多维循环迭代空间的仿射变换及循环上下界不等式的矩阵表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的博客

 

当然要注意这里的0就是公式中说的向量

至此我们就把循环迭代空间的矩阵化的主要方法讲解完了。进一步的还可以将这个表达式“齐次化”扩展成:

循环体并行优化(二) ——多维循环迭代空间的仿射变换及循环上下界不等式的矩阵表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的博客

 

这样我们看到左右的矩阵分别多了一列和一行,大家可以通过矩阵乘法验证这个矩阵乘积的结果和之前的矩阵乘积结果以及之前得到的不等式组是一样的。

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

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

for(i=0;i<=6;i++)

    for(j=i;j>=0;j--)


来说,虽然我们要先进行一下仿射变换才能变成我们要求的“连续”空间,但是此处我们发现j的增量绝对值为1,所以不变换也行。因此“观察”后得到:

循环体并行优化(二) ——多维循环迭代空间的仿射变换及循环上下界不等式的矩阵表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的博客

 

再观察一个例子:

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

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


变换后得到:

循环体并行优化(二) ——多维循环迭代空间的仿射变换及循环上下界不等式的矩阵表示法 - Gamebaby Rock Sun - Gamebaby Rock Sun的博客

 

其中规律就是,对于d层循环,就先把d层循环的每个下标变量列为列向量,矩阵B的大小就是有2*d行、d列、因为每重循环有一个上界和一个下界,所以有2*d行。对于矩阵的每一行来说,每个列的位置就对应该重循环下标变量的系数,如果循环下标变量出现在>=或者<=符号尖端所指的一侧那么系数就是-1(更正确的说法应该是变量系数*-1),比如i<=10尖端指向i,对应系数就是-1。否则开口一侧的变量系数就取1(或原系数)。当然未出现的变量系数就是0。对于i+j0,则对应的ij的系数都是-1。若是i-j0i的系数是-1,而j的系数就是-1*-1=1

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

至此我们可以先小结一下,将循环迭代空间的分离不等式“综合”为一个矩阵表达式的好处如下:

首先对于数学上的向量空间处理来说,矩阵化带来的好处就是能够整体的去看待和分析一个空间,更进一步可以通过连续的矩阵乘法及加法将一个向量空间变换到另一个等价的向量空间,也就是常说的线性变换(因为有常数列,其实应该是仿射变换更准确)。

其次对于此系列文章的主题——循环体并行优化来说,将循环迭代空间进行矩阵表示之后,就可以与之后的数据空间、处理器空间(均做类似的矩阵化表示)等,通过解线性方程组(丢番图方程)的方式来分析得到数据访问冲突,也就是之前说的读/写冲突,通常就是丢番图方程组的解集,或者处理器分配、数据连续放置高速缓存等问题的答案,从而最终可以高效高速的将循环体并行化执行,以最大化利用多处理器系统的优势。这也就是为什么这系列文章最终绕到矩阵表达方式的根本原因。

  评论这张
 
阅读(328)| 评论(3)

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018