連立線形方程式と逆行列

/****************************/
/* 逆行列の計算             */
/*      例： 2 1            */
/*           1 1            */
/*      coded by Y.Suganuma */
/****************************/
#include <stdio.h>

int gauss(double **, int, int, double);

int main ()
{
	double **w, eps;
	int i1, ind, n = 2, m = 2;
/*
          データの設定
*/
	w = new double * [n];
	for (i1 = 0; i1 < n; i1++)
		w[i1] = new double [n+m];

	eps     = 1.0e-10;
	w[0][0] = 2.0;
	w[0][1] = 1.0;
	w[0][2] = 1.0;
	w[0][3] = 0.0;
	w[1][0] = 1.0;
	w[1][1] = 1.0;
	w[1][2] = 0.0;
	w[1][3] = 1.0;
/*
          実行と出力
*/
	ind = gauss(w, n, m, eps);

	printf("ind = %d\n", ind);
	printf("     %f %f\n", w[0][2], w[0][3]);
	printf("     %f %f\n", w[1][2], w[1][3]);

	return 0;
}

/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める）              */
/*      w : 方程式の左辺及び右辺                       */
/*      n : 方程式の数                                 */
/*      m : 方程式の右辺の列の数                       */
/*      eps : 逆行列の存在を判定する規準               */
/*      return : =0 : 正常                             */
/*               =1 : 逆行列が存在しない               */
/*******************************************************/
#include <math.h>

int gauss(double **w, int n, int m, double eps)
{
	double y1, y2;
	int ind = 0, nm, m1, m2, i1, i2, i3;

	nm = n + m;

	for (i1 = 0; i1 < n && ind == 0; i1++) {

		y1 = .0;
		m1 = i1 + 1;
		m2 = 0;
					// ピボット要素の選択
		for (i2 = i1; i2 < n; i2++) {
			y2 = fabs(w[i2][i1]);
			if (y1 < y2) {
				y1 = y2;
				m2 = i2;
			}
		}
					// 逆行列が存在しない
		if (y1 < eps)
			ind = 1;
					// 逆行列が存在する
		else {
						// 行の入れ替え
			for (i2 = i1; i2 < nm; i2++) {
				y1        = w[i1][i2];
				w[i1][i2] = w[m2][i2];
				w[m2][i2] = y1;
			}
						// 掃き出し操作
			y1 = 1.0 / w[i1][i1];

			for (i2 = m1; i2 < nm; i2++)
				w[i1][i2] *= y1;

			for (i2 = 0; i2 < n; i2++) {
				if (i2 != i1) {
					for (i3 = m1; i3 < nm; i3++)
						w[i2][i3] -= w[i2][i1] * w[i1][i3];
				}
			}
		}
	}

	return(ind);
}