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

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

#************************************************************/
# 実対称行列の固有値・固有ベクトル（ヤコビ法）              */
#      n : 次数                                             */
#      ct : 最大繰り返し回数                                */
#      eps : 収束判定条件                                   */
#      a : 対象とする行列                                   */
#      a1, a2 : 作業域（nxnの行列），a1の対角要素が固有値   */
#      x1, x2 : 作業域（nxnの行列），x1の各列が固有ベクトル */
#      return : =0 : 正常                                   */
#               =1 : 収束せず                               */
#      coded by Y.Suganuma                                  */
#************************************************************/
def Jacobi(n, ct, eps, a, a1, a2, x1, x2)
					# 初期設定
	k   = 0
	ind = 1
	p   = 0
	q   = 0
	for i1 in 0 ... n
		for i2 in 0 ... n
			a1[i1][i2] = a[i1][i2]
			x1[i1][i2] = 0.0
		end
		x1[i1][i1] = 1.0
	end
					# 計算
	while ind > 0 && k < ct
						# 最大要素の探索
		max = 0.0
		for i1 in 0 ... n
			for i2 in 0 ... n
				if i2 != i1
					if a1[i1][i2].abs() > max
						max = a1[i1][i2].abs()
						p   = i1
						q   = i2
					end
				end
			end
		end
						# 収束判定
							# 収束した
		if max < eps
			ind = 0
							# 収束しない
		else
								# 準備
			s  = -a1[p][q]
			t  = 0.5 * (a1[p][p] - a1[q][q])
			v  = t.abs() / Math.sqrt(s * s + t * t)
			sn = Math.sqrt(0.5 * (1.0 - v))
			if s*t < 0.0
				sn = -sn
			end
			cs = Math.sqrt(1.0 - sn * sn)
								# akの計算
			for i1 in 0 ... n
				if i1 == p
					for i2 in 0 ... n
						if i2 == p
							a2[p][p] = a1[p][p] * cs * cs + a1[q][q] * sn * sn -
                                       2.0 * a1[p][q] * sn * cs
						elsif i2 == q
							a2[p][q] = 0.0
						else
							a2[p][i2] = a1[p][i2] * cs - a1[q][i2] * sn
						end
					end
				elsif i1 == q
					for i2 in 0 ... n
						if (i2 == q)
							a2[q][q] = a1[p][p] * sn * sn + a1[q][q] * cs * cs +
                                       2.0 * a1[p][q] * sn * cs
						elsif i2 == p
							a2[q][p] = 0.0
						else
							a2[q][i2] = a1[q][i2] * cs + a1[p][i2] * sn
						end
					end
				else
					for i2 in 0 ... n
						if i2 == p
							a2[i1][p] = a1[i1][p] * cs - a1[i1][q] * sn
						elsif i2 == q
							a2[i1][q] = a1[i1][q] * cs + a1[i1][p] * sn
						else
							a2[i1][i2] = a1[i1][i2]
						end
					end
				end
			end
								# xkの計算
			for i1 in 0 ... n
				for i2 in 0 ... n
					if i2 == p
						x2[i1][p] = x1[i1][p] * cs - x1[i1][q] * sn
					elsif i2 == q
						x2[i1][q] = x1[i1][q] * cs + x1[i1][p] * sn
					else
						x2[i1][i2] = x1[i1][i2]
					end
				end
			end
								# 次のステップへ
			k += 1
			for i1 in 0 ... n
				for i2 in 0 ... n
					a1[i1][i2] = a2[i1][i2]
					x1[i1][i2] = x2[i1][i2]
				end
			end
		end
	end

	return ind
end

					# データの設定
ct   = 1000
eps  = 1.0e-10
n    = 3
a    = Array.new(n)
a1   = Array.new(n)
a2   = Array.new(n)
x1   = Array.new(n)
x2   = Array.new(n)
for i1 in 0 ... n
	a[i1]  = Array.new(n)
	a1[i1] = Array.new(n)
	a2[i1] = Array.new(n)
	x1[i1] = Array.new(n)
	x2[i1] = Array.new(n)
end

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 in 0 ... n
		printf(" %f", a1[i1][i1])
	end
	printf("\n固有ベクトル（各列が固有ベクトル）\n")
	for i1 in 0 ... n
		for i2 in 0 ... n
			printf(" %f", x1[i1][i2])
		end
		printf("\n")
	end
end