連立線形方程式と逆行列

/****************************/
/* 逆行列の計算             */
/*      例： 2 1            */
/*           1 1            */
/*      coded by Y.Suganuma */
/****************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double w[][], eps;
		int ind, n = 2, m = 2;
	/*
	          データの設定
	*/
		w = new double [n][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);

		System.out.println("ind = " + ind);
		System.out.println("     " + w[0][2] + " " + w[0][3]);
		System.out.println("     " + w[1][2] + " " + w[1][3]);
	}

	/*******************************************************/
	/* 線形連立方程式を解く（逆行列を求める）              */
	/*      w : 方程式の左辺及び右辺                       */
	/*      n : 方程式の数                                 */
	/*      m : 方程式の右辺の列の数                       */
	/*      eps : 逆行列の存在を判定する規準               */
	/*      return : =0 : 正常                             */
	/*               =1 : 逆行列が存在しない               */
	/*******************************************************/
	static 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 = Math.abs(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);
	}
}