MASS 高性能庫(kù)
MASS 指的是數(shù)學(xué)加速子系統(tǒng)（Mathematical Acceleration Subsystem）。它由數(shù)學(xué)函數(shù)組成，這些數(shù)學(xué)函數(shù)是為在各種 IBM 計(jì)算平臺(tái)上優(yōu)化性能所設(shè)定的。MASS 最初是由 IBM 公司在 1995 年啟動(dòng)的，并在隨后的發(fā)展中繼續(xù)得到改善，一直到現(xiàn)如今仍然在改進(jìn)。
現(xiàn)在所有的 IBM® POWER™ 處理器都有相應(yīng)的各種版本的 MASS，運(yùn)行 AIX® 或者 Linux® 操作系統(tǒng)。還有其他版本的 IBM System BlueGene®/L 和 IBM System BlueGene®/P 超級(jí)電腦，以及 IBM Cell Broadband Engine™ （Cell/B.E.™）方案。庫(kù)包含有元素函數(shù)的加速實(shí)施方案，例如三角函數(shù)以及雙曲線函數(shù)以及它們的倒數(shù)、乘方、對(duì)數(shù)、指數(shù)、錯(cuò)誤函數(shù)以及其他函數(shù)。包含函數(shù)的完整列表可以在 IBM Mathematical Acceleration Subsystem 頁(yè)面中找到。
有標(biāo)量的庫(kù)也有向量的庫(kù)，而對(duì)于 Cell/B.E. 和 POWER7 來(lái)說(shuō)，還有單個(gè)指示的多個(gè)數(shù)據(jù)（SIMD）庫(kù)。注意精確性與例外情況的處理在 MASS 函數(shù)與系統(tǒng)庫(kù)函數(shù)中可能是不一樣的。對(duì)于目標(biāo)硬件的其他匯編器（例如 gcc）的用戶來(lái)說(shuō)，MASS 庫(kù)與 IBM XL C/C++ 還有 XL Fortran 匯編器封裝到一起，并且可以通過 MASS Web 網(wǎng)站來(lái)獲得。
可以通過 C、 C++ 或者 Fortran 源程序來(lái)訪問庫(kù)。IBM XL C/C++ 與 IBM XL Fortran 匯編器可以識(shí)別機(jī)會(huì)以使用 MASS 來(lái)加速源程序，并自動(dòng)激活它而不用更改源程序。本文向您介紹了怎樣實(shí)施一項(xiàng)技術(shù)幫助您的公司更好地使用這些強(qiáng)大的技術(shù)。
什么程序可以獲益？
任何包含有對(duì)數(shù)學(xué)庫(kù)函數(shù)（例如 exp、 log、sin、cos 等等）調(diào)用的 C、C++ 或者 Fortran 程序，潛在意義上都會(huì)從本文中所描述的技術(shù)中受益。
什么是自動(dòng)化向量？
自動(dòng)化向量技術(shù)是一種過程，在這個(gè)過程中 IBM XL C/C++ 或者 Fortran 匯編器會(huì)識(shí)別一個(gè)機(jī)會(huì)，去改善匯編過程中程序的性能，方法就是將對(duì)一次循環(huán)中一個(gè)標(biāo)準(zhǔn)庫(kù)（C/C++ 庫(kù)或者 Fortran 本質(zhì)）的訪問替換為對(duì)相應(yīng) MASS 向量函數(shù)的訪問。因?yàn)?MASS 向量函數(shù)要比對(duì)一個(gè)標(biāo)準(zhǔn)庫(kù)函數(shù)的重復(fù)訪問快很多（倍數(shù)接近 30 倍），所以最后得到的性能改善效果將會(huì)是驚人的。
一個(gè)簡(jiǎn)單的例子就是為多個(gè)論斷計(jì)算特定函數(shù)的循環(huán)，例如接下來(lái)的 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
有了適當(dāng)?shù)膮R編器選項(xiàng)，匯編器就會(huì)意識(shí)到機(jī)會(huì)去給程序加速，方法就是將對(duì) exp() 的重復(fù)訪問替換為相應(yīng)的 MASS 向量函數(shù) vexp()，結(jié)果會(huì)產(chǎn)生一個(gè)程序，好像最開始是這樣寫成的這樣：
include 'massv.include'
subroutine sub (y, x, n)
real*8 y(*), x(*)
integer n
call vexp (y, x, n) ! MASS vector function
end subroutine
這只是一個(gè)簡(jiǎn)單的范例，演示了自動(dòng)向量化背后的基本思想。XL 匯編器實(shí)際上能夠識(shí)別更加復(fù)雜的機(jī)會(huì)，并在需要的條件下重新安排源程序中的指南，以創(chuàng)建自動(dòng)向量化的機(jī)會(huì)。
在本文中的范例研究部分中，會(huì)檢查一個(gè)更加復(fù)雜和實(shí)際的范例。
自動(dòng)向量化的匯編器選項(xiàng)
您可以使用以下的幾個(gè)選項(xiàng)來(lái)匯編程序：
-qhot -qnostrict (for Fortran)
-qhot -qnostrict –qignerrno (for C/C++)
-qhot -O3
-O4
-O5
當(dāng)您在使用這些選項(xiàng)集中的一個(gè)時(shí)，通過調(diào)用等價(jià) MASS 向量函數(shù)（除了對(duì)以下函數(shù)的訪問除外：vatan2、vsatan2、 vdnint、 vdint、 vcosisin、vscosisin、vqdrt、vsqdrt、vrqdrt、vsrqdrt、vpopcnt4、vpopcnt8、vexp2、 vexp2m1、vsexp2、 vsexp2m1、vlog2、 vlog21p、 vslog2 和 vslog21p），匯編器會(huì)自動(dòng)嘗試對(duì)系統(tǒng)數(shù)學(xué)函數(shù)的訪問向量化。如果匯編器不能對(duì)程序進(jìn)行向量化，它會(huì)自動(dòng)試著調(diào)用等價(jià) MASS 標(biāo)量函數(shù)。對(duì)于自動(dòng)化的標(biāo)量或者向量，匯編器會(huì)使用匯編器庫(kù) libxlopt.a 中包含的 MASS 函數(shù)的版本。您不需要向代碼中的 MASS 函數(shù)添加任何特意的調(diào)用，或者鏈接 xlopt 庫(kù)。
除了一系列的選項(xiàng)之外，當(dāng) -qipa 選項(xiàng)處于可用狀態(tài)時(shí)，如果匯編器不能進(jìn)行向量化，那么它會(huì)試著在決定調(diào)用它們之前去內(nèi)聯(lián) MASS 標(biāo)量函數(shù)。
如果您想要取消自動(dòng)向量化的激活，那么您可以添加選項(xiàng) –qhot=novector。
用例研究
接下來(lái)的部分是一個(gè)實(shí)際程序的范例 — 一個(gè)離散的 Fourier 轉(zhuǎn)變（DFT） — 顯示了在匯編不同匯編器選項(xiàng)時(shí)的改善結(jié)果。程序已經(jīng)足夠簡(jiǎn)單以方便演示，然后又足夠的復(fù)雜以提供非瑣細(xì)的優(yōu)化機(jī)會(huì)。
兩個(gè)程序的計(jì)時(shí)都是在附錄 3 中給出的驅(qū)動(dòng)器程序完成的，運(yùn)行的環(huán)境是在 4.704 GHz 下運(yùn)行的 POWER6 電腦。
附錄 1 顯示了 Fortran DFT 源程序。它包含了一個(gè)嵌套的循環(huán)，該循環(huán)會(huì)調(diào)用 exp()、cos() 以及 sin()，接下來(lái)是一個(gè)調(diào)用 sin() 和 sqrt() 的循環(huán)。程序會(huì)使用 -O3（它并不能進(jìn)行自動(dòng)向量化） 并使用 –O4 （它能使用自動(dòng)向量化）。
注意自動(dòng)向量化帶來(lái)的好處會(huì)隨著問題規(guī)模的增加而增加，最終當(dāng)問題的規(guī)模達(dá)到 2000 時(shí)加速的程度會(huì)達(dá)到 8.94x 。
附錄 2 顯示了附錄 1 中 Fortran DFT 程序的 C 版本（它包含了一個(gè)虛 consume() 路徑，這樣匯編器的內(nèi)部程序化分析[IPA]就不能看到，計(jì)算的結(jié)果實(shí)際上在演示范例中并沒有用得上，并因此可以改善整個(gè)的程序）。
程序?qū)?huì)使用 -O3（它并不會(huì)提供自動(dòng)向量化） ，使用 -O4 （它提供自動(dòng)向量化），使用 –O5 （它提供自動(dòng)向量化并提供 IPA）。
正如在 Fortran 范例中演示的那樣，自動(dòng)向量化帶來(lái)的好處隨著問題規(guī)模的增加而增加，最后當(dāng) n=2000 的時(shí)候達(dá)到了。另外，IPA 在 -O5 處提供的活化能夠提供一個(gè)額外的 1.22x 加速，因?yàn)樗梢詻Q定輸入與輸出沒有別名（這就是說(shuō)，它沒有在內(nèi)存中重疊），允許它去向量化進(jìn)行極坐標(biāo)的轉(zhuǎn)變。-O5 在 –O3 的基礎(chǔ)上加速的程度是 7.33x 。
結(jié)論
本文向您提供了對(duì) IBM MASS 庫(kù)以及 IBM XL C/C++ 和 XL Fortran 匯編器的自動(dòng)向量化功能的描述。另外，本文演示了對(duì)范例程序（離散 Fourier 轉(zhuǎn)變）使用各種匯編器選項(xiàng)的操作，向您展示了通過使用 MASS 自動(dòng)向量化的自動(dòng)調(diào)用功能，使得與以前版本相比速度提高了 8.94 倍的效果。
這種演示想要通過一種程序來(lái)鼓勵(lì)用戶，這種程序會(huì)訪問數(shù)學(xué)函數(shù)以驗(yàn)證可用的匯編器選項(xiàng)，并從 IBM XL C/C++ 或者 XL Fortran 匯編器的自動(dòng)向量化加速中獲益。
附錄 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
#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 y_re[i]=y_im[i]=0;
}
for(k=0;k {
// compute y(k), k-th DFT output
for(i=0;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 {
// 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 {
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[])
{
}