最大固有値と固有ベクトル(べき乗法)

/****************************************/
/* 最大固有値と固有ベクトル(べき乗法) */
/*      coded by Y.Suganuma             */
/****************************************/
import java.awt.*;
import java.awt.event.*;
import java.applet.*;
import java.util.StringTokenizer;

public class Power extends Applet implements ActionListener {
	TextField order, mod, max, ini;
	TextArea ta1, ta2;
	Button bt;
	Font f = new Font("TimesRoman", Font.BOLD, 20);

	/************/
	/* 初期設定 */
	/************/
	public void init()
	{
					// レイアウト,背景色,フォント
		setLayout(new BorderLayout(5, 5));
		setBackground(new Color(225, 255, 225));
		setFont(f);
					// 上のパネル
		Panel pn1 = new Panel();
		add(pn1, BorderLayout.NORTH);

		pn1.add(new Label("次数:"));
		order = new TextField("3", 3);
		pn1.add(order);

		pn1.add(new Label("修正:"));
		mod = new TextField("15", 3);
		pn1.add(mod);

		pn1.add(new Label("最大繰返し:"));
		max = new TextField("200", 5);
		pn1.add(max);

		pn1.add(new Label("初期値:"));
		ini = new TextField("1 1 1", 20);
		pn1.add(ini);

		pn1.add(new Label(" "));
		bt = new Button("OK");
		bt.setBackground(Color.pink);
		bt.addActionListener(this);
		pn1.add(bt);
					// 中央のパネル(問題)
		Panel pn2 = new Panel();
		add(pn2, BorderLayout.CENTER);

		ta1 = new TextArea("11 6 -2\n-2 18 1\n-12 24 13", 10, 30);
		pn2.add(ta1);
					// 下のパネル(答え)
		Panel pn3 = new Panel();
		add(pn3, BorderLayout.SOUTH);

		ta2 = new TextArea(10, 70);
		pn3.add(ta2);
	}

	/******************************/
	/* ボタンが押されたときの処理 */
	/******************************/
	public void actionPerformed(ActionEvent e)
	{
		if (e.getSource() == bt) {
			double A[][], B[][], C[][], r[], v[][], v0[], v1[], v2[], eps;
			int ind, ct, m, n;
					// 問題の設定
			eps = 1.0e-10;
			n   = Integer.parseInt(order.getText());
			m   = Integer.parseInt(mod.getText());
			ct  = Integer.parseInt(max.getText());
			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];

			String s = ini.getText();
			StringTokenizer str = new StringTokenizer(s, " \n");
			int k = 0;
			while (str.hasMoreTokens()) {
				v0[k] = Double.parseDouble(str.nextToken());
				k++;
			}

			s   = ta1.getText();
			str = new StringTokenizer(s, " \n");
			k   = 0;
			while (str.hasMoreTokens()) {
				A[k/n][k%n] = Double.parseDouble(str.nextToken());
				k++;
			}

			for (int i1 = 0; i1 < n; i1++) {
				for (int i2 = 0; i2 < n; i2++)
					B[i1][i2] = 0.0;
				B[i1][i1] = 1.0;
			}
					// 計算
			ind = power(0, n, m, ct, eps, A, B, C, r, v, v0, v1, v2);
					// 結果
			if (ind == 0) {
				ta2.setForeground(Color.red);
				ta2.setText(" 解を求めることができません");
			}
			else {
				ta2.setForeground(Color.black);
				ta2.setText("");
				for (int i1 = 0; i1 < ind; i1++) {
					ta2.append("固有値: " + String.format("%.5f",r[i1]) + "\n");
					ta2.append("  固有ベクトル:");
					for (int i2 = 0; i2 < n; i2++)
						ta2.append(" " + String.format("%.5f",v[i1][i2]));
					ta2.append("\n");
				}
			}
		}
	}

	/****************************************/
	/* 最大固有値と固有ベクトル(べき乗法) */
	/*      i : 何番目の固有値かを示す      */
	/*      n : 次数                        */
	/*      m : 丸め誤差の修正回数          */
	/*      ct : 最大繰り返し回数           */
	/*      eps : 収束判定条件              */
	/*      A : 対象とする行列              */
	/*      B : nxnの行列,最初は,単位行列 */
	/*      C : 作業域,nxnの行列           */
	/*      r : 固有値                      */
	/*      v : 各行が固有ベクトル(nxn)   */
	/*      v0 : 固有ベクトルの初期値       */
	/*      v1,v2 : 作業域(n次元ベクトル) */
	/*      return : 求まった固有値の数     */
	/*      coded by Y.Suganuma             */
	/****************************************/
	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;
	}
}