数理解析研究所講究録
915 巻 1995 年 28-33

28

スカラー計算機とベクトル計算機における固有値コードの性能評価

別府良孝 (Yoshitaka BEPPU)

聖徳学園女子短大・情報系

51. はじめに
実対称三重対角行列$\mathrm{T}$に関する標準固有値問題を解く$\mathrm{T}\mathrm{v}_{\mathrm{i}}=\lambda \mathrm{i}\mathrm{V}_{\mathrm{i}}(\mathrm{i}--1,2, \cdots, \mathrm{m}, \cdot\cdot \mathrm{n})$)種々の固有値コードが開発されているが、 それらのコードを速度と精度の両面から性能評価したので報告する。 固有値コードとしては、 分割統治コード・ホモトピーコ$-$ ト ‘IUM AIAA Aを取り上げ、 スカラー計算機とベクトル計算機における測定結果を報告する。
ために

$1$

$\backslash$$\mathrm{N}$$\mathrm{C}\mathrm{E}\mathrm{R}\cdot \mathrm{N}$$\mathrm{P}$$\mathrm{C}\cdot \mathrm{E}$$\mathrm{S}\mathrm{P}$$\mathrm{C}\mathrm{K}\cdot \mathrm{L}$$\mathrm{C}\mathrm{K}$$\mathrm{P}$

52. テスト行列について以下の行列をテスト行列として用いた。 これらは、$\mathrm{H}$$\mathrm{J}_{\text{、}}$.

C. G. L. W. F.

行列と略記する。

Jahn
Constant
GradedLaguerre.

$\mathrm{d}_{\mathrm{i}}=\mathrm{m}+1,\mathrm{s}_{2}=\mathrm{k}\sqrt{\mathrm{i}},\mathrm{s}_{2\mathrm{i}}-1=\mathrm{k}\sqrt{\mathrm{m}+\mathrm{i}}$$\mathrm{d}_{\mathrm{i}}=4,\mathrm{s}_{\mathrm{i}}=1$$\mathrm{d}_{\mathrm{i}}=\mathrm{i},\mathrm{s}_{\mathrm{i}}=1$$\mathrm{d}_{\mathrm{i}}=\mathfrak{U}-1,\mathrm{s}_{\mathrm{i}}=\mathrm{i}$Wilkinson $\mathrm{d}_{\mathrm{i}}=$$\mathrm{F}\mathrm{r}\mathrm{A}\mathrm{c}\mathrm{k}$$\mathrm{A}_{\mathrm{i}\mathrm{j}}=\mathrm{n}+1$

Block$\mathrm{s}_{\mathrm{i}}=-1$-max(i,j)$\mathrm{A}_{\mathrm{i}\mathrm{j}}=4\delta(\mathrm{i},\mathrm{j})-\sum_{\mathrm{k}=1\mathrm{m}}^{1}\sum_{=0}\mathrm{B}\mathrm{r}-\mathrm{r}-2$

where$\mathrm{r}=\sqrt{\mathrm{n}}$$\mathrm{B}=\delta(\mathrm{i},\mathrm{k}+1+\mathrm{m}\mathrm{l})\cdot\delta \mathrm{O}^{\mathrm{i}},+\Gamma-1)$.

$+\delta(\mathrm{i},\mathrm{k}+\mathrm{n}\mathrm{n})\cdot\delta\zeta|,\mathrm{i}+\mathrm{r}+1)$$+\mathrm{t}\mathrm{e}\mathrm{x}\mathrm{c}\mathrm{h}_{\mathrm{A}\mathrm{g}}\mathrm{e}$RandomHuckel$\mathrm{a}$$\mathrm{i}$and }$-32768<\mathrm{A}_{\mathrm{i}\mathrm{j}}<+32768$$\mathrm{A}_{\mathrm{i}\mathrm{j}}=-7.2\delta(\mathrm{i},\mathrm{j})-\frac{3.0}{(\mathrm{i}-\mathrm{j})^{2}}\{1-\delta(\mathrm{i},\mathrm{j})\}$$\mathrm{j}$$\mathrm{B}_{\text{、}}\mathrm{R}_{\text{、}}$

2953.

$\mathrm{Q}\mathrm{R}$法におけるシフト戦略について 2)

ループに入る前法” の収束を加速するためのシフト戦略として、“法” で求める ループ前究極シフト法” や QR ループに全ての固有値を 無平方根
“平方根付き$\mathrm{Q}\mathrm{R}$$\mathrm{Q}\mathrm{R}$“$\mathrm{Q}\mathrm{R}$の内で固有値を–個つつ求める “ループ内高次シフト法” が提案されてきた。 “ループ内高次シフト法” の場合、与えられた がデフレーションの結果、 次の左上三重対角行列次の右下対角行列とに分解されるので、 シフト量 としては、 次の左上三と重対角行列の右下 L 次の小行列の固有値を用いればよい。 小行列の固有値は、 “無平方根ラゲール法 3) のいずれかを用いて見積もる事ができる
ニュートン法”
法”$\mathrm{k}$$\mathrm{T}$$\mathrm{k}$$(\mathrm{n}-\mathrm{k})$$\mathrm{Q}\mathrm{R}$.

$\mu$.

“”“

既発表のデータ 2) および新規データに基づいて、 以下のごとく結論した。
(a) ループ前究極シフト法について
このシフト法は、 近接・縮退固有値が苦手で、 2$\mathrm{n}$語の作業領域が必要逆反復法に劣る場合が多い。 だからこそ、Parlettb [The fastest way to compute a full set of eigenvectors is
to execute QL algorithm with the fast Give..ns rotation using theである。 速度も精度も、$\mathrm{L}\mathrm{U}$previously computed eigenvalues as origin shifts. This combinationof techniques seems to have been overlooked so far.

$\rfloor$と述べているのであろう。

(b) ループ内高次シフト法について$\mathrm{b}1$$\mathrm{b}2$. ニュートン法やラゲール法は、 近接・縮退固有値が苦手で、 あふれの心配が有る。 高速化したければ、 L 語の作業領域が必要である。
法は、 近接・縮退固有値にも強く、 あふれの心配は無平方根.

$\mathrm{Q}\mathrm{R}$無いが,$\mathrm{b}3$$\mathrm{b}4$2

語の作業領域が必要である。
固有値の精度は、 を大きくしても、 余り改善されない。

.

$\mathrm{L}$$\mathrm{L}$. 固有ベクトルの精度は、$\mathrm{L}=0$. 10$\mathrm{n}$と$\mathrm{L}=\mathrm{n}$の時に改善される場合が多い。
$\mathrm{b}5$.

$\mathrm{n}=100_{\text{、}}$

ベストな$\mathrm{L}$$\mathrm{n}=625$$\mathrm{b}6$.

$\mathrm{F}$$\mathrm{F}$$\tau$

AA

は、
なら$\mathrm{C}\mathrm{O}\mathrm{M}$$\mathrm{n}$にも依存する。

$\mathrm{L}=0$–V$\mathrm{P}$.

$03\mathrm{n}$–20$625$ でのテスト結果によれば、$\mathrm{n}=100$

ならにすれば良い 4). 20$\mathrm{L}=0$. 20$\mathrm{n}$に、$\text{。}$$0_{\text{、}}\mathrm{N}\mathrm{E}\mathrm{C}-\mathrm{S}\mathrm{X}2_{\text{、}}\mathrm{H}$$\mathrm{C}\mathrm{O}\mathrm{M}-\mathrm{v}\mathrm{p}-2600_{\text{、}}\mathrm{N}\mathrm{E}\mathrm{C}-\mathrm{S}\mathrm{X}3$ば、 ベストな$\mathrm{L}=0$$265_{\text{、}}400_{\backslash }$I$\mathrm{T}$A$\mathrm{C}-\mathrm{S}820$

でのテスト結果によれ$\mathrm{L}$は、 ベクトル加速率にも依存する。 V $\mathrm{P}-200$ なら$\mathrm{n}$に、$\mathrm{S}\mathrm{X}3$なら$\mathrm{L}=0$.

$03\mathrm{n}$にすれば良い 4)$\circ$30$\mathrm{b}7$. 収束に必要な平面回転の数や$\mathrm{L}$$\mathrm{Q}\mathrm{R}$に関して振動する場合があるが、変換の数は、$\mathrm{L}\propto \mathrm{n}$とすれば$\mathrm{L}\propto \mathrm{k}$$\mathrm{L}$とするととともに漸減する。

34. 固有値のみを求める場合の性能評価
表 1 には使用コード名が、 表 2 には 900 次元の求めるのに要した$\mathrm{C}\mathrm{P}\mathrm{U}$機で求めるのに要した$\mathrm{T}$のすべての固有値をスカラー機で時間が、 表 3 には 900 次元の$\mathrm{C}\mathrm{P}\mathrm{U}$$\mathrm{T}$のすべての固有値をベクトル時間が、 表 4 には得られた固有値の絶対誤差が、示されている。 これらのデータに基づいて、 以下のごとく結論した。

(a). 速度に関しては、 NUMPAC の中の無平方根(b). 精度に関しては、. 分割併合ホモトピー法の$\mathrm{Q}\mathrm{R}$Symmetric Tridiagonal)5.
(C).

$\mathrm{Q}\mathrm{R}$

法は、$\mathrm{n}$$\cdot$6)

コードが最高速である。
(Homotopy for Real$\mathrm{H}\mathrm{R}\mathrm{S}\mathrm{T}$が最高精度である。

が 1000 を超えるとデフレーションエラーのために精度が低下する。 だからこそ、 Parlett は [The$\mathrm{Q}\mathrm{R}$algorithm is the mosteffective way of finding all the eigenvalues of a small matrix.
述べている。

algorithm

code

coderyear

二分割法

Lapack多分割法

DSTEBZHOBSVN

Ninomiya

1992
1984

平方根付 QRHOQRVW

Ninomiya

1984

無平方根 QR

Lapack無平方根 QRDSTERFNSHOUD

Beppu

1992
1982

分割併合ホモトピー法

HRST

K. Li, T. Y. Li & Z. Zeng1994

表 1. テストに使用した固有値コード。

$\lrcorner \text{と}$

31

code

$\mathrm{F}$DSTEBZDSPNGDSTERFNSHOUD
HOQRVNHOBSVN$\mathrm{B}$$\mathrm{H}$$\mathrm{R}$$\mathrm{C}$$\mathrm{J}$$\mathrm{G}$6236 8179 8469
3207 1042 3072
1277 790 1172
420 313 484
864 553 887
7915 7919 7919

6112 1807 9548 6681 8178
2492 4294 3130 3177 1458
951 947 1355 1197 1040
294 436 519 397 346
466 753 923 702 619
7851 7911 7903 7908 7923

表 2. Cpu-time (milli sec. ) for finding all the eigenvalues when
was used on March, 1994.

$\mathrm{n}--900$$\mathrm{F}\mathrm{A}\mathrm{C}\mathrm{O}\mathrm{M}-\mathrm{N}-18\mathrm{o}\mathrm{o}$code$\mathrm{B}$$\mathrm{F}$$\mathrm{H}$$\mathrm{R}$$\mathrm{G}$$\mathrm{C}$$\mathrm{J}$$\mathrm{L}$6089 1797 9500 6637 8131 6199 8133 8419
2398 4141 3019 3061 1407 3088 1004 2951
485 478 689 607 531 651 401 596
290 427 512 391 340 414 308 476
462 745 911 693 611 855 546 876
132 132 131 131 131 131 131 132

DSTEBZDSPMGDSTERFNSHOUD
HOQRVNHOBSVN

表 3. Cpu-time (milli sec. ) for finding all the eigenvalues when$\mathrm{F}\mathrm{A}\mathrm{C}\mathrm{O}\mathrm{M}^{-}\mathrm{V}\mathrm{P}-26\mathrm{o}\mathrm{o}$code

DSTEBZDSPNGDSTERFNSHOUD
HOQRVNHOBSVN$\mathrm{F}$2. lD-ll
1.
4.
2.
3.
1.

$3\mathrm{D}-14$$3\mathrm{D}-13$$9\mathrm{D}-07$$2\mathrm{D}-12$$6\mathrm{D}-07$

was used on March, 1994.
$\mathrm{R}$$\mathrm{B}$1.
1. lD-14
1. OD-131.
1.
2.
$8\mathrm{D}-11$$2\mathrm{D}-13$$5\mathrm{D}-13$$5\mathrm{D}-10$$\mathrm{n}^{--}900$$\mathrm{H}$$\mathrm{L}$$\mathrm{J}$1.
1.
2OD-08 1 OD-13 1.
3.
7.
5.
2. lD-ll
5. lD-08 1.
5.
2.
3.
5.
2.

$8\mathrm{D}-10$$6\mathrm{D}-11$

1.
1.

$3\mathrm{D}-11$$3\mathrm{D}-11$$9\mathrm{D}-16$$9\mathrm{D}-14$

1.
2.
6.
6.
3. lD-ll
1. OD-07$2\mathrm{D}-11$$6\mathrm{D}-14$$2\mathrm{D}-11$$3\mathrm{D}-08$$9\mathrm{D}-13$$6\mathrm{D}-12$$7\mathrm{D}-12$$5\mathrm{D}-12$$4\mathrm{D}-13$$3\mathrm{D}-05$

表 4. Eigenvalue error$\mathrm{F}\mathrm{A}\mathrm{C}\mathrm{O}\mathrm{M}^{-}\mathrm{V}\mathrm{p}-2600$$2\mathrm{D}-10$$7\mathrm{D}-08$

for the case of
was used on March, 1994.
$\triangle\lambda$$\mathrm{n}--9\mathrm{o}\mathrm{o}$.

$\mathrm{W}$1.
3. lD-14
1. ID-II6.
5.
1.
$2\mathrm{D}-11$$5\mathrm{D}-11$$3\mathrm{D}-11$$3\mathrm{D}-08$.

32

95. 固有対を求める場合の性能評価
すべての固有値と固有ベクトルを求めるために、 表 5 に示した固有値コードを使用し$\text{、}$2)

以下の結論を得た。

(a). スカラー機で固有対を高速に求めたければ、$\mathrm{Q}\mathrm{R}$. 逆反復法の$\mathrm{Q}\mathrm{R}$

コード 7) を用いるとよい。

(b). ベクトル機で固有対を高速に求めたければ、 中次シフト法の$\mathrm{Q}10$

コードを用いるとよい。

(c). 固有ベクトルを精度よく求めたければ、$\mathrm{H}\mathrm{O}\mathrm{M}\mathrm{O}$$\mathrm{H}\mathrm{R}\mathrm{S}\mathrm{T}$II コードまたは

コード 8) を用いるとよい。

(d). ホモトピーコードの高速性は確認できなかった。

year

codecoder

平方根付 QR(D2 shift)

D2

Ninomiya

平方根付 QR(Q10 shift)

Q10

Beppu(D2 shift)

TQL2Eispack平方根付 QR/QR (D2 shift)

DSTEV

Lapack1984
1990
1976 (1992)
1992

分割統治 9)TREEQR

Dongarra

1987

QRII

HRSTII

Beppu-Ninomiya
Li, Li &Zeng

1982
1994

HOMO

Li &Rhee1989

algorithm

平方根付無平方根$\mathrm{Q}\mathrm{L}$$\mathrm{Q}\mathrm{R}+$分割併合
ホモトピー$+$

逆反復

逆反復

表 5. テストに使用した固有値コード。

II33

96. おわりに.

今後も固有値コードの性能評価は機会有るごとになされるであろうが、コードを含めての性能評価が望まれる。
II コード.

コード$\mathrm{Q}10$$\mathrm{Q}\mathrm{R}$謝$\mathrm{H}\mathrm{R}\mathrm{S}\mathrm{T}$辞: 電算機を使用させていただいた、 名大京大東大分子研の大型計算機セン

ターに深謝いたします。

参考文献:$1)\mathrm{B}$$2)\mathrm{Y}$$3)\mathrm{H}$. Parlett: The Symmetric Eigenvalue Problem (1980) Prentice-Hall.
. Beppu, “Origin-Shift Straregies in the QR method”, PCG94 (1994)251.
. Isaka, Y. Beppu and K. Takeuchi: 法の原点移動方法の比較、 数値解析$\mathrm{Q}\mathrm{R}$シンポジウム (1991)$4)\mathrm{S}$. Sasaki, H. Isaka and . Yorozu:$\mathrm{J}$$\mathrm{S}\mathrm{X}-3$

における$\mathrm{Q}\mathrm{R}$法の原点移動方法の比較、 情報処理学会第 43 回全国大会 (1991)
5). Li & Z. Zeng: Laguerre’

Iteration in Solving the Symmetric
Tridiagonal Eigenploblem , SIAM J. Sci. Comp. 15 (1994)1145.
. Li, “A Fully Parallel Method for Tridiagonal Eigenvalue Problem”,
I . J. Math. Math. Sci. , (1994)
st iagonali zation Subroutine”,. Beppu &I.
$\mathrm{T}\mathrm{Y}$$\mathrm{s}$$6)\mathrm{K}$$\mathrm{n}\mathrm{t}$‘$7)\mathrm{Y}$$\mathrm{H}\mathrm{Q}\mathrm{R}\mathrm{I}\mathrm{I}:\mathrm{A}\mathrm{F}\dot{\mathrm{a}}$$\dot{\mathrm{D}}$$\mathrm{N}\mathrm{i}\mathrm{n}\mathrm{o}\mathrm{m}\mathrm{i}\mathrm{y}\mathrm{a},$Computers & Chemistry, Vol.6(1982)87.
and N.Rhee: Homotopy Algorithms for Symmetric Eigenvalue$8)\mathrm{T}$.

$\mathrm{L}\mathrm{i}$Problems, Numer. Nath. 55(1989)265.
$9)\mathrm{J}$. Dongarra and D. Sorensen: A Ful ly Paral lel Algorithm for theSymmetric Eigenvalue Problem, SIAM J. Sci. Stat. Comput. 8(1987)139.