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

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

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

	return ind
end

				# データの設定
ct  = 200
eps = 1.0e-10
n   = 3
m   = 15

a   = Array.new(n)
b   = Array.new(n)
c   = Array.new(n)
v   = Array.new(n)
for i1 in 0 ... n
	a[i1] = Array.new(n)
	b[i1] = Array.new(n)
	c[i1] = Array.new(n)
	v[i1] = Array.new(n)
end

r  = Array.new(n)
v0 = Array.new(n)
v1 = Array.new(n)
v2 = Array.new(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 in 0 ... n
	for i2 in 0 ... n
		b[i1][i2] = 0.0
	end
	b[i1][i1] = 1.0
	v0[i1]    = 1.0
end
				# 計算
ind = power(0, n, m, ct, eps, a, b, c, r, v, v0, v1, v2)
				# 出力
if ind == 0
	printf("固有値が求まりませんでした！\n")
else
	for i1 in 0 ... ind
		printf("固有値 %f", r[i1])
		printf(" 固有ベクトル ")
		for i2 in 0 ... n
			printf(" %f", v[i1][i2])
		end
		printf("\n")
	end
end