5.1検証方法について
今回のLAPACKルーチンの移植では、実行速度を重視し演算結果に対し影響を及ぼす最適化レベル-Oe用い、演算結果を個々検証して、その影響を最適化制御行の挿入で抑止するという方針を採った。ここでは、その検証方法と最適化の影響について述べる。
最適化の影響の検証には、netlibから配布されるLAPACKパッケージに含まれる、演算結果および性能評価のためのテストセットを使用した。
具体的には、次の手順で行った。
- LAPACKルーチンを演算結果に影響を与える最適化レベル-Oeと基本的な最適化レベル-Obを指定してコンパイルしたLAPACKライブラリをそれぞれ用意する。
- テストセット本体は最適化の影響を受けないように-Obでコンパイルする。
- それぞれの最適化レベルのLAPACKライブラリを組み込んだテストセットを二種類作成し、同じテストデータを用いてそれぞれ実行する。
さらに、最適化の影響部分を特定するためには、テストセットおよびLAPACKルーチンを抽出して、WRITE文を挿入し、中間データでの演算結果の比較も必要であった。
5.2最適化の影響と対処
今回の検証では、二種類のテストセットの結果を比較すると演算結果が異なるものが幾つかあったが、許容誤差の範囲では問題にならないものであると思われた。 ここでは、最適化の影響を受けテストセットの実行で演算例外が発生したものについて紹介する。
5.2.1ループ内でのスカラに対するIF文がある
-Obで作成したテストセットは正常に実行されるが、-Oeで作成したテストセットの三重対角行列の連立一次方程式を解くテストセットでゼロ割りのエラーが発生した。
実行時のエラーメッセージからエラーの発生箇所は、図15に示すルーチンDGTTRS内のループ(DO30)であった。
このDGTTRSは、DGTTRFルーチンで求めたLU分解を用いて三重対角行列Aをもつ連立一次方程式
を解くルーチンである。
 |
図15. DGTTRSのベクトル化リスト
コーディングを眺めると上三角行列を解く処理において、ループ(DO 30)内にIF文が置かれ、変数Nが1より大きい場合だけ、
| ( N-1, J ) = ( B( N-1, J )-DU( N-1 )*B( N, J ) ) / D(N-1) |
を実行するように制御されている。
このようなコーディングがある場合に、最適化処理としては、まず、この式に現れるD(N-1)による除算をD(N-1)の逆数の乗算に置換える。さらに、逆数を求める演算は不変式となるので、ループ(DO 30)外に移動する。このような場合、結果としてIF文の制御からはずれてしまい、変数Nが1の場合にD(0)に対する参照が起こり、本来起こらないゼロ割のエラーが発生する。
-Oeでコンパイルした場合の最適化メッセージを確認すると、
| jwd8201i-i "dgttrs.f", line 136: Invariant expression was moved. jwd8206i-i "dgttrs.f", line 136: Division was changed to multiplication. |
という二つのメッセージが当該IF文について出力されていた。
このような最適化の影響を回避するために、ループ(DO 30)の前に最適化制御行!ocl noevalを挿入し、このループ内での計算式の評価順序の変更を伴う最適化を抑止した。
なお、他の三重対角行列をもつ連立一次方程式の解を求めるルーチンにも同様なコーディングが存在する可能性があるので調べた結果、DGTSVルーチンにも同じコーディングがあったので、同様に最適化を抑止した。
5.2.2情報落ちによるゼロ割りが発生する
-Oeで作成したdivide and conquer法を用いて固有値、固有ベクトルを求めるLAPACKルーチンDSTEDCのテストセットでゼロ割りのエラーが発生した。
実行時のエラーメッセージからエラーの発生箇所は、図16に示すDSTEDCのスレーブルーチンDLAED4のループ(DO 40)内で起こっている。すなわち、配列DELTAの要素にゼロが含まれてしまうことが原因である。コーディングを調べると配列DELTAの値は、同じルーチン内で以前に実行されるループ(DO 30)で求められているので、このループに対する最適化の影響が考えられた。
 |
図16. DLAED4のベクトル化リスト
ベクトル化リストを見る限り、このループ(DO 30)は制御変数Jでベクトル化されており、ループ内のコーディングの式通りに実行される限り、配列DELTAはJ = Iの場合一旦ゼロが求まるが後の-TAUがあるのでゼロが求まることはあり得ない。
デバッグとしてループ(DO 30)の後に、配列DELTAの値がゼロとなった場合の関連する変数の値を打出すループを挿入すると、
・配列DELTAがゼロになるのはI=J=28の場合である。
・また、D(28)は5.6617e-04、TAUは8.1315e-25
という結果になった。
ループ(DO 30)の式は括弧により演算の順序を明示してあるが、これを無視した最適化が行われ、D(I)-TAUのスカラ演算を先に行い、その結果のスカラとD(J)のベクトルとの引算が行われるとD(I)とTAUの差が倍精度実数の有効桁(10進15桁)を越えているので、情報落ちが発生しD(I)-TAU = D(I)となりゼロが求まってしまうことが考えられた。
Fortranの言語規格では、演算式の括弧を保証することを求めている。この件に関してメーカに調査を依頼したが、「ベクトル演算では、演算式の括弧を無視した最適化を行う」との回答であった。このループ(DO 30)の式で括弧を保証した場合の命令は、
| Vreg = Vdj - Sdi Vdelta = Vreg - Stau |
という二つのベクトルとスカラの引算で処理される。これをレベル-Oeで最適化すると、
| Sreg = Sdi + Stmp Vdelta = Vdj - Sreg |
というように先にスカラ同士の演算とその結果とベクトルの引算となり、一つのベクトル演算を無くすことができるので高速化が可能となる。
DSTEDCのスレーブルーチンDLAED4でのゼロ割りの発生は、最適化の影響であることは明らかである。これを回避するには、ループ(DO 30)の前に最適化制御行!ocl noevalを前に挿入し、このループ内での計算式の評価順序の変更を伴う最適化を抑止する方法が考えられるが、このルーチンのコーディングには、同様なループが複数あったので、プログラムの先頭に最適化指示行を挿入して、サブルーチン全体に対してこの最適化を抑止した。