連立線形方程式と逆行列

<?php

/****************************/
/* 逆行列の計算             */
/*      例： 2 1            */
/*           1 1            */
/*      coded by Y.Suganuma */
/****************************/
/*
          データの設定
*/

	$eps  = 1.0e-10;
	$w    = array(2);
	$w[0] = array(2.0, 1.0, 1.0, 0.0);
	$w[1] = array(1.0, 1.0, 0.0, 1.0);
/*
          実行と出力
*/
	$ind = gauss($w, 2, 2, $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]);

/******************************************/
/* 線形連立方程式を解く（逆行列を求める） */
/*      $w : 方程式の左辺及び右辺         */
/*      $n : 方程式の数                   */
/*      $m : 方程式の右辺の列の数         */
/*      $eps : 逆行列の存在を判定する規準 */
/*      return : =0 : 正常                */
/*               =1 : 逆行列が存在しない  */
/*******************************************/
function gauss(&$w, $n, $m, $eps)
{
	$ind = 0;
	$nm  = $n + $m;

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

		$y1 = 0.0;
		$m1 = $i1 + 1;
		$m2 = 0;
				// ピボット要素の選択
		for ($i2 = $i1; $i2 < $n; $i2++) {
			$y2 = 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;
}

?>