4 数据并行性

BACKWARD


基本物理现象,例如热量产生和消失特性和电子信号速度,在理论和实践上限制了单处理机系统的计算速度。虽然这些限制目前在大约"gigafiop"(每秒十亿个数值操作)范围内,但是一些计算科学的当代应用需要更高的速度,例如多数"grand challenge"中的应用问题. 通过增加处理器的个数来提高计算能力比增加单处理机的速度更可行。将来有可能所有涉及大量计算的应用将使并行技术得到极大应用。

应用程序可以使用两种(或者其中一种)基本形式,在这里被称作“数据并行”和“处理并行”。数据并行指同时对许多数据对象执行相似的计算。这种并行的原型(尤其对于计算科学应用)是对一个数组的所有元素同时进行操作——例如,用一个给定值去除数组的每个元素(比如在矩阵归约时将主元行基准化)

针对特定目的,这里数据并行将意味着对数组元素的并发操作。处理并行指并行地执行不同处理过程,在这里过程是任意顺序的操作。例如,一个子例程是一个过程,就象任意一个语句块一样。处理并行的一个越来越重要的形式是同时计算两个表达式——例如过程调用中的两个实参。已经开发了几个指定处理并行的途径,处理并行的标准也将出现[3].这种并行将在章的以后发表中讨论。

虽然处理并行在计算科学的许多应用中都很重要,但是由于在这些应用中数组经常被使用,因此在大多应用中数据并行可能更有用,甚至在一些应用中成为关键所在。所以这一节集中讨论并行数组操作,重点在于Fortran 90中这种操作的丰富的集合。

4.1数组操作

APL可能是第一个把数组当作完整的数据对象而不仅仅是一个数据对象的简单集合的程序设计语言. 操作可以对整个数组进行,如

在这里A,B,C都可以是数组,比如A,B,C可能都是200×300的二维数组。这样一个操作的含义是对所有200个i值和300个j值计算 Ci,j=Ai,j+Bi,j,在这种情况下,共有60,000个涉及数组元素的独立计算。因此APL模型是对数组所有元素的并发计算。(可是很少有APL的实现利用这种概念的并行性)。

除了支持与上面所示的加法操作相似的一般数学操作,APL还定义了其它以数组为单位的操作.这些操作包括归约操作(例如,+/A表示求数组A的所有元素和), 构造操作(例如, i4构造一个四个元素的向量, 值为1,2,3,4), 和查询操作(例如, 如果ρB返回数组B的状态——一个向量,它的大小是B的秩,它的元素值是B的相应维的大小——那么*/ρB计算数组B中元素的个数)。所有这样的操作可被用来组成更复杂的数组值表达式(例如,对一个一维数组Q,*/Q+i*/ρQ从右向左计算(这是APL!),产生一个大小为Q的向量1,2,3,…,并与Q相加,结果元素相乘;如果Q是向量(3,2,4),那*/Q+i*/ΦρQ=*/((3,2,4)+(1,2,3))=*(4,4,7)=56)。在数据并行之前引入APL变得完全实用( )

Fortran90的数组操作实际上提供了所有的APL的逐个元素(elenent-by-element)数据并行特征,并且摒弃了APL的缺点。这些操作使用类拟Fortran的表达式计算规则,是对标量运算、函数、表达式的自然扩展,归约,构造和查询操作通过增加一些富有意义命名的内部函数来提供。为了保证这些操作能够在当代并行处理系统上有效地执行,进行了仔细的考虑

一般说来,除了在少数几个上下文中表达式限制为标量外,任意Fortran90表达式都可以含有数组操作,结果可以是数组值。Fortran77仅允许标量表达式;Fortran 90中几乎所有这样的表达式都可以是数据并行数组值。(在控制上下文中,如IF语句控制条件(标量逻辑表达式),DO循环索引表达式中,以及I/O说明符(如单元数,文件名,开放语句说明符等等)中,都需要标量表达式)。数组操作的例子如下所示(变量中的任意一个或全部都可以是数组):

注意,在这些情况中,数组表达式与标量表达式无法区分——用户需要从其它上下文获知这些变量已被定义为数组——但是每个表达式可能代表上百万次并行计算。如果A,B,C,P,Q和R是二维数组,Z和V是一维数组,那么这三个语句可以写成如下的等价形式(证明确定议了数组操作):

函数可以被定义成数组值,因此可以作为数组值表达式的操作数,这以及与数组相关的内部函数在4.4节加以说明。

因此,用数组代替标量作为操作数,那么element-by-element数据并行的数组值的表达式是标量表达式的直接自然的扩展/综合。

一致性要求,建立一个数组表达式的基本要求是操作数的一致性,这意味着数组操作的每个操作数必须具有相同的维数,并且每一维具有相同个数的元素——也就是说一致的数组具有完全一样的形式。当然,这种操作的结果与操作数一致,数组结果的每个元素的值是数组操作的相应元素的标量计算。

这样,如果A和B是2×3的数组,如下所示

则 A + B 的结果为:

A X B的结果为:

如果一个表达式中的运算不只一个,那么第一个子表达表的(赋值的数组)结果是第二个运算的操作数,以此类推. 例如,对于上述的 A 和 B,在表达式A+B×A中, B×A的结果被加到A上,这样A+B×A的结果是

注意,一个3×2的数组与2×3的数组不一致——虽然它们有相同的秩和元素个数,但相应维的大小不同——因此这样两上数组不能作为同一个数组运算的操作数,这个基本的一致性规则的唯一例外是操作数中的一个是标量,在这种情况下,标量被“广播”到与另外一个操作数相一致的数组中,这个广播数组的每个元素值是就是标量的值,例如,B+2是一个有效的数组操作(假设B为前面定义的数组),B+2的值是

在数组操作中,(广播)标量通常用在数组初始化和换算时:

最后这个例子说明了Fortran 90 数组操作的一个关键方面:在一个数组值的赋值语句中,在赋值之前,右边的数组值已被算出,否则右边的数组值在运算未完成之前可能会受到影响(虽然在这个简单的例子中不会这样),因此,Fortran 90的概念模型是:先对右边数组值的所有元素并行计算(或以任意顺序)然后再赋值,任何能保证这种行为特性的实现是允许的。

体现这种规则重要性的一个例子是4.5节的选主元过程。在这个例子中,主元行用如下的数组操作加以规格化:

这个操作使用一个数组片段G(P,:);数组片将在4.2节描述,在这个例子中G(P,:)是矩阵G的第P行(主元行),G(P,K)是主元素(K是主元列)。规格化操作定标主元行,使得主元的值是1. 注意,如果这个元素的值在右边的计算完成之前已经变成了1,那么这一行就没有被适当地规格化(在顺序编程中这是一个典型的错误). 因此,数组操作不能被认为是对数组元素的“循环”操作,这种"循环操作"暗示着顺序的操作;通常当涉及到赋值的时候,把数组操作当作循环就会产生错误结果. 数组操作应当被看作是积分/并行计算。

数组构造子(constructor)

数组值可以用数组构造器来显式构造,如果需要的合成数组维数大于1,那么用RESHAPE内部函数来构造;数组构造器构造一个一维数组。数组构造器只不过是结果元素值的一个列表, 这个列表封闭在在定界符(/... /)中,各元素值之间用逗号分开。前述的APL例子可以用Fortran 90的数组构造器写成

(尽管有更一般的方法来指定这样一个“iota”顺序)。这个例子中,每个元素都是一个简单标量常量. 通常,数组构造器中的任一元素可以是任意的标量表达式。然而如果它们全是常量,那么这样一个构造器(可能结合RESHAPE函数一见后文)表示一个数组常量,可以在PARAMETER定义中出现. 因此,数组构造器(结合RESHAPE)是Fortran 90中表示数组常量的一种方法。

如果每个元素都须显式列出,那么对于定义一个非常大的数组值来说,数组构造器就不是那么实用了。因此,除了标量表达式外,Fortran 90还提供了两种方式来定义灵敏组构造器的表项,即隐含do结构和数组表达式,第一种具有如下形式

索引变量是一个标量整型变量, 用来作为迭代索引,与Fortran 77 的输入/输出说语句中的隐含do循环方式一样. 数组构造器中隐含do结构的一个简单应用是产生任意的“iota”顺序。例如,数组构造器

产生元素值为1,2,3,4,5...n的向量;如果n为4,则结果是(/1,2,3,4/)。再例如,一个具有一百万个1和0交叉序列的向量,

(/1,0,1,0,1,0,1,$\1dots$/),可用(/(1,0,j=1,500000)/)来定义。隐含do循环仅仅以指定的次数来复制这个列表,如果索引变量是表达式列表一个表达式的操作数,则用相应的索引变量值来对那一项进行复制。在隐含do表达式列表中的项可以是数组构造器本身允许的三种形式之一——即标量表达式,隐含do结构,和数组表达式。上面所举的两个例子在隐含do列表中仅使用了简单的标量表达式。

一个任意维数的数组表达式可以出现在数组构造器中。例如,如果A是一个1000×1000的数组,那么

是一个含有100万个元素的数组构造器,每个元素值比相应的数组A的元素值大1.3。A+1.3的所有元素按照与Fortran相似的“以列为主”的顺序被放到数组构造器中,即一列接一列地运行完第一维,然后再运行第二维,这样(/A+1.3/)等价于

隐含do结构可以用来实现不同用顺序. 例如,如果需要A+1.3的逐个行向量,而不是列向量,则下述数组构造器可以实现:

最后,RESHAPE内部函数的简单形式可以用来将(一维)数组构造器的结果重定形为需要的形状. 这种简单定义形式为

这里形状向量的每个元素代表所需要的数组形状的每一维,元素的值代表该维的元素个数。例如, 可用如下形式来定义一个名为 Ident_1000的1000×1000的实型单位矩阵

因而,对于构造数组值(包括数组常量)来说, 数组构造器结合RESHAPE内部函数是一个极其有力的工具。

屏蔽数组赋值

一个“屏蔽(mask)”是一个逻辑类型的数组。一个屏蔽数组操作是这样的一个操作,在这个操作中,用与操作结果相一致的屏蔽来指明只有并行元素操作的一个子集被执行。这种功能在一些内部函数中和数组赋值中都是可用的。在后一种情况下,数组的赋值是在WHERE语句的屏蔽码控制下执行,一般形式为:

WHERE屏敝码必须与赋值语句左边的数组(与赋值号右边的表达式相一致)相一致。当屏蔽码的元素值为真(.TRUE.)时,才对相应的元素执行赋值操作,否则为假(.FALSE.)时,不执行赋值操作,使用屏蔽数组赋值的典型例子是:

在这个例子中,对于数组C中那些值为0(或负数)的元素,不进行除操作和赋值操作. 根据一致性原则,数组A,B和C都相容,因此(数组值)逻辑表达式C.gt.0是一个与这些数组相一致的屏蔽码。通常,象所举的例子一样,一个WHERE屏蔽是一个涉及一个或多个赋值操作数的逻辑表达式。

一个屏蔽码可以控制多个赋值语句,这种情况下WHERE是一个语句块的形式:

用这种方法可将任意多个数组赋值放在一起; 当然,它们都必须与屏蔽码相一致.

上面所述的WHERE形式使赋值语句左边的数组中的一些元素未被赋值。WHERE的块形式的一个扩展,ELSEWHERE操作,允许在屏蔽码为假(.FALSE.)时,将左边相应的数组元素赋值,这种形式定义为:

这种形式的简单例子是:

这里当数组C中元素值小于或等于零时,将数组B的相应元素值赋给H中的对应元素。这是WHERE的一个重要形式,因为它产生了一个完全定义的数组H,这个数组可在后续数组操作中使用。如果没有ELSEWHERE语句,那么数组H可能未被完全意义,这样它就不能在其它数组表达式中使用。

虚拟形状的哑参

Fortran77中允许(无下标)数组名出现的一个地方是作为过程参数出现。在 Fortran 77 的局限性的数组概念下, 这对于向过程传递数组对象是足够的. 但是,由于一个Fortran 77数组总是占用一块连续的存贮区,所以只需将这个存贮区的位置传给过程,过程可以把它当为具有相同形状的一个(连续)数组的位置,或不同形状的一个数组, 甚至是当作一个标量。然而这对于数组是完全独立的对象,采用一致性规则以及数组对象不必占用连续空间(比如具有许多数组片段——见4.2节)的Fortran 90语言来说,这就显得不足够了。

在Fortran 90中,任意数组表达式都可用来作为过程调用的实参,被调用过程必须适当地将这些参数处理为数组值对象,如果仅仅传递一个简单的位置参数,就不一定能处理得当。虚拟形状的哑参解决了这个问题。这些哑参调节传递的数组“描述信息”,描述信息除了包括数组的位置外,还包括数组的其它信息。这些附加的信息包括被传递的数组对象的秩(维数),每个元素的类型和大小,每一维元素的个数以及每一维的“跨距”;跨距表示同一维的元素间的距离,因此说明了和相连位置的偏移. 所以,可将任意一个数组表达式传递给虚拟形状的哑参。(也可传递给一个“老式”的哑参,但那将可能导致昂贵的现场后的连续临时存贮区的存入和提取)。

一个虚拟形状的哑参用冒号定义每一维,如下例所示,T是标量,U是一个二维虚拟形状的数组,V是一个一维虚拟形状的数组。

调用CALC3时,任意二维数组表达式(实型)可被传递给U,任意一维数组表达式可被传递给V;相反,一个二维实型数组必须传递给U,一个一维实型数组必须传递给V。U和V的定义中的冒号告诉CALC3接收由调用程序提供的描述信息。然后U和V就正确地表示了实参表中相应的数组对象,并且可在CALC3过程体的数组表达式中使用。如果CALCS是调用程序的一个内部过程,或是被调用程序使用的模块中的一个模块过程,那么虚拟形状的哑参与相应实参之间的适当联系将被透明地实现。但是, 如果CALC3是一个外部过程,那么是调用程序中必须为CALC3提供一个接口块,这样它才能知道虚拟形状的哑参是接收者,从而传递有效的描述信息;否则调用程序不能假定哑参是虚拟形状的,因此必须提供一个连续的实参,必要的话需将数组聚集或提取。CALC3的一个适当接口块为:

4.2 数组片段

含有多于一个元素的数组的一部分叫一个数组片段。在很多情况下一个数组操作只需在一个数组片段而不是整个数组上进行. 前面所举的将一个矩阵的主元行规格化的例子就是一个典型。在这个例子中,计算只牵涉到矩阵的一行而不是整个数组. 在这种情况下,数组片段是二维数组中的一行;一般地,一个数组的任意直线型的子集都能作为一个数组片段,因此可以作为一个对象在数组操作中使用. 一个数组片段可以是任意维数,这个维不超过定义该片段的数组的维数。

一个标量数组元素,很自然地由数组名和每一维的一个下标值来指定。它的一般形式是我们已经熟知的:

这里下标个数是数组的维数,每个下标是一个标量的整型表达式,或者就是一个标量下标。一个数组片段是通过用至少一个“片段下标”代替标量下标来定义。一个片段下标是一系列代表那一维的标量下标值,因此,一个片段下标可被认为是(或构造成)一个下标值的一维数组,叫做向量下标。如果(仅当)一个标量下标由一个向量下标代替,那么结果是一个一维数组片段;如果两个标量下标被向量下标代替,则结果是一个二维数组片段,等等。一个数组片段的维数与它所具有的向量下标个数相等。

现在举个例子,考虑如下所示的大小为5×6的数组Q,Q的三个片段用粗体字表示:整个第二列(一个一维片段),子疑辖堑�2×2矩阵(一个二维片段),第五行的后半部分(一个一维片段)。

     Q=[13  11  25  2   1   9  ]
       |9   3   31  14  52  27 |
       |16  45  54  36  15  20  |               (8)
       |7   20  18  19  8   19  |
       [37  56  54  66 77  90 ]
Q((/1,2,3,4,5/),2) = (/11,3,45,20,56)  ! 第二列
Q((/1,2/),(/5,6/)) = |1  9  |          ! 右上角
                     |52 27 |
Q(5,(/4,5,6/)) = (/66,77,90/)          ! 第五行的最后一部分

注意这些例子中的所有向量下标都能写成隐含do结构的形式:

隐含do形式更易扩充,而且,对于大的片段,它比显式列表更紧凑. 隐含do结构对空间规则但非连续的向量下标也很有用。例如:

隐含do形式如此常用,因此对于索引do控制三元组,又提供了一种更可读的简写表示法,叫做“三元组下标”。

一个三元组下标只是用逗号代替冒号将隐含do的控制值分开,最后一个控制值(增量或步长)是可选的。因此,上面四个例子用三元组简记法可以写成更清楚的形式:

更进一步的简化形式是,如果三元组的第一个值被省略,则那一维的下界为假定值; 如果第二个值被省略,则那一维的上界为假定值。这样,表示这四个片段的最紧凑的方法是:

Q(:,:)是一个片段,实际上是整个数组Q,这解释了本节前面所举的一个例子。(注意Q(:,:)也具有虚拟形状的哑参定义形式;但是在这里不存在二义性,因为如果这种简记法在一个数组表达式中出现,它总是代表一个数组片断)

虽然上面例子使用了数组构造器,但是再回到向量下标的更一般形式,它允许使用任何一维数组表达式。唯一的要求是向量下标的每一个元素值对那一维来说应该是一个有效下标值。向量下标的普通形式是一个一维整型数组名(或片段),其元素值先前已被建立。这种形式对于间接访问极其有用,例如对表格内部的索引;表格元素可以通过用一个包含想要的表格索引值的数组给表格数组加上下标来检索(或设置)。

最后两个例子将结束对数组片段的介绍. 为了明确起见,又用到了数组构造器,但实际上更可能用到的是更简洁的一维数组名或片段,首先,一个数组片断不一定和上面的例子一样容易用图形表述出来. 对上例中定义的数组Q, 片段

Q((/2,5,3/),(/6,4/))

它代表组片段

       [Q2,6, Q2,4]   [27  14]
       |Q5,6, Q5,4| = |90  66|                  (9)
       [Q3,6, Q3,4]   [20  36]

这个片段可在任意一个数组表达式中使用,在该表达式中,3×2的数组对象是有效的. 它也可以出现在数组赋值语句的左边,这时右边表达式结果的第(1,1)个元素值被赋给Q2,6,右边第(3,2)个元素值被赋给Q3,4,以此类推。

一个向量下标包含的元素个数可以大于那个数组中那个维的大小,这种情况下,由于所有值必须在数组维数范围内,所以有重复值。实际上,既使向量的大小比数组维数小一个向量下标的下标值也可重复(唯一的要求是下标值必须在范围之内). 这两种情况在下例中都将说明,这个例子从Q的元素中定义了一个7×4的片段。

                                           | Q4,1  Q4,4  Q4,4  Q4,3 |
                                           | Q1,1  Q1,4  Q1,4  Q1,3 |
                                           | Q2,1  Q2,4  Q2,4  Q2,3 |
        Q((/4,1,2,3,4,2,5/),(/1,4,4,3)) =  | Q3,1  Q3,4  Q3,4  Q3,3 |
                                           | Q4,1  Q4,4  Q4,4  Q4,3 |
                                           | Q2,1  Q2,4  Q2,4  Q2,3 |
                                           | Q5,1  Q5,4  Q5,4  Q5,3 |

注意,这个片段的第一行与第五行相同,第三行与第六行相同。因此在这个数组片段中Q的许多元素出现两次,Q2,4和Q4,4出现四次。在数组表达式中,具有多次出现的给定父数组元素的数组片段是完全合法的,但这样的数组片段一定不能出现在数组赋值句的左边。

4.3动态数组

Fortran 90有三种不同的动态数组。这三种数组都允许在运行时创建,其大小由计算(或输入)值决定。这三种动态数组为:

  1. 自动数组(automatic arrays)
  2. 可分配数组(allocatable arrays)
  3. 指针数组(pointer arrays)

自动数组

自动数组是局部数组,它的大小依赖于与哑参相联的值。自动数组在进入过程时自动创建(分配),在退出过程时自动释放。自动数组的大小随过程的激活点不同而不同,如下是自动数组的一个例子:

function F18(A,N)
integer N                             ! 一个标量
real A(:,:)                           ! 一个假定形状的数组
real F18(size(A,1))                   ! 函数结果本身是一个动态数组
complex Local_1(N, 2*N+3)             ! Local_1是一个自动数组
                                      ! 其大小以 N 为基础
real Local_2(size(A,1),size(A,2))     ! Local_2是一个自动数组
                                      ! 其大小与A完全相同
real Local_3(4*size(A,2))             ! Local_3是一个一维数组
                                      ! 其大小是A 的第二维的4倍
..
end function F18

注意内部查寻函数在定义自动数组中的重要作用, 例如SIZE(它返回在指定的维中参数数组的大小)。Fortran 90提供了很多允许在定义中出现的查寻函数。数组界限和大小,字符串长度以及类型kind值,都可以用涉及这些查询函数的表达式来指定。粗略地说,一个指定表达式当被调用时是一个标量整型表达式,它的操作数值在进入过程入口时是可确定的。这样的操作数包括常数,指向内部过程的引用,以及通过哑参、模块、common和(在模块和内部过程的情况下)宿主过程(host procedure)可以访问的变量。

可分配数组

可分配数组是那些用ALLOATABLE明确定义的数组。一个可分配数组可以对一个过程是局部的,或者可以放在一个模块中,从而对应用程序的所有过程来说都是全局的。一个可分配的数组由语句ALLOCATE明确定义,由语句DEALLOCATE显式释放,或者如果它是没有指明SAVE的局部数组,可在过程的出口自动释放。(如果指明了SAVE,那么局部可分配数组可在过程的一个执行和下一个执行之间一直存在——这些数组必须用DEALLOCATE语句显式释放)。一个全局可分配数组在它被显式释放之前一直存在,它可以在分配该数组的过程以外的过程中出现。当数组大小依赖于计算值而不是哑参或模块,公共体以及主过程的变量时,使用可分配(或指针)数组,一个可分配数组的分配状态(分配或未分配)可用ALLOCATED内部函数来测试。如下是可分配数组的例子。

例子

subroutine Peach
  use Recipe                   ! 访问全局可分配数组 Jam.
  real, allocable :: Pie(:,:)  ! Pie 是一个 2-维可分配数组
  ...
  allocate  (Pie(N,2*N))       ! 分配一个局部可分配数组
  if (.not.allocated(Jam)) allocate (Jam(4*M))
                               ! 如果一个全局可分配数组还没有分配
                               ! 则分配它 
  .......deallocate(Pie)
  ....
end subroutine Peach
module Recipe                  ! Jam 是一个全局可分配数组, 它可以
  real, allocable :: Jam(:)    ! 在任何过程中使用该模块分配和释放
  ...
end module Recipe

注意可分配数组的界限的定义只是冒号。表示这些值将在以后分配时提供。这使得可分配数组的定义看起来好象与虚拟形状哑参的定义类似,大概是因为它们在维的大小的“延迟”特性在概念上相似。

指针数组

指针数组也用ALLOCATE语句明确定义,并且有任意可计算的大小以及用DELOCATE语句显式释放, 这一点与可分配数组类似。在第一节的"指针"小节中给出了指针数组的例子。这些例子也说明了目标数组和指针赋值的使用,后者不能使用可分配数组。指针数组的另外的简单例子就是在前面所举的可分配数组的例子中把“allocatable”用“pointer”来代替。

另外,指针数组可被用来作为(“指向”)其他数组和数组片段的替换入口;用指针赋值语句来建立这样的替换入口,与指针相联目标(当这种替换入口被调用时)可以是其它显式分配的数组,或者是明确指定为指针允许目标的静态或自动数组。一个指针数组的相联状态可用ASSOCIATED内部函数来测试。最后,指针数组可以是哑参和结构成分(structure components),这两种情况对可分配数组来说都不允许. 给定外表上相似的可分配数组和指针数组,这两种动态数组的根本区别是什么? 什么时期用可分配数组而不用指针数组? 指针数组包含了可分配数组的所有功能,从这一点来看,尽管永远不需要可分配数组——指针数组就足够了,但指针数组有一个效率问题。尽管指针数组必须总是指向明确的目标,只有这样才使优化变得可行,可是指针赋值使指针数组的优化比可分配数组困难得多。由于可分配数组的更多的限制特性和功能,它们比指针数组更简单、更有效。

因此, 当所有的需求是数组简单的动态分配和释放, 而且使用自动数组是不够的, 则使用可分配数组. 这种情况的一个常见的例子是当一个"工作数组"的大小取决于一次局部计算的结果的时候. 另一方面, 如果算法调用一个动态替换(alias), 例如一个宿主数组(host array)的"移动"部分, 则可能需要使用一个指针数组.

4.4 返回赋值数组的函数(array-valued functions)

Fortran 90函数的返回结果可以是赋值数组。很多内部函数总是返回数组值, 大多数本征函数可以返回数组值。另外,用户写的函数也可以返回数组值。返回数组值的函数提供两个(相关)巨大好处。第一,返回赋值数组的函数可以在数组赋值表达式中用作操作数, 允许以最自然的形式表示数组并行计算。第二,这使得用数组值子表达式来构成可并行求值的计算更方便。因此返回赋值数组的函数为将数据并行与处理并行用自然的方式结合起来提供了更多的机会。

例如,g和u是两个大的相互一致的二维数组,g+cshift(u,1,2)-cshift(u,-1,2)是一个表达式,类似这样的表达式在地震模型计算中经常出现(CSHIFT是一个内部函数,它按第二个参数所指定的次数,循环移动第一个参数(一个数组)中由第三个参数指定的那一维)。一旦考虑到数据并行,且提供给编译器并行计算任何子表达式的机会, 则这个表达式非常清楚地表示了计算特征。这也许尤其恰当,例如,如果g本身是个函数引用,那么并行计算子表达式cshift(u,1,23)-chift(u,-1,2)的值和 g 的值也许会更有利。

返回赋值数组的内部函数有两类,一类叫作“变换(transformational)”函数,一类叫做“基本(elemental)”函数。一个变换函数接收一个输入数组,产生一个不同的数组作为结果--它将输入数组“转换”成其它形式,甚至可能是不同形状的数组,或者是标量。(变换函数甚至能将标量转换成数组)。CSHIFT是变换函数的一个简单例子,其结果与(第一个)参数相一致。内部函数MATMUL(矩阵相乘)也是变换函数的一个例子,它返回一个与(任一个)参数形状不同的数组结果。归约函数,SUM,PRODUCT,COUNT,等等都是变换函数的例子,它们将数组参数归约成标量结果。Fortran 90的数组变换内部函数(共42个)如表4所列:

变换内部函数 说明
环境查询函数(9) 见 3.5 节
数组函数(21) 见下文
ASSOCIATED 检验指针的相连状态
BIT-SIZE 一个整型数中位的个数
DOT-PRODUCT 两个向量的数学点积
KIND 见3.1节和3.2节
LEN 字符串的长度
MATMUL 数学矩阵乘积
PRESENT 检验可选参数的存在
REPEAT 复制一个字符串
SELECTED_INT_KIND 见3.1节和3.2节
SELECTED_REAL_KIND 见3.1节和3.2节
TRIM 移走串中的结尾空白符
TRANSFER 将位模式转换为不同类型

表 4

基本内部函数多数是用标量哑参定义的. 这些函数可用数组实参调用,返回一个与实参相一致的结果数组。结果数组的每个元素值就是如果函数用实参的相应(标量)元素来调用时所获得的值。因此一个基本函数被自动(或在概念上并行地)应用到实参的每个元素中。任意一个普通的计算内部函数可被直接调用,例如,在

中,X可以是标量,这时COS函数返回一个标量结果; 或者X是一个数组(任意维),此时COS函数返回一个与X相一致的数组值结果,如果把前面所举的类似地震模型的例子修改为

那么表达式中的每一项成为一个内部函数调用,该函数返回与 g 和 u 相一致的结果. 第一项是基本调用,其他两项是变换调用。Fortran 90的所有这108个内部函数中,除了上面所列的42个变换函数外,都可直接调用。

注意,基本函数调用可被认为是一些独立的标量函数调用,一个变换函数调用被认为是在完整独立的"integral selfcontained"计算中,“全部一起, 同时”交付结果。

Fortran 90提供了21个内部数组函数,其中有些函数(如SIZE)是用来查询数组特性的,其他函数或者是用来构造数组,或者用来从数组中提取信息。这些函数都是变换函数,如表5所列。

数组内部函数 说明
ALL 当所有元素值为真时, 其值为真
ANY 当有一个元素值为真时, 其值为真
ALLOCATED 检验数组是否被分配
COUNT 值为真的元素的个数
CSHIFT 循环移动数组的某一维
EOSHIFT "去尾(end-off)"移动数组的某一维
LBOUND 数组下界
MAXLOC 数组最大元素的位置
MAXVAL 数组最大元素的值
MERGE 在屏蔽码控制下, 合并两个数组
MINLOC 数组最小元素的位置
MINVAL 数组最小元素的值
PACK 在屏蔽码控制下, 把一个数组收集到一个向量中
PRODUCT 数组所有元素的乘积
SHAPE 数组形状
SIZE 整个数组的大小
SPREAD 通过增加一维来扩展一个数组
SUM 数组所有元素的和
TRANSPOSE 一个二维数组的矩阵转置
UBOUND 数组上界
UNPACK 在屏蔽码的控制下, 将一个向量分散到一个数组中

表 5

用户可以自己定义返回赋值数组的函数,所有这样的函数都是变换函数. 前一节所(21)举的例子中,函数F18是用定义返回数组质的函数的一个例子。在这种情况下,返回的数组形状,由参数(如数组参数的形状特性)动态决定;这也许是返回数组值的函数的最有用的形式。同样见4.5节中的例子。

注意,函数结果用一般定义语句来定义为赋值数组,就好象函数名是一个普通变量(和它在函数体内)一样. 虽然自动数组也许是用户定义返回赋值数组函数的最有用形式,但是其他形式的数组也是有效的: 如确定形状的数组,可分配数组,指针数组。这些函数也可在过程中定义和使用,就象函数名是另外一个变量一样. 主要的附加要求是数组值在从函数的一个执行返回之前必须被全定义。另一方面,返回赋值数组的函数的接口在使用该函数的地方必须显示说明,这样调用程序才能知道它在处理一个返回赋值数组的函数。

我们用一个返回赋值数组函数的简单例子来结束本节. 假设一个数组表达式中需要一个含有n个元素的一维数组的部分和——也就是说,需要的第k个值是sum(p(1:k))。一个返回赋值数组的函数对于交付所需集合的值是比较理想的(尽管在这个简单例子中,在表达式中使用数组构造器的效果和调用函数Partial-sums的效果几乎一样):

function Partial_sums(P)
  real P(:)                     !  假定形状的哑元数组
  real Partial_sums(size({))    !  要返回的部分和
  integer k
  Partial_sumns = (/(sum(P(1:k),k=1,size(P))/)
                     ! 在功能上, 它等价于:
                     ! do k=1, size(P)
                     !  Partial_sums(k) = sum(P(1:k))
                     ! end do
                     ! 但是, do 循环用来说明一个串型计算,
                     ! 而不是并行计算
end function Partial_sums

下面的更复杂的数据并行计算的例子也用用户定义的返回赋值数组的函数来交付结果。

4.5 例子:高斯消元法

为了说明数据并行的实际应用,这个例子提出两种形式的传统高斯消元算法来解决线性方程组。之所以选择这个特殊的例子,是因为我们普遍熟悉高斯消元,这样就可以将大部分注意力集中于数据并行技术,而尽可得少地将注意力分散到对这个问题的熟悉上。这个例子的一种形式叫做简单高斯消去法(Simple_Gauss), 沿着矩阵的主对角线(在下面的程序代码中称为Grid 或者 G)向下消主元;另一种形式叫做主元高斯法(Pivot-Gauss),实现更复杂但更健壮的高斯最大主元消去策略。

Simple-Gauss和Pivot-Gauss都有两种版本——一个标量顺序版本和一个数据并行版本。如例所示,数据并行版本能够编译和运行;顺序版本在行首用“!!”注明为注释。如果把这些注释符去掉,把以“!!!!”结尾的行注释掉,则这个顺序版本可以编译和运行。

顺序版本不包括数组操作(除了对G的初始化外),以所熟悉的标量do循环对矩阵操作为特征。对这些循环的数据并行版本的替换紧随其后,以便顺序版本和并行版本能方便地进行比较. 为了使比较更容易更有用,从这些例子中去掉了一些(不太多)自由性. 同样,为了使这种比较更容易些,对整个矩阵按主元进行了归约, 而不是仅对那些需要归约的列,这样每个算法的每个版本涉及的操作大约是实际需要的(标量元素)操作的二倍;这些算法可以很容易地进行调整, 来相应地限制操作的个数, 不过这要损失一些算法的清晰性。

对这些算法的一个分析表明, 简单高斯消去算法的顺序版本大约有4N3+7N2+9N+1个标量操作,并行版本在10个并行操作中大约有5N3+8N2+4N+1个标量操作. 主元高斯法的顺序版本大约有N4+7N3+4N2+5N+1个标量操作,并行版本在LnN+14个并行操作中大约有(lnN)N4+8N3+9N2+5N+1个标量操作。对相当大的N值,其结果总结在表6中. 表的最后一列(理想化地)假设一个数据并行操作与一个标量操作时间相等。

标量操作个数 并行操作个数 执行时间
简单高斯方法, 顺序 4N3 - 4N3
主元高斯法, 顺序 (N+7)N3 - (N+7)N3
简单高斯方法, 并行 5N3 10 10
主元高斯法, 并行 (NlnN+8)N3 lnN + 14 lnN + 14

表 6

因此,两种算法的并行版本包含的标量运算都比顺序版本多,但相比而言其并行操作却相当少。就标量操作而言,一个并行操作的有效耗费目前随系统不同而不同,但趋势看起来好象是逼近标量操作耗费。从这些方面来看,高斯消元法的数据并行版本的确更具吸引力。

最后,在简单高斯法和主元高斯法的基本归约操作中都使用了Fortran 90的一个内部函数SPREAD。 SPREAD函数复制(广播)一个标量到一个一维数组中,或将一个 n 维数组复制到一个n+1维数组中。这里使用的是标量到一维数组的形式,就是将标量操作G(i,j)=G(i,L)*G(L,j)转化成一个对数组G的整个数组操作。在这个表达式中,在对i和j的循环中L是“常量”,因此必须将它广播以填满数组,从而可以对整个数组操作。理解这一点是理解这个算法的数据并行版本的关键, 也是最困难的部分。SPREAD函数有三个参数:第一个是标量或被广播的数组,第二个是发生广播的维数(如果广播一个标量,则必须为1),第三个是复制的个数(在这些例中是N或N+1)

function Simple.Gauss(Grid)         ! 高斯消元 - 非最大主元
   real :: Grid(:,:)               ! 要归约的矩阵
   real ::  Simple_Gauss(size(Grid,1))  ! 返回结果向量
   real ::  G(size(Grid,1), size(Grid,2))  ! G 是一个局部工作向量
   logical :: Not_pivot_row(size(Grid,1), size(Grid,2 )) ! 主元行屏蔽
 if (size(Grid,2).ne.size(Grid,1)+1 stop \"bad Grid shape\"
 N = size(Grid,1)
 G = Grid     ! 在G 上操作, 而不是Grid上
 do L=1,N     ! G(L,L) 是下一个主元元素
  if (abs(G(L,L)).lt.1E-4) stop \"zero encounterd in pivot\"
!! G_pivot = G(L,L)    !!
!! do j=1,N+1      !! 主元行规格化
!!  (G(L,j) = G(L,j)/G_pivot    !!
!! enddo     !!
  G(L,1) = G(L,1)/G(L,L)     !!!!  数据并行版本
!! do i = 1,N   !!
!!   do j = 1, n+1    !!
!! if ( i.ne.L.and.j.ne.L )   then    !! 然后用该主元归约矩阵
!! G(i,j) = G(i,j) - G(i,L)*G(L,j)    !!
!! end if !!
!!  end do !!
!! end do !!
  Not_pivot_row = .true.; Not_pivot_row(L,:) = .false.
  where ( NOt_pivot_row) &      !!!!  数据并行版本
   G = G - G(:,spread(L,1,N+1))*G(spread(L,1,N),:)    !!!!  数据并行版本
  end do            ! 对所有主元重复
!! do i = 1, N         !! 最后, 提取结果向量形式
!!  Simple_Gauss(i) = G(i, N+1)    !!
!! end do           !!   G 的最后一列
  Simple_Gauss = G(:, N+1)    !!!! 数据并行版本
end function Simple_Gauss
function Pivot_Gauss(Grid)       ! 高斯消元, 最大主元
  real :: Grid(:,:)   ! 要归约的矩阵
  real :: Pivot_Gauss(size(Grid,1)) ! 返回结果向量
  real :: G(size(Grid,size(Grid,2))  ! G 是一个局部工作数组
  integer :: P(size(Grid,1),2)  ! P 是一个主元数组
  logical :: Not_pivot_row(size(Grid,1),size(Grid,2))
                                ! 只屏蔽当前主元行
  logical :: Not_pivot_row(size(Grid,1),size(Grid,1))
                                ! 屏蔽所有以前的主元行和列
  if (size(Grid,2).ne.size(Grid,1)+1) stop \"bad Grid shape\"
  N = size(Grid,1)
  G = Grid       ! Work on G, not Grid.
  do L = 1, N    !! L 是下一个主元数
!!  G_pivot = -1    !!
!!  do i = 1,N   !! 首先, 寻找下一个主元
!!    Inner_pivot_search :  & !!
!!    do j=1,N   !!
!!     do k=1,L-1   !!
!!       if(i.eq.P(k,1).or.j.eq.P(k,2)) cycle Inner_pivot_search
!!     end do  !! 跳过这个元素; 它在前一个主元行或列中
!!     if (abs(G(i,j)).gt.G_pivot) then   !!
!!      G_pivot = abs(G(i,j))   !!
!!      P(L,1) = i !!
!!      P(L,2) = j !!
!!     end if   !!
!!    end do   Inner_pivot_search   !!
!!   end do     !!
   Not_pivot_rows_or_cols = .true.     !!!! 数据并行版本
   Not_pivot_rows_cr_cols(P(1:L-1,1),:) = .false.  !!!! 数据并行版本
   Not_pivot_rows_or_cols(:,P(1:L-1,2)) = .false.  !!!! 数据并行版本
 P(L,:) = maxloc(abs(G(:,1:N)),mask=Not_pivot_rows_or_cols)
                                                   !!!! 数据并行版本
  if (abs(G(P(L,1),P(L,2))).lt.1E-4) stop \"ill-conditioned matrix\"
!!  G_pivot = G(P(L,1),P(L,2))   !!
!!  do j=1,N+1     !!  然后规格化主元行,
!!   G(P(L,1),j) = g(P(L,1),j)/G_pivot     !! 建立主元行标识
!!  end do  !!
  G(P(L,1),:) = G(P(L,1),:)/G(P(L,1),P(L,2))  !!!! 数据并行版本
!!  do i = 1,N    !! 然后用这个主元元素归约矩阵
!!   do j = 1,N +1  !! 
!!    if ( i.ne.P(L,1).and.j.ne.P(L,2)) then  !!
!!    G(i,j) = G(i,j)-G(i,P(L,2))*G(P(L,1),j)  !!
!!    end if  !!
!!   end do  !!
!!  end do  !!
  Not_pivot_rov = .trur.: Not_pivot_row(P(L,1),:) = .false. 
                                         !!!!  数据并行版本
  where ( Not_pivot_row ) &              !!!!  数据并行版本
  G = G-G(:,spread(P(L,2),1,N+1))*G(spread(P(L,1),1,N),:)
                                         !!!!  数据并行版本
end do   !  对所有主元重复
!! do i=1,N                    !! 最后, 从 G 的最后一列中恢复解向量
!! Pivot_Gauss(P(i,2)) = G(P(i,1),N+1)  !! 
!! end do  !!
Pivot_Gauss(P(:,2)) = G(P(:,1),N=1)    !!!! 数据并行版本
end function Pivot_Gauss

参考文献

[1] Fortran 90 Handbook, Adams, Brainerd, Martin, Smith, Wagener, McGraw-Hill, 1992

[2] Programming Language Fortran, ANSI standard X3.198-1992

[3] DRAFT - Process Parallelism Standard, ANSI committee X3H5


Copyright: NPACT BACKWARD