1 準備

1.1 関数を使おう

ガウス・ジョルダン法は、連立1次方程式の解、あるいは逆行列を求める汎用的な手法で、あらゆるこの種の問題に適用できる。良いプログラムを書けば、再利用するのである。このように、再利用できそうなルーチンは、メイン関数には書かないで、関数という機能でプログラムに組み込む²。そうすると、ほかの場所で連立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)のxと yを実引数と呼ぶ。そして、その処理を書いている関数add(double xin, double yin)のxinとyinを仮引数という。呼ぶ方の変数を実引数、呼ばれる方の変数を仮引数と呼ぶのである。

	#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}$ は一意に決まらないのは、線形代数の言うところである。その場合、ガスス・ジョルダン法で計算する関数は、呼び出し側へ警告を出さなくてはならない。呼び出し側へゼロを返すことで実現できる。それは、

行列が特異な場合、
return(0);
行列が正則な場合、
return(1);

と書けば良い。

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日