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

package main

import "fmt"
import "math"

const N = 4

// LR分解で固有値を求める
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 l [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 u [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++ {
// LU分解
decomp(a, &l, &u)
// 行列の積
multiply(u, l, &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.Fabs(a[i][j])
}
}
if (e < 0.00000000001) {
break
}
}

fmt.Println("\neigenvector")
disp_eigenvalue(a)
}
// LU分解
func decomp(a[N][N]float64, l *[N][N]float64, u *[N][N]float64) {
for i := 0; i < N; i++ {
for j := 0; j < N; j++ {
l[i][j] = 0.0
u[i][j] = 0.0
}
}

l[0][0] = 1.0
for j := 0; j < N; j++ {
u[0][j] = a[0][j]
}

for i := 1; i < N; i++ {
u[i][0] = 0.0
l[0][i] = 0.0
l[i][0] = a[i][0] / u[0][0]
}
for i := 1; i < N; i++ {
l[i][i] = 1.0
var t float64 = a[i][i]
for k := 0; k <= i; k++ {
t -= l[i][k] * u[k][i]
}
u[i][i] = t
for j := i + 1; j < N; j++ {
u[j][i] = 0.0
l[i][j] = 0.0
t       = a[j][i]
for k := 0; k <= i; k++ {
t -= l[j][k] * u[k][i]
}
l[j][i] = t / u[i][i]
t       = a[i][j]
for k := 0; k <= i; k++ {
t -= l[i][k] * u[k][j]
}
u[i][j] = t
}
}
}
// 行列の積
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 GO1103.go

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

z:\go>GO1103
1       8.6000000000    1.8444444444    4.6143790850    2.9411764706
2       9.6046511628    1.1012311902    4.8997514499    2.3943661972
3       9.8377723971    1.0107124514    4.9934604403    2.1580547112
4       9.9219788334    1.0010980897    5.0142271522    2.0626959248
5       9.9611291643    1.0001111458    5.0138771204    2.0248825695
6       9.9805335651    1.0000111820    5.0095550207    2.0099002322
7       9.9902522897    1.0000011216    5.0057991328    2.0039474559
8       9.9951218389    1.0000001123    5.0033019159    2.0015761328
9       9.9975597740    1.0000000112    5.0018103821    2.0006298327
10       9.9987795937    1.0000000011    5.0009686042    2.0002518010
(中略)
140      10.0000000000    5.0000000000    1.9999999973    1.0000000027
141      10.0000000000    5.0000000000    1.9999999986    1.0000000014
142      10.0000000000    5.0000000000    1.9999999993    1.0000000007
143      10.0000000000    5.0000000000    1.9999999997    1.0000000003
144      10.0000000000    5.0000000000    1.9999999998    1.0000000002
145      10.0000000000    5.0000000000    1.9999999999    1.0000000001
146      10.0000000000    5.0000000000    2.0000000000    1.0000000000
147      10.0000000000    5.0000000000    2.0000000000    1.0000000000
148      10.0000000000    5.0000000000    2.0000000000    1.0000000000
149      10.0000000000    5.0000000000    2.0000000000    1.0000000000

eigenvector
10.0000000000    5.0000000000    2.0000000000    1.0000000000