区間演算法のスツルム二分法への適用実験

全文

(1)Vol.2013-HPC-139 No.16 2013/5/30. 情報処理学会研究報告 IPSJ SIG Technical Report. 区間演算法のスツルム二分法への適用実験村上弘1,a). 概要：二分法では数値の符号判定を行う必要があるが，通常の数値計算では数値誤差により符号の判定が正しくない可能性がある．数値の範囲を上下から押さえる数値区間演算法を導入することで，使用している数値の精度で符号を正しく判定できるかあるいは出来ないかを調べることができる．例として，スツルムの二分法に対して区間演算を用いることで，実対称三重対角行列の各固有値を正しく包含する区間を求めることを試みてみたい．キーワード：スツルム法，二分法，対称三重対角行列，固有値，バンチの方法. Application of Interval Arithmetics to Sturm’s Bisection Method Hiroshi Murakami1,a). Abstract: In the bisection method, we have to determine signs of numbers. However, those determined signs obtained by the usual numerical calculation under the presence of numerical errors may be wrong. If the numerical interval arithmetic method is used by which both lower and upper bounds for all numbers are calculated, within the precision of calculation in use we can determine whether the signs can be decided correctly or they can not be. As an example, we would like to try to obtain true enclosing intervals of eigenvalues of a real symmetric tridiagonal matrix by the application of the interval arithmetic to Sturm’s bisection method. Keywords: Sturm’s method, bisection, symmetric tridiagonal matrix, eigenvalue, Bunch’s method. 1. はじめに. できる．またスツルム法に類似した方法に，実対称行列や複素. 実対称あるいは複素エルミートの三重対角行列の固有値. エルミート行列に対してシルベスタの慣性律（Sylvester’s. は全て実数であり，その近似値を求める方法としては，ス. law of inertia）に基づく修正コレスキ分解やバンチの分解. ツルムの定理に基づく二分法，QR 反復法，最近のもので. を用いる方法もある．. は分割統治法が良く知られている．QR 反復法や分割統治法は，全固有値の近似値を求めるときには能率的とされる．他方で，スツルムの方法は，固有値の全部ではなくて一. 2. スツルム列を用いたカウント実対称三重対角行列 T の主対角要素を a1 , a2 , . . . , aN ，. 部だけを求める場合に特に能率的であって例えば，指定さ. 副対角要素を b1 , b2 , . . . , bN −1 として，T − xI の k 次首. れた区間内にある固有値を全て求めること，固有値に値順. 位行列の行列式を fk (x) と置くとき，f0 (x) = 1 と置. に番号を付けて指定範囲の番号の固有値を全て求めるこ. くとき，k = 1 では f1 (x) = (a1 − x)，k > 1 では. と，指定した値より大きい（あるいは小さい）固有値を値. fk (x) = (ak − x)fk−1 (x) − b2k−1 fk−2 (x) となり，T が既. の順に指定個数を越えない範囲で全て求めること，などが. 約（副対角要素に零が無い）である場合に多項式の列. 1. f0 (x), f1 (x), . . . , fN (x) がスツルム列となることを用いる．. a). 首都大学東京・数理情報科学専攻 Department of Mathematics and Information Sciences, Tokyo Metropolitan University [email protected]. ⓒ 2013 Information Processing Society of Japan. スツルム列中の隣り合う多項式の符号の変化数が x よりも小さい固有値の個数に一致する（但し零は飛ばして数え. 1.

(2) Vol.2013-HPC-139 No.16 2013/5/30. 情報処理学会研究報告 IPSJ SIG Technical Report. る）．実対称三重対角行列では固有値は全て実であることと，既約であるならば重複する固有値が無いこと，などが，一般的な実多項式の実根を求める場合のスツルム法に比べて事情を簡単にしている．スツルム法（Sturm’s method）のより詳しい原理については，参考文献（[4], [5], [6], [7]）などに譲り，説明を省略する．. 2.1 スツルムのカウントを行うプログラムスツルムの方法を Fortran90 で書いたサブルーチンを以下の図 1 に示す（RK は実数型の種別を表す整数値で用いるシステム依存の定数である．例えば Intel Fortan や. GNU Fortran など多くのシステムでは IEEE 754 の 64bit 倍精度実数型に対する RK は 8 である．（後で揚げる実行時間の計測では，DO ループ中の共通部分式 A(I)-X をループ冒頭で変数 AIX と置き，また I = 1 の場合をループの外に出して無駄な判定回数を減らすなどの書き換えをしている）．このプログラムは実対称三重対角行列 T の次数を N に，主対角要素 a1 , a2 , . . . , aN を配列 A(1:N) に，副対角要素 b1 , b2 , . . . , bN −1 を 2 乗した値を配列 BSQ(1:N-1) に入れて呼び出すと，行列 T の実数 X の値よりも小さい固有値の個数を，スツルムの定理に基づいてスツルム列の符号変化数を用いてカウントして，整数 ICNT に入れて返す．（行. SUBROUTINE STURM_CNT1(N,A,BSQ,X,ICNT) INTEGER,INTENT(IN)::N REAL(RK),INTENT(IN)::A(N),BSQ(N-1),X INTEGER,INTENT(OUT)::ICNT !-------------------------------------INTEGER I; REAL(RK) QOLD,Q,QNEW ICNT=0 DO I=1,N IF(I==1)THEN Q=1.0; QNEW=A(I)-X ELSEIF(BSQ(I-1)==0.0)THEN Q=1.0; QNEW=A(I)-X ELSE QOLD=Q; Q=QNEW QNEW=(A(I)-X)*Q-BSQ(I-1)*QOLD ENDIF IF(Q==0.0.OR.QNEW*Q<0.0)ICNT=ICNT+1 #define RENORMALIZE #ifdef RENORMALIZE ! 毎回あるいはたまに再規格化を行う． IF(Q==0.0)THEN; QNEW=1.0 ELSE; QNEW=QNEW/Q; Q=1.0; ENDIF #endif ENDDO END SUBROUTINE STURM_CNT1 図 1 スツルムの定理による固有値の計数（その 1）. 済ませることも可能であろう．. 列 T が複素エルミート三重対角の場合には，A はエルミート行列の主対角要素なので実数値を入れ，BSQ には実数値である副対角要素の絶対値の 2 乗を入れ，固有値のシフト. X も実数値として同じルーチンを用いて，カウントが計算できる．）ICNT 以外の引数である N，配列 A や配列 BSQ および X の値は変更されない．三重対角行列の副対角に零要素があって既約でない場合も正しく扱えるようになっている．マクロ#ifdef から#endif に挟まれた部分は，スツルム列の値がオーバーフローやアンダーフローを発生させて精度がなるべく低下しないように反復中に再規格化を行うためのもので，必ずしも反復ごとに行う必要はない（むしろ丸め誤差の導入をなるべく避けるために頻度を抑えた方が良いかもしれない）．また，このプログラム中で読み易さのために定数を 1.0 や 0.0 などと書いているが，Fortran の規格ではこれらは単精度の定数であるから，1.0 RK あるいは倍精度に限定するならば 1.0D0 などと書いた方がより安全である．隣り合うスツルム列の値の比を新たに Q と置いて. STURM CNT1 の算法を変形すると，新たに下記のサブルーチン STURM CNT2 が得られる（図 2）．中間変数が Q 一個で. SUBROUTINE STURM_CNT2(N,A,BSQ,X,ICNT) INTEGER,INTENT(IN)::N REAL(RK),INTENT(IN)::A(N),BSQ(N-1),X INTEGER,INTENT(OUT)::ICNT !-------------------------------------INTEGER I; REAL(RK) Q; LOGICAL QINF ICNT=0 DO I=1,N IF(I==1)THEN Q=A(I)-X; QINF=.FALSE. ELSEIF(BSQ(I-1)==0.0)THEN Q=A(I)-X; QINF=.FALSE. ELSEIF(QINF)THEN Q=A(I)-X; QINF=.FALSE. ELSEIF(Q==0.0)THEN QINF=.TRUE. ELSE Q=(A(I)-X)-BSQ(I-1)/Q ENDIF IF(QINF.OR.Q<0.0)ICNT=ICNT+1 ENDDO END SUBROUTINE STURM_CNT2 図 2 スツルムの定理による固有値の計数（その 2）. あるこの方式が通常よく用いられている算法に近い．論理変数 QINF は状態フラグで Q の値が無限大であることを真値で表している．浮動小数点数が IEEE754 であって，非. 2.2 二分法を適用する場合の注意. 数である正負の無限大とその符号の判定を Fortran で適切. 二分法として用いる場合，シフトの値 x0 ，x1（x0 < x1 ）. に利用できるならば QINF のような状態フラグを設けずに. におけるそれぞれのカウントの値を m0 ，m1 とするとき，. ⓒ 2013 Information Processing Society of Japan. 2.

(3) Vol.2013-HPC-139 No.16 2013/5/30. 情報処理学会研究報告 IPSJ SIG Technical Report. 数学的にはカウントは常に非減少 m0 ≤ m1 で，固有値は. ターの慣性律）．いま x を実数とするとき，A の固有値で. 区間 [x0 , x1 ) 内には m1 − m0 個ある．もしも m1 − m0 = 0. 「x より大きい」，「x に等しい」，「x より小さい」ものの個. ならば，区間 [x0 , x1 ) は固有値を含まないので，それ以上. 数はそれぞれ，x だけシフトされた行列 B = A − xI の正，. の分割を中止する．またもしも m1 − m0 < 0 となったら，. 零，負である固有値の個数に等しい．そこでシルベスター. 誤差の影響によりカウントの値が既に正しくない．. の慣性律を使うにはまず B を C に変換してから C の固有. x1 − x0 の値が予め与えた許容値（tolerance）TOL より. 値の各符号の個数を数えれば良い．B の C への変換と C. も大きければ再分割を試み，そうでなければその区間は更. の符号の計数が容易に実行できれば，スツルムの方法と同. なる再分割をせずに両端の値（あるいは区間の中点の値と. 様にして二分法で存在区間を任意に縮小することで固有値. 区間幅）とその区間中の固有値の個数 m1 − m0 を報告す. の近似値を求めることができる．. る．（またもしも与えた TOL の値が小さすぎて，再分割の. もっとも単純な適用は B = A − xI の修正コレスキ分解. ために作った中点 xc = (x0 + x1 )/2 の値が浮動小数点数と. B = LDLT で，M を単位下三角行列 L に，D を対角行列. して x0 = xc あるいは xc = x1 となる場合には，使用して. にした場合である．実対称行列 A が三重対角などの幅の. いる数値の精度では中間点がもはや存在しないので，それ. 狭い帯行列では計算量と記憶量の両方の面で適用は容易で. 以上の再分割を中止する）．また，中点 xc でのカウントの. あるが，行列 B が正則であっても不定値の場合には（シ. 値を mc とするとき，条件「m0 ≤ mc かつ mc ≤ m1 」が. フトが A の固有値分布の間にあると B は不定値になる）. 成立しない場合にも，誤差の影響によりカウントの非減少. 計算の途中で対角要素が零になって分解が破綻する可能. 性が崩れているものとして再分割を中止する．. 性がある．ピボット交換を対称性を保つように入れた分解. 全ての固有値を求めるのには，例えばゲルシュゴリン（Gershgorine）の定理で求めた固有値の存在範囲の上下限の両方を外側に向けて僅かに拡大した区間を初期区間. [x0 , x1 ) として採用して二分法を適用することができる．. P BP T = LDLT （P は置換の行列）でも破綻する場合が 0 1 でああり，その最も簡単な例は正則な 2x2 行列 1 0 る．さらに計算の安定性を狙ってピボット交換を行うと，その代償として行列の帯幅が広がるという困難もある．. 2.3 性能測定例 Intel Core i7-3770K（3.5GHz，HT オフ，TB オフ）の PC で，Intel Fortran Compiler v13.0 で最適化オプション. 修正コレスキ分解は途中でピボットが零とならなければ一応計算は可能だが，それだけでは分解結果が数値的には不安定になりうる．（但し，対角行列 D の符号分布だけが. -fast によりコンパイルした並列化無しの 1 スレッドでの. 必要な場合には，途中でピボットが零にならなければ計算. 実行例を示す．用いたテスト行列は主対角が全て 2，副対. を続けても，一応意味のある結果が得られるようにも思わ. 角が全て 1 のもので，全ての固有値は良く分離していて，固. れる．これはスツルム法のプログラム STURM CNT1 に於い. 有値の上限と下限はゲルシュゴリンの定理を用いるとそれ. て途中でスツルム列を再規格化しながら計算を進めてい. ぞれ 0 と 4 になる．このテスト行列の全ての固有値に対してその包含する区間幅を 10−12 以下まで 2 分法で狭めるのに掛かった合計の経過時間は，STURM CNT1 と STURM CNT2. ることや，特にスツルム法の STURM CNT2 に於いて途中で. Q=0 が起きない場合と同様である．もしもそうであるならば，ある対角シフトの値 x を選んで分解を行って，途中で. とを用いた場合にそれぞれ，次数 N = 1, 000 のときは 0.27. ピボットが零になって計算が破綻した場合でも，x の値を. 秒および 0.22 秒で，次数 N = 3, 000 のときは 2.36 秒お. 僅かにずらして再計算を行えば困難を実質的に避られるこ. よび 1.88 秒で，次数 N = 10, 000 のときは 24.8 秒および. とになりそうで，さらにピボット交換を省くこともできる. 19.8 秒で，次数 N = 30, 000 のときは 211.1 秒および 168.2. かもしれない．これらのことが実際にどうであるかについ. 秒で，次数 N = 100, 000 のときは 2,199 秒および 1,751 秒. ては調査研究してみるべきであろう．）. であった．（表 1 および図 1）．全ての固有値 N 個を求める合計時間が，概ね N の二乗に比例していることがわかる．. 3.1 バンチの方法. また，STURM CNT1 に比べて STURM CNT2 は約 1.25 倍程度. 数値計算では算法はなるべくならば安定であるのが望ま. 高速であった（注：STURM CNT1 のループ内で毎回ごとに. しい．不定値の実対称行列に対して対称性を有効に利用で. 再規格化をしている場合）．. きる数値安定な分解法として，バンチ（Bunch）の方法（文. 3. シルベスターの慣性律を用いるカウント法. いて，P が置換の行列で，L を単位下三角行列とするとこ. 実対称行列 B が実正則行列 M を用いて B = M CM T. ろまでは修正コレスキ法と同様であるが，行列 D は 1 次. と変換されるときは，C も実対称行列で，B と C でそれぞれ正，零，負である固有値の個数は一致する（シルベス. 献:[1], [2], [3]）がある．それは分解 P BP T = LDLT に於. あるいは 2 次の対称行列からなるブロック対角行列に拡張されている．D の正，零，負である固有値の個数は簡単に求められる（バンチの方法で生じる 2 次のブロック対角要. ⓒ 2013 Information Processing Society of Japan. 3.

(4) Vol.2013-HPC-139 No.16 2013/5/30. 情報処理学会研究報告 IPSJ SIG Technical Report. 素は正と負の固有値を各 1 個ずつ持つ）．そこで，慣性律に基づく固有値のカウンティングにバンチの分解が利用されている．一般の実対称な帯行列にバンチの分解法を適用するとき，対称性を保つピボット交換のたびに L の帯幅が次第に広がってしまう性質を持つが，行列 B が特に三重対角および五重対角である場合には，L の帯幅がピボット交換で広がらない方法が示されている．特に B が三重対角の場合にはピボット操作が不要で，分解で得られる単位下三角行列 L の下帯幅は 2 以下（(i, j) 要素は i > j + 2 では零）となる [1]．近似固有値を得たあとに固有ベクトルを逆反復法で求めたりあるいは改良を行う場合に，このバンチの方法で得られた分解を用いることもできるが，二分法で固有値を分離する計算では D の固有値の各符号のカウントだけが必要で，分解で得られる L や D の要素を記憶に格納する必要はない．シフト行列の正および負の固有値の個数をバンチの方法を用いて数える手続きを Fortran90 で書いたサブルーチン BUNCH CNT を以下の図 3 に示す．実対称三重対角行列. T の次数を整数 N とし，主対角要素 a1 , a2 , . . . , aN を配列 A(1:N) に，副対角要素 b1 , b2 , . . . , bN −1 の 2 乗の値を配列 BSQ(1:N-1) に入れておく．このルーチンを呼ぶと，整数 CNTP にはバンチの分解により求めた実数 X よりも大きい固有値の個数が，整数 CNTN には X よりも小さい固有値の個数が，それぞれ結果として返される．（X に等しい固有値の個数は N-CNTP-CNTN になる．）もちろん，これらは浮動小数点数を用いた計算で求めた個数なので，誤差の影響により数学的に厳密な値とは異なり得る．二分法を用いて区間を縮小する目的のためには，CNTP と CNTN のどちらか一方（例えば X よりも小さい固有値の個数をカウントするのであれば，CNTN）だけが必要である．. 3.2 バンチの方法を用いた性能測定例 Intel Core i7-3770K（3.5GHz，HT オフ，TB オフ) の PC で，Intel Fortran Compiler v13.0 で最適化オプション -fast によりコンパイルした並列化無しの 1 スレッドでの実行例を示す．用いたテスト行列は主対角が全て 2，副対. !================================================== ! COUNT THE NUMBERS OF POSITIVE AND NEGATIVE ! EIGENVALUES OF A REAL SYMMETRIC TRIDIAGONAL MATRIX ! WITH A DIAGONAL SHIFT BY THE BUNCH’S METHOD. !================================================== SUBROUTINE BUNCH_CNT(N,A,BSQ,X,CNTP,CNTN) INTEGER,INTENT(IN)::N ! THE SIZE OF MATRIX. REAL(RK),INTENT(IN)::A(N) ! MAIN-DIAGONALS. REAL(RK),INTENT(IN)::BSQ(N-1)! SQUARES OF SUB-DIAG. REAL(RK),INTENT(IN)::X ! DIAGONAL SHIFT. INTEGER,INTENT(OUT)::CNTP,CNTN ! COUNTS OF EIG. !-------------------------------------------------REAL(RK) COEF,ALP,BET,Q; INTEGER K COEF=(SQRT(5.0)-1)/2 ! NEED NOT BE SO ACCURATE. BET=MAX(MAXVAL(ABS(A(1:N)-X)), & SQRT(MAXVAL(BSQ(1:N-1)))) ALP=COEF/BET CNTP=0; CNTN=0 K=1; Q=(A(1)-X) LOOP:DO WHILE(K<=N) IF(K==N)THEN ! CASE THE LAST 1X1 PIVOT. IF (Q>0.0)THEN; CNTP=CNTP+1 ELSEIF(Q<0.0)THEN; CNTN=CNTN+1; ENDIF K=K+1 ELSEIF(BSQ(K)==0.0)THEN ! CASE BSQ(K)==0.0 (1X1 PIVOT). IF (Q>0.0)THEN; CNTP=CNTP+1 ELSEIF(Q<0.0)THEN; CNTN=CNTN+1; ENDIF Q=(A(K+1)-X); K=K+1 ELSEIF(ABS(Q)>=ALP*BSQ(K))THEN ! CASE 1X1 PIVOT. IF (Q>0.0)THEN; CNTP=CNTP+1 ELSEIF(Q<0.0)THEN; CNTN=CNTN+1; ENDIF Q=(A(K+1)-X)-BSQ(K)/Q; K=K+1 ELSE ! CASE 2X2 PIVOT. CNTP=CNTP+1; CNTN=CNTN+1 IF(K<=N-2)THEN Q=(A(K+2)-X)-BSQ(K+1)*Q/(Q*(A(K+1)-X)-BSQ(K)) ENDIF K=K+2 ENDIF ENDDO LOOP END SUBROUTINE BUNCH_CNT 図 3 バンチの方法による固有値のカウント. 角が全て 1 のもので，全ての固有値は良く分離していて，ゲルシュゴリンの定理から固有値の上限と下限はそれぞれ 0 と 4 である．このテスト行列の全ての固有値に対してその. 表 1 全固有値を二分法で求める経過時間（秒）次数 N STURM CNT1 STURM CNT2 BUNCH CNT. 1,000. 0.27. 0.22. 3,000. 2.36. 1.88. 1.78. ントによる 2 分法で狭めるのに掛かった合計の経過時間は，. 10,000. 24.8. 19.8. 18.8. 次数 N = 1, 000 では 0.20 秒で，次数 N = 3, 000 では 1.78. 30,000. 211.1. 168.2. 161.1. 秒で，次数 N = 10, 000 では 18.8 秒で，次数 N = 30, 000. 100,000. 2,199. 1,752. 1,683. 包含する区間幅を 10−12 以下まで BUNCH CNT を用いたカウ. 0.20. では 161.1 秒で，次数 N = 100, 000 では 1,683 秒であった．（表 1 および図 1）．全ての固有値 N 個を求める合計時間が，概ね N の二乗に比例していることがわかる．また，. STRUM CNT2 と計算時間がほぼ同じであることもわかる． ⓒ 2013 Information Processing Society of Japan. 4. 区間演算の導入の試み計算に現れる数値の絶対確実な上限と下限を常に求めな. 4.

(5) Vol.2013-HPC-139 No.16 2013/5/30. 情報処理学会研究報告 IPSJ SIG Technical Report. ELAPSE TIME (SEC). 1000. 点とする区間による区間演算が実現できる（文献 [8], [9]）．. STURM 1 STURM 2 BUNCH O(N2) curve. 但し現状の CPU の構成では演算モードを変更する際には多大なオーバーヘッドを伴う．. 100. 10. 1. 0.1. 1000. 10000. 100000. MATRIX SIZE N. 図 4 全固有値を二分法で求める経過時間（秒）. がら計算を行うことにより，例えば符号の判定が保証付きで出来る場合とそうでない場合を区別する．それにより（用いた精度の計算で），結論の正当性を保証できるかまたは保証ができないことが示せる．実数の演算を実数の区間に対するものに拡張して，区間同士の演算（例えば加減乗除）を導入する．区間の両端を浮動小数点数（FP 数）で置き換えて，区間同士の演算の際に行う FP 数の演算による丸め誤差の混入を，丸めを行うモードを制御することにより，区間同士の演算結果が数学的な正しさを失わないようにする．つまり FP により実現した区間演算の結果の区間は丸め誤差の影響分だけ余分に広がる．通常は演算を繰り返すたびに，丸め誤差の影響が蓄積して区間は広がっていくので，確実な判断が導けないことが示されることもある．それでは区間演算を実対称三重対角行列に対するスツルム法あるいはそれの類似算法に適用して固有値の個数の計数を行うとどのようになるだろうか，それが調べたいことである．. 4.1 FP 演算の丸め方式の制御図 5 の C プログラムは，Intel や AMD 製の x86，x86 64 系 CPU の FPU や SSE に対する，IEEE754 の FP 演算の際の丸め方式を制御するためのもので，名前の最後にアンダースコアが付いている関数は Fortran からリンクし. /* For x86, x86_64 FPU. */ #include <float.h> #include <fpu_control.h> void fpu_near() {short cw= _FPU_DEFAULT; _FPU_SETCW(cw); } void fpu_chop() {short cw= _FPU_RC_ZERO | _FPU_DEFAULT; _FPU_SETCW(cw); } void fpu_up() {short cw= _FPU_RC_UP | _FPU_DEFAULT; _FPU_SETCW(cw); } void fpu_down() {short cw= _FPU_RC_DOWN | _FPU_DEFAULT; _FPU_SETCW(cw); } short fpu_getcw() {short cw; _FPU_GETCW(cw); return cw; } void fpu_setcw(short cw) { _FPU_SETCW(cw); } /* For MMX/SSE/SSE2. */ #define _SSE_DEFAULT 0x00001f80 #define _SSE_RC_DOWN 0x00002000 #define _SSE_RC_UP 0x00004000 #define _SSE_RC_ZERO 0x00006000 #define _SSE_GETMXCSR(cw) \ __asm__ ("stmxcsr %0" : "=m" (*&cw)) #define _SSE_SETMXCSR(cw) \ __asm__ ("ldmxcsr %0" : : "m" (*&cw)) void sse_near() {int csr= _SSE_DEFAULT; _SSE_SETMXCSR(csr); } void sse_chop() {int csr= _SSE_RC_ZERO | _SSE_DEFAULT; _SSE_SETMXCSR(csr); } void sse_up() {int csr= _SSE_RC_UP | _SSE_DEFAULT; _SSE_SETMXCSR(csr); } void sse_down() {int csr= _SSE_RC_DOWN | _SSE_DEFAULT; _SSE_SETMXCSR(csr); } int sse_getmxcsr() {int csr; _SSE_GETMXCSR(csr); return csr; } void sse_setmxcsr(int csr) {_SSE_SETMXCSR(csr); } /* For FPU and MMX/SSE/SSE2. */ void set_round_near_() {fpu_near(); sse_near(); void set_round_chop_() {fpu_chop(); sse_chop(); void set_round_up_() {fpu_up(); sse_up(); void set_round_down_() {fpu_down(); sse_down();. } } } }. て呼ぶためのものである．Fortran からサブルーチンとして SET ROUND NEAR を呼ぶと FPU/SSE の演算機構の状態が演算結果の値を最も近い FP 数の値に丸める（標準の）. 図 5 x86，x86 64 系 CPU の FPU/SSE の丸めモード制御用 C プログラム. モードに設定され，SET ROUND CHOP を呼ぶと演算結果の値を 0 の方向で最も近い FP 数に切り捨てるモードに，. SET ROUND UP を呼ぶと演算結果の値を正の無限大の方向で最も近い FP 数に丸めるモードに，SET ROUND DOWN を呼ぶと演算結果の値を負の無限大の方向で最も近い FP 数に丸めるモードに，それぞれ設定される．これらを使って. FP 演算の丸めモードを適切に切り替えて使うことにより，丸めを完全に意識した演算が可能になり，特に FP 数を端 ⓒ 2013 Information Processing Society of Japan. 5.

(6) 情報処理学会研究報告 IPSJ SIG Technical Report. Vol.2013-HPC-139 No.16 2013/5/30. 5. おわりに実対称行列の固有値の個数をカウントするスツルムの定理（あるいはシルベスタの慣性律）に基づく方法は通常，数値と演算に浮動小数点（FP）を用いて近似計算されている．精度を固定した近似計算であっても数値と演算に（FP 数を両端とする）区間数を用いたものを使うことにより，求めた固有値の個数のカウントが正しいことを保証できる場合とそうでない場合とをはっきりと区別できる．以上のことまでを実証する予定でしたが，本予稿の作成時点に於いては，具体的実験例を記述するまでに至りませんでした．現在，行列が 100 次程度の小規模の場合にだいたい動作するプログラムができたところで，まだ動作の不良などがあり確認をしているところです．研究会当日までには，倍精度の FP 区間数による演算による多少の実例（うまくいく場合も行かない場合も含めて）をお見せできるかと思います。その他，本予稿の内容の記述にも至らぬ点が多いことを，お詫び申し上げます．参考文献 [1]. [2]. [3]. [4]. [5] [6] [7] [8] [9]. James R.Bunch:”Partial Pivoting Strategies for Symmetric Matrices”, SIAM J.Numer.Anal., Vol.11, No.3, (June,1974), pp.521–528. 特に p.526 からの §4.”Another Approach”の記述． James R.Bunch, Landa Kaufman and Beresford N.Parlett:”Decomposition of a Symmetric Matrix”, Numer.Math., Bd.27, (1976), pp.95–109. James R.Bunch and Linda Kaufman:”Some Stable Methods for Calculating Inertia and Solving Symmetric Linear Systems”, Math.Comp., Vol.11, No.137, (Jan.,1977), pp.163–179. Gene H.Golub and Charles F.van Loan: ”Matrix Computations”, 3rd Edition, The Johns Hopkins University Press, Baltimore and London, 1996. 特に §8.5.1 の”Eigenvalues by Bisection”と §8.5.2 の”Sturm Sequence Methods”．高木貞治:「代数学講義」, 改訂新版, 共立出版,1965 年. 特に第 3 章「スツルムの問題，根の計算」の記述．森正武:「数値解析法」, 朝倉現代物理学講座 7, 朝倉書店,1984 年. 森正武:「数値解析」, 第 2 版, 共立数学講座 12, 共立出版,2002 年. 大石進一：「精度保証付き数値計算」, コロナ社, 1999. 大石進一：「Linux 数値計算ツール」, コロナ社, 2000.. ⓒ 2013 Information Processing Society of Japan. 6.

(7)