超高速行列演算チップの作成

1996

年度 卒業論文

超高速行列演算チップの開発

9320028

木村・齋藤研 中島 瑞樹

電気通信大学 電子工学科 電子デバイス工学講座

指導教官 齋藤 理一郎 助教授

提出日 平成

9

年

2

月

5

日

概要

Ax = x

を満たす数



とベクトル

x

をそれぞれ固有値・固有ベクトルといい科

学技術計算などの計算において、必要不可欠になっているが現在

,

行列の計算

(

固有

値・固有ベクトル

)

は時間がかかる。そこで計算を早くさせるためソフトウエア又は

ハードウエアによる解決方法があり、ソフトウエアではプログラムのアルゴリズム

や並列計算機を使用した並列化プログラムの作成があるが、本研究ではハードウエ

アによる方法を用いた。そのため、固有値固有ベクトルを求める為だけの専用プロ

セッサーを仮想し

(

これを専用チップと定義する

)

そのプロセッサー上でどのように

したら固有値固有ベクトルを求めることができるかを考えるのが本研究の目的であ

る。そしてその過程が正しいか確かめるために

C

言語を用いプログラムでシミュレー

ションすることにした。これにより固有値固有ベクトルを求める際の過程がわかり、

どの部分が重要であるか、またこの過程は時間がかかるため専用チップ上で行うよ

うにして速く計算させるなど専用チップの制作へスムーズに移行しやすくなる。

ここではまず固有値固有ベクトルを求める方法としてハウスホルダー変換

, 2

分法

,

そして逆反復法について説明し次に専用チップの概要の説明、最後に専用チップ上

でこれらの方法によって固有値・固有ベクトルが求まるよう実現できるかを考える。

その結果専用チップ上で求まることが確認できた。

目次

1

序論

1 2

ハウスホルダー変換による行列の三重対角化

1 2.1

二分法

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 6 2.2

逆反復法

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 8 2.3

連立１次方程式と

LU

分解

: : : : : : : : : : : : : : : : : : : : : : : 10 2.3.1

ガウス消去法と

LU

分解

: : : : : : : : : : : : : : : : : : : : : 10 2.3.2 LU

分解の行列表現

: : : : : : : : : : : : : : : : : : : : : : : 13 2.3.3 LU

分解による連立１次方程式の解

: : : : : : : : : : : : : : : 15 3

専用チップについて

17 3.1

専用チップを用いた計算例

: : : : : : : : : : : : : : : : : : : : : : : 18 3.2

計算時間について

(

実際の数値

) : : : : : : : : : : : : : : : : : : : : : 20 4

過程について

20 4.1

アルファベットの意味

: : : : : : : : : : : : : : : : : : : : : : : : : : 20 4.2

実際の例

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 21 4.3

特殊なもの

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 21 4.4

機能一覧表

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 22 4.4.1

三重対角行列を求める際

(

ハウスホルダー変換

)

の機能。

: : : 22 4.4.2

２分法

(

固有値を求める

)

での機能

: : : : : : : : : : : : : : : 22 4.4.3

逆反復法での機能

: : : : : : : : : : : : : : : : : : : : : : : : 23 5

専用チップでの過程

24 5.1

三重対角化までの過程

: : : : : : : : : : : : : : : : : : : : : : : : : : 24 5.2

二分法の過程

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : 32 5.3

逆反復法の過程

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : 36 6

考察

40

7

結論

41

1

序論

1 1

序論

本研究では行列の固有値

,

固有ベクトルを求める方法としてハウスホルダー変換

, 2

分法

,

そして逆反復法によって求めている。この章ではまずハウスホルダー変換

, 2

分

法

,

そして逆反復法について説明する。

対称行列の固有値を計算するとき，通常は行列を三重対角行列に変換し，その三重

対角行列の固有値を計算する。その理由は，与えられた行列の固有値を一度に求め

るよりも，このように三重対角化を中間におく方が，全体として手間が少なくなる

からである。ここでは，三重対角化の方法としてハウスホルダ一変換，そしてその

固有値を求める方法として二分法を使うプログラムについて説明する。対称行列の

固有ベクトルの計算については次章で述べる。

2

ハウスホルダー変換による行列の三重対角化

n2n

行列

A

に関する固有値問題（

elgenvalueproblem

）

Ax =x (2:1)

を考える。本章では，行列

A

は対称であるとしておく。固有値を求めるために，ま

ず

A

を相似変換

~ A=P 01 AP (2:2)

によって三重対角行列に変換する。

~ A= 0 B B B B B B B B B B B B B B @  1  1  1  2  2 0  2  3 . . . . . . . . . . . . 0 . . .  n01  n01  n01  n 1 C C C C C C C C C C C C C C A (2:3)

三重対角行列（

tridiagonal matrix

）とは，

(2.3)

のように，対角成分およびそれ

に隣接する副対角成分のみが

0

でない行列のことである．相似変換によって固有値

2

ハウスホルダー変換による行列の三重対角化

2

は変わらないから，

A

の固有値の代わりにより簡単な計算によって

~ A = P 01 AP

の

固有値を求める，というのがここで述べる方法の考え方である。

三重対角化を行う全体の変換行列

P

を構成する前に，まず適当な

n

次元ベクトル

w

を使って

Q=I0cww T ;c= 2 w T w (2:4)

なる

n 2n

行列

Q

を定義する。この行列は対称な直交行列であること、すなわち

Q 01 = Q 0T = Q

をみたすことは容易に確かめられる。この型の行列

Q

による相似

変換を，ハウスホルダー（

Householder

）変換という。行列

Q

による

A

のハウスホ

ルダー変換を考え，これを

B

とおく。

B =Q 01 AQ (2:5)

三重対角化の第

1

のステップは，

B

の第

1

列の第

3

行目以下の成分がすべて

0

に

なるようにすることである。そのようにすると，

B

の第

1

行の第

3

列目より右の成

分も自動的に

0

になる。なぜならば，

A

は対称であるから，

(2.5)

より

B

も対称に

なるからである。つまり，

Q

すなわち

w

をうまく選んで，

B

を次の形にもってゆく

わけである。

B = 0 B B B B B B B B B B B @  1  1 0 ::: 0  1 3 3 ::: 3 0 3 3 ::: 3 . . . . . . . . . . . . 0 3 3 ::: 3 1 C C C C C C C C C C C A (2:6) *

は一般には

0

でない成分を表す。

2

ハウスホルダー変換による行列の三重対角化

3

いま，

A

の第

1

列のベクトルを次のように書く。

a = 0 B B B B B B B B B B B @ a 1 a 2 a 3 . . . a n 1 C C C C C C C C C C C A (2:7)

そして，ハウスホルダー変換を定義する

(2.4)

のベクトル

w

を

w = 0 B B B B B B B B B B B @ 0 a 2 +s a 3 . . . a n 1 C C C C C C C C C C C A ; s 2 = n X j=2 a 2 j (2:8)

に選んでみる。

s

の符号は

a 2 +s

の計算で桁落ちの生じないように

a 2

と同符号にと

る。このとき

,(2.4)

の

c

は次のようになる。

c= 1 s 2 +jsa 2 j (2:9)

一方，

w

を

(2.8)

のようにとると，

(2.4)

の

Q

の第

1

列は第

1

成分が

1

である単位

ベクトルとなり，

Q 01 A

の右から

0

を乗ずるときその結果の

Q 01 AQ

の第

1

列は

Q 01 A

の第

1

列と変わらない．したがって，

B = Q 01 AQ

の第

1

列のベクトルを

b

とおく

と，

(I0cww T )a=b (2:10)

が成り立つ。これから直ちに

b = 0 B B B B B B B B B B B @ a 1 0s 0 . . . 0 1 C C C C C C C C C C C A (2:11)

となり，

B

が

(2.8)

の形をもつこと，および

8 > < > :  1 =a 1  =0s (2.12)

2

ハウスホルダー変換による行列の三重対角化

4

となることがわかった。実際に

(2.5)

の変換を計算するときには，次の関係を利用

する。

B =(I0cww T )A(I0cww T )=A0(wq T +qw T ) (2:13)

ただし，最後の式は，まず

P =cAw (2:14)

を計算してから

q =P 0 c 2 wP T w (2:15)

を求め，この

q

を使って計算する。ハウスホルダー変換によって行列が

(2.6)

の形に

なったから，次に

(2.6)

の第

1

行と第

1

列を除いた残りの 小行列の部分について同様

の変換を行う。この操作を

n-2

回くり返せば，もとの行列

A

は最後に

(2.3)

の形をも

つ三重対角行列になるわけである。つまり，第

1

行から第

k-1

行までおよび第

1

列

から第

K-1

列まで三重対角化の終了した行列を

A (k )

とするとき

(2.8)

と同しように

選んだベクトル

w (k)

から行列

Q (k ) = I0c k w (k) w (k)T

を構成し，この

Q (k)

による相

似変換

A (k +1) =(Q (k) ) 01 A (k) Q (k ) (2:16)

を行う。これを

n-2

回くり返して

~ A=A n01 =p 01 AP (2:17) P =Q n02 Q n01

・・・

Q (1) (2:18)

とすればよい。以上の手順をまとめると，次のようになる。ただし，

A (k)

の第

ij

成分を

a i j k

と書く。

A (1) =A

とおき

,k =1;2;3;

・・・

;n02

について次の手順をくり返す。

(i) s k = v u u t n X a (k)2 jk (2:19)

2

ハウスホルダー変換による行列の三重対角化

5

もしも

a (k) k+1;k <0

ならば

s k

の符号を負にする。

もしも

s k =0

ならば

(ii)(iii)

を実行せずに次の

k

へ進む。

(ii) c= 1 s 2 k +a (k) k+1;k 2s k (2:20) w (k) = 0 B B B B B B B B B B B B B B @ 0 . . . a (k) k+1;k +s k a k+2;k (k ) . . . a (k) n;k 1 C C C C C C C C C C C C C C A (2:21) (iii) P (k) =c k A (k) w (k) (2:22) q (k) =P (k ) 0 c k 2 w (k) P (k)T w (k) (2:23) A (k+1) =A k 0(w (k) q (k)T +q (k ) w (k)T ) (2:24)

2

ハウスホルダー変換による行列の三重対角化

6 2.1

二分法

固有値を求める次のステップは、もとの行列

A

の相似変換である

(2.3)

の三重対角

行列

~ A

の固有値を計算することである。この計算方法として、ここではスツルムの

定理に基づく二分法を採用する。

三重対角行列

(2.3)

に対応して行列

I 0 ~ A

を考え、その第

k

主小行列式を

p k ()

とおく。

0 B B B B B B B B B B B @ 0 1 0 1 0 1 0 2 . . . 0 . . . . . . . . . 0 . . . 0 k01 0 k01 0 k 01 0 k 1 C C C C C C C C C C C A (2:25)

これを最後の行について展開すると

p k ()=(0 k )p k01 ()0 k01 2 p k 02 +() (2:26)

が得られる。これは



の多項式列

fp k ()g

に関する漸化式であるが、これが

k = 2

の時にも成立するように便宜的に

p 0 ()=1 (2:27)

と定義しておく。

k =n

のとき

p n ()=jI0 ~ Aj (2:28)

となるが、この

n

次多項式

p n ()

の根が

~ A

の固有値、すなわち

A

の固有値に等しい。

ここでは証明しないが、上の多項式の列

p n ();p n01 ();111;p 1 ()p 0 () (2:29)

の符号の変化の回数に関して、いわゆるスツルム

(Sturm)

の定理から次の事実が導

かれる。すなわち、



を固定して

(2.29)

を左から右へ見ていったときの符号の変化

の回数を

N()

とすると、

N()

は、



より大きい固有値の個数に等しい。

行列

A

および

~ A

は対称であるから、その固有値

 i ;i = 1;2;111

はすべて実数で

ある。固有値は大きい方から

 1 ; 2 ;111

の順に並べるものとする。もしも、二つの数

2

ハウスホルダー変換による行列の三重対角化

7 a j01 ;b j01

に関して

N(a j01 )  k; N(b j01 ) < k;

が成立しているとすると、上に述べ

たことから、大きい方から第

k

番目の固有値

 k

は

a j01 <  k < b j01

の間に存在す

る。この範囲を十分に狭くすることができたならば、

 k

の近似値が求められたこと

になる。区間の幅を半分にしながら、次第に

 k

の存在する範囲を狭くしてゆく方法

が、二分法

(bisection method)

である。

2

ハウスホルダー変換による行列の三重対角化

8 2.2

逆反復法

B

を

n 2n

行列として、その固有値を

 1 ; 2 ;111; n

とする。これらの固有値は絶

対値の大きい方から順にならんでいて、絶対値最大の固有値には縮退はないものと

する。つまり

j 1 j > j 2 j  j 3 j  111  j n j (2:30)

をみたしているものとする。このとき、適当な初期ベクトル

(0)

から出発して、

(l ) =B (l 01) l =1;2111 (2:31)

を計算してゆくと、よく知られているように、

(l )

は次第に

B

の絶対値最大の固有

値

 1

に対応する固有ベクトルに近づいてゆく。これがいわゆるべき乗法

(powermethod)

の原理である。

さて与えられた

n2n

行列

~ A

の

k

番目の固有値を

 k

、それに対応する固有ベクト

ルを

とする。いま、



を

 k

の適当な近似値として、

(I0 ~ A) 01 (2:32)

なる行列を考えると、

(I0 ~ A) 01 = 1 j0 k j (2:33)

が成り立つ。近似値



は

 k

に十分近く、他の固有値

 j ;j 6=k

とは

1 j0 k j > 1 j0 j j j 6=k (2:34)

なる関係になっているとすると、

1=(0 k )

は行列

(I0 ~ A) 01

の絶対値最大の固有

値になる。したがって、適当な初期値

(l )

から出発して

(I 0 ~ A) 01

にべき情報を適

用し、次々に

(l ) =(I 0 ~ A) 01(l 01) (2:35)

を計算してゆくと、

(l )

は次第に固有ベクトル

に近づいてゆくことになる。これ

が逆反復法

(inverse iteration)

の原理である。

実際の計算では、

(2.35)

の反復は連立１次方程式

~ (l) (l 01)

2

ハウスホルダー変換による行列の三重対角化

9

を

(l )

について解くことによって進めることが出来る。



が

 k

に十分近いと

(2.35)

の収束は速くなるが、一方



が

 k

ひじょうに近いと、方程式

(2.36)

の解はいわゆる

不定の状態に近くなる。しかし、

(2.36)

を解いて得た解は固有ベクトルとしてはあ

くまで正しいものになっている。

もしも、固有値

 k

に縮退があると、

(2.35)

の反復によって

(k )

は

 k

に対応する

複数個の固有ベクトルの１次結合に近づく。したがって、縮退がある場合には、

 k

に対応する２番目の固有ベクトルを求める反復からは、すでに求めてある固有ベク

トルと直交するような初期ベクトルをとって反復を開始する必要がある。

方程式

(2.36)

を解くときには、

LU

分解を使うとよい。すなわち、

I0 ~ A

をまず

I 0 ~ A=LU (2:37)

のように

L

と

U

の積に

LU

分解する。このとき、

(2.36)

は、

LU ( ) = (l 01) (2:38)

と書くことができる。この方程式は次のように、

L = (l01) (2.39) U (l ) = (2.40)

に分けて解けばよい。

以上の議論は一般の

n2n

行列

~ A

に対するものであるが、以後、

~ A

として前に得

た三重対角行列を考えることにする。したがって、以下では対称行列

A

の固有ベク

トルを求めることになる。

逆反復法の手順

(2.38)

において

l=1

とおくと

U (l ) =L 01(0) (2:41)

となるが、初期値

(0)

が任意にとれることから実は

L 01(0)

を任意にとれることに

なる。したがって、最初の反復では

(2.40)

の前進代入の部分

L 01(0)

は省略するこ

とができる。

 k

の近似値



が十分に近いとすると２回目の反復

LU (2) = (1)

で

十分精度の高い固有ベクトルを得ることができる。ただし、

 k

に縮退がある場合に

は、縮退度が増すにつれて反復回数を多くする必要がある

2

ハウスホルダー変換による行列の三重対角化

10

これまでの段階で計算された固有ベクトル

k

は、三重対角行列

~ A

の固有ベクトル

である。

~ A

は

A

を相似変換

~ A = P 01 AP

を行って得たものであるから、

~ A

と

A

の

固有値は同じである。そこで、固有値

 k

に対応する

A

の固有ベクトルを

k

とする

と

~ A(P 01 k )= k (P 01 k ) (2:42)

なる関係が成立するから、

k =P 01 k

であることがわかる。したがって、

k

は

k =P k (2:43)

によって正しい固有値

k

へ変換しなければならない。この変換は

(2.4)

と

(2.18)

の

関係を通じて

(k)

および

c k

を使って実行することができる。

2.3

連立１次方程式と

LU

分解

2.3.1

ガウス消去法と

LU

分解

解くべき連立１次方程式は

A = (2:44)

であるとする。

A

は

n2n

行列

A = 0 B B B B B B B @ a 11 a 12 ::: a 1n a 21 a 22 ::: a 2n . . . . . . . . . a n1 a n2 ::: a nn 1 C C C C C C C A (2:45)

で、 および はそれぞれ次のようなベクトルである。

= 0 B B B B B B B @ b 1 b 2 b n 1 C C C C C C C A (2:46) = 0 B B B B B B B @ x 1 x 2 x 1 C C C C C C C A (2:47)

2

ハウスホルダー変換による行列の三重対角化

11

方程式

(2.44)

を最も標準的な方法はガウス消去法であるが、ガウス消去法によっ

て

(2.44)

を直接解く代わりに、係数行列を左下三角行列

L

と右上三角行列

U

の積と

して

A=LU (2:48)

の形に

LU

分解すると便利である。分解のためのアルゴリズムはガウス消去法その

ものであり、この分解を経由して方程式

(2.44)

を解いても、必要な手間は同じであ

る。

そこで、まずガウス消去法（

Gauss elimination

）について説明しよう。係数行列

A

と右辺のベクトル の成分を

a ij (1) =a ij ; b i (1) =b i (2:49)

と書くと、方程式

(2.44)

は次のようになる。

8 > > > > > > > > > > > < > > > > > > > > > > > : a 11 (1) x 1 +a 12 (1) x 2 +111+a 1n (1) x n =b 1 (1) [1] a 21 (1) x 1 +a 22 (1) x 2 +111+a 2n (1) x n =b 2 (1) [2] a 31 (1) x 1 +a 32 (1) x 2 +111+a 3n (1) x n =b 3 (1) [3] 111111111 a n1 (1) x 1 +a n2 (1) x 2 +111+a nn (1) x n =b n (1) [n] (2:50)

上記添字

(1)

は、これから第１段目の消去が行われることを示す。消去の第１段で

は、

(2.50)

の

[2]

以下の式から

x 1

を含む項を消去する。そのために、まず

d 1 01 = 1 a 11 (1) (2:51)

とおいて、

[1]

に

m 21 =a 21 (1) 2d 1 01 (2:52)

を乗じて

[2]

から引くと、

[2]

の

x 1

含む項が消えて

(a 22 (1) 0m 21 a 12 (1) )x 2 +111+(a 2n (1) 0m 21 a 1n (1) )x n =b 2 (1) 0m 21 b 1 (1)

となる。次に、

[1]

に

m =a (1) 2d 01

2

ハウスホルダー変換による行列の三重対角化

12

を乗じて

[3]

から引くと、

[3]

の

x 1

含む項が消える。この操作を繰り返し、

[n]

の

x 1

含む項までを消去すると、

(2.50)

は次の形になる。

8 > > > > > > > > > > > < > > > > > > > > > > > : a 11 (1) x 1 +a 12 (1) x 2 +111+a 1n (1) x n =b 1 (1) [1] 0 a 22 (1) x 2 +111+a 2n (1) x n =b 2 (1) [2] 0 a 32 (1) x 2 +111+a 3n (1) x n =b 3 (1) [3] 0 111111111 a n2 (1) x 2 +111+a nn (1) x n =b n (1) [n] 0 (2.53)

ここで、

a ij (2) ;b i (2)

は次式で定義される。

8 > < > : a ij (2) =a ij (1) 0m i1 a 1j (1) ; i;j =2;3;111;n b i (2) =b i (1) 0m i1 b 1 (1) ; i=2;3;111;n

次に、消去の第２段では、

(2.53)

の

[3] 0

以下の式から

x 2

を含む項を消去する。こ

のための手順は、

[2] 0 ;111;[n] 0

を新たに与えられた方程式と考えれば、第１段におけ

る消去の手順と全く同じである。以下同様にして、

x 3 ;x 4 ;111;x n01

を含む項の消去

を続けることができる。

以上の消去の手順をまとめると、次のようになる。

2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 k =1;2;111;n01 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 d k (01) =1=a kk (k) i=k+1;k+2;111;n 2 6 6 6 6 6 6 6 6 4 m i;k =a ik (k) 2d k 01 j =k+1;k+2;111;n  a ij (k+1) =a ij (k ) 0m ik 2a kj (k) b i (k+1) =b i (k) 0m ik 2b k (k) d n 01 =1=a nn (n) (2.54) m ik

を乗数（

multiplier

）という。

d k

は対角成分

a kk (k )

を意味するが、これは以後の

計算では常に逆数の形でのみ現れるので、あらかじめ

d k 01

を計算しておくわけであ

る。

この操作の結果、もとの方程式は次のように右上三角行列

U

を係数にもつ形に変

形される。

(n)

2

ハウスホルダー変換による行列の三重対角化

13 U = 0 B B B B B B B @ a 11 (1) a 12 (1) ::: a 1n (1) a 22 (2) ::: a 2n (2) . . . . . . 0 a nn (n) 1 C C C C C C C A ;u ij =a ij (i) (ij) (2.56) (n) = 0 B B B B B B B @ b 1 (1) b 2 (2) . . . b n (n) 1 C C C C C C C A (2.57)

ここまでの消去の操作を、ガウス消去法における前進消去（

forwardelimination

）と

いう。

2.3.2 LU

分解の行列表現

ところで、すぐに確かめられるように、消去の第

k

段目の操作

(2.54)

は、単位行

列の第

k

列目の対角線より下に

m ik ;i=k+1;111;n

の符号を変えたものを挿入した

行列

M k = 0 B B B B B B B B B B B B B B B B B B @ 1 1 0 1 0m k +1;k 1 0 0m k +2;k . . . 0 . . . 0m n;k 1 1 C C C C C C C C C C C C C C C C C C A (2:58)

を使って、

A (k+1) =M A (k) (2:59)

2

ハウスホルダー変換による行列の三重対角化

14

と書くことができる。ここで、

A (k)

はこれから第

k

段の消去を行う直前の係数行列

である。

A (k) = 0 B B B B B B B B B B B B B B B B B B @ a 11 (1) a 12 (1) 111 a 1k (1) 111 a 1n (1) a 22 (2) 111 a 2k (2) 111 a 2n (2) . . . . . . . . . a k;k (k) a k;n (k) 0 a k+1;k (k) a k+1;n (k) . . . . . . . . . a n;k (k) 111 a n;n (k) 1 C C C C C C C C C C C C C C C C C C A (2:60)

一方、

A (1) =A;A (n) =U

であるから、

(2.59)

より

U =M n01 M n02 111M 1 A

すなわち

A = (M n01 M n02 111M 1 ) 01 U = M n01 01 M n02 01 111M 1 01 U (2.61)

が成り立つ。ところが、

k = 0 B B B B B B B B B B B B B B B B B B @ 0 . . . 0 m k+1;k m k+2;k . . . m n;k 1 C C C C C C C C C C C C C C C C C C A ::::::k+1 k = 0 B B B B B B B B B B B B B B B B B B @ 0 . . . 0 1 0 . . . 0 1 C C C C C C C C C C C C C C C C C C A ::::::k

なるベクトルを定義すると、行列

M k

は

M k =I0 kk T

と書くことができる。

I

は単位行列で、上付添字

T

は一般に転置を表すものとする。

一方、

k T k =0

であるから

(I0 T )(I+ T )=I0 ( T ) T =I

2

ハウスホルダー変換による行列の三重対角化

15

が成り立ち

M k 01 =I kk T = 0 B B B B B B B B B B B B B B B B B B @ 1 . . . 0 1 0m k+1;k 1 0 0m k+2;k . . . 0 . . . 0m n;k 1 1 C C C C C C C C C C C C C C C C C C A

であることがわかる。これから、直接計算することにより

M 1 01 M 2 01 111M n01 01 = 0 B B B B B B B B B B B @ 1 m 21 1 m 31 m 32 1 . . . . . . . . . m n1 m n2 111 m n;n01 1 1 C C C C C C C C C C C A ; m ij = a ij (j) a jj (j) (2:62)

となることを確かめることができる。すなわち、この行列は左下三角行列である。

したがって、この行列を

L=M 1 01 M 2 01 111M n01 01 (2:63)

とおけば、

(2.61)

によって、もとの行列

A

は

A=LU (2:64)

の形に、左下三角行列

L

と右上三角行列

U

の積に

LU

分解されたことになる。

2.3.3 LU

分解による連立１次方程式の解

行列

A

の

LU

分解を実行しておくと、方程式

(2.44)

を解くことは簡単である。方

程式

(2.44)

は

LU = (2:65)

と書くことができるから、これを解くためには、まず

2

ハウスホルダー変換による行列の三重対角化

16

の解 を求め、次にこの解 を右辺に持つ方程式

U = (2:67)

を解いて を求めればよい。

L

と

U

は共に三角行列であるから、解を求める操作は

きわめて簡単である。

まず、

(2.66)

は、

(2.62)

と

(2.63)

より、

L

の対角成分が１であることに注意すれ

ば、

8 > > > < > > > : y 1 =b 1 y i =b i 0 i01 X j=1 m ij y j ; i=2;3;111;n (2:68)

によって解くことができる。この部分を、前進代入（

forwardsubstitution

）という。

方程式

(2.67)

と

(2.57)

とを比較すれば

= (n)

が成り立ってることがわかる。

次に、

(2.67)

は、

(2.56)

より、後ろから順に

8 > > > > < > > > > : x n = 1 a nn (n) y n x i = 1 a ii (i) 0 @ y i 0 n X j=i+1 a ij (i) x j 1 A ; i=n01;n02;111;1 (2:69)

によって解くことができる。この部分を後進代入（

backwardsubstitution

）という。

この段階で、ガウス消去法は完了したことになる。

係数行列

A

が同じで右辺のベクトル が異なる連立方程式を幾組も解く場合には、

行列

A

をあらかじめ

LU

分解しておくと、各々の に対して

(2.68)

と

(2.69)

を計算

するだけできわめて効率良く解を求めることができる。

3

専用チップについて

17 3

専用チップについて

専用チップの形については図

3.1

に示したものを用いる。

s2

R2

M1

M2

R1

s0

s1

C2

C1

ector

matrix

scalar

_vector

図

3.1:

専用チップの形

1 S 0 ;S 1 ;S 2

にはスカラーが入り

C 1 ;C 2 ;R 1 ;R 2

にはベクトルがそして

M 1 ;M 2

には行列

が入る。

以上これらのレジスターに値をいれて

,

いろいろな過程

(

例えば

C 1

から

C 2

に コ

ピーするなど

)

を駆使してハウスホルダー変換、

2

分法、逆反復法によって固有値・

固有ベクトルを求める。実際の

1

つ

1

つの過程についてはあとの章で述べる。ただ

し、現段階でプログラム上で新しいレジスターを定義して実行しているのもあるの

でその点については今後検討しなければならない。

1

3

専用チップについて

18 3.1

専用チップを用いた計算例

例えば

C =AB

の場合の計算例を以下に示す。

ただし

A= 0 B @ 1 2 3 4 1 C AB = 0 B @ 5 6 7 8 1 C AC = 0 B @ 19 22 43 50 1 C A

として考えます。

例えば

222

行列のかけ算の場合、普通に計算すると以下のようになる。

0 B @ 1 2 3 4 1 C A 0 B @ 5 6 7 8 1 C A= 0 B @ 5+14 6+16 15+28 18+32 1 C A= 0 B @ 19 22 43 50 1 C A

これを専用チップを用いて計算するには以下のような過程で計算する。

4*8*2 bit

data bus

n

A

C

B

,B,C Matrix Register

図

3.2:

専用チップの配置例

2

まず図

3.2

のように専用チップを

2

個用いて掛算を実行する。その結果を下の右側

のチップに入れている。それぞれのチップは相互につながっておりデータを伝送さ

せることができる。

次に

C=AB

の場合の計算を専用チップ上で実現するには以下のことを行う。

2

3

専用チップについて

19

6

7

8

5

7

5

19

50

19

43

22

43

1

3 4

2

7

5

A

C

B

図

3.3:

専用チップの計算例

3 (

注

)

ただし図

3.3

は図

3.1

と少し形が違うが

C 1

の位置が移動しているだけなので問

題はない。

まず

A

と

B

に 行列を入れ結果を

C

に入れる。

1. RegisterA

の

1

列目をとりだします。そして

C1

に格納する。

2. RegisterA

の

C 1

から

Register B

の

R 1

に移動する。

ここからは

RegisterB

での仕事になり、

3. R 1

と

M 1

の行ごとに掛ける。つまりここでは

N 2N

行列の場合普通は

N 2

回

必要になるが専用チップを用いることにより

N

回で済むことになる。これが

このチップの特徴であります。

4. M 1

の各行ごとに足し合わせる。

5.

足し合わせた結果を

C 1

に入れる。

6.

これを

C

の

C 1

に入れる。

(

または

C

の

M 1

にそのまま入れる。

)

そして

1

から繰り返すが 今度は

2

列目からとりだし

1

から繰り返す。

ここでは６つの過程を用いたがそれぞれ１つ１つの過程を機能を定義する。つまり

機能を組み合わせることによりハウスホルダー変換などが実現でき固有値固有ベク

トルが求まることとなる。

3

4

過程について

20 3.2

計算時間について

(

実際の数値

) 

行列の足し算の場合

この場合は行ベクトルを格納するレジスター同士がそれぞれ１成分どうしつな

がっていれば １回の動作で実行できる。



行列のかけ算の場合

普通

n2n

行列の場合、

n 3

に比例する時間がかかる。つまり

10000210000

の

行列なら

10 12

かかることになるがこの専用チップでは

n 2

に比例する時間がか

かり、大行列ほど時間の短縮になる。

10000210000

の行列なら

10 8

かかリ

, 1 10000

の時間の短縮となる。

以上のようにこの専用チップはかけ算が速く計算できるので過程の中にかけ算

が多いと計算速度向上が期待できる。

4

過程について

ここでは実際に固有値固有ベクトルを求める際の過程について述べる。次の機能

を用いて実現している。

プログラムと関数の名前の見方

(

過程のファイル

) 4.1

アルファベットの意味



始めのアルファベットの意味

コピー

,

移動

c : copy

、

m : move

ただし

move

も

copy

として実行している。

足し算

a: add

かけ算

t: times 

次の

m 1 (

レジスター

)

などの意味

次の

m 1

や

r 1 ,s 0

などは レジスタの位置を表し、足すものなどががそのあとに

くる。

4

過程について

21

又、格納場所は

(

アンダーバー

)

の後ろの場所になる。

最後の括弧は場所

(

何行または何列

)

を表す。

4.2

実際の例

例えば

(1) am1 m2

は

am1 m2

は

m 1

と

m 2

を足し 結果を

m 2

に格納する。

(2) mc2i s1(i ) i

は行

これは

C 2 i

を

S 1

に入れるのだが

(i )

は

C 2

の

i

行目の成分を

S 1

に入れるというこ

とを意味している。

(3)mc1i m1 ( i,k,l )

mc1i m1 ( i,k,l)

は

c1i

から

m

1

に移動するがそのあとの

( i,k,l)

の意味は

c 1

の

i

行目を

m 1

の

k

行

l

列に入れるという意味である。

ちなみに

i;k

は

i

行

,k

行を意味し

j;l

は

j

列

,l

列を表す。

(i

行

;j

列 または

k

行

;l

列

と用いることが多い。

) 4.3

特殊なもの

特殊な操作は、操作自体の名前になる。

(

例

)

平方根 は

sqrt

逆数を求めるときは

gyaku

値を入れるときは

Input M 1

の対角成分を

C 1

に入れる。

C 2

から

M 2

の副対角成分に入れる。

などがある。

(

細かくは実際の過程に述べる。

)

4

過程について

22 4.4

機能一覧表

また

i;j

などは値を入れる場所

(

行や列など

)

で、次の括弧の中の数値はプログラ

ムの順番です

4.4.1

三重対角行列を求める際

(

ハウスホルダー変換

)

の機能。

mm1 c1i ( j ) (1) j

は列

ms0 c1i ( i )(3) i

は行

mc1i c2i (4) tc2 c1(5) ac1i s0(6) (13) sqrt s0(7) mc2i s1(i )(8) i

は行

as1 s0 (9) ms0 c2i( i )(10) mc1i c2i (11-1) mc2i r1j (11-2) tc2i c1i (12) ts1 s0 (16) ts0 r1j (17) tr1j m1 (18) (38) am1 c1i( i )(19) am1 c1i 2(i )(30) clr m1 (20) mc1i m1( j ) (21) (29) (31)(37)(51) ts0 c1i (22) (36) ac1i s0 2(24) mc2i c1i (11-1)(27) (32) (50) mc1i r1j (33) (40) mm1 c1i( j ) (34) am1 m2 (35) mm2 m1 (end) mm2 m3 mm3 m2 mr1j m2 Q

の計算

tc1i m2 ms0 c1i( i ) ts2 m2 ar2j m2 (i ) 4.4.2

２分法

(

固有値を求める

)

での機能

mm1 m2( )(0) ftrace c2i( ) (1) tc2i c2i( ) (2) ttei c2i( ) (3) max c1i () (6) Trace m2( ) (7) ms0 r3j ( ) ms0 r4j ( )

4

過程について

23 mr4j r5j( ) ts1 r5j ( ) mr3j r4j( ) ttei r4j( ) ttei r5j ( )(8-0) mc1i m1 ( i,k,l) (8) mm1 r1j( k,l,j ) (9) gyaku r1j( )(11) tc2i r1j( i,j )(12) ar1j m1 (i,k,l ) (13) ar5j m1 (j,k,l) (14) mm1 m1 (j,l) (15) count (l ) (16) compare () (17) mr5j r6j ( j ) last 4.4.3

逆反復法での機能

ms0 r1j ( )(2) tr1j m1 () (3) and (19) Trace m1 ( ) (4) as0 c1i ( ) (5-2) gTrace m1 ( )(6) mc2i m1 ( l) (7-0) mm1 c1i ( l) (8) mm1 r1j( k ) (9) ts0 c1i ( ) (10) mr1j s0 ( j )(11) gyaku s0( ) (12) mr1j m1 ( i ) (14) ms0 c1i ( i )(16) mm1 m2 ( )(17) and (21) mc1i m1 ( l) (18) am2 m1 ( )(20) mm1 c2i ( j ) (22) mc2i c1i ( )(23) tc1i c2i ( ) (24) ac2i s2( ) (25) gyaku s2 () (26) sqrt s2( )(27)

5

専用チップでの過程

24 5

専用チップでの過程

ここからは実際に専用チップ上で固有値固有ベクトルを求める際の過程を説明す

る。

左側には途中過程を専用チップ上で説明しており、右側には機能を表す。なお右の

機能の意味

(

例えば

CM 1 (k

列

)C 1 )

などは

4.2

章で述べたのと同じである。

5.1

三重対角化までの過程

A

M1,M2

に行列

A

を入れる。

0

0

x

x

_nk

kk

α

α

β

β

0

0

0

1

nn

1

1

2

0 x

_nk

x

0

x

_kk

x

_kn

0

0

(2)

(1)

β

_k-1

(1)CM 1 (k

列

)C 1 M 1

上の

k

列目を

C 1

に入れる。

(2)IS 0 (0) S 0

に

0

を入れる。

5

専用チップでの過程

25

0

α

α

β

β

0

0

0

1

nn

1

1

2

0 x

_nk

x

0

x

_kk

x

_kn

0

0

0

x

_{k+1 k}

x

_nk

(3)

(3)CS 0 C 1 (i

行

) (C 1

の

i

行目成分を

0

にする

)

0

α

α

β

β

0

0

0

1

nn

1

1

2

0 x

_nk

x

0

x

_kk

x

_kn

0

0

0

x

_{k+1 k}

x

_nk

0

0

x

_{k+1 k}

x

_nk

(4)

(4)CC 1 C 2 C 1

成分を

C 2

へ

copy

0

α

α

β

β

0

0

0

1

nn

1

1

2

0 x

_nk

x

0

x

_kk

x

_kn

0

0

x

_{k+1 k}

x

_nk

0

0

x

_{k+1 k}

x

_nk

(5)

(6)

s

2

(5)TC 2 C 1 C 2

成分と

C 1

成分を掛け合わせ

C 1

に格納する。

(6)AC 1 S 0 C 1

上の成分を足し 結果を

S 0

へ入

れる。

(S 2

ができる

)

5

専用チップでの過程

26

0

α

α

β

β

0

0

0

1

nn

1

1

2

0 x

_nk

x

0

x

_kk

x

_kn

0

0

x

_{k+1 k}

x

_nk

0

x

_{k+1 k}

x

_nk

s

0

x

(8)

(9)

(7)

(7)SQR TS 0 S 0

の平方根をとる。

(8)CC 2 (i+1)S 1 C 2

の

i+1

行目成分を

S 1

へ入れ

る。

(9)AS 1 S 0 S 1

と

S 0

の 足 し 算 を

S 0

上 で 行 な

う。

α

α

β

β

0

0

0

1

nn

1

1

2

0 x

_nk

x

0

x

_kk

x

_kn

0

0

x

_nk

0

x

+s

k+1 k

x

x

k+1 k

+s

x

_{k+2 k}

(10)

0

0

nk

x

x

_{k+2 k}

x

k+1 k

+s

(11)

(12)

0

0

k+1 k

x +s

x

nk

(10)CS 0 C 2 (i+1) W

の完成

(11-1)CC 2 C 1 (11-2)CC 2 R 1 W

を

C 1 R 1

へ入れる。

w

w

w

w

w

w

(11-4)

M2

(11-3)

(11-3)CR 1 M 2 (11-4)TC 1 M 2

これは逆反復法で用いるため

5

専用チップでの過程

27

α

α

β

β

0

0

0

1

nn

1

1

2

0 x

_nk

x

0

x

_kk

x

_kn

0

w

w

w w

(k)

(k)

T

w w

(12)

(13)

T

x

(12)TC 2 C 1

これより

W T W

の計算

(13)AC 1 S 0 W T W

の計算。

w

α

α

β

β

0

0

0

1

nn

1

1

2

0 x

_nk

x

0

x

_kk

x

_kn

0

w w

T

w

(k)

(k)

2

2

(14)

₍₁₅₎

(16)

(14)IS 1 (2) S 1

に

2

を入れる。

(15)GS 0 S 0

の逆数をとる。

(16)TS 1 S 0 (12)

∼

(16)

で

c= 2 w T w

の計算

c

c

c

c

c

c

c

c

c

c

-1

1

2

n

n

n

1

2

2

1

(16-1)

(16-3)

(16-2)

これもあとで逆反復法で用いる。

(16-1)CS 0 C 1 (16-2)TC 1 M 2 (16-3)InputS 2 (01)

5

専用チップでの過程

28

c

c

c

c

c

c

c

c

c

c

-1

1

2

n

n

n

1

2

2

1

(16-4)

これもあとで逆反復法で用いる。

(16-4)TS 2 M 2 (16-5)InputR 2 (1) (16-6)AR 2 M 2

w

α

α

β

β

0

0

0

1

nn

1

1

2

0 x

_nk

x

0

x

_kk

x

_kn

0

w

(k)

c

(k)

(17)

(19)

(18)

c

p

(17)TS 0 R 1 CW

の計算

(18)TR 1 M 1 (19)AM 1 C 1 P =cAaw

の計算終了

w

(k)

c

p

p

(22)

(21)

(20)CLR M 1 M 1

をすべて

0

にする。

(21)CC 1 M 1 (1

列目

) P

を

M 1

上の任意の位置へ。

(22)TS 0 C 1 CP T

の計算

5

専用チップでの過程

29

c

p

w

(k)

c

p

(23)

(24)

(23)TC 2 C 1 CP T w

の計算

(

実はこれはスカラー

です。

) (24)AC 1 S 0 CP T w

の計算終了

w

(k)

- 1

2

cp w

T

p

(26)

(25)

(27)

(25) IS 1 (0 1 2 ) (26) TS 1 S 0 0 1 2 cp T w

をつくる。

(27) CC 2 C 1

w

(k)

cp w

T

w

(k)

1

2

-1

2

-p

(28)

(29)

cwp w

T

(28)TS 0 C 1 (29)CC 1 M 1 (2)

固定する。

5

専用チップでの過程

30

w

(k)

1

2

-

cwp w

T

p

q

(30)

(30)AM 1 C 1 q = p0 c 2 wP T w

の計

算終了。

w

(k)

w

(k)

(34)

(32)

(31)

(33)

q

q

(30-2)CLR M 1 (31)CC 1 M 1 (1

列目

) q

を

M 1

上へ入れる

(32)CC 2 C 1 (33)CC 1 R 1 (34)CM 1 (1)C 1

w

(k)

w

(k)

(38)

-1

- q

(35)

(37)

(36)

(40)

(35)IS 0 (01) (36)TS 0 C 1 -q

できる

(37)CC 1 M 1 (38)TR 1 M 1 0qw T

の計算終了。

(39)AM 1 M 2 A 0 qw T

の計算終

了。

5

専用チップでの過程

31

w

(k)

w

(k)

(36)

-q

(43)

(41)

₍₄₂₎

(40)

(44)

(40)CC 1 R 1 (41)CC 2 C 1 (42)CC 1 M 1 (all) (43)TR 1 M 1 0wq T

の計算終了。

(44)AM 1 M 2

以上を

N 02

回繰り返すことにより

M 2

の行列が対角化され三重対角化は完成され

る。

4 4 /home2/students/mizuki/xg/no1.eps

∼

no18.eps

5

専用チップでの過程

32 5.2

二分法の過程

ここで

k

は

1

から

25

まで動き

kk

は １から

N

まで動く。なぜ

k

が

25

なのかとい

うと ２分法では１回ごとに幅が

1 2

ずつ狭くなるのでこの回数が精度に影響する。い

まは

25

回としているが

50

回以上

(

できれば

64

回くらい

)

がよい。

kk

は

N

回でこれは行列の行数

(

または列数

)

である。

~

A

(0)M 1 ;M 2

に三重対角化された

~ A

が入る。

α

₁

α

_n

β

n-1

α

_n-1

β

_n-1

1

β

α

2

β

2

β

₂

α

₃

1

β

β

₂

β

n-1

1

β

0

2

2

2

-0

0

(1)

(2)(3)

(4)

(1)(FM 1 )C 2

副対角成分を

C 2

に

入れる。

(2)TC 2 C 2 (3)T(01)C 2 01

を

C 2

に掛け結果

を

C 2

に入れる。

(4)ZM 1 M 1

の絶対値をとる。

5

専用チップでの過程

33

α

₁

α

_n

β

_n-1

α

_n-1

β

n-1

1

β

α

2

β

₂

β

₂

α

₃

1

β

0

0

0

0

(5)

r

r

max

(6)

(5)AM 1 C 1 (6-1)Max(C 1 )S 0 C 1

の最大値をと

り

S 0

に入れる。

α

₁

α

2

α

₃

α

n

α

_n-1

β

_n-1

β

_n-1

1

β

α

n

α

_n-1

α

2

β

₂

β

₂

α

₃

1

β

α

₁

0

0

(6-2)

(6-2)(TM 1 )C 1

ここ

(7)

からは仮に

R 3 ;R 4 ;R 5

を

用いている。これは

2

分法の上限

下限を表している。

(7-1)CS 0 R 4 (7-2)InputS 1 (01) (7-3)TS 1 R 4 (7-4)MS 0 R 3 (7-5)AR 3 R 4

α

₁

α

2

α

₃

α

_n

α

_n-1

β

_n-1

β

_n-1

1

β

α

_n

α

_n-1

α

2

β

₂

β

₂

α

₃

1

β

α

₁

0

0

(7-6)MR 4 R 5 (7-7)InputS 1 (0:5) (7-8)TS 1 R 5 (7-9)MR 3 R 4 (7-10)T(01)R 4 012R 4

の計算

以 上 で は じ め の 範 囲

(0r

か ら

r)

が確定する。

5

専用チップでの過程

34

α

₁

α

2

α

₃

α

_n

α

_n-1

β

_n-1

β

n-1

1

β

α

_n

α

_n-1

α

2

β

₂

β

₂

α

₃

1

β

α

₁

0

0

(8-0)T(01)R 5 (all) (8-1)MC 1 i(1)M 1 (1;k)k

は移動す

る。

(8-2)AR 5 (1)M 1 (1;kk) kk

は移動

する。

(9)MM 1 (k01;kk)R 1 (kk) (10)MC 1 (k)M 1 (k;kk) (11)GR 1 (kk)R 1

の逆数をとる

(12)TC 2 (k-1)R 1 (kk)kk

は列です

(13)AR 5 (1)M 1 (1;kk) kk

は列で 移動する。

(14)MM 1 (k-1,kk)R 1 (kk)M 1

の

(k-1

行

kk

列

)

の値を

R 1

の

kk

列に入れる。

なお

(9)

∼

(14)

までは

N

回繰り返す。

(

これは

1

列目から

N

列目までといること

である。

) (15)MC 1 (k)M 1 (k;kk) (16)Count

これは

M 1

のそれぞれの列ベクトルの値が正である数を数えて結果を

R 1

に入れる。

(17)Compare R 1

と

R 2

の比較をする。

R 2

には左から

1

、

2

、

3

、と順に値が入って

いる。この結果

r1j[j]>= r2j[j]

なら

r4j[j]=r5j[j]

そして

r5j[j]= r4[j]+r3[j] 2

となり、

そうでないなら

r3j[j]=r5j[j]

そして

r5j[j]= r4j[j]+r3j[j] 2

となる。

5

専用チップでの過程

35

このあと

(18)MR 5 R 6 (kk) kk

は移動する。

(19)InputR 4 (0) (20)MR 6 R 3 (21) R 3 (kk+1)+R 4 (kk+1) 2

を

R 5 (kk+1)

へ

(kk

は列を表す。

)

このあと

(8-0)

へいき

kk

回繰り返す。

そしてそれぞれの列で計算するので

N

回繰り返す。

以上の過程で

R 6

に固有値が入る。

(

ただし

2

分法は

R 6

まであるのでさらに改良する可能性がある。

) 5 5 /home2/students/mizuki/xg2/no1.eps

∼

no6.eps

5

専用チップでの過程

36 5.3

逆反復法の過程

µ

(0)M 1

に行列を入れる。この行列

を

A

とする。

µ

-1

-1

-1

(2)

(1)

(3)

-A

(1)IS 0 (01) (2)CS 0 R 1 (3)TR 1 M 1 M 1

には

0A

がはいる

nn

-a

-a

-a

-a

11

22

33

-a

-a

-a

-a

11

22

33

nn

µ

-1

-1

(4)

µ

(5-1)

(4)MTM 1 C 1 M 1

の対角成分を

C 1

に入れる。

(5-1)CS 1 S 0

5

専用チップでの過程

37

11

22

-a

µ

-a

µ

33

-a

nn

-a

µ

µ

µ

-a

₂₂

33

-a

µ

nn

-a

µ

µ

-a

µ

₁₁

(6)

(5-2)

(5-2)AS 0 C 1 (6)C 1

を

M 1

の対角成分に入れる。

11

X

n1

X

1n

X

nn

X

1

x

n

x

2

x

N

N+1

ここで

X =I 0A

とおき成分を

X 11

で表す。

また

x= 0 B B B B B B B @ x 1 x 2 . . . x n 1 C C C C C C C A

を任意の初期ベクトルとする。

nn

X

11

X

n1

X

x

_n

2

x

1n

X

x

1

X

11

N

N+1

-X

₁₁

-X

_n1

-1

X

1n

x

1

X

₁₁

(10)

(8)

(9)

(7)

(7)IS 0 (01) (8)CM 1 (j

列

)C 1 (9)CM 1 (i

行

)R 1 (10)TS 0 C 1

5

専用チップでの過程

38

nn

X

x

n

2

x

1n

X

x

1

X

11

n1

X

N

N+1

-X

₁₁

-X

_n1

X

1n

x

1

X

11

1

X

11

(12)

(11)

(13)

(11)CR 1 (j

列

)S 0 (12)GS 0

nn

X

x

n

11

1

_X

X

12

11

x

1

X

X

1n

X

11

2

x

11

1

_X

X

12

11

x

1

X

X

1n

X

₁₁

n1

X

-X

_n1

1

X

11

X

22

X

21

0

X

₂₁

(14)

(17) M

₂

(15)

(16)

(13)

(18)

(13)TS 0 R 1 (14)CR 1 M 1 (j

列

) (15)IS 0 (0) (16)CS 0 C 1 (i

行

) (17)CM 1 M 2

11

1

_X

X

12

11

x

1

X

X

1n

X

11

X

21

-

X

21

-

X

X

12

11

X

21

-X

11

x

1

-X

_n1

0

X

₂₁

(18)

(19)

0

0

0

0

(18)CC 1 M 1 (19)TR 1 M 1 (20)AM 2 M 1

5

専用チップでの過程

39

U

_k

1

1

1

1

0

0

(7)

∼

(21)

を

N

回繰り返す。この

結果

n+1

列目にできるベクトル

U k

が固有ベクトルとなる。

結果

n + 1

列目にできるベクトル

U k

が固有ベクトルとなり、これをとりだせば良

い。

6 6

以上の図は

/home2/students/mizuki/xg3/no1.eps

から

no10.eps

6

考察

40 6

考察

前章までで述べたように多くの機能を現段階では用いねばならないが実現できる

ことがわかった。しかしまだ問題点が多くあり改良しなければならない。

まずシミュレーションプログラムについてはさらなる検討の余地がある。特に無

駄な部分を取りシミュレーションプログラム自体の時間の短縮が必要である。

実際に計算させたところハウスホルダー変換では

1282 128

行列で

7.12

秒なのに

対し 逆行列での計算

(

逆反復法

)

がかなりの時間

(1282128

行列で

284

秒

)

を要し、

改良しなければならない。

(

しかしこれはファイルアクセスの時間も含んでいる。

)

また

2

分法でのプログラムはレジスターの数が多く定義

(R 1

から

R 6

の

6

つ

)

して

いる。これは途中結果を１回１回格納しなければならないためだか、なるべくレジ

スターを有効に利用し少なくなるようにしなければならない。あまりに多過ぎると

チップをつくるのが困難になってしまう可能性がある。

7

結論

41 7

結論

この専用チップを用いた結果、現段階では多くの過程が必要だが専用チップを用

いて固有値固有ベクトルを求めることが可能である。

そしてチップのさらなる効率化を考えハード化することにより固有値固有ベクト

ルを求める際専用チップを用い計算時間が短縮すると考えられる。

7

結論

42

謝辞

本研究及び論文作成に当たり、終始御懇切なる御指導、御鞭撻を賜わりました指

導教官である齋藤理一郎助教授に衷心より御礼の言葉を申しあげます。

また、本研究を進めるにあたり、熱心な御指導をいただくとともに種々の御高配を

賜わりました木村忠正教授、湯郷成美助教授に深謝の意を表します。

また、研究活動をともにし、多くの援助をいただいた八木将志氏に深謝いたします。

そして、数々の御援助、御助言をしていただいた中平政男氏、竹谷隆夫氏、はじめ

木村

1

齋藤研究室の大学院生、卒研生の方々、高田さんをはじめ

(

株

)

画像技研の方々

に感謝します。

参考文献

43

参考文献

[1]

数値計算の手順と実際 高田勝 春海佳三朗 コロナ社

[2] FORTRAN77

数値計算プログラミング 森正武 岩波書店

[3] C

言語によるプログラミング

(

基礎編・応用編

)

内田智史 オーム社

[4] C

言語による最新アルゴリズム辞典 奥村春彦 技術評論社

A

付録・シミュレーションプログラム

44 A

付録・シミュレーションプログラム

このプログラムは

hhc.c , 2bun.c ,gyaku.c

の

3

つのプログラムからな

り

,hhc.c

はハウスホルダー変換の実行

2bun.c

は

2

分法の実行

gyaku.c

は逆反復法を実行する。

なお

hhc.c

の実行ファイルを実行すると

al.dat ,be.dat

というファイルができる。

これは三重対角行列の対角成分

(al.dat)

副対角成分

(be.dat)

が入っている。

2bun.c

の実行ファイルを実行すると

al.dat

と

be.dat

を読み込み

koyuu.dat

という固有値がはいったファイルができる。

最後に

gyaku.c

の実行ファイルを実行すると

koyuu.dat

を読んで

koyuuvec.dat

という固有ベクトルが入ったファイルができる。

/*---ハウスホルダー変換プログラム

hhc.c 1997 1 23 Mizuki Nakajima ---*/ /*---include ---*/ #include<stdio.h> #include<math.h> #include<time.h> #define N 4 /*

行列

N

行

N

列 の配列を定義する

*/ #include"register.h" /*

専用チップのレジスタの定義

*/ /* #define DBGPRN */ /*

途中経過を出力するか？

#define DBGPRN

をコメントアウトすると出力しない

*/ /*---*/ int i,j ; /* i

行

j

列

*/ int k ; int s = 1; /* s = start */

/*---プロトタイプ宣言

---*/

void input ( void );/*

この関数を用いる

*/ void input_m2 ( void );

void input_test ( void ); /* test */

/*---以下は

input_c1i

などの関数は

c1i

に値

0

を入れる関数

print_m1

などの関数は

m1

に値

0

を入れる関数

_

の後には格納場所が入る

これは後の２分法、逆反復法も同様である

---*/

void input_c1i( void );

void input_c2i( void );

void input_r1j( void );

void input_r2j( void );

void print_m1( void );

void print_m2( void );

void print_c1i( void );

void print_c2i( void );

void print_r1j( void );

void print_r2j( void );

/*---以下は

ハウスホルダー変換により三重対角化するため

一つの過程を一つの関数として表示し

関数を組み合わせてハウスホルダー変換を実現している

---*

void m1_c1i ( int j ); /* (1) j

は列

*/ void s0_c1i ( int i ); /* (3) i

は行

*/ void c1i_c2i( void ); /* (4) */

void tc2_c1 ( void ); /* (5) */

double ac1i_s0 ( void );/* (6) (13) */

void sqrt_s0( void ); /* (7) */

void mc2i_s1 ( int i ); /* (8) i

は行

*/ void as1_s0 ( void ); /* (9) */

void ms0_c2i( int i ); /* (10) */

void mc1i_c2i( void ); /* (11-1) */

void mc2i_r1j( void ); /* (11-2) */

void tc2i_c1i( void ); /* (12) */

void ts1_s0 ( void ); /* (16) */

void ts0_r1j ( void ); /* (17) */

void tr1j_m1 ( void ); /* (18) (38) */

double am1_c1i ( int i ); /* (19) */

double am1_c1i_2 ( int i ); /* (30) */

void clr_m1 ( void ); /* (20) */

void mc1i_m1 ( int j ); /* (21) (29) (31)(37)(51)*/

void ts0_c1i ( void ); /* (22) (36) */