怎样提高调用数学函数的程序的性能

利用 IBM XL C/C++ 或 XL Fortran 汇编器的自动向量化功能

本文向您介绍了 IBM MASS 高性能数学库的内容,并展示了怎样通过 IBM® XL C/C++ 与 XL Fortran 汇编器的自动向量化功能去使用它们,而不用需要源程序方面的更改。在介绍自动向量化的概念以及相应的汇编器选项以后,本文会提供一个离散 Fourier 转变程序作为自动向量化的真实范例。计时的结果演示了通过 MASS 自动向量化的自动激活功能,从而对范例程序运行汇编器使得速度增加了 8.94 倍。

Robert F. Enenkel, 软件开发人员, IBM

Robert Enenkel 照片Robert Enenkel 工作于 IBM 多伦多实验室的 Optimizing Compiler Group,他之前是 IBM Centre for Advanced Studies (CAS) 的一名研究会员。Enenkel 拥有多伦多大学的学士、硕士和博士学位,其论文题目是不同等式的并行数字方法领域。他目前主要研究和开发与编译器和操作系统相关的数字计算,包括浮点算法,数学函数库,算法性能调优。Enenkel 博士得到了多个 IBM 发明成就奖和创造认可奖。更多的信息可以在他的 Web 主页 https://www-927.ibm.com/ibm/cas/toronto/people/members/renekel.shtml 上找到。



Daniel M. Zabawa, 软件开发人员, IBM

Daniel Zabawa 照片Daniel Zabawa 在 IBM 多伦多研究室的 Optimizing Compiler Group 工作。他拥有多伦多大学的计算机科学学士和硕士学位。Zabawa 目前研究和开发的主要方面在超标量体系架构的循环优化和调度、浮点算法库,以及浮点算数函数的性能调优。



2010 年 7 月 12 日

MASS 高性能库

MASS 指的是数学加速子系统(Mathematical Acceleration Subsystem)。它由数学函数组成,这些数学函数是为在各种 IBM 计算平台上优化性能所设定的。MASS 最初是由 IBM 公司在 1995 年启动的,并在随后的发展中继续得到改善,一直到现如今仍然在改进。

现在所有的 IBM® POWER™ 处理器都有相应的各种版本的 MASS,运行 AIX® 或者 Linux® 操作系统。还有其他版本的 IBM System BlueGene®/L 和 IBM System BlueGene®/P 超级电脑,以及 IBM Cell Broadband Engine™ (Cell/B.E.™)方案。库包含有元素函数的加速实施方案,例如 三角函数以及双曲线函数以及它们的倒数、乘方、对数、指数、错误函数以及其他函数。包含函数的完整列表可以在 IBM Mathematical Acceleration Subsystem 页面中找到。

有标量的库也有向量的库,而对于 Cell/B.E. 和 POWER7 来说,还有单个指示的多个数据(SIMD)库。注意精确性与例外情况的处理在 MASS 函数与系统库函数中可能是不一样的。对于目标硬件的其他汇编器(例如 gcc)的用户来说,MASS 库与 IBM XL C/C++ 还有 XL Fortran 汇编器封装到一起,并且可以通过 MASS Web 网站来获得。

可以通过 C、 C++ 或者 Fortran 源程序来访问库。IBM XL C/C++ 与 IBM XL Fortran 汇编器可以识别机会以使用 MASS 来加速源程序,并自动激活它而不用更改源程序。本文向您介绍了怎样实施一项技术帮助您的公司更好地使用这些强大的技术。

什么程序可以获益?

任何包含有对数学库函数(例如 exp、 log、sin、cos 等等)调用的 C、C++ 或者 Fortran 程序,潜在意义上都会从本文中所描述的技术中受益。


什么是自动化向量?

自动化向量技术是一种过程,在这个过程中 IBM XL C/C++ 或者 Fortran 汇编器会识别一个机会,去改善汇编过程中程序的性能,方法就是将对一次循环中一个标准库(C/C++ 库或者 Fortran 本质)的访问替换为对相应 MASS 向量函数的访问。因为 MASS 向量函数要比对一个标准库函数的重复访问快很多(倍数接近 30 倍),所以最后得到的性能改善效果将会是惊人的。

一个简单的例子就是为多个论断计算特定函数的循环,例如接下来的 Fortran 程序。

subroutine sub (y, x, n)
    real*8 y(*), x(*)
    integer n, i
    intrinsic exp
    do i=1,n
        y(i) = exp(x(i))  ! standard intrinsic
    end do
end subroutine

有了适当的汇编器选项,汇编器就会意识到机会去给程序加速,方法就是将对 exp() 的重复访问替换为相应的 MASS 向量函数 vexp(),结果会产生一个程序,好像最开始是这样写成的这样:

include 'massv.include'
subroutine sub (y, x, n)
    real*8 y(*), x(*)
    integer n
    call vexp (y, x, n)  ! MASS vector function
end subroutine

这只是一个简单的范例,演示了自动向量化背后的基本思想。XL 汇编器实际上能够识别更加复杂的机会,并在需要的条件下重新安排源程序中的指南,以创建自动向量化的机会。

在本文中的范例研究部分中,会检查一个更加复杂和实际的范例。


自动向量化的汇编器选项

您可以使用以下的几个选项来汇编程序:

  1. -qhot -qnostrict (for Fortran)
  2. -qhot -qnostrict –qignerrno (for C/C++)
  3. -qhot -O3
  4. -O4
  5. -O5

当您在使用这些选项集中的一个时,通过调用等价 MASS 向量函数(除了对以下函数的访问除外:vatan2、vsatan2、 vdnint、 vdint、 vcosisin、vscosisin、vqdrt、vsqdrt、vrqdrt、vsrqdrt、vpopcnt4、vpopcnt8、vexp2、vexp2m1、vsexp2、 vsexp2m1、vlog2、 vlog21p、 vslog2 和 vslog21p),汇编器会自动尝试对系统数学函数的访问向量化。如果汇编器不能对程序进行向量化,它会自动试着调用等价 MASS 标量函数。对于自动化的标量或者向量,汇编器会使用汇编器库 libxlopt.a 中包含的 MASS 函数的版本。您不需要向代码中的 MASS 函数添加任何特意的调用,或者链接 xlopt 库。

除了一系列的选项之外,当 -qipa 选项处于可用状态时,如果汇编器不能进行向量化,那么它会试着在决定调用它们之前去内联 MASS 标量函数。

如果您想要取消自动向量化的激活,那么您可以添加选项 –qhot=novector。


用例研究

接下来的部分是一个实际程序的范例 — 一个离散的 Fourier 转变(DFT) — 显示了在汇编不同汇编器选项时的改善结果。程序已经足够简单以方便演示,然后又足够的复杂以提供非琐细的优化机会。

两个程序的计时都是在附录 3 中给出的驱动器程序完成的,运行的环境是在 4.704 GHz 下运行的 POWER6 电脑。

附录 1 显示了 Fortran DFT 源程序。它包含了一个嵌套的循环,该循环会调用 exp()、cos() 以及 sin(),接下来是一个调用 sin() 和 sqrt() 的循环。程序会使用 -O3(它并不能进行自动向量化) 并使用 –O4 (它能使用自动向量化)。结果如图 1 所示 。

注意自动向量化带来的好处会随着问题规模的增加而增加,最终当问题的规模达到 2000 时加速的程度会达到 8.94x 。

图 1:DFT Fortran 性能与汇编选项 -O3 和 -O4 在各种规格问题上的比较
显示各种问题规模下元素的图

附录 2 显示了附录 1 中 Fortran DFT 程序的 C 版本(它包含了一个虚 consume() 路径,这样汇编器的内部程序化分析[IPA]就不能看到,计算的结果实际上在演示范例中并没有用得上,并因此可以改善整个的程序)。

程序将会使用 -O3(它并不会提供自动向量化) ,使用 -O4 (它提供自动向量化),使用 –O5 (它提供自动向量化并提供 IPA)。结果如图 2 所示。

正如在 Fortran 范例中演示的那样,自动向量化带来的好处随着问题规模的增加而增加,最后当 n=2000 的时候达到了。另外,IPA 在 -O5 处提供的活化能够提供一个额外的 1.22x 加速,因为它可以决定输入与输出没有别名(这就是说,它没有在内存中重叠),允许它去向量化进行极坐标的转变。-O5 在 –O3 的基础上加速的程度是 7.33x 。

图 2:DFT C 性能与汇编选项 -O3, -O4 与 -O5 对于不同规格问题的比较
显示不同规格问题的图

结论

本文向您提供了对 IBM MASS 库以及 IBM XL C/C++ 和 XL Fortran 汇编器的自动向量化功能的描述。另外,本文演示了对范例程序(离散 Fourier 转变)使用各种汇编器选项的操作,向您展示了通过使用 MASS 自动向量化的自动调用功能,使得与以前版本相比速度提高了 8.94 倍的效果。

这种演示想要通过一种程序来鼓励用户,这种程序会访问数学函数以验证可用的汇编器选项,并从 IBM XL C/C++ 或者 XL Fortran 汇编器的自动向量化加速中获益。


附录 1 – Fortran DFT 源程序

  subroutine dft (x, a, phi , n)
    real*8 x(n), a(n), phi(n)
    integer n

    ! Compute discrete Fourier transform of real inputs         
    ! x(i) and convert to polar form.

    real*8, parameter :: pi=3.1415926535897932384d0
    real*8 y_re(n), y_im(n), t, term_re, term_im
    intrinsic exp, cos, sin, sqrt, atan

    y_re(1:n) = 0.d0
    y_im(1:n) = 0.d0
    do k=1,n
      ! compute y(k), k-th DFT output
      do i=1,n
        ! compute i-th term of y(k):
        ! x(k)*exp(-2*pi*I*(k-1)*(i-1)/n)
        ! compute real and imaginary parts of i-th term
        ! using exp(I*t)=exp(t)*(cos(t)+I*sin(t))

        t = -2.d0*pi*(k-1)*(i-1)/n
        term_re = x(i) * cos(t) * exp(t)
        term_im = x(i) * sin(t) * exp(t)

        ! add term to sum
        y_re(k) = y_re(k) + term_re
        y_im(k) = y_im(k) + term_im
      end do
    end do

    ! transform y to polar coordinates
    do k=1,n
      ! compute amplitude of y(k)
      a(k) = sqrt (y_re(k)**2 + y_im(k)**2)
      ! compute phase of y(k)
      phi(k) = atan (y_im(k) / y_reim(k))
    end do

  end subroutine

  ! initialize input data
  subroutine init (a, n)
    real*8 a(n)
    integer n
    intrinsic sin,sqrt
    
    do j=1,n
      a(j)=sin(1.d0/sqrt(real(j,8)))
    end do
  end subroutine

附录 2 – C DFT 源程序

#include <math.h>
#define PI 3.1415926535897932384
void dft(double x[],double a[],double phi[],int *m)
{
  double y_re[NMAX], y_im[NMAX], t, s, term_re, term_im;
  int i,j,k,n=*m;
  for(i=0;i<n;++i) {
    y_re[i]=y_im[i]=0;
  }
  for(k=0;k<n;++k) 
  {
    // compute y(k), k-th DFT output
    for(i=0;i<n;++i)
    {
      // compute i-th term of y(k):
      // x(k)*exp(-2*pi*I*(k-1)*(i-1)/n)
      // compute real and imaginary parts of i-th term
      // using exp(I*t)=exp(t)*(cos(t)+I*sin(t))

      t=-2.*PI*k*i/(double)n;
      term_re=x[i]*exp(t)*cos(t);
      term_im=x[i]*exp(t)*sin(t);

      // add term to sum
      y_re[k]+=term_re;
      y_im[k]+=term_im;
    }
  }

  // transform y to polar coordinates
  for(k=0;k<n;++k)
  {
    // compute amplitude of y(k)
    a[k]=sqrt(y_re[k]*y_re[k]+y_im[k]*y_im[k]);
    // compute phase of y(k)
    phi[k]=atan2(y_im[k],y_re[k]);
  }
}

// initialize input data
void init(double a[],int *m)
{
  int j,n=*m;
  for(j=0;j<n;++j) 
  {
    a[j]=sin(1./sqrt((double)j+1.0));
  }
}

// Dummy function to use result, preventing compiler from
// optimizing away the computation.
void consume(double a[],double b[],double c[])
{
}

附录 3 – 驱动器程序

下面是一个用于计时 DFT 代码的 Fortran 驱动器程序。

  program main
    interface
      subroutine dft(x,a,phi,n)
      real*8 x(n),a(n),phi(n)
      integer n
      end subroutine
      subroutine init(a,n)
      real*8 a(n)
      integer n
      end subroutine
      subroutine consume(a,b,c)
      real*8 a(*),b(*),c(*)
      end subroutine
    end interface

    ! Parameters:
    !   nmax is the problem size.
    !   nrep is the number of repetitions of the
    !        problem. This should be chosen so that
    !        the elapsed time is long enough to give
    !        sufficient timing resolution.
    !   cyc is the clock frequency in Hz for the
    !       processor that the program is to be run on.
    !       (Can be found from AIX command pmcycles.)
    
    integer, parameter :: nmax=1000
    integer, parameter :: nrep=100
    real*8, parameter :: cyc=4704000000.d0

    real*8 x(nmax), a(nmax), phi(nmax)
    real*8 tx, ty, accum, del(4)
    intrinsic sin, sqrt
    real*8 rtc

    del(4)=0.d0
    acc = 0.d0
    do k=1,nrep
      tx=rtc()
      call init(x,nmax)
      call dft(x,a,phi,nmax)
      ty=rtc()
      call consume(x,a,phi)
      do j=1,nmax
        acc = acc + a(j) + phi(j)
      end do
      del(4) = del(4) + (ty-tx)
    end do
    del(1) = del(4)/real(nmax,8)
    del(2) = del(1)/real(nrep,8)
    del(3) = cyc*del(2)
    print *,'acc=',acc/real(nrep,8),' n=',nmax,
 &    ' r=',nrep, ' a=',del(1),' b=',del(2),
 &    ' c=',del(3),' w=',del(4)
  end program

声明

获得的精确性可能会随着使用处理器的模型及其配置而发生变化,同时依赖的因素还有所用汇编器的版本。因此,在不同的操作条件下您将获得不同的使用经验 。

参考资料

学习

获得产品和技术

  • IBM 试用版软件:使用从 developerWorks 中直接下载的软件,构建您的下一个开发项目。

讨论

条评论

developerWorks: 登录

标有星(*)号的字段是必填字段。


需要一个 IBM ID?
忘记 IBM ID?


忘记密码?
更改您的密码

单击提交则表示您同意developerWorks 的条款和条件。 查看条款和条件

 


在您首次登录 developerWorks 时,会为您创建一份个人概要。您的个人概要中的信息(您的姓名、国家/地区,以及公司名称)是公开显示的,而且会随着您发布的任何内容一起显示,除非您选择隐藏您的公司名称。您可以随时更新您的 IBM 帐户。

所有提交的信息确保安全。

选择您的昵称



当您初次登录到 developerWorks 时,将会为您创建一份概要信息,您需要指定一个昵称。您的昵称将和您在 developerWorks 发布的内容显示在一起。

昵称长度在 3 至 31 个字符之间。 您的昵称在 developerWorks 社区中必须是唯一的,并且出于隐私保护的原因,不能是您的电子邮件地址。

标有星(*)号的字段是必填字段。

(昵称长度在 3 至 31 个字符之间)

单击提交则表示您同意developerWorks 的条款和条件。 查看条款和条件.

 


所有提交的信息确保安全。


static.content.url=http://www.ibm.com/developerworks/js/artrating/
SITE_ID=10
Zone=Rational, Multicore acceleration
ArticleID=500109
ArticleTitle=怎样提高调用数学函数的程序的性能
publish-date=07122010