実対称行列の固有値・固有ベクトル（ヤコビ法）

<?php

/************************************************/
/* 実対称行列の固有値・固有ベクトル（ヤコビ法） */
/*      coded by Y.Suganuma                     */
/************************************************/

					// データの設定
	$ct   = 1000;
	$eps  = 1.0e-10;
	$n    = 3;
	$A    = array($n);
	$A1   = array($n);
	$A2   = array($n);
	$X1   = array($n);
	$X2   = array($n);
	for ($i1 = 0; $i1 < $n; $i1++) {
		$A[$i1]  = array($n);
		$A1[$i1] = array($n);
		$A2[$i1] = array($n);
		$X1[$i1] = array($n);
		$X2[$i1] = array($n);
	}

	$A[0][0] = 1.0;
	$A[0][1] = 0.0;
	$A[0][2] = -1.0;

	$A[1][0] = 0.0;
	$A[1][1] = 1.0;
	$A[1][2] = -1.0;

	$A[2][0] = -1.0;
	$A[2][1] = -1.0;
	$A[2][2] = 2.0;
					// 計算
	$ind = Jacobi($n, $ct, $eps, $A, $A1, $A2, $X1, $X2);
					// 出力
	if ($ind > 0)
		printf("収束しませんでした！\n");
	else {
		printf("固有値\n");
		for ($i1 = 0; $i1 < $n; $i1++)
			printf(" %f", $A1[$i1][$i1]);
		printf("\n固有ベクトル（各列が固有ベクトル）\n");
		for ($i1 = 0; $i1 < $n; $i1++) {
			for ($i2 = 0; $i2 < $n; $i2++)
				printf(" %f", $X1[$i1][$i2]);
			printf("\n");
		}
	}

/*************************************************************/
/* 実対称行列の固有値・固有ベクトル（ヤコビ法）              */
/*      n : 次数                                             */
/*      ct : 最大繰り返し回数                                */
/*      eps : 収束判定条件                                   */
/*      A : 対象とする行列                                   */
/*      A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値   */
/*      X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル */
/*      return : =0 : 正常                                   */
/*               =1 : 収束せず                               */
/*      coded by Y.Suganuma                                  */
/*************************************************************/
function Jacobi($n, $ct, $eps, $A, &$A1, $A2, &$X1, $X2)
{
					// 初期設定
	$k   = 0;
	$ind = 1;
	$p   = 0;
	$q   = 0;

	for ($i1 = 0; $i1 < $n; $i1++) {
		for ($i2 = 0; $i2 < $n; $i2++) {
			$A1[$i1][$i2] = $A[$i1][$i2];
			$X1[$i1][$i2] = 0.0;
		}
		$X1[$i1][$i1] = 1.0; 
	}
					// 計算
	while ($ind > 0 && $k < $ct) {
						// 最大要素の探索
		$max = 0.0;
		for ($i1 = 0; $i1 < $n; $i1++) {
			for ($i2 = 0; $i2 < $n; $i2++) {
				if ($i2 != $i1) {
					if (abs($A1[$i1][$i2]) > $max) {
						$max = abs($A1[$i1][$i2]);
						$p   = $i1;
						$q   = $i2;
					}
				}
			}
		}
						// 収束判定
							// 収束した
		if ($max < $eps)
			$ind = 0;
							// 収束しない
		else {
								// 準備
			$s  = -$A1[$p][$q];
			$t  = 0.5 * ($A1[$p][$p] - $A1[$q][$q]);
			$v  = abs($t) / sqrt($s * $s + $t * $t);
			$sn = sqrt(0.5 * (1.0 - $v));
			if ($s*$t < 0.0)
				$sn = -$sn;
			$cs = sqrt(1.0 - $sn * $sn);
								// Akの計算
			for ($i1 = 0; $i1 < $n; $i1++) {
				if ($i1 == $p) {
					for ($i2 = 0; $i2 < $n; $i2++) {
						if ($i2 == $p)
							$A2[$p][$p] = $A1[$p][$p] * $cs * $cs + $A1[$q][$q] * $sn * $sn -
                                          2.0 * $A1[$p][$q] * $sn * $cs;
						else if ($i2 == $q)
							$A2[$p][$q] = 0.0;
						else
							$A2[$p][$i2] = $A1[$p][$i2] * $cs - $A1[$q][$i2] * $sn;
					}
				}
				else if ($i1 == $q) {
					for ($i2 = 0; $i2 < $n; $i2++) {
						if ($i2 == $q)
							$A2[$q][$q] = $A1[$p][$p] * $sn * $sn + $A1[$q][$q] * $cs * $cs +
                                          2.0 * $A1[$p][$q] * $sn * $cs;
						else if ($i2 == $p)
							$A2[$q][$p] = 0.0;
						else
							$A2[$q][$i2] = $A1[$q][$i2] * $cs + $A1[$p][$i2] * $sn;
					}
				}
				else {
					for ($i2 = 0; $i2 < $n; $i2++) {
						if ($i2 == $p)
							$A2[$i1][$p] = $A1[$i1][$p] * $cs - $A1[$i1][$q] * $sn;
						else if ($i2 == $q)
							$A2[$i1][$q] = $A1[$i1][$q] * $cs + $A1[$i1][$p] * $sn;
						else
							$A2[$i1][$i2] = $A1[$i1][$i2];
					}
				}
			}
								// Xkの計算
			for ($i1 = 0; $i1 < $n; $i1++) {
				for ($i2 = 0; $i2 < $n; $i2++) {
					if ($i2 == $p)
						$X2[$i1][$p] = $X1[$i1][$p] * $cs - $X1[$i1][$q] * $sn;
					else if ($i2 == $q)
						$X2[$i1][$q] = $X1[$i1][$q] * $cs + $X1[$i1][$p] * $sn;
					else
						$X2[$i1][$i2] = $X1[$i1][$i2];
				}
			}
								// 次のステップへ
			$k++;
			for ($i1 = 0; $i1 < $n; $i1++) {
				for ($i2 = 0; $i2 < $n; $i2++) {
					$A1[$i1][$i2] = $A2[$i1][$i2];
					$X1[$i1][$i2] = $X2[$i1][$i2];
				}
			}
		}
	}

	return $ind;
}

?>