1 準備

1.1 関数を使おう

ガウス・ジョルダン法のプログラムはメイン関数には書かないで，C言語の関数にすべき である²．この計算手法は，連立1次方程式の解，あるい は逆行列を求める場合，広範囲に利用できる．したがって，いつでも簡単に利用できる形 態として関数を利用するのである．良いプログラムを一度書けば，この種の問題にいつで も適用(再利用)することができるようになる．他の場所で連立1次方程式を解きたい場合， その関数をコールするで済む．プログラム作成の手間が大幅に省ける．また，ほかのプロ グラムを書く場合でも，それをコピーして，再利用できるので便利である．さらに，プロ グラムは分かりやすくなるメリットもある．C言語は関数の独立性が高いので，関数にし て再利用することが比較的簡単にできる．

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

図 1: ガウスジョルダン法を計算するC言語の関数のブラックボックス

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

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

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

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

いかなるプログラム言語でも，C言語の関数(mainではない)に対応するものがある．同じ ような処理がある場合，1つの独立した処理のブロックとしてまとめ，どこからでもコー ルすることができるようにすると便利である．このような，機能のブロックをサブルーチ ン，C言語では関数と呼ぶ．例えば，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関数 =================*/
	int main(void)
	{
	   double x, y, wa;

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

	   いろいろな処理

	   return 0;
	}

	/* ========== 足し算の関数 =================*/
	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関数 ==============================================*/
	int main(void){
	   double a[100][100], b[100];       // a[][]係数行列   b[]非同次項
	   double inv_a[100][100], x[100];   // inv_a[][]逆行列 x[]解

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

	   いろいろな処理

	   return 0;
	}

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

	いろいろな処理

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

	}

処理する関数の方は，一番最初のサイズを除いて書く必要がある．これは，例え ばz[10][20][30]の大きさの配列のz[3][5][7]というデータに アクセスする場合を考えれば分かる．このデータがあるアドレスは，

$$= \times(3\times20\times30+5\times30+7)$$ (1)

になる．これから，最初の配列の添え字のサイズは不要で，2番目の添え字からアドレス 計算に必要となることが分かる．2番目以降の配列の添え字のサイズを処理する関数側で 明示する必要が生じる．

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

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
著者: 山本昌志