割り算は遅い

出典: トータル・ディスクロージャ・サイト(事実をありのままに)

割り算は遅い。現行の汎用CPUでは普遍的に成り立つこの事実は、HPC向けに限らず、どんなコンピュータプログラムを書く際にも覚えておいたほうが良い。今回の記事は基本的にx86をターゲットにしているが、CUDAをはじめとするGPGPUでも原則は同じである。

目次

どれだけ遅いのか

Intelが公開しているIntel® 64 and IA-32 Architectures Optimization Reference Manualという資料のAppendix C Instruction Latency and Throughputという章に、x86の各命令に必要なクロック数のデータが載っている。Latencyとは命令の実行の開始から完了までにかかるクロック数であるが、他の命令で計算結果を利用しない限り、命令の実行を続けながら後続の命令をどんどん処理することができる命令も多い。CPUのパイプラインを占有するために後続の命令をブロックしてしまうクロック数はThroughputで示されている。

現行のXeon E5シリーズは表の06_2AH,06_2DHの欄である(具体的には、Linuxの/proc/cpuinfoファイルに書いてあるcpu familymodelの数字を16進数に置き換えれば良い)が、C-6ページに書いてあるAVXのベクトル倍精度除算命令VDIVPDは、Latencyが45、Throughputが44とある。これは、命令の開始から完了までに45クロックかかり、そのうち最初の44クロックはパイプラインを占有するのでパイプラインに後続の命令を投入できず、残りの1クロックは、この命令の結果を使わない命令であれば同じパイプラインを使う命令でも同時実行できる、ということを示している。

というのが表の見方であるが、周りに書いてある数字が1や2ばかりの中、このように45といったやたら大きな数字が目を引くのがまさに除算命令なのである。ちなみにVADDPD(AVXベクトル倍精度加算命令)やVSUBPD(AVXベクトル倍精度減算命令)はLatency 3、Throughput 1VMULPD(AVXベクトル倍精度乗算命令)でさえLatency 5、Throughput 1である。加減乗除だけを見比べても、除算は圧倒的に遅いのである。

逆数命令

このように除算命令は遅いのだが、実はVRCPPSという逆数の近似値を求めるベクトル命令が在り、これはLatency 7、Throughput 1と、VDIVPDと比べると圧倒的な速さで動く。ただしあくまでも近似値であり、Intelが公開している命令の仕様には比の誤差が最大で1.5 × 2-12ほど出るとある。それでもここから数値計算で精度を改善してもお釣りの来る速さ(むしろ、1ステップごとに精度が2倍となるNewton法を使ったときに、1ステップで単精度の精度23ビットにギリギリとなるように仕様を決めたのだろう)なのか、Intel Compiler(12.1でも13.0でも)で-no-prec-div (Linux, Mac OS X)や/Qprec-div- (Windows)オプションを指定すると、単精度除算を「VRCPPSで除数の逆数を求め、Newton法を1ステップ実行して23ビットの精度を確保した後に、この逆数と被序数との乗算を行う」という方法で計算するようになる。この命令は精度が低いために単精度命令しか存在しないためか、倍精度計算だとこのオプションでもVRCPPSを使ってくれないが、一度単精度に変換して逆数を求め、倍精度に変換し直した後に手でNewton法を実装することで、変換のオーバーヘッドにも関わらず単なるVDIVPDを使ったプログラムより高速化することが多い。以下にIntel Fortranを使うことを前提とした例を示す。

逆数命令とNewton法を使う例

      PROGRAM RCP
      IMPLICIT NONE
      INTEGER I,N,ITER,NITER
      PARAMETER(N=1024,NITER=1000000)
      INTEGER ISTART,IEND,IRATE,IMAX
      REAL*8 OA,OB,A,B
      REAL*8 BI,X
      COMMON/ORIG/OA(N),OB(N)
      COMMON/DATA/A(N),B(N)
C
C     初期化
      DO 10 I=1,N
         CALL RANDOM_NUMBER(OA(I))
         CALL RANDOM_NUMBER(OB(I))
 10   CONTINUE
      DO 20 I=1,N
         A(I)=OA(I)
         B(I)=OB(I)
 20   CONTINUE
C
C     単純な除算、おそらくVDIVPDにコンパイルされる
      CALL SYSTEM_CLOCK(ISTART,IRATE,IMAX)
      DO 30 ITER=1,NITER
         DO 30 I=1,N
            B(I)=A(I)/B(I)
 30   CONTINUE
      CALL SYSTEM_CLOCK(IEND,IRATE,IMAX)
      WRITE(*,*)DBLE(IEND-ISTART)/DBLE(IRATE),SUM(B)
C
C     再初期化
      DO 40 I=1,N
         A(I)=OA(I)
         B(I)=OB(I)
 40   CONTINUE
C
C     単精度除算 + Newton法、単精度除算はおそらくVRCPPSにコンパイルされる
      CALL SYSTEM_CLOCK(ISTART,IRATE,IMAX)
      DO 50 ITER=1,NITER
         DO 50 I=1,N
            BI=B(I)
            X=DBLE(1.0/REAL(BI))
            X=X*(2.0D0-BI*X)
            X=X*(2.0D0-BI*X)
            B(I)=A(I)*X
 50   CONTINUE
      CALL SYSTEM_CLOCK(IEND,IRATE,IMAX)
      WRITE(*,*)DBLE(IEND-ISTART)/DBLE(IRATE),SUM(B)
C
C     ちなみに加算だとここまで速い
      CALL SYSTEM_CLOCK(ISTART,IRATE,IMAX)
      DO 60 ITER=1,NITER
         DO 60 I=1,N
            B(I)=A(I)+B(I)
 60   CONTINUE
      CALL SYSTEM_CLOCK(IEND,IRATE,IMAX)
      WRITE(*,*)DBLE(IEND-ISTART)/DBLE(IRATE),SUM(B)
      STOP
      END

X=DBLE(1.0/REAL(BI))と、倍精度変数BIを単精度の値に変換して単精度定数1.0に対して除算を行い、その結果を倍精度に変換して倍精度変数Xに代入する。ここが最初のステップであるが、-xAVX -no-prec-divないし/QxAVX /prec-div-のオプションが在ればVRCPPSと1ステップのNewton法を使った計算としてコンパイルされる。あとは、Newton法を使ってX=X*(2.0D0-BI*X)の漸化式で精度を上げていく。Xの初期値が23ビット精度であるところから、理論上は1ステップで46ビット精度、2ステップで92ビット精度となり、52ビット精度の倍精度としては十分な精度となる。最後にこうして求まった逆数Xを乗算して、一連の計算が終了する。

例の実行結果

上記の例を、以下の環境で実行してみる。

CPU: Core i7 2600K 3.4 GHz
OS: Ubuntu 11.10
Compiler: Intel Fortran Composer XE 2013 Update 1 (13.0.1.117)
   3.21520000000000        102400.000000000     
  0.992200000000000        102400.000000000     
  0.159100000000000        1024102400.00000     

一番上が単純に割り算した場合で、2番目がIntel Fortranの最適化でVRCPPSを使うように書き換え、Newton法で精度を高めて計算を行ったプログラムである。3倍以上の高速化は特筆すべきポイントであろう。ただし、それでも3番目の加算に比べれば6倍遅いので、できればプログラム中に出てくる除算の回数を極力減らすようなプログラミングを心掛けたい。たとえば、同じ数による割り算が複数箇所に登場する場合は、逆数の計算1回とそれを使った複数回の乗算に置き換えたり、配列から値を読んでその値で割り算を行うような場合には、配列の値のほうを逆数にしておく、などの方法を採ることで、除算の回数を減らすことができる。

平方根についての余談

平方根の計算にも実は似た性質がある。

  • 平方根を直接計算するベクトル命令を持つが、非常に遅い(VSQRTPD: Latency 45, Throughput 44)。
  • 平方根の逆数の近似値を求める高速なベクトル単精度命令(VRSQRTPS: Latency 7, Throughput 1)を持つ。
  • 平方根の逆数Y=1.0/SQRT(X)のNewton法による漸化式はY=Y*(1.5-0.5*X*Y*Y)であり、割り算を含まないので高速に計算できる。
  • 平方根SQRT(X)=X*(1.0/SQRT(X))であるので、平方根の逆数が求まれば乗算だけで計算できる。
  • Intel Compilerは単精度計算で上記のことを自動的に行ってくれる-no-prec-sqrt, /Qprec-sqrt-というオプションを持つ。

割り算の場合は128ビットのSSEでもRCPPSとNewton法を使ったほうがDIVPDを使ったプログラムより速くなることが多いのであるが、倍精度の53ビット精度を確保する平方根計算だと実は最近のCPU(WestmereとSandy Bridgeで確認)のSSEではSQRTPDの直接使用のほうが速い。一方AVXであれば平方根そのものの計算でもVRSQRTPSとNewton法を使うほうが速くなることが多い(が、割り算ほど劇的には速くならない)。どちらにせよ、平方根の逆数を求めるのであればSSEでもAVXでも(V)RSQRTPSとNewton法を使った計算は、遅い命令である(V)SQRTPDと(V)DIVPDの2命令の使用を避けることができるため非常に速くなる。頭の片隅に覚えておくと良いだろう。

Xeon Phiについて

Xeon Phiにも単精度逆数命令(VRCP23PS)や単精度逆数平方根命令(VRSQRT23PS)は存在するが、逆に除算や平方根のベクトル命令は存在しない。x87命令としては盲腸のように残っている(FDIV, FSQRT)が、意識的にx87命令をターゲットにしたコンパイルオプションを指定したり、x87命令に対応するアセンブリ命令を直接使ったりしない限り、除算や平方根のx87命令が使われることは無いだろう。除算や平方根の計算を素直に書けば倍精度でも自動的にこれらの命令が使われるようである。


この記事へのコメントをお寄せください

  • サイトへの書き込みに差し支えございましたら トータルディスクロージャーサイトサポート係へメールをお送りください
  • トータル・ディスクロージャ・サイトに投稿された文章と画像は、すべてその著作権がHPCシステムズ株式会社に帰属し、HPCシステムズ株式会社が著作権を所有することに同意してください。
  • あなたの文章が他人によって自由に編集、配布されることを望まない場合は、投稿を控えてください。
  • コメントを書き込む場合は名前にひらがなを織り交ぜてください。
  • あなたの投稿する文章と画像はあなた自身によって書かれたものであるか、パブリック・ドメインかそれに類する自由なリソースからの複製であることを約束してください。あなたが著作権を保持していない作品を許諾なしに投稿してはいけません!

<comments hideform="false" />


Comments

ノート:割り算は遅い

個人用ツール