まず始めに、オリジナルの FFTE パッケージをそのまま PGI コンパイル・オプションを指定して実行モジュールを作成してみます。実行するプログラムは、アーカイブを展開した中にある tests ディレクトリ下の speed3d.f（三次元FFT）です。この初期性能を測定します。なお、使用したシステムは、以下の表にある、AMD のプロセッサとインテル(R) のデュアル・プロセッサを有する二つのシステムです。

Makefile定義するコンパイル・オプション並びに パッケージで使用されているヘッダーファイル param.h は、以下の通りです。なお、64ビットアドレッシングに対応するため、コンパイル・オプションには -mcmodel=medium を付加しています。32bit システムの場合は必要ありません。また、param.h の中のパラメータ NP は、プロセッサ特性で NP の値が変わりますが、プロセッサの如何に拘らず、4 あるいは 8 を指定して予め性能を評価することをお勧めします。

【シリアル計算時のオプション】
  F77    = pgf90
  FFLAGS = -fastsse -Minfo -Mvect=nosizelimit -mcmodel=medium -Mprefetch=distance:8,nta -I..
  LDFLAGS= 

【OpenMP 並列計算時のオプション】
  F77    = pgf90
  FFLAGS = -fastsse -Minfo -Mvect=nosizelimit -mcmodel=medium -Mprefetch=distance:8,nta -mp -I..
  LDFLAGS=

【param.hファイル】
  AMD Athlon64x2 プロセッサの場合　        ： PARAMETER (NBLK=16), PARAMETER (NP=4)
  Intel(R) Pentium(R) D プロセッサの場合 　： PARAMETER (NBLK=16), PARAMETER (NP=8)
  Cache Size of L2 cache
      PARAMETER (L2SIZE=1048576)　キャッシュサイズの指定

【speed3d.fの時間関数変更:オリジナルは ETIME 関数】
      REAL*8 second
     ---------------
　    TIME1=second()
      DO I=1,LOOP
      CALL ZFFT3D(A,NX,NY,NZ,-1)
      END DO
      TIME2=second()
      TIME0=TIME2-TIME1

上記の表の中で問題となる点は、シリアル実行用にコンパイルした実行モジュールの性能とOpenMP 並列モジュールの1スレッドの性能差が大きいことです。これは、一概には言えませんが何らかのコンパイル最適化等の問題があることをほのめかす事実として、調査を行ってみました。

表 1 の性能は、このままでも良い性能であることは確かですが、上記で述べた問題は、何らかの最適化処理によるものと考え、プログラム内の精査を行いました。結論から述べると、この問題は PGI コンパイラは、より高度な最適化処理を施してしまったため、オリジナルソースレベルで考慮されている手動での最適化技術を無駄なものにしてしまっていることでした。本来は、コンパイラによる最適化処理を行った方が良いのが一般的な考え方ですが、この FFTE パッケージの著者である高橋氏は、ソースレベルでキャッシュ最適化（キャッシュタイリング、ストリップマイニング）を行っているため、このソースレベルにコンパイラがさらに最適化を施した場合、著作者が考えているソースレベル最適化のイメージを壊すことになってしまいます。この FFTE パッケージは、このようなコーディングレベルでのキャッシュ活用最適化を施しているため、性能の良い、さらにポータビリティの良い FFT パッケージになっています。今回の性能の問題は、結果的にこのような箇所に PGI コンパイラがより高度な最適化を施してしまったために生じたものであり、本来のコンパイラ機能とは全く逆の作用を及ぼしてしまったと言うことになります。言ってみれば、コンパイラが余計なお世話をしたことになります。

以下に、zfft3d.f のソース内の該当する最適化部分をリストします。また、リストの下部に「コンパイル最適化メッセージ」を表示しています。ここでは、行番号 85、95 の周辺部の最適化により、コンパイラがループ構造を再構成した旨のメッセージが現れています。これは、本来、著作者が意図した最適化イメージではありません。従って、コンパイラ・ディレクティブを挿入して、以下の該当部分のループの最適化を抑止します。

C     (C) COPYRIGHT SOFTWARE, 2000-2004, ALL RIGHTS RESERVED
C                BY
C         DAISUKE TAKAHASHI

zfft3d0:　オリジナルソース（zfft3d.f)

(   80) !$OMP DO
(   81)       DO 80 J=1,NY
(   82)         DO 70 II=1,NX,NBLK
(   83)           DO 30 KK=1,NZ,NBLK
(   84)             DO 20 I=II,MIN0(II+NBLK-1,NX)
(   85)               DO 10 K=KK,MIN0(KK+NBLK-1,NZ)
(   86)                 BZ(K,I-II+1)=A(I,J,K)
(   87)    10         CONTINUE
(   88)    20       CONTINUE
(   89)    30     CONTINUE
(   90)           DO 40 I=II,MIN0(II+NBLK-1,NX)
(   91)             CALL FFT235(BZ(1,I-II+1),C,WZ,NZ,LNZ)
(   92)    40     CONTINUE
(   93)           DO 60 K=1,NZ
(   94)             DO 50 I=II,MIN0(II+NBLK-1,NX)
(   95)               A(I,J,K)=BZ(K,I-II+1)
(   96)    50       CONTINUE
(   97)    60     CONTINUE
(   98)    70   CONTINUE
(   99)    80 CONTINUE

zfft3d0: コンパイルメッセージ
　　　　　　　　　
    84, Distributing loop; 2 new loops                 （この最適化は本来行って欲しくない）
        Interchange produces reordered loop nest: 85, 84
        Loop unrolled 4 times
    93, Interchange produces reordered loop nest: 94, 93
        Loop unrolled 4 times

● PGI コンパイラ・ディレクティブを挿入してコンパイラによる高度最適化を抑止する

コンパイラによる（ベクトル）最適化を抑止するためのディレクティブとして、cpgi$ novector が用意されています。このディレクティブを以下に示すように該当する部分に挿入します。PGI では、その他各種のディレクティブが用意されておりますが、この詳細は、 PGI USER'S GUIDE をご参照ください。なお、 cpgi$l の　l　(エル) は、次の行の「ループ(loop)」のみに適用すると言う意味となります。このディレクティブを挿入してコンパイルすると、以下のように、コンパイルメッセージは、「86, Loop unrolled 4 times」のみとなります。

zfft3d0:　オリジナルソース（zfft3d.f) にベクトル最適化抑止ディレクティブを挿入する

(   80) !$OMP DO
(   81)       DO 80 J=1,NY
(   82)         DO 70 II=1,NX,NBLK
(   83)           DO 30 KK=1,NZ,NBLK
(   84) cpgi$l novector
(   85)             DO 20 I=II,MIN0(II+NBLK-1,NX)
(   86)               DO 10 K=KK,MIN0(KK+NBLK-1,NZ)
(   87)                 BZ(K,I-II+1)=A(I,J,K)
(   88)    10         CONTINUE
(   89)    20       CONTINUE
(   90)    30     CONTINUE
(   91)           DO 40 I=II,MIN0(II+NBLK-1,NX)
(   92)             CALL FFT235(BZ(1,I-II+1),C,WZ,NZ,LNZ)
(   93)    40     CONTINUE
(   94) cpgi$l novector
(   95)           DO 60 K=1,NZ
(   96)             DO 50 I=II,MIN0(II+NBLK-1,NX)
(   97)               A(I,J,K)=BZ(K,I-II+1)
(   98)    50       CONTINUE
(   99)    60     CONTINUE
(  100)    70   CONTINUE
(  101)    80 CONTINUE
(  102) !$OMP DO

zfft3d0:  コンパイルメッセージ

    86, Loop unrolled 4 times
    96, Loop unrolled 4 times

zfft3d.f の中には、さらに同じようにディレクティブの挿入を行う箇所があります。

(  103)       DO 170 K=1,NZ
(  104)         DO 150 II=1,NX,NBLK
(  105)           DO 110 JJ=1,NY,NBLK
(  106) cpgi$l novector
(  107)             DO 100 I=II,MIN0(II+NBLK-1,NX)
(  108)               DO 90 J=JJ,MIN0(JJ+NBLK-1,NY)
(  109)                 BY(J,I-II+1)=A(I,J,K)
(  110)    90         CONTINUE
(  111)   100       CONTINUE
(  112)   110     CONTINUE
(  113)           DO 120 I=II,MIN0(II+NBLK-1,NX)
(  114)             CALL FFT235(BY(1,I-II+1),C,WY,NY,LNY)
(  115)   120     CONTINUE
(  116) cpgi$l novector
(  117)           DO 140 J=1,NY
(  118)             DO 130 I=II,MIN0(II+NBLK-1,NX)
(  119)               A(I,J,K)=BY(J,I-II+1)
(  120)   130       CONTINUE
(  121)   140     CONTINUE
(  122)   150   CONTINUE
(  123)         DO 160 J=1,NY
(  124)           CALL FFT235(A(1,J,K),C,WX,NX,LNX)
(  125)   160   CONTINUE
(  126)   170 CONTINUE

● zfft1d.f 並びに zfft2d.f に関しても同じようにディレクティブを挿入

【zfft1d.fの挿入箇所】

(  105) !$OMP DO PRIVATE(IJ,IJ0,IR,J,TEMP)
(  106)       DO 110 II=1,N1,NBLK
(  107)         DO 30 JJ=1,N2,NBLK
(  108) cpgi$l novector
(  109)           DO 20 I=II,MIN0(II+NBLK-1,N1)
(  110)             DO 10 J=JJ,MIN0(JJ+NBLK-1,N2)
(  111)               C(J,I-II+1)=A1(I,J)
(  112)    10       CONTINUE
(  113)    20     CONTINUE
(  114)    30   CONTINUE
（中略）
(  146) !$OMP DO
(  147)       DO 150 JJ=1,N2,NBLK
(  148)         DO 120 J=JJ,MIN0(JJ+NBLK-1,N2)
(  149)           CALL FFT235(B(1,J),C,W1,N1,IP1)
(  150)   120   CONTINUE
(  151) cpgi$l novector
(  152)         DO 140 I=1,N1
(  153)           DO 130 J=JJ,MIN0(JJ+NBLK-1,N2)
(  154)             A2(J,I)=B(I,J)
(  155)   130     CONTINUE
(  156)   140   CONTINUE
(  157)   150 CONTINUE

【zfft2d.fの挿入箇所】

(   77) !$OMP DO
(   78)       DO 70 II=1,NX,NBLK
(   79)         DO 30 JJ=1,NY,NBLK
(   80) cpgi$l novector
(   81)           DO 20 I=II,MIN0(II+NBLK-1,NX)
(   82)             DO 10 J=JJ,MIN0(JJ+NBLK-1,NY)
(   83)               B(J,I-II+1)=A(I,J)
(   84)    10       CONTINUE
(   85)    20     CONTINUE
(   86)    30   CONTINUE
(   87)         DO 40 I=II,MIN0(II+NBLK-1,NX)
(   88)           CALL FFT235(B(1,I-II+1),C,WY,NY,LNY)
(   89)    40   CONTINUE
(   90) cpgi$l novector
(   91)         DO 60 J=1,NY
(   92)           DO 50 I=II,MIN0(II+NBLK-1,NX)
(   93)             A(I,J)=B(J,I-II+1)
(   94)    50     CONTINUE
(   95)    60   CONTINUE
(   96)    70 CONTINUE
(   97) !$OMP DO

FFTE パッケージをコンパイルする際の性能最適なコンパイル・オプションを以下に示します。
以下のオプションの中で、-Mvect=nosizelimit は必ず付加してください。これは、ループ内のステートメントが多い箇所があり、PGI のデフォルトのベクトル最適化の対象とする閾値を超えるものがあり、このオプションを指定することにより、全てのループを最適化の対象として扱います。このオプションを付けるだけでも性能が大きく向上します。また、3次元FFT で、データサイズが大きくなり 64ビットアドレッシングが必要な場合は、必ず、-mcmodel=medium オプションを付加してください。さらに、OpenMP による並列化モジュールの作成を行うためには、 -mp オプションを付加してください。（PGI 6.2より、-Mvect=nosizelimit がデフォルトとなりましたので、PGI 6.2 以降では指定する必要がありません）

【FFTE の最適なオプション】

FFLAGS = -fastsse -Minfo -Mvect=nosizelimit -Mprefetch=distance:8,nta -Minline=size:50

(OpenMP並列の場合)

FFLAGS = -fastsse -Minfo -Mvect=nosizelimit -Mprefetch=distance:8,nta -Minline=size:50 -mp

上記で述べた最適化を施した後、 FFTE の 3 次元 FFT （256**3) の性能は、以下のようにオリジナルをそのままコンパイルして得た性能（表1）に較べて、OpenMP の性能が大きく向上しました。特に、2スレッド並列性能は、オリジナルのものに較べて最大 1.57 倍向上したことが分かります。

以下の表では、三次元、二次元 FFT のデータサイズを変更して性能測定をした結果を示しました。プロセッサは、インテル(R)とAMDのプロセッサの二つのケースを示しました。これらのプロセッサ（L2 chache size 1MB)においては、128データサイズの場合が最適高速な性能を得られることが分かります。3次元データのFFTでは、プロセッサのクロック値がその性能に反映される傾向があります。一方、二次元 FFT では、、クロックの低い Athlon64x2 の性能は、Pentium(R) Dの性能とほぼ同じと言う結果が得られています。このように、扱うデータサイズによっても絶対性能が異なるため、一概にプロセッサの種類、クロック値で、性能が速い・遅いを議論できないことがご理解いただけることでしょう。
以下の結果は、一つのメモリ帯域を共有する「デュアルコア」プロセッサ上でのスレッド並列性能ですので、以下のような並列効率が得られるソースベースのライブラリは驚きに値するものです。従って、FFT を使用しているプログラムにおいては、OpenMP対応の FFTE パッケージを使用することで、大きな性能向上を得ることができます。