#逆行列を求める(C言語)

##方法
逆行列の求め方は様々な方法がありますが、
今回はGauss-Jordan法を使って求めていきます。

もし間違っていたら、コメントにて指摘していただけると幸いです。

###枢軸選択（ピボット選択）
プログラム上で数値を扱う上で気を付けなければならないことはたくさんありますが、その一つに0除算があります。

今回はそれを回避するために、ピボット選択を行います。

ピボット選択とは簡単に説明すると、
行列の対角成分から0をなくそうというものです。
またピボットの絶対値を大きなものすることで
割り算の誤差を小さくする効果もあります。
（※ピボットとは対角成分のことです。）

#プログラム
まずこれがプログラムの全文になります。

ROWSとCOLUMNSは事前にマクロで定義したもので、
行と列の要素数が入っています。

またaという逆行列を求めたい行列と
bという単位行列を作成しています。

static float a[ROWS][COLUMNS] = { {3,1,1},{5,1,3},{2,0,1} };
static float b[ROWS][COLUMNS] = { {1,0,0},{0,1,0},{0,0,1} };

printf("\na行列\n");
for (int i = 0; i < ROWS; i++) {
	for (int k = 0; k < COLUMNS; k++) {
		printf("%f ", a[i][k]);
	}
	printf("\n");
}
printf("\nb行列\n");
for (int i = 0; i < ROWS; i++) {
	for (int k = 0; k < COLUMNS; k++) {
		printf("%f ", b[i][k]);
	}
	printf("\n");
}


float max = a[0][0];
float tmp;

for (int i = 0; i < ROWS; i++) {
	
	//部分ピボット選択
	for (int k = i; k < ROWS; k++) {
		if (fabs(max) < fabs(a[k][i])) {
			for(int l = 0; l < ROWS;l++){
				tmp = a[i][l];
				a[i][l] = a[k][l];
				a[k][l] = tmp;

				tmp = b[i][l];
				b[i][l] = b[k][l];
				b[k][l] = tmp;

				max = a[k][i];
			}
		}
	}
	if (max == 0) {
		printf("解けません");
	}

	//行正規化
	float p = a[i][i];
	for (int l = i; l < COLUMNS; l++) 
		a[i][l] /= p;
	
	for (int l = 0; l < COLUMNS; l++) 
		b[i][l] /= p;


	for (int k = 0; k < COLUMNS; k++) {
		if (i != k) {
			float d = a[k][i];

			for (int l = i; l < ROWS; l++) 
				a[k][l] += -d * a[i][l];

			for (int l = 0; l < ROWS; l++)
				b[k][l] += -d * b[i][l]; 				
		}
	}
		
}

printf("\n変換後\na行列\n");
for (int i = 0; i < ROWS; i++) {
	for (int k = 0; k < COLUMNS; k++) {
		printf("%f ", a[i][k]);
	}
	printf("\n");
}
printf("\nb行列(逆行列)\n");
for (int i = 0; i < ROWS; i++) {
	for (int k = 0; k < COLUMNS; k++) {
		printf("%f ", b[i][k]);
	}
	printf("\n");
}

###ピボット選択実装部分

float max = a[0][0];
	float tmp;
	
	for (int i = 0; i < ROWS; i++) {
		
		//部分ピボット選択
		for (int k = i; k < ROWS; k++) {
			if (fabs(max) < fabs(a[k][i])) {
				for(int l = 0; l < ROWS;l++){
					tmp = a[i][l];
					a[i][l] = a[k][l];
					a[k][l] = tmp;

					tmp = b[i][l];
					b[i][l] = b[k][l];
					b[k][l] = tmp;

					max = a[k][i];
				}
			}
		}
		if (max == 0) {
			printf("解けません");
		}
                         :
                         :
                         :

先ほどの説明通り絶対値を比較して、
条件に合ったときのみ、行全体を入れ替えています。
（※厳密には部分ピボット選択といいます）

###行正規化

for (int i = 0; i < ROWS; i++) {
                :
                :
                :

float p = a[i][i];
		for (int l = i; l < COLUMNS; l++) 
			a[i][l] /= p;
		
		for (int l = 0; l < COLUMNS; l++) 
			b[i][l] /= p;
                :
                :
                :

ここでは単位行列になるように、
対角成分で行の正規化を行います。

そして次のコード部分で
a行列を単位行列に、bの行列には同じ操作を加えて逆行列を完成させていきます。

for (int k = 0; k < COLUMNS; k++) {
			if (i != k) {
				float d = a[k][i];

				for (int l = i; l < ROWS; l++) 
					a[k][l] += -d * a[i][l];

				for (int l = 0; l < ROWS; l++)
					b[k][l] += -d * b[i][l]; 				
			}
		}

結果的に出力されたbの行列が、
aの逆行列になっているはずです。

##まとめ
今回はGauss-Jordan法で逆行列を求めましたが、
似た方法としてガウス消去法やLU分解などでも求めることができます。
（ガウス消去法の方が効率が良らしい...)

また行列を多次元配列として定義すると、
関数の引数に渡すときに面倒だと感じたので、
1次元配列で定義することをお勧めします。

###参考文献
ピボット選択を行ったGauss-Jordan法