最大固有値と固有ベクトル（べき乗法）

/****************************************/
/* 最大固有値と固有ベクトル（べき乗法） */
/*      coded by Y.Suganuma             */
/****************************************/
import java.io.*;

public class Test {
	public static void main(String args[]) throws IOException
	{
		double A[][], B[][], C[][], r[], v[][], v0[], v1[], v2[], eps;
		int i1, i2, ind, ct, m, n;
					// データの設定
		ct  = 200;
		eps = 1.0e-10;
		n   = 3;
		m   = 15;

		A   = new double [n][n];
		B   = new double [n][n];
		C   = new double [n][n];
		v   = new double [n][n];

		r  = new double [n];
		v0 = new double [n];
		v1 = new double [n];
		v2 = new double [n];

		A[0][0] = 11.0;
		A[0][1] = 6.0;
		A[0][2] = -2.0;

		A[1][0] = -2.0;
		A[1][1] = 18.0;
		A[1][2] = 1.0;

		A[2][0] = -12.0;
		A[2][1] = 24.0;
		A[2][2] = 13.0;

		for (i1 = 0; i1 < n; i1++) {
			for (i2 = 0; i2 < n; i2++)
				B[i1][i2] = 0.0;
			B[i1][i1] = 1.0;
			v0[i1]    = 1.0;
		}
					// 計算
		ind = power(0, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
					// 出力
		if (ind == 0)
			System.out.println("固有値が求まりませんでした！");
		else {
			for (i1 = 0; i1 < ind; i1++) {
				System.out.print("固有値 " + r[i1]);
				System.out.print(" 固有ベクトル ");
				for (i2 = 0; i2 < n; i2++)
					System.out.print(" " + v[i1][i2]);
				System.out.println();
			}
		}
	}

	/****************************************/
	/* 最大固有値と固有ベクトル（べき乗法） */
	/*      i : 何番目の固有値かを示す      */
	/*      n : 次数                        */
	/*      m : 丸め誤差の修正回数          */
	/*      ct : 最大繰り返し回数           */
	/*      eps : 収束判定条件              */
	/*      A : 対象とする行列              */
	/*      B : nxnの行列，最初は，単位行列 */
	/*      C : 作業域，nxnの行列           */
	/*      r : 固有値                      */
	/*      v : 各行が固有ベクトル（nxn）   */
	/*      v0 : 固有ベクトルの初期値       */
	/*      v1,v2 : 作業域（n次元ベクトル） */
	/*      return : 求まった固有値の数     */
	/*      coded by Y.Suganuma             */
	/****************************************/
	static int power(int i, int n, int m, int ct, double eps, double A[][], double B[][],
	                 double C[][], double r[], double v[][], double v0[], double v1[], double v2[])
	{
		double l1, l2, x;
		int i1, i2, i3, ind, k, k1;
					// 初期設定
		ind = i;
		k   = 0;
		l1  = 0.0;
		for (i1 = 0; i1 < n; i1++) {
			l1     += v0[i1] * v0[i1];
			v1[i1]  = v0[i1];
		}
		l1 = Math.sqrt(l1);
					// 繰り返し計算
		while (k < ct) {
						// 丸め誤差の修正
			if (k%m == 0) {
				l2 = 0.0;
				for (i1 = 0; i1 < n; i1++) {
					v2[i1] = 0.0;
					for (i2 = 0; i2 < n; i2++)
						v2[i1] += B[i1][i2] * v1[i2];
					l2 += v2[i1] * v2[i1];
				}
				l2 = Math.sqrt(l2);
				for (i1 = 0; i1 < n; i1++)
					v1[i1] = v2[i1] / l2;
			}
						// 次の近似
			l2 = 0.0;
			for (i1 = 0; i1 < n; i1++) {
				v2[i1] = 0.0;
				for (i2 = 0; i2 < n; i2++)
					v2[i1] += A[i1][i2] * v1[i2];
				l2 += v2[i1] * v2[i1];
			}
			l2 = Math.sqrt(l2);
			for (i1 = 0; i1 < n; i1++)
				v2[i1] /= l2;
						// 収束判定
							// 収束した場合
			if (Math.abs((l2-l1)/l1) < eps) {
				k1 = -1;
				for (i1 = 0; i1 < n && k1 < 0; i1++) {
					if (Math.abs(v2[i1]) > 0.001) {
						k1 = i1;
						if (v2[k1]*v1[k1] < 0.0)
							l2 = -l2;
					}
				}
				k    = ct;
				r[i] = l2;
				for (i1 = 0; i1 < n; i1++)
					v[i][i1] = v2[i1];
				if (i == n-1)
					ind = i + 1;
				else {
					for (i1 = 0; i1 < n; i1++) {
						for (i2 = 0; i2 < n; i2++) {
							C[i1][i2] = 0.0;
							for (i3 = 0; i3 < n; i3++) {
								x          = (i1 == i3) ? A[i1][i3] - l2 : A[i1][i3];
								C[i1][i2] += x * B[i3][i2];
							}
						}
					}
					for (i1 = 0; i1 < n; i1++) {
						for (i2 = 0; i2 < n; i2++)
							B[i1][i2] = C[i1][i2];
					}
					for (i1 = 0; i1 < n; i1++) {
						v1[i1] = 0.0;
						for (i2 = 0; i2 < n; i2++)
							v1[i1] += B[i1][i2] * v0[i2];
					}
					for (i1 = 0; i1 < n; i1++)
						v0[i1] = v1[i1];
					ind = power(i+1, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
				}
			}
							// 収束しない場合
			else {
				for (i1 = 0; i1 < n; i1++)
					v1[i1] = v2[i1];
				l1 = l2;
				k++;
			}
		}

		return ind;
	}
}