1 準備

1.1 関数を使おう

ガウス・ジョルダン法は、連立1次方程式の解、あるいは逆行列を求める汎用的な手法で、 あらゆるこの種の問題に適用できる。良いプログラムを書けば、再利用するのである。こ のように、再利用できそうなルーチンは、メイン関数には書かないで、関数という機能で プログラムに組み込む2。そうすると、ほかの場所で連立1 次方程式を解きたい場合は、その関数をコールするだけで済むのでプログラム作成の手間 が省ける。また、ほかのプログラムを書く場合でも、それをコピーして、再利用できるの で便利である。そういうわけで、メイン関数ではなく、ガウス・ジョルダン法の専用の関 数を作ることが望ましい。その方がプログラムは分かりやすくなる。。C言語は関数の独 立性が高いので、その辺は比較的簡単にできる。

ガウス・ジョルダン法の専用の関数を作る場合、どうするか?。まず、その入 出力と機能を考える。入力は、係数行列の $ \boldsymbol{A}$と非同次項の $ \boldsymbol{b}$が必要である。 そして、出力は逆行列の $ \boldsymbol{A}^{-1}$と解のベクトル $ \boldsymbol{x}$となる。 要するに、図[*]のような機 能の関数をつくるわけである。

図: ガウスジョルダン法を計算するC言語の関数のブラックボックス
\includegraphics[keepaspectratio, scale=1.0]{figure/GaussJordan_BlackBox.eps}

この関数ができると、問題を解く時に必要な係数行列と非同次項を入力さえす れば、逆行列と解を計算してくれる。ガウス・ジョルダン 法の手続きを、関数で実現する方法について、次節以降で丁寧に説明する。

1.2 関数へのデータの受け渡し

1.2.1 値の渡し方を思い出そう

機能は決まったので、つぎに、関数へのデータの受け渡しを考えなくてはならない。C言 語のデータの渡し方は、少し複雑なので、復習をする。

いかなるプログラム言語でも、C言語の関数(mainではない)に対応するものが有り、それ は、サブルーチンと呼ばれている。同じような処理がある場合、1つの独立した処理のブ ロックとしてまとめ、どこからでもコールすることができるようにすると便利である。こ のような、機能のブロックをサブルーチンと呼ぶ。例えば、sin(x)などである。 $ \sin x$の計算が必要な都度、その処理を書いていたのではたまらないので、独立した関 数としてその処理が書くのである。これは、ライブラリーとなっているので、その処理内 容は通常は分からない。

この関数に、データを与える変数のことを引数と言う。先ほどの例で言うと、 sin(x)xが引数である。プログラムの中では、sinと言う 名前がついている処理にxを与え、それを処理する。

引数には、2種類(実引数と仮引数)が有る。それについて、次のプログラムで説明する。 この場合、main関数でコールするときの文、add(x,y)xyを実引数と呼ぶ。そして、その処理を書いている関数add(double xin, double yin)xinyinを仮引数という。呼ぶ方の変数を実引 数、呼ばれる方の変数を仮引数と呼ぶのである。

	#include <stdio.h>
	double add(double xin, double yin);

	/* ========== main関数 =================*/
	main(){
	   double x, y, wa;

	   いろいろな処理
		
	   wa=add(x,y);

	   いろいろな処理

	}

	/* ========== 足し算の関数 =================*/
	double add(double xin, double yin){
	   double zout;

	   zout = xin+yin;
	   return(zout);
	}

実引数から仮引数に値を送る方法は、C言語では2通りの方法がある。以前学習した通りで、

値渡し(Call by value)
コールした(呼び出した)関数と処理する関数 では、別々のメモリー領域を用意する。そして、コールしたと きに呼び出し側のメモリー(実引数)の値を処理する関数のメモリー (仮引数)にコピーする。従って、処理する関数が仮引数の値を 変えても、呼び出し側の実引数の値が変わることはない。
アドレス渡し(Call by reference)
処理する関数では値を 格納するメモリー領域を用意しない。実引数は、コールした関数の 実引数の値が格納されているメモリーのアドレスとなる。そ のアドレスが仮引数に渡される。従って、処理する関数はそ のアドレスを格納するメモリー領域を用意するだけである。処理す る関数は、実引数のアドレスに格納されている値を処理することに なる。この場合は、呼び出された関数が仮引数の値を変える と、呼び出し側の実引数の値も変わります。
である。

C言語の場合、通常の変数(配列でない)の場合、値渡しである。これは良くできた仕様で ある。処理する関数が呼び出し側のデータを変えることが無いので、プログラミングの時、 余計な気を使わないですむ。関数の独立性が高いといわれる所以である。実際の例では、 先のaddという関数は、main関数から呼び出されており、main の実引数の値(x,y)の値が、addの仮引数(xin,yin)にコピー される。そこでの処理の結果は、戻り値(返却値)zoutに入れられて、元の関数 に戻す。元の関数のwaに、zoutの値がコピーされるのである。

一方、配列を処理する関数に渡す場合は、アドレス渡しになる。一般に配列のデータは、 通常の変数よりもかなり大きく、それをいちいちコピーしていたら不経済ということが理 由と言われている。

ということで、今回の場合、配列を渡すためアドレス渡しになる。処理する関数でその配 列の値を変えると、コールした関数のその値も変わる。しかし、これは便利なこともある。 いちいち戻り値を与える必要が無く、気軽に呼び出した関数に結果を戻せる(FORTRANと同 じ)。

従って、図[*]のような入出 力の関数を実現するための引数の書き方は、

	#include <stdio.h>
	void gaussjordan(double a[][100], double b[100],
	                 double inv_a[][100], double x[100]);

	/* ========== main関数 ==============================================*/
	void main(){
	   double a[100][100], b[100];
	   double inv_a[100][100], x[100];

	   いろいろな処理
		
	   gaussjordan(a,b,inv_a,x);

	   いろいろな処理
	}

	/* ========== ガウスジョルダン法の関数 =============================*/
	void gaussjordan(double a[][100], double b[100],
	                 double inv_a[][100], double x[100]){

	いろいろな処理

	inv_a[i][j] = いろいろな計算
	b[i] = これも計算
	
	いろいろな処理

	}

となる。処理する関数の方は、一番最初のサイズを除いて書く必要がある。これは、例え ばz[100][200][300]の大きさの配列のz[23][73][36]というデータに アクセスする場合を考えれば分かる。このデータがあるアドレスは、[Zのアド レス+配列1個のデータサイズ $ \times(36+73\times300+23\times200\times300)$]になる。 したがって、最初の配列のサイズを除いて、配列のサイズがアドレス計算に必要となり、 処理する関数側で明示する必要がある。

実際のプログラムでは、もう少し効率よく配列を使うが、大筋はこの通り である。これで、関数に値を与える復習は終わり。

1.2.2 行列が特異な場合の警告

もし係数行列 $ \boldsymbol{A}$が特異だと、解 $ \boldsymbol{x}$は一意に決まらないのは、線形代数の言う ところである。その場合、ガスス・ジョルダン法で計算する関数は、呼び出し側へ警告を 出さなくてはならない。呼び出し側へゼロを返すことで実現できる。それは、 と書けば良い。

1.2.3 行列のサイズはどうするの

本当にこれだけでよいのか?。先の例で配列を値を、処理するプログラムに知らせること ができるが、これで全て計算の準備が整ったわけではない。行列やベクトルのサイズを関 数に知らせる必要がある。これは、大きな配列を用意しておいて、その一部分に係数や非 同次項の値を入れるため、処理するときに行列やベクトルの大きさが必要となる。このよ うな理由から、行列やベクトルのサイズを渡さなくてはならない。そこを考慮すると、ガ ウス・ジョルダン法の関数のプロトタイプ宣言は、次のようになる。

	int gaussjordan(int n, double a[][100], double b[100],
	                double inv_a[][100], double x[100]);




ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成16年11月17日


no counter