cwp w T

5 専用チップでの過程 ²⁹

c ^p w ^(k)

c

p (23)

(24)

(23)TC

2 C

CP T

wの計算⁽実はこれはスカラーです。⁾

(24)AC

1 S

CP T

wの計算終了

w ^(k)

- 1 2 cp w ^T

p (25) (26)

(27)

(25) IS

1 (0

2 )

(26) TS

1 S

0 0

2 cp

wをつくる。

(27) CC

2 C

w ^(k)

cp w ^T

w ^(k)

1 - 2

1 - 2 p (28)

(29)

5 専用チップでの過程 ³⁰

w ^(k) ^q ^p ^- ¹ ₂ ^{cwp w} ^T

(30)

(30)AM

1 C

q = p0 c

2 wP

w の計

算終了。

w ^(k)

(34)

(32) (31)

(33)

q q

(30-2)CLR M

(31)CC

1 M

(1列目⁾

qを^M¹上へ入れる

(32)CC

2 C

(33)CC

1 R

(34)CM

1 (1)C

w ^(k)

-1 (38)

- q (35)

(37) (36)

(40)

(35)IS

0 (01)

(36)TS

0 C

-q できる

(37)CC

1 M

(38)TR

1 M

1 0qw

T の計算終了。

(39)AM

1 M

A 0 qw

T の計算終

了。

5 専用チップでの過程 ³¹

w ^(k) w ^(k)

(36)

-q (43)

(41) (42)

(40)

(44)

(40)CC

1 R

(41)CC

2 C

(42)CC

1 M

1 (all)

(43)TR

1 M

1 0wq

T の計算終了。

(44)AM

1 M

以上を^N ⁰²回繰り返すことにより^M2 の行列が対角化され三重対角化は完成される。⁴

/home2/students/mizuki/xg/no1.eps〜^no18.eps

5 専用チップでの過程 ³²

5.2 二分法の過程

ここで^kは¹から²⁵ まで動き^kkは１から^N まで動く。なぜ^kが²⁵なのかというと２分法では１回ごとに幅が ¹

ずつ狭くなるのでこの回数が精度に影響する。いまは ²⁵回としているが⁵⁰回以上⁽できれば⁶⁴回くらい⁾がよい。

kkは^N 回でこれは行列の行数⁽または列数⁾である。

A ~

(0)M

2に三重対角化された^A^~ が入る。

α ₁

α _n β n-1

α _n-1 β _n-1 β 1

α ₂ β 2

β ₂ α ₃ β 1

β ₂

β n-1

β 1

0 2 2

- 2

-- 0

0 (1) (2)(3)

(4)

(1)(FM

1 )C

2 副対角成分を^C²に入れる。

(2)TC

2 C

(3)T(01)C

01を^C2に掛け結果を^C2に入れる。

(4)ZM

1 M

1の絶対値をとる。

5 専用チップでの過程 ³³

α ₁

α _n β _n-1 α _n-1

β n-1

β 1

α ₂ β ₂

β ₂ α ₃

β 1 0

0 0 0

(5)

r r max (6)

(5)AM

1 C

(6-1)Max(C

1 )S

0 C

1の最大値をとり^S0に入れる。

α ₁ α ₂ α ₃

α n

α _n-1 β _n-1

β _n-1 β 1

α n

α _n-1 α ₂

β ₂ β ₂ α ₃ β 1

α ₁

0 0

(6-2)

(6-2)(TM

1 )C

ここ⁽⁷⁾からは仮に^R3

5を用いている。これは ²分法の上限下限を表している。

(7-1)CS

0 R

(7-2)InputS

1 (01)

(7-3)TS

1 R

(7-4)MS

0 R

(7-5)AR

3 R

α ₁ α ₂ α ₃

α _n

α _n-1 β _n-1

β _n-1 β 1

α _n α _n-1 α ₂

β ₂ β ₂ α ₃ β 1

α ₁

0 0

(7-6)MR

4 R

(7-7)InputS

1 (0:5)

(7-8)TS

1 R

(7-9)MR

3 R

(7-10)T(01)R

012R

4の計算

以上ではじめの範囲^(0rから^r) が確定する。

5 専用チップでの過程 ³⁴

α ₁ α ₂ α ₃

α _n

α _n-1 β _n-1

β n-1

β 1

α _n α _n-1 α ₂

β ₂ β ₂ α ₃ β 1

α ₁

0 0

(8-0)T(01)R

5 (all)

(8-1)MC

1 i(1)M

(1;k)k は移動する。

(8-2)AR

5 (1)M

(1;kk) kk は移動する。

(9)MM

(k01;kk)R

1 (kk)

(10)MC

1 (k)M

(k;kk)

(11)GR

(kk)R

1の逆数をとる

(12)TC

(k-1)R

(kk)kkは列です

(13)AR

5 (1)M

(1;kk) kkは列で移動する。

(14)MM

(k-1,kk)R

(kk)M

1の^(k-1 行^kk列⁾の値を^R1の^kk列に入れる。

なお⁽⁹⁾〜⁽¹⁴⁾までは^N 回繰り返す。⁽これは¹列目から^N 列目までといることである。⁾

(15)MC

1 (k)M

(k;kk)

(16)Countこれは^M1のそれぞれの列ベクトルの値が正である数を数えて結果を^R1

に入れる。

(17)Compare R

1と^R2の比較をする。^R2 には左から¹、²、³、と順に値が入っている。この結果

r1j[j]>= r2j[j]なら

r4j[j]=r5j[j]そして^r5j^[j]⁼ ^r4[j]⁺^r3[j]

となり、

そうでないなら^r3j[j]⁼^r5j^[j]そして ^r5j[j]⁼ ^r4j[j]⁺^r3j[j]

となる。

5 専用チップでの過程 ³⁵ このあと

(18)MR

5 R

(kk) kkは移動する。

(19)InputR

4 (0)

(20)MR

6 R

(21) R

(kk+1)+R

(kk+1)

を^R5

(kk+1)へ^(kkは列を表す。⁾

このあと^(8-0)へいき^kk回繰り返す。

そしてそれぞれの列で計算するので^N 回繰り返す。

以上の過程で^R6に固有値が入る。

(ただし²分法は^R6まであるのでさらに改良する可能性がある。⁾ ⁵

/home2/students/mizuki/xg2/no1.eps〜^no6.eps

5 専用チップでの過程 ³⁶

5.3 逆反復法の過程

µ

(0)M

1に行列を入れる。この行列を^Aとする。

µ -1 -1 -1

(1) (2)

(3)

-A

(1)IS

0 (01)

(2)CS

0 R

(3)TR

1 M

1 には^0Aがはいる

nn

-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¹ に入れる。

(5-1)CS

1 S

5 専用チップでの過程 ³⁷

11 22

µ -a µ -a

33 -a

-a nn

µ

µ-a ₂₂

33 µ -a

-a nn

µ µ

µ -a ₁₁

(6) (5-2)

(5-2)AS

0 C

(6)C

1 を^M1 の対角成分に入れる。

X 11

X n1

X 1n

X nn

x 1

x n

x 2

N

N+1

ここで^X ⁼^I ⁰^Aとおき成分を

11で表す。

また

x= 0

@ x

n 1

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

X nn

X 11

X n1 x _n

x 2

X 1n x 1

X 11

N

N+1 -X ₁₁

-X _n1

-1 X

1n x X ₁₁ 1

(10)

(8)

(9) (7)

(7)IS

0 (01)

(8)CM

(j列^)C1 (9)CM

(i行^)R1 (10)TS

0 C

5 専用チップでの過程 ³⁸

X nn x _n x 2

X 1n x 1

X

11 X n1

N

N+1 -X ₁₁

-X _n1

X 1n x X 1

11 1

X

11 (12) (11)

(13)

(11)CR

(j列^)S⁰

(12)GS

X nn x _n

11 1 _X ^X ¹²

11 x 1

X X _1n X 11

x 2 ¹¹

1 _X ^X ¹²

11 x 1

X X _1n X ₁₁

X n1

-X _n1 1 X

11 X ₂₂ X 21

0 X ₂₁

(14)

(17) M ₂ (15)

(16)

(13)

(18)

(13)TS

0 R

(14)CR

1 M

1 (j列⁾

(15)IS

0 (0)

(16)CS

0 C

1 (i行⁾

(17)CM

1 M

11 1 _X ^X ¹²

11 x 1

X X _1n X 11

X 21

- X

- 21 X X

12 11

X 21

- X ₁₁ x ₁

-X _n1 0 X ₂₁

(18)

(19)

0 0 0 0

(18)CC

1 M

(19)TR

1 M

(20)AM

2 M

5 専用チップでの過程 ³⁹

U _k 1

1 1 1

0 0

(7)〜⁽²¹⁾を^N 回繰り返す。この結果ⁿ⁺¹ 列目にできるベクトル

kが固有ベクトルとなる。

結果ⁿ ⁺ ¹列目にできるベクトル^Ukが固有ベクトルとなり、これをとりだせば良い。⁶

6以上の図は/home2/students/mizuki/xg3/no1.epsから^no10.eps

6 考察 ⁴⁰

考察

前章までで述べたように多くの機能を現段階では用いねばならないが実現できることがわかった。しかしまだ問題点が多くあり改良しなければならない。

まずシミュレーションプログラムについてはさらなる検討の余地がある。特に無駄な部分を取りシミュレーションプログラム自体の時間の短縮が必要である。

実際に計算させたところハウスホルダー変換では¹²⁸² ¹²⁸行列で^7.12秒なのに対し逆行列での計算⁽逆反復法⁾がかなりの時間⁽¹²⁸²¹²⁸行列で ²⁸⁴秒⁾を要し、

改良しなければならない。⁽しかしこれはファイルアクセスの時間も含んでいる。⁾ また²分法でのプログラムはレジスターの数が多く定義^(R¹から^R⁶の⁶つ⁾している。これは途中結果を１回１回格納しなければならないためだか、なるべくレジスターを有効に利用し少なくなるようにしなければならない。あまりに多過ぎるとチップをつくるのが困難になってしまう可能性がある。

7 結論 ⁴¹

結論

この専用チップを用いた結果、現段階では多くの過程が必要だが専用チップを用いて固有値固有ベクトルを求めることが可能である。

そしてチップのさらなる効率化を考えハード化することにより固有値固有ベクトルを求める際専用チップを用い計算時間が短縮すると考えられる。

7 結論 ⁴²

謝辞

本研究及び論文作成に当たり、終始御懇切なる御指導、御鞭撻を賜わりました指導教官である齋藤理一郎助教授に衷心より御礼の言葉を申しあげます。

また、本研究を進めるにあたり、熱心な御指導をいただくとともに種々の御高配を賜わりました木村忠正教授、湯郷成美助教授に深謝の意を表します。

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

そして、数々の御援助、御助言をしていただいた中平政男氏、竹谷隆夫氏、はじめ木村¹齋藤研究室の大学院生、卒研生の方々、高田さんをはじめ⁽株⁾画像技研の方々に感謝します。

参考文献 ⁴³

参考文献

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

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

[3] C言語によるプログラミング⁽基礎編・応用編⁾ 内田智史オーム社

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

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

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

このプログラムは ^hhc.c ^, ^2bun.c ^,gyaku.c の ³つのプログラムからなり^,hhc.c はハウスホルダー変換の実行 ^2bun.c は ²分法の実行

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に値⁰を入れる関数

print_m1 などの関数は^m1に値⁰を入れる関数

_ の後には格納場所が入る

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

---*/

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) */

/* (23)はある ^*/

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

void mc2i_c1i ( void ); /* (11-1)(27)(32) (50)*/

void mc1i_r1j( void ); /* (33) (40)*/

void mm1_c1i ( int j ); /* (34) */

void am1_m2 ( void ); /* (35) */

void mm2_m1( void ); /* end */

void mm2_m3( void );

void mm3_m2( void );

void mr1j_m2 ( void ); /* Qの計算 ^*/

void tc1i_m2 ( void );

void ms0_c1i( int i );

void ts2_m2 ( void ) ;

void ar2j_m2 ( int a );

/*---main ( )

---*/

void main( void )

{

FILE *q ;

FILE *fp ;

FILE *fp1;

#ifdef DBGPRN

time_t t;

t = time( NULL );

printf("%s\n",ctime( &t ));

#endif

/* 専用レジスタ ^clear ^*/

s0 = 0;

s1 = 0;

s2 = 0;

/* 専用レジスタに数値代入 ^*/

input( ); /* m1 m2 に数値代入 ^*/

input_c1i( );

input_c2i( );

input_r1j( );

input_r2j( );

#ifdef DBGPRN

printf("もとの行列表示^\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

/*---３重対角化 ^start

---*/

/* kを動かす ^k=1 から ^k を用いるのは ^1,3,です ^*/

for ( k = 1 ; k <= (N-2) ; k++){

m1_c1i ( k ); /* (1) m1_c1移動 ^kは列^*/

s0 = 0 ; /* (2) clr_s0 */

s0_c1i( k ); /* (3)-1 */

#ifdef DBGPRN /* 確認 ^*/

printf("(3)終了後 ^\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

c1i_c2i( ); /* (4) mc1i_c2i */

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

s0 = ac1i_s0( );/* (6) ac1_s0 */

#ifdef DBGPRN /* 確認 ^*/

printf("(6)終了後 ^\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

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

print_r1j( );

print_r2j( );

#endif

sqrt_s0( ); /* (7) */

#ifdef DBGPRN /* 確認 ^*/

printf("(7)\n");

printf("s0 = %f (sの答え)\n",s0);

printf("(8)start\n");

printf("k= %d \n",k);

#endif

mc2i_s1( k+1 ); /* (8) */

as1_s0( ); /* (9) */

#ifdef DBGPRN

printf("(9) \n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

#ifdef DBGPRN

printf("符号確認 ^\n");

printf("s0 = %f\n",s0);

#endif

ms0_c2i( k+1 ); /* (10) */

#ifdef DBGPRN

printf(" w の値^\n");

print_c2i( );

#endif

mc2i_c1i( ); /* (11-1) */

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

#ifdef DBGPRN /* 確認 ^*/

printf("(11-2)\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

/*---Q の計算 ^(No ¹⁾ ---*/

mr1j_m2( ); /* Q の計算 ⁽¹⁾ ^*/

tc1i_m2( ); /* Q の計算 ⁽²⁾ ^w ^* ^{w^T} の計算終了 ^*/

/*---Q の計算 ^end ^(No ¹⁾ ---*/

tc2i_c1i( ); /* (12-1) */

s0 = ac1i_s0( );/* (13) w^T * w の計算終了 ^*/

#ifdef DBGPRN /* 確認 ^*/

printf("(13)\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( ) ;

print_m2( ) ;

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

s1 = 2 ; /* (14) */

if ( s0 == 0 ){

printf("s0 = 0 です。⁰で割れません^!! ^\n");

}

else{

s0 = 1 / s0 ; /* (15) */

}

#ifdef DBGPRN /* 確認 ^*/

printf("(15)\n");

printf("s0 = %f\n",s0);

#endif

ts1_s0( ); /* (16) c の計算終了 ^*/

/*---Q の計算 ^(No ²⁾ ---*/

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

for ( i = 1 ; i <= N ; i++){

ms0_c1i( i ); /* (3) */

}

tc1i_m2( ); /* (4) */

s2 = -1 ; /* (5) */

ts2_m2( ); /* (6)*/

for ( i = 1 ; i <= N ; i++ ){

r2j[i] = 1 ; /* (7) */

ar2j_m2( i ); /* (8) */

#ifdef DBGPRN

printf(" (7) Input r2j(%d) \n ",i);

printf(" (8) ar2j(%d)_m2(%d) \n ",i,i);

print_r2j( );

#endif

for( j = 1 ; j <= N ; j++){

r2j[j] = 0 ;

}

/* printf(" result \n");

* print_m1( );

* print_m2( );

/*---Q の計算 ^(No ²⁾ end---*/

if((q=fopen("q.dat","a")) == NULL){

printf("Can not open file q.dat!! \n");

exit(1);

}

fprintf (q," k = %d\n",k);

for ( i = s; i <= N ; i++ ){

for ( j = s ; j <= N ; j++ ){

fprintf( q, "%10.6lf",m2[i][j] );

}

fprintf(q,"\n");

}

fclose(q);

ts0_r1j( ); /* (17) */

#ifdef DBGPRN /* 確認 ^*/

printf("(17)\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

tr1j_m1( ); /* (18) */

for ( i = s; i <= N; i++){

c1i[i] = am1_c1i( i ); /* (19) P=cAW の計算終了 ^*/

/* printf(" c1i[%d]= %f \n",i,c1i[i]); */

}

#ifdef DBGPRN /* 確認 ^*/

printf("(19)\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

clr_m1( ); /* (20) */

j = s; /* jは固定 ⁽一番左端の列⁾ ^*/

mc1i_m1 ( j ); /* (21 ) */

ts0_c1i ( ); /* (22) */

#ifdef DBGPRN /* 確認 ^*/

printf("(22)\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

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

print_c2i( );

print_r1j( );

print_r2j( );

#endif

tc2i_c1i ( ); /* (23) */

s0 = ac1i_s0_2( ) ;

#ifdef DBGPRN /* 確認 ^*/

printf("(24)\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

s1 = -0.5 ; /* (25) */

ts1_s0( ) ; /* (26) */

mc2i_c1i ( ); /* (27) */

#ifdef DBGPRN /* 確認 ^*/

printf("(27)\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

ts0_c1i( ) ; /* (28) */

#ifdef DBGPRN /* 確認 ^*/

printf("(28)終了後^\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

j=s+1; /* jは任意 ⁽ただし^jは⁽²¹⁾でいれた値⁺¹です^!!)(一番左端の列⁾ ^*/

mc1i_m1( j ); /* (29) */ /* 注意ｊの値 ^*/

#ifdef DBGPRN /* 確認 ^*/

printf("(29)終了後^\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

for( i= s; i <=N ; i++){

c1i[i] = am1_c1i_2(i ); /* (30-1) */

/* printf(" c1i[%d] = %f \n",i,c1i[i]); */

}

#ifdef DBGPRN /* 確認 ^*/

printf("(30-1)終了後^\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

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

clr_m1( );/* (30-2) */

j = s ; /* jは固定 ⁽一番左端の列⁾ ^*/

mc1i_m1( j ); /* (31) */

mc2i_c1i( );/* (32) */

mc1i_r1j( );/* (33) */

#ifdef DBGPRN /* 確認 ^*/

printf("(33)end \n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

j = s ;

mm1_c1i( j ); /* (34) j 固定 ^*/

#ifdef DBGPRN /* 確認 ^*/

printf("(34)\n");

print_m1( );

print_c1i( );

#endif

s0 = -1; /* (35) */

ts0_c1i ( ); /* (36) */

#ifdef DBGPRN /* 確認 ^*/

printf("(36)\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

for ( j =s; j <=N ; j++){

mc1i_m1( j ); /* (37) */

}

#ifdef DBGPRN /* 確認 ^*/

printf("(37)\n");

print_c1i( );

print_m1( );

#endif

tr1j_m1( ) ;/* (38) */

#ifdef DBGPRN /* 確認 ^*/

printf("(38)\n");

print_r1j( );

print_m1( );

#endif

if ( k == 1 ) {

input_m2( );

}

else {

mm3_m2( );

}

am1_m2 ( ) ; /* (39) */

#ifdef DBGPRN /* 確認 ^*/

printf("(39)\n");

print_m1( );

print_m2( );

#endif

#ifdef DBGPRN /* 確認 ^*/

printf("(39)\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

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

mc1i_r1j( ); /* (40) */

#ifdef DBGPRN /* 確認 ^*/

printf("(40)\n");

print_c1i( );

print_r1j( );

#endif

mc2i_c1i( ); /* (50) or (41)*/

#ifdef DBGPRN /* 確認 ^*/

printf("(50) -->(41)\n");

print_c2i( );

print_c1i( );

#endif

for ( j =s; j <=N ; j++){

mc1i_m1 ( j ); /* (51) 本当は^(42)*/

}

#ifdef DBGPRN /* 確認 ^*/

printf("(51) -->(42)\n");

print_c1i( );

print_m1 ( );

#endif

tr1j_m1( ); /* (52) */

#ifdef DBGPRN /* 確認 ^*/

printf("(52) -->(43)\n");

print_r1j( );

print_m1 ( );

#endif

am1_m2 ( ); /* (53)本当は⁽⁴⁴⁾ ^*/

#ifdef DBGPRN /* 確認 ^*/

printf("(53)-->(44)\n");

printf("s0 = %f\n",s0);

printf("s1 = %f\n",s1);

printf("s2 = %f\n",s2);

print_m1( );

print_m2( );

print_c1i( );

print_c2i( );

print_r1j( );

print_r2j( );

#endif

#ifdef DBGPRN

printf("結果^\n");

print_m2( );

#endif

for ( i = s; i<=N ;i++){

al[i] = m2[i][i] ;

/* printf ("al[%d]= %lf \n",i,al[i]); */

}

for ( i =s; i<=N-1 ;i++){

be[i] = m2[i+1][i] ;

/* printf ("be[%d]= %lf \n",i,be[i]); */

}

mm2_m3( );

mm2_m1( ); /* m2 m1 copy */

#ifdef DBGPRN /* 確認 ^m2 ^and ^m1 ^*/

print_m2( );

print_m1( );

#endif

s = s + 1 ;

/* 専用レジスタ ^clear ^*/

s0 = 0;

s1 = 0;

s2 = 0;

input_c1i( );

input_c2i( );

input_r1j( );

input_r2j( );

} /* for end */

if((fp=fopen("al.dat","w")) == NULL){

printf("Can not open flie al.dat!! \n");

exit(1);

}

if((fp1=fopen("be.dat","w")) == NULL){

printf("Can not open file be.dat!! \n");

exit(1);

}

for ( i = 1; i <=N ; i++ ){

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

fprintf(fp1,"%f \n ",be[i]);

}

fclose(fp);

fclose(fp1);

} /* main end */

/*********************************************************

* 行列 ^input 関数

*********************************************************/

void input ( void )

{

m1[1][1] = 4; m1[1][2] = 3; m1[1][3] = 2; m1[1][4] = 1;

m1[2][1] = 3; m1[2][2] = 3; m1[2][3] = 2; m1[2][4] = 1;

m1[3][1] = 2; m1[3][2] = 2; m1[3][3] = 2; m1[3][4] = 1;

m1[4][1] = 1; m1[4][2] = 1; m1[4][3] = 1; m1[4][4] = 1;

m2[1][1] = 4; m2[1][2] = 3; m2[1][3] = 2; m2[1][4] = 1;

m2[2][1] = 3; m2[2][2] = 3; m2[2][3] = 2; m2[2][4] = 1;

m2[3][1] = 2; m2[3][2] = 2; m2[3][3] = 2; m2[3][4] = 1;

m2[4][1] = 1; m2[4][2] = 1; m2[4][3] = 1; m2[4][4] = 1;

}

/*********************************************************

* 行列 ^input_m2 関数

*********************************************************/

void input_m2 ( void )

{

m2[1][1] = 4; m2[1][2] = 3; m2[1][3] = 2; m2[1][4] = 1;

m2[2][1] = 3; m2[2][2] = 3; m2[2][3] = 2; m2[2][4] = 1;

m2[3][1] = 2; m2[3][2] = 2; m2[3][3] = 2; m2[3][4] = 1;

m2[4][1] = 1; m2[4][2] = 1; m2[4][3] = 1; m2[4][4] = 1;

}

/*********************************************************

* 行列 ^input 関数

*********************************************************/

void input_test ( void )

{

for ( i=1;i<=N ;i++){

for(j=i;j<=N;j++){

m1[i][j]=(i+j)*(i+j);

m2[i][j]=(i+j)*(i+j);

}

for ( i=2;i<=N;i++){

for (j =1; j <=i-1;j++){

m1[i][j]=m1[j][i];

m2[i][j]=m2[j][i];

}

#ifdef DBGPRN

for ( i=1;i<=N ;i++){

for(j=i;j<=N;j++){

print_m1( );

print_m2( );

}

#endif

}

/*********************************************************

* < 行列 ^r1j ^input 関数 ^> ^*

*********************************************************/

void input_r1j ( void )

{

int a,b;

for(a=1;a<=N;a++){

for(b=1;b<=N;b++){

r1j[a]= 0;

}

/*********************************************************

* < 行列 ^r2j ^input 関数 ^> ^*

*********************************************************/

void input_r2j ( void )

{

int a,b;

for(a=1;a<=N;a++){

for(b=1;b<=N;b++){

r2j[a]= 0;

}

/*********************************************************

< 行列 ^c1i ^input 関数 ^>

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

*********************************************************/

void input_c1i ( void )

{

int a,b;

for(a=1;a<=N;a++){

for(b=1;b<=N;b++){

c1i[a]=0;

}

/*********************************************************

* < 行列 ^c2i ^input 関数 ^> ^*

*********************************************************/

void input_c2i ( void )

{

int a,b;

for(a=1;a<=N;a++){

for(b=1;b<=N;b++){

c2i[a]=0;

}

/*************************************************************

* < print_r1j 関数 ^> ^*

**************************************************************/

void print_r1j ( void )

{

int a;

for(a=s;a<=N;a++){

printf(" r1j[%d]=%f ",a,r1j[a]);

}

printf("\n\n");

}

/*************************************************************

* < print_r2j 関数 ^> ^*

**************************************************************/

void print_r2j ( void )

{

int a;

for(a=s;a<=N;a++){

printf(" r2j[%d]=%f ",a,r2j[a]);

}

printf("\n\n");

}

/*************************************************************

* < print_c1i 関数 ^> ^*

**************************************************************/

void print_c1i ( void )

{

int a;

for(a=s;a<=N;a++){

printf(" c1i[%d]=%f \n",a,c1i[a]);

}

printf("\n");

}

/*************************************************************

* < print_c2i 関数 ^> ^*

**************************************************************/

void print_c2i ( void )

{

int a;

for(a=s;a<=N;a++){

printf(" c2i[%d]=%f \n",a,c2i[a]);

}

printf("\n");

}

/*************************************************************

* < print_m1 関数 ^> ^*

**************************************************************/

void print_m1 ( void )

{

int a,b;

for(a=s;a<=N;a++){

for(b=s;b<=N;b++){

printf("m1[%d][%d]=%f ",a,b,m1[a][b]);

}

printf("\n");

}

printf("\n");

}

/*************************************************************

* < print_m2 関数 ^> ^*

**************************************************************/

void print_m2 ( void )

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

int a,b;

for(a=s;a<=N;a++){

for(b=s;b<=N;b++){

printf("m2[%d][%d]=%f ",a,b,m2[a][b]);

}

printf("\n");

}

printf("\n");

}

/*---関数---*/

/******************************************************

(1) m1_c1i ( int b ) m1から^c1iへ移す^(copy)

******************************************************/

void m1_c1i ( int b )

{

int a;

for(a=s;a<=N;a++){

c1i[a]=m1[a][b];

}

/*****************************************************

(3) s0_c1i ( int a ) 移動^(copy)

/*****************************************************/

void s0_c1i ( int a )

{

c1i[a] = s0;

}

/*****************************************************

(4) c1i_c2i ( ) 移動

*****************************************************/

void c1i_c2i ( void )

{

int a;

for ( a = s; a <=N ; a++ ){

c2i[a] = c1i[a];

}

/***********************************************

(5) tc2_c1( ) (times( )) 掛け算関数

***********************************************/

void tc2_c1( )

{

int a; /* aは行 ^*/

for( a = s ; a <=N ; a++ ){

c1i[a] = c2i[a] * c1i[a] ;

}

/*************************************************************

< c1iの¹行目から^N行目まで足しあわせる。^>

(6) double ac1i_s0( ) ( sum_gyou( ) )

**************************************************************/

double ac1i_s0( )

{

int a;

double t = 0 ; /* t = total */

for ( a = s ; a <=N ; a++){

t = t + c1i[a];

}

return t;

}

/******************************************************

(7) sqrt_s0( )

******************************************************/

void sqrt_s0 ( void )

{

s0 = sqrt( s0 );

}

/*****************************************************

(8) mc2i_s1 ( int a ) 移動 ^(copy)

*****************************************************/

void mc2i_s1 ( int a )

{

#ifdef DBGPRN

printf(" (8) \n");

printf("a = %d \n",a);

#endif

s1 = c2i[a];

if ( s1 < 0 ){

s0 = -s0 ;

}

#ifdef DBGPRN

printf("s1 = %f\n",s1);

#endif

}

ドキュメント内超高速行列演算チップの作成 (ページ 33-81)

c p w (k)

c

p (23)

(24)

w (k)

- 1 2 cp w T

p (25) (26)

(27)

w (k)

cp w T

w (k)

1 - 2

1 - 2 p (28)

(29)

w (k) q p - 1 2 cwp w T

(30)

w (k)

w (k)

(34)

(32) (31)

(33)

q q

w (k)

w (k)

-1 (38)

- q (35)

(37) (36)

(40)

w (k) w (k)

(36)

-q (43)

(41) (42)

(40)

(44)

A ~

α 1

α n β n-1

α n-1 β n-1 β 1

α 2 β 2

β 2 α 3 β 1

β 2

β n-1

β 1

0

2 2

- 2

-- 0

0

(1) (2)(3)

(4)

α 1

α n β n-1 α n-1

β n-1

β 1

α 2 β 2

β 2 α 3

β 1 0

0

0 0

(5)

r r max (6)

α 1 α 2 α 3

α n

α n-1 β n-1

β n-1 β 1

α n

α n-1 α 2

β 2 β 2 α 3 β 1

α 1

0 0

(6-2)

α 1 α 2 α 3

α n

α n-1 β n-1

β n-1 β 1

α n α n-1 α 2

β 2 β 2 α 3 β 1

α 1

0 0

c ^p w ^(k)

w ^(k)

- 1 2 cp w ^T

w ^(k)

cp w ^T

w ^(k)

w ^(k) ^q ^p ^- ¹ ₂ ^{cwp w} ^T

w ^(k)

w ^(k)

w ^(k)

w ^(k)

w ^(k) w ^(k)

α ₁

α _n β n-1

α _n-1 β _n-1 β 1

α ₂ β 2

β ₂ α ₃ β 1

β ₂

α ₁

α _n β _n-1 α _n-1

α ₂ β ₂

β ₂ α ₃

α ₁ α ₂ α ₃

α _n-1 β _n-1

β _n-1 β 1

α _n-1 α ₂

β ₂ β ₂ α ₃ β 1

α ₁

α ₁ α ₂ α ₃

α _n

α _n-1 β _n-1

β _n-1 β 1

α _n α _n-1 α ₂

β ₂ β ₂ α ₃ β 1

α ₁

α ₁ α ₂ α ₃

α _n

α _n-1 β _n-1

α _n α _n-1 α ₂

β ₂ β ₂ α ₃ β 1

α ₁

µ-a ₂₂

µ -a ₁₁

X n1 x _n

N+1 -X ₁₁

-X _n1

1n x X ₁₁ 1

X nn x _n x 2