1 準備

1.1 関数を使おう

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

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

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

この関数ができると、問題を解く時に必要な係数行列と非同次項を入力さえすれば、逆行列と解を計算してくれます。この関数、実際はガウス・ジョルダン法の手続きをC言語でいかにして実現するかは、次節以降に述べます。

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

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

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

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

この関数にデータを与える変数のことを引数と言います。先ほどの例で言うと、 C言語の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 = x+y;
	   return(zout);
	}

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

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

です。

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

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

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

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

	#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
平成19年8月21日