行数 : 58 , 配列 : real*8 h(1023,1023) (8.0MB)

倍精度実数の1023×1023次元 正方行列をベクタに掛ける操作を 6000回 繰り 返す。

ここで、6000=NSOL × ITER であり、NSOL は9行目の、ITER は10行目の代入 文で、値がそれぞれ 3 と 2000 に設定されている。計算時間が、測定の便宜 上、長すぎるときや短すぎるときには、10行目の代入文「ITER=2000」の右辺 の値を変えると、その値に比例して行列×ベクタの反復回数が変わる。例えば、 ITER=500として計算した結果得られたCPU時間を 2000/500 倍すれば、我々の 示す表の数値と直接比較できるCPU時間の値が得られる。

このベンチマーク・テスト第２版では、すべてのコードにおいて、計算が行わ れたことを確認するために、なるべく意味のある量を計算して出力するように している。また、複数のCPUを持ち並列実行が可能なマシンの場合、大きなか たまり(matvec.fの場合は、行列×ベクタというかたまり)が複数個、同時に実 行されないように、逐次的に計算しなければ求まらない量を出力するようにして いる。

このコードで出力されるのは、冪乗法によりもとめた行列の固有値で、絶対値 の大きい順に３個出力される。ただし、このコードはベンチマークテストのた めを第一に考えて作ったコードであって、行列の固有値をもとめるために設計 したものではないので以下の☆で示したような不十分な事項があり、行列の 固有値を求める目的には他のコードを使用するのが賢明である。
☆ベンチマーク・テストは計算量が一定であることが必要なの で、マシンに依存しうる数値表現や演算の誤差の微小な差によって反復回数が 変わらないようにするため、収束状況を見て反復回数を決めるのでなく、指定 された回数だけ反復動作をするアルゴリズムに変更してある。このため結果の 精度は保証されない。
☆固有値を求めることの効率の点からは、冪乗法は非効率的であり Lanczos法等の他の手法によるべきである。
☆行列は対称行列であるが、行列×ベクタの計算では、その対称性 は利用していない。

行数 : 66 , 配列 : complex*16 h(723,723) (8.0MB)

matvec.f の倍精度複素数版である。次元は723×723であるが、これは、 配列のサイズがmatvec.fとほぼ同じになるように設定したのである。 倍精度複素数の正方行列を倍精度複素数のベクタに掛ける操作を 3000 回繰り返す。

ここで、3000=NSOL × ITER であり、NSOL は16行目の、ITER は17行目の代入 文で、値がそれぞれ 3 と 1000 に設定されている。計算時間が、測定の便宜 上、長すぎるときや短すぎるときには、17行目の代入文「ITER=1000」の右辺 の値を変えると、その値に比例して行列×ベクタの反復回数が変わる。例えば、 ITER=300として計算した結果得られたCPU時間を 1000/300 倍すれば、我々の 示す表の数値と直接比較できるCPU時間の値が得られる。

行列はエルミートであるが、その対称性は、行列×ベクタの操作では利用していない。 計算内容は、その行列の固有値(エルミートであるため実数である)を絶対値の大きい 順に3個求めることである。

行数 : 102 , 配列 :real*8 h(1001,1001) (7.6MB)

1001次元の連立一次方程式をガウス法で解くことを 5 回繰り返す。前回の解が 次の方程式に影響を与えるため、解に応じて次の方程式が変更されるため、 解は逐時的に求めなければならず、同時に並列に求めることはできない。

7行目の「CALL LEQS2(1.0D0, 1.0D0, 1.0D0, 5)」の最後の引数 5 が反復回数 をあらわしている。この回数に 概ね比例して CPU時間が変わる。

このテストや、その配列サイズを変えただけのテストの leqs4h では、高級ワー ク・ステーションの IBM RS6000 や HP9000などの一部のマシンで、他のテス トの結果から予測されるより、ずっと長い計算時間がかかるという結果が出た。 その原因は詳しくは未検討であるが、メモリ・アクセスの番地の飛びかたが不 得意な場合があるとか、配列添字の計算が苦手だとか、いろいろ仮説は浮かぶ。

計算結果は特別な意味のある量ではないが、計算を行わなければもとまらない 量を書かせることで、計算が行われたことを確認する意味がある。

ガウス法による連立一次方程式を解くサブルーチン SLEQS は、ベンチマーク テスト専用に若干の修正がされているので、一般の用途には、他のサブルーチ ンを利用するのが望ましい。

なお、このサブルーチンを始めベンチマークの全てのコードは田嶋が full scratch から作ったものだけを用いているので、すべて、使用、修正、ご自分 のコードへの組み込み等、ご自由になさっていただいて、かまいません。

行数 : 102 , 配列 : real*8 h(401,401) (1.2MB)

401次元の連立一次方程式を標準的な方法であるガウス法で解くことを78回繰 り返す。leqs1k.f との相違は次元が 1001 から 401 に減っていることと、解 を求める回数が 5 から 78 に増えていることである。回数を増やしたのは、 leqs1k.fと同じくらいのCPU時間になるように調節するためであり、 「78 = 5*(1001/401)^3」 として決めたものである。

行数 : 168 , 配列 : real*8 h(101,101), v(101,101) (159KB), NB: h0(101,101) is rarely accessed

101 次元の正方対称行列の固有値をヤコビ法で求めるこ とを 200 回反復する。より正確には、ヤコビ変換の特別巡回を 1000 (= NCYCLE × LMAX = 5 × 200 )回 繰り返す。

なお、「Pentium-Pro 200MHz cache 256KB, linux, g77」というシステムの場 合について計測した結果によると、サブルーチン JACOBI(ヤコビ変換などを行う) 以外の部分で CPU 時間の 15% が消費されている。

機械に依存しうる計算誤差に左右されずに計算量が一定であることを保証する ため、収束状況によらず特別巡回を ちょうど 5 回だけ繰り返すように設定し てある。このため、このヤコビ法のサブルーチン JACOBI はこのベンチマーク テスト以外の用途に流用しないのが賢明である。

行数 : 168 , 配列 : real*8 h(11,11),v(11,11) (1.9KB), NB:h0(11,11) is rarely accessed.

11次元の正方対称行列の固有値をヤコビ法で求めることを繰り返す。より正確 には、ヤコビ変換の特別巡回を1,000,000 (= NCYCLE × LMAX = 5 × 200000 ) 回 繰り返す。

行数 : 78 , 配列 : なし

４変数の連立一次微分方程式を４次のRunge-Kutta法(1/6公式)で 200,000,000 ステップ解く。このステップ数は、16 行目のパラメータ文中の 「N=200000000」で与えられている。CPU時間は、N に比例する。

内容は、平面内の質点に対する Newton の運動方程式の解である。ポテンシャ ルエネルギーは LOG_E(R), R=SQRT(X**2+Y**2)である。計算に必要なのはポテ ンシャルでなく、その gradient であるが、その計算には LOG 演算、SQRT 演算 は不要である。

行数 : 38 , 配列 : なし

反復 200回のループを４重にネストしたものの中で、積分計算を行う。 12行目のパラメータ文中の「N=200」でループ長の200が設定されている。 CPU時間は、Nの4乗に比例する。

コードでは、正方形内に一様に分布する独立な２点の位置座標について、それ らの４座標を変数とするある被積分関数を中点則で積分している。ただし、そ の被積分関数は特に意味のある関数ではない。

このテストはベクトル化に好適なコードであるが、配列を使用しないので、ベ クターのロード・ストア命令をもたないワーク・ステーションでも、スーパー・ コンピュータに匹敵できるのではないかと期待していたが、結果を見ると、スー パー・コンピュータのS3800/480(500MHz)がUltraSparc2(296MHz)の 10.5 倍の 速度を出していた。クロック比の補正をしても前者が6.3倍効率がよい。パイ プライン数(加算・乗算８本、除算１本) と スーパー・スカラーの並列実行可能 数(不詳) の差がでているのだろう。あるいは、ワークステーションでは、除 算パイプラインへのフィードサイクルが１クロック以上かかる(？)というよう な事情があるのかもしれない。

行数 : 23 , 配列 : なし

Fortranの内部関数の計算速度を見るため、 SIN(X)+SQRT(COS(X)) を 100,000,000 回計算する。すぐあとで示すように、SIN, COS, SQRT は 内部関数のなかでも計算時間がかからない部類に属している。 なお、このコードは並列実行が可能である。

出力されるのは、上記の関数の x が 0.1 から 0.9 までの積分値である。

参考に、各内部関数について、10,000,000回計算するのに必要なCPU時間を計測 した結果を下記に示す。system A は「Pentium Pro 200MHz, Linux, g77」、 system Bは 「Pentium II 300MHz, Linux, g77」である。

       (system A)        (system B)
            (sec)             (sec)
  TANH      16.80    SINH     10.51
  SINH      16.73    TANH     10.48
  ACOS      15.01    ACOS      9.96
  ASIN      13.61    ASIN      8.87
  COSH      13.11    COSH      8.34
  EXP       10.53    EXP       7.22
  ATAN2      9.31    ATAN2     6.26
  TAN        8.96    TAN       6.13
  ATAN       8.42    ATAN      5.84
  LOG10      7.62    LOG10     5.19
  LOG        6.58    LOG       4.49
  COS        5.90    COS       3.99
  SIN        5.38    SIN       3.89
  SQRT       3.47    SQRT      3.04
  NO FUNC    0.52    NO FUNC   0.77  ←loopを回すのに必要なCPU時間

行数 : 23 , 配列 : なし

Fortranの内部関数の計算速度を見るため、 ATAN2(EXP(X),LOG(X+1.0D0)) を 100,000,000 回計算する。先に示したように、EXP, LOG, ATAN2 は 内部関数のなかでは計算時間がかかる部類に属している。 なお、このコードは並列実行が可能である。

出力されるのは、上記の関数の x が 0.1 から 0.9 までの積分値である。

行数 : 52 , 配列 : なし

整数のみをあつかうあまり単純でない計算である。配列は使用しない。３重ルー プのなかに、IF文による飛び出し、加算、乗算、IF BLOCK(中に加減乗除算を 含む)がある。IF BLOCK は4バイト整数の溢れ対策のためにあるので、 中身が実行される確率は低い。

計算結果は、円周率πの近似値である。 半径 N の球の内部で、N**2 - R**2 を積分した結果をもちいて近似する。

サブルーチン INTPI に渡される引数 PARM1 はダミーである。引数が全く ないと、コンパイルの段階で計算を行ってしまうということもありうると考えた ので つけたものである。ループ長を引数にしなかったのは、そうするとコンパ イラによる最適化の効率が落ちることがあると思ったからである。

行数 : 21 , 配列 : なし

「K=(205*J)/103-I ; I=J ; J=K」という整数だけの操作を 1,000,000,000 回反復する。

出力されるのは、反復終了後の I と J の値である。

逐次的に実行するとすれば、ほとんど、クロックに比例するはずだが、 いろいろのマシンをためしてみると、ループ・アンローリングをして 加速していると解釈したくなるような、早いマシンがある。あるいは、 コンパイルのときに途中まで計算しているということはないか。

ベンチマーク・テスト(第一版)の intosc.f との違いは、反復回数が約 15.7 倍になったことである。なお、 20907、20908番目の I が、1、2番目と同じ値 (それぞれ、0、 1045)に戻ってしまうため、以降は周期20906のループに入っ てしまう。したがって、NMAX=1000000000回計算を繰り返す必要は無く、 MOD(NMAX,20906)回の計算で最終の出力を出すには十分なのだが、コンパイラ が自動的にそれを考慮する能力をもつことは近い将来にはないだろうから、こ のままでベンチマークテストの一つとして継続して使うことにした。

行数 : 51 , 配列 : integer*4 p(1024) (4KB)

1024個のinteger*4変数からなる配列から２要素を選び値を交換する操作を 1024*262144 = 268435456 = 2^28 回 反復する。配列のサイズは4KBであるか ら、CPU内のキャッシュ・メモリーにおさまる大きさである。したがって、こ のテスト・コードの実行速度は、主に CPU cache へのアクセスの速さに依存 する。あるいは、合同乗積乱数の生成などにより多くのCPU時間を消費してい るかもしれない。

要素の選び方は、疑似的にランダムになるように、値を交換すべき 第一の要素と第２の要素のインデックスを別々の合同乗積乱数系列で生成する。 合同乗積乱数系列は I = MOD(I*A+B,M) として逐次求めて行くが、Mが2の冪乗 で、MOD(A,4)=1、MOD(B,2)=1 のときには、0,...,M-1のすべての値を１度ずつ とる周期 M の系列になることが知られている。このコードでは、そのような、 周期 M の合同乗積乱数系列をもちいている。

行数 : 51 , 配列 : integer*4 p(32768) (128KB)

32768個のinteger*4変数からなる配列から２要素を選び値を交換する操作を 32768*8192 = 2684354562 = 2^28 回 反復する。配列のサイズは128KBであるか ら、level-2 キャッシュ・メモリーにおさまる大きさである。したがって、こ のテスト・コードの速度は、主に cache memory へのアクセスの速さを反映す る。

行数 : 51 , 配列 : integer*4 p(1048576) (4MB)

1048576個のinteger*4変数からなる配列から２要素を選び値を交換する操作を 1048576*256 = 2684354562 = 2^28 回 反復する。配列のサイズは 4MB であり、 キャッシュ・メモリーにはおさまらない。したがって、このテスト・コードの 実行に要する時間は、メモリー本体へのアクセスの速度に依存するはずである。