Option Explicit On
Option Strict On
Public Class Form1

    Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim n As Integer
        Dim a(,) As Double
        Dim ev() As Double
        Dim vec(,) As Double
        Dim dat As String
        Dim ar(,) As Double
        Dim c(,) As Double

        'n：固有方程式の次数(配列宣言はn-1)
        'ev(k)：固有値(大きい順に出力)
        'vec(i=0〜n-1,k)：固有ベクトル

        n = 3
        ReDim a(n - 1, n - 1)
        ReDim ev(n - 1)
        ReDim vec(n - 1, n - 1)
        ReDim ar(n - 1, n - 1)
        ReDim c(n - 1, n - 1)


        'Case-1
        a(0, 0) = 4 : a(0, 1) = 2 : a(0, 2) = 3
        a(1, 0) = 2 : a(1, 1) = 6 : a(1, 2) = 1
        a(2, 0) = 3 : a(2, 1) = 1 : a(2, 2) = 5

        'a(0, 0) = 3 : a(0, 1) = 2 : a(0, 2) = -2
        'a(1, 0) = 2 : a(1, 1) = 3 : a(1, 2) = -2
        'a(2, 0) = -2 : a(2, 1) = -2 : a(2, 2) = 8

        'a(0, 0) = 5 : a(0, 1) = 4 : a(0, 2) = 3 : a(0, 3) = 2 : a(0, 4) = 1
        'a(2, 0) = 4 : a(2, 1) = 4 : a(2, 2) = 3 : a(2, 3) = 2 : a(2, 4) = 1
        'a(3, 0) = 3 : a(3, 1) = 3 : a(3, 2) = 3 : a(3, 3) = 2 : a(3, 4) = 1
        'a(4, 0) = 2 : a(4, 1) = 2 : a(4, 2) = 2 : a(4, 3) = 2 : a(4, 4) = 1
        'a(5, 0) = 1 : a(5, 1) = 1 : a(5, 2) = 1 : a(5, 3) = 1 : a(5, 4) = 1

        For i = 0 To n - 1
            For j = 0 To n - 1
                ar(i, j) = a(i, j)
            Next j
        Next i


        Call JACOBI(n, a, ev, vec)

        Console.WriteLine("*** Eigen value ***")
        dat = String.Format("{0,16:0.000000e+000}", ev(0))
        For j = 1 To n - 1
            dat = dat & String.Format("{0,16:0.000000e+000}", ev(j))
        Next j
        Console.WriteLine(dat)
        Console.WriteLine("*** Eigen vector ***")
        For i = 0 To n - 1
            dat = String.Format("{0,16:0.000000e+000}", vec(i, 0))
            For j = 1 To n - 1
                dat = dat & String.Format("{0,16:0.000000e+000}", vec(i, j))
            Next j
            Console.WriteLine(dat)
        Next i

        'Case-1 出力結果
        '*** Eigen value ***
        '9.000000e+000   4.732051e+000   1.267949e+000
        '*** Eigen vector ***
        '5.773503e-001   2.113249e-001   7.886751e-001
        '5.773503e-001  -7.886751e-001  -2.113249e-001
        '5.773503e-001   5.773503e-001  -5.773503e-001

        '検算
        For i = 0 To n - 1
            For j = 0 To n - 1
                c(i, j) = 0.0
                For k = 0 To n - 1
                    c(i, j) = c(i, j) + ar(i, k) * vec(k, j)
                Next k
            Next j
        Next i
        Console.WriteLine("*** 検算 ***")
        For i = 0 To n - 1
            dat = String.Format("{0,16:0.000000e+000}", c(i, 0))
            For j = 1 To n - 1
                dat = dat & String.Format("{0,16:0.000000e+000}", c(i, j))
            Next j
            dat = dat & "  "
            dat = dat & String.Format("{0,16:0.000000e+000}", vec(i, 0) * ev(0))
            For j = 1 To n - 1
                dat = dat & String.Format("{0,16:0.000000e+000}", vec(i, j) * ev(j))
            Next j
            Console.WriteLine(dat)
        Next i

    End Sub

    Private Sub JACOBI(ByVal n As Integer, ByRef a(,) As Double, ByRef ev() As Double, ByRef vec(,) As Double)
        'Jacobi法による固有値・固有ベクトル計算（奥村先生）
        '固有値ev(k),固有ベクトルvec(k,i)
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim s As Double
        Dim c As Double
        Dim aoff As Double
        Dim eps As Double
        Dim t As Double
        Dim aa As Double
        Dim aj As Double
        Dim ak As Double
        Dim vj As Double
        Dim vk As Double
        Dim vmax As Double
        Dim im As Integer
        Dim bb As Double

        For i = 0 To n - 1
            For j = 0 To n - 1
                vec(i, j) = 0.0
            Next j
            vec(i, i) = 1.0
        Next i

        aoff = 0.0
        For j = 0 To n - 2
            For k = j + 1 To n - 1
                aoff = aoff + a(j, k) * a(j, k)
            Next k
        Next j
        eps = aoff * 2.0
        For j = 0 To n - 1
            eps = eps + a(j, j) * a(j, j)
        Next j
        eps = eps * 0.00000000005

        Do While aoff > eps
            For j = 0 To n - 2
                For k = j + 1 To n - 1
                    If a(j, k) = 0.0 Then
                        t = 0.0
                    Else
                        aa = (a(k, k) - a(j, j)) / (2.0 * a(j, k))
                    End If
                    If aa >= 0 Then
                        t = 1.0 / (aa + Math.Sqrt(1.0 + aa * aa))
                    Else
                        t = 1.0 / (aa - Math.Sqrt(1.0 + aa * aa))
                    End If
                    c = 1.0 / Math.Sqrt(1.0 + t * t)
                    s = t * c
                    a(j, j) = a(j, j) - t * a(j, k)
                    a(k, k) = a(k, k) + t * a(j, k)
                    a(j, k) = 0.0
                    For i = 0 To j - 1
                        aj = a(i, j)
                        ak = a(i, k)
                        a(i, j) = aj * c - ak * s
                        a(i, k) = aj * s + ak * c
                    Next i

                    For i = j + 1 To k - 1
                        aj = a(j, i)
                        ak = a(i, k)
                        a(j, i) = aj * c - ak * s
                        a(i, k) = aj * s + ak * c
                    Next i

                    For i = k + 1 To n - 1
                        aj = a(j, i)
                        ak = a(k, i)
                        a(j, i) = aj * c - ak * s
                        a(k, i) = aj * s + ak * c
                    Next i

                    For i = 0 To n - 1
                        vj = vec(i, j)
                        vk = vec(i, k)
                        vec(i, j) = vj * c - vk * s
                        vec(i, k) = vj * s + vk * c
                    Next i

                Next k
            Next j

            aoff = 0.0
            For j = 0 To n - 2
                For k = j + 1 To n - 1
                    aoff = aoff + a(j, k) * a(j, k)
                Next k
            Next j
        Loop

        For k = 0 To n - 1
            ev(k) = a(k, k)
        Next k

        For k = 0 To n - 2
            vmax = -1.0E+30
            For i = k To n - 1
                If ev(i) >= vmax Then
                    im = i : vmax = ev(i)
                End If
            Next i
            aa = ev(k) : ev(k) = ev(im) : ev(im) = aa
            For i = 0 To n - 1
                bb = vec(i, k)
                vec(i, k) = vec(i, im)
                vec(i, im) = bb
            Next i
        Next k

    End Sub

End Class