Public Class Form1

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

        'a(1, 1) = 0.35355339 : a(1, 2) = -0.35355339 : a(1, 3) = 0 : a(1, 4) = 0
        'a(2, 1) = -0.35355339 : a(2, 2) = 1.35355339 : a(2, 3) = 0 : a(2, 4) = -1
        'a(3, 1) = 0 : a(3, 2) = 0 : a(3, 3) = 1.35355339 : a(3, 4) = 0.35355339
        'a(4, 1) = 0 : a(4, 2) = -1 : a(4, 3) = 0.35355339 : a(4, 4) = 1.35355339

        Call MASSMAT(a) 'n=6

        For i = 1 To n
            For j = 1 To n
                b(i, j) = a(i, j)
            Next j
        Next i

        Call MATINV(n, a, er)

        Console.WriteLine("error=" & er.ToString("0"))
        Console.WriteLine("逆行列")
        For i = 1 To n
            dat = i.ToString("0")
            For j = 1 To n
                dat = dat & "  " & a(i, j).ToString("e")
            Next j
            Console.WriteLine(dat)
        Next i

        For i = 1 To n
            For j = 1 To n
                c(i, j) = 0.0
                For k = 1 To n
                    c(i, j) = c(i, j) + a(i, k) * b(k, j)
                Next k
            Next j
        Next i
        Console.WriteLine("検算結果")
        For i = 1 To n
            dat = i.ToString("0")
            For j = 1 To n
                dat = dat & "  " & c(i, j).ToString("0")
            Next j
            Console.WriteLine(dat)
        Next i

    End Sub

    Private Sub MATINV(ByVal n As Integer, ByRef a(,) As Double, ByRef er As Integer)
        '*************************************
        ' 逆行列計算
        ' 山田嘉昭著：マトリックス法材料力学
        '*************************************
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim l As Integer
        Dim l1 As Integer
        Dim irow As Integer
        Dim icol As Integer
        Dim jrow As Integer
        Dim jcol As Integer
        Dim amax As Double
        Dim s As Double
        Dim t As Double
        Dim determ As Double = 1.0
        Dim factor As Double = 1.0E-30
        Dim ipivot(n) As Double
        Dim pivot(n) As Double
        Dim index(n, 2) As Integer

        For j = 1 To n
            ipivot(j) = 0
        Next j

        For i = 1 To n

            amax = 0.0
            For j = 1 To n
                If ipivot(j) - 1 <> 0 Then
                    For k = 1 To n
                        If 0 < ipivot(k) - 1 Then GoTo owari
                        If ipivot(k) - 1 < 0 Then
                            If Math.Abs(amax) - Math.Abs(a(j, k)) < 0 Then
                                irow = j
                                icol = k
                                amax = a(j, k)
                            End If
                        End If
                    Next k
                End If
            Next j
            ipivot(icol) = ipivot(icol) + 1

            If (irow + icol) Mod 2 <> 0 Then
                determ = -determ
            End If
            If irow - icol <> 0 Then
                For l = 1 To n
                    s = a(irow, l)
                    a(irow, l) = a(icol, l)
                    a(icol, l) = s
                Next l
            End If
            index(i, 1) = irow
            index(i, 2) = icol
            pivot(i) = a(icol, icol)
            If Math.Abs(pivot(i)) < factor Then GoTo okashii
            If determ <= 1.0E+30 Then
                determ = determ * pivot(i)
            End If
            a(icol, icol) = 1.0
            For l = 1 To n
                a(icol, l) = a(icol, l) / pivot(i)
            Next l

            For l1 = 1 To n
                If l1 - icol <> 0 Then
                    t = a(l1, icol)
                    a(l1, icol) = 0.0
                    For l = 1 To n
                        a(l1, l) = a(l1, l) - a(icol, l) * t
                    Next l
                End If
            Next l1

        Next i

        For i = 1 To n
            l = n + 1 - i
            If index(l, 1) - index(l, 2) <> 0 Then
                jrow = index(l, 1)
                jcol = index(l, 2)
                For k = 1 To n
                    s = a(k, jrow)
                    a(k, jrow) = a(k, jcol)
                    a(k, jcol) = s
                Next k
            End If
        Next i

owari:
        er = 0 '正常終了
        Return
okashii:
        er = 1 '異常終了

    End Sub

    Private Sub MASSMAT(ByRef am(,) As Double)
        '****************************
        '平面骨組要素質量行列
        '****************************
        Dim i As Integer
        Dim j As Integer
        Dim AL As Double

        AL = 1.0

        For i = 1 To 6
            For j = 1 To 6
                am(i, j) = 0.0
            Next j
        Next i

        am(1, 1) = 1.0 / 3.0
        am(2, 1) = 0.0
        am(3, 1) = 0.0
        am(4, 1) = 1.0 / 6.0
        am(5, 1) = 0.0
        am(6, 1) = 0.0
        am(2, 2) = 13.0 / 35.0
        am(3, 2) = 11.0 * AL / 210.0
        am(4, 2) = 0.0
        am(5, 2) = 9.0 / 70.0
        am(6, 2) = -13.0 * AL / 420.0
        am(3, 3) = AL * AL / 105.0
        am(4, 3) = 0.0
        am(5, 3) = 13.0 * AL / 420.0
        am(6, 3) = -AL * AL / 140.0
        am(4, 4) = 1.0 / 3.0
        am(5, 4) = 0.0
        am(6, 4) = 0.0
        am(5, 5) = 13.0 / 35.0
        am(6, 5) = -11.0 * AL / 210.0
        am(6, 6) = AL * AL / 105.0

        '対称項セット
        For i = 1 To 5
            For j = i + 1 To 6
                am(i, j) = am(j, i)
            Next j
        Next i

    End Sub

End Class