1 準備

1.1 関数を使おう

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

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

**図 1:** ガウスジョルダン法を計算する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と同じ)．

従って，図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
2005-11-18