# Go 言語で固有値を求める (QR分解)

package main

import "fmt"
import "math"

const N = 4

// QR分解で固有値を求める
func main() {
var a [N][N]float64 = [N][N]float64{{5.0, 4.0, 1.0, 1.0},
{4.0, 5.0, 1.0, 1.0},
{1.0, 1.0, 4.0, 2.0},
{1.0, 1.0, 2.0, 4.0}}
var q [N][N]float64 = [N][N]float64{{0.0 ,0.0 ,0.0 ,0.0},
{0.0 ,0.0 ,0.0 ,0.0},
{0.0 ,0.0 ,0.0 ,0.0},
{0.0 ,0.0 ,0.0 ,0.0}}
var r [N][N]float64 = [N][N]float64{{0.0 ,0.0 ,0.0 ,0.0},
{0.0 ,0.0 ,0.0 ,0.0},
{0.0 ,0.0 ,0.0 ,0.0},
{0.0 ,0.0 ,0.0 ,0.0}}

for k := 1; k < 200; k++ {
// QR分解
decomp(a, &q, &r)
// 行列の積
multiply(r, q, &a)
// 対角要素を表示
fmt.Printf("%3d\t", k)
disp_eigenvalue(a)

// 収束判定
var e float64 = 0.0
for i := 1; i < N; i++ {
for j := 0; j < i; j++ {
e += math.Abs(a[i][j])
}
}
if (e < 0.00000000001) {
break
}
}

fmt.Println("\neigenvalue")
disp_eigenvalue(a)
}
// QR分解
func decomp(a[N][N]float64, q *[N][N]float64, r *[N][N]float64) {
var x = []float64 {0.0 ,0.0 ,0.0 ,0.0}

for k := 0; k < N; k++ {
for i := 0; i < N; i++ {
x[i] = a[i][k]
}

for j := 0; j < k; j++ {
var t float64 = 0.0
for i := 0; i < N; i++ {
t += a[i][k] * q[i][j]
}
r[j][k] = t
r[k][j] = 0.0
for i := 0; i < N; i++ {
x[i] -= t * q[i][j]
}
}

var s float64 = 0.0
for i := 0; i < N; i++ {
s += x[i] * x[i]
}
r[k][k] = math.Sqrt(s)
for i := 0; i < N; i++ {
q[i][k] = x[i] / r[k][k]
}
}
}
// 行列の積
func multiply(a[N][N]float64, b[N][N]float64, c *[N][N]float64) {
for i := 0; i < N; i++ {
for j := 0; j < N; j++ {
var s float64 = 0.0
for k := 0; k < N; k++ {
s += a[i][k] * b[k][j]
}
c[i][j] = s
}
}
}
// 対角要素を表示
func disp_eigenvalue(a[N][N]float64) {
for i := 0; i < N; i++ {
fmt.Printf("%14.10f\t", a[i][i])
}
fmt.Println("")
}
z:\go>8g GO1104.go

z:\go>8l -o GO1104.exe GO1104.8

z:\go>GO1104
1       9.6046511628    1.1012311902    4.8997514499    2.3943661972
2       9.9219788334    1.0010980897    5.0142271522    2.0626959248
3       9.9805335651    1.0000111820    5.0095550207    2.0099002322
4       9.9951218389    1.0000001123    5.0033019159    2.0015761328
5       9.9987795937    1.0000000011    5.0009686042    2.0002518010
6       9.9996948428    1.0000000000    5.0002648858    2.0000402713
7       9.9999237072    1.0000000000    5.0000698501    2.0000064427
8       9.9999809266    1.0000000000    5.0000180426    2.0000010308
9       9.9999952316    1.0000000000    5.0000046034    2.0000001649
10       9.9999988079    1.0000000000    5.0000011657    2.0000000264
(中略)
85      10.0000000000    5.0000000000    2.0000000000    1.0000000000
86      10.0000000000    5.0000000000    2.0000000000    1.0000000000
87      10.0000000000    5.0000000000    2.0000000000    1.0000000000
88      10.0000000000    5.0000000000    2.0000000000    1.0000000000
89      10.0000000000    5.0000000000    2.0000000000    1.0000000000
90      10.0000000000    5.0000000000    2.0000000000    1.0000000000
91      10.0000000000    5.0000000000    2.0000000000    1.0000000000
92      10.0000000000    5.0000000000    2.0000000000    1.0000000000
93      10.0000000000    5.0000000000    2.0000000000    1.0000000000
94      10.0000000000    5.0000000000    2.0000000000    1.0000000000

eigenvalue
10.0000000000    5.0000000000    2.0000000000    1.0000000000