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

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

public class Test {
	public static void main(String args[]) throws IOException
	{
		double A[][], A1[][], A2[][], X1[][], X2[][], eps;
		int i1, i2, ind, ct, n;
					// データの設定
		ct   = 1000;
		eps  = 1.0e-10;
		n    = 3;
		A    = new double [n][n];
		A1   = new double [n][n];
		A2   = new double [n][n];
		X1   = new double [n][n];
		X2   = new double [n][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)
			System.out.println("収束しませんでした！");
		else {
			System.out.println("固有値");
			for (i1 = 0; i1 < n; i1++)
				System.out.print(" " + A1[i1][i1]);
			System.out.println("\n固有ベクトル（各列が固有ベクトル）");
			for (i1 = 0; i1 < n; i1++) {
				for (i2 = 0; i2 < n; i2++)
					System.out.print(" " + X1[i1][i2]);
				System.out.println();
			}
		}
	}

	/*************************************************************/
	/* 実対称行列の固有値・固有ベクトル（ヤコビ法）              */
	/*      n : 次数                                             */
	/*      ct : 最大繰り返し回数                                */
	/*      eps : 収束判定条件                                   */
	/*      A : 対象とする行列                                   */
	/*      A1, A2 : 作業域（nxnの行列），A1の対角要素が固有値   */
	/*      X1, X2 : 作業域（nxnの行列），X1の各列が固有ベクトル */
	/*      return : =0 : 正常                                   */
	/*               =1 : 収束せず                               */
	/*      coded by Y.Suganuma                                  */
	/*************************************************************/
	static int Jacobi(int n, int ct, double eps, double A[][], double A1[][], double A2[][],
                  double X1[][], double X2[][])
	{
		double max, s, t, v, sn, cs;
		int i1, i2, 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 (Math.abs(A1[i1][i2]) > max) {
							max = Math.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  = Math.abs(t) / Math.sqrt(s * s + t * t);
				sn = Math.sqrt(0.5 * (1.0 - v));
				if (s*t < 0.0)
					sn = -sn;
				cs = Math.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;
	}
}