Option Strict On
Option Explicit On
Public Class Form1


    Private Sub ToolStripButton1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles ToolStripButton1.Click

        Dim sr As System.IO.StreamReader
        Dim dat As String
        Dim sbuf() As String
        Dim delim() As Char = {" "c}
        Dim fnameL As String = "" 'ファイルリスト
        Dim fnameR As String = "" '入力ファイル名
        Dim fnameW As String = "" '出力ファイル名
        Dim kprog As Integer   '第一中間点設定法
        Dim coefAB As Double   'スケーリングパラメータ算定用係数
        Dim deltaS0 As Double  '弧長
        Dim coefS As Double    '弧長に乗じる係数
        Dim maxload As Integer '載荷ステップ最大数

        '-----------------------------------------------------
        '入出力データファイルリスト指定
        '-----------------------------------------------------
        OpenFileDialog1.InitialDirectory = My.Computer.FileSystem.CurrentDirectory
        If OpenFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fnameL = OpenFileDialog1.FileName
        '-----------------------------------------------------
        '読み込み・実行
        '-----------------------------------------------------
        sr = New System.IO.StreamReader(fnameL, System.Text.Encoding.Default)
        Do Until sr.EndOfStream
            dat = sr.ReadLine()
            sbuf = dat.Split(delim)
            fnameR = System.IO.Path.GetDirectoryName(fnameL) & "\" & sbuf(0)
            fnameW = System.IO.Path.GetDirectoryName(fnameL) & "\" & sbuf(1)
            kprog = CInt(sbuf(2))
            coefAB = CDbl(sbuf(3))
            deltaS0 = CDbl(sbuf(4))
            coefS = CDbl(sbuf(5))
            maxload = CInt(sbuf(6))
            Call main_part(fnameR, fnameW, kprog, coefAB, deltaS0, coefS, maxload)
        Loop
        sr.Close()

        MessageBox.Show("計算完了", "通知")
    End Sub


    Private Sub main_part(ByVal fnameR As String, ByVal fnameW As String, ByVal kprog As Integer, ByVal coefAB As Double, _
                          ByVal deltaS0 As Double, ByVal coefS As Double, ByVal maxload As Integer)
        '平面骨組構造解析
        Dim i As Integer : Dim j As Integer
        Dim ne As Integer : Dim n As Integer
        Dim ia As Integer : Dim ja As Integer
        Dim sr As System.IO.StreamReader
        Dim sw As System.IO.StreamWriter
        Dim dat As String
        Dim sbuf() As String
        Dim delim() As Char = {","c}
        Dim strcom As String

        Dim NODT As Integer       '節点総数
        Dim NELT As Integer       '要素総数
        Dim MATEL As Integer      '材料種別数
        Dim KOX As Integer        'x方向変位既知節点数
        Dim KOY As Integer        'y方向変位既知節点数
        Dim KOZ As Integer        '回転変位既知節点数
        Dim NF As Integer         '載荷節点数
        Dim nt As Integer         '全自由度
        Dim mm As Integer         '変位既知節点処理語の自由度
        Dim kakom(,) As Integer   '要素構成接点番号
        Dim Em() As Double        '要素弾性係数
        Dim AA() As Double        '要素断面積
        Dim AI() As Double        '要素断面二次モーメント
        Dim matno() As Integer    '材料種別No
        Dim x() As Double         '節点x座標
        Dim y() As Double         '節点y座標
        Dim nokx() As Integer     'x方向変位既知節点番号
        Dim noky() As Integer     'y方向変位既知節点番号
        Dim nokz() As Integer     '回転変位既知節点番号
        Dim ff() As Double        '全体節点外力ベクトル
        Dim df() As Double        '全体節点外力基準増分ベクトル（入力値）
        Dim ftvec() As Double     '全体節点力ベクトル
        Dim index() As Integer    '変位既知節点処理のためのインデックス
        Dim sm(6, 6) As Double    '要素剛性マトリックス
        Dim tsm(,) As Double      '全体剛性マトリックス
        Dim dis() As Double       '全体節点変位増分
        Dim dist() As Double      '全体節点変位
        Dim dsav() As Double      '全体節点変位(釣合点記憶)
        Dim dis0() As Double      '全体節点変位増分(基準外力増分相当)
        Dim disr() As Double      '全体節点変位増分(不平衡力相当)
        Dim reac() As Double      '作業用ベクトル
        Dim work() As Double      '作業用ベクトル
        Dim hen(6, 6) As Double   '座標変換行列
        Dim fevec(6) As Double    '要素節点力ベクトル
        Dim wm(,) As Double       '材料種別作業用
        Dim sttime As Date        '計算開始時刻
        Dim entime As Date        '計算終了時刻
        Dim tspan As TimeSpan     '計算時間
        Dim ib As Integer         'バンド幅
        Dim fact0 As Double = 1.0E-30 '0判定
        Dim reforce(,) As Double  '断面力保存
        Dim refsave(,) As Double  '断面力保存
        Dim iload As Integer      '載荷ステップインクリメント
        Dim itaration As Integer  '収束過程インクリメント
        Dim maxita As Integer = 50 '繰り返し回数最大値
        Dim nocon As Integer       '収束・非収束記憶フラグ
        Dim lam As Double          '弧長増分パラメータ
        Dim dlam As Double         '弧長増分パラメータ
        Dim dtest0 As Double       '変位増分収束判定用変数
        Dim ftest0 As Double       '不平衡力収束判定用変数
        Dim dtest As Double        '変位増分収束判定用変数
        Dim ftest As Double        '不平衡力収束判定用変数
        Dim jdg As Integer         'コレスキー法使用可否
        Dim deltaS As Double       '弧長増分パラメータ
        Dim Spara As Double        'スケーリングパラメータ
        Dim fjdg0 As Double        '弧長変更用
        Dim fjdg As Double         '弧長変更用
        Dim fend0 As Double        '荷重変化監視用
        Dim fend1 As Double        '荷重変化監視用
        Dim dend0 As Double        '変位変化監視用
        Dim dend1 As Double        '変位変化監視用
        Dim jflagc As Integer      '計算終了状態
        Dim jdginv As Integer      '掃出法特異行列

        '-----------------------------------------------------
        '時間計測開始
        '-----------------------------------------------------
        sttime = DateTime.Now
        '*****************************************************
        'データ入力
        '*****************************************************
        sr = New System.IO.StreamReader(fnameR, System.Text.Encoding.Default)
        '--------------------------------------------------------------------------------------
        '書出用コメント入力
        '--------------------------------------------------------------------------------------
        dat = sr.ReadLine : sbuf = dat.Split(delim) : strcom = sbuf(0)
        '--------------------------------------------------------------------------------------
        '節点数，要素数，材料種別数,X方向変位既知節点数，Y方向変位既知節点数，回転変位既知節点数，載荷点数
        '--------------------------------------------------------------------------------------
        dat = sr.ReadLine : sbuf = dat.Split(delim)
        NODT = CInt(sbuf(0))
        NELT = CInt(sbuf(1))
        MATEL = CInt(sbuf(2))
        KOX = CInt(sbuf(3))
        KOY = CInt(sbuf(4))
        KOZ = CInt(sbuf(5))
        NF = CInt(sbuf(6))
        nt = NODT * 3
        '--------------------------------------------------------------------------------------
        '配列寸法宣言
        '--------------------------------------------------------------------------------------
        ReDim kakom(NELT, 2)
        ReDim Em(NELT)
        ReDim AA(NELT)
        ReDim AI(NELT)
        ReDim matno(NELT)
        ReDim wm(MATEL, 3)
        ReDim x(NODT)
        ReDim y(NODT)
        ReDim nokx(KOX)
        ReDim noky(KOY)
        ReDim nokz(KOZ)
        ReDim index(nt)
        ReDim tsm(nt, nt)
        ReDim ff(nt)
        ReDim df(nt)
        ReDim ftvec(nt)
        ReDim dis(nt)
        ReDim dist(nt)
        ReDim dsav(nt)
        ReDim dis0(nt)
        ReDim disr(nt)
        ReDim reac(nt)
        ReDim work(nt)
        ReDim reforce(NELT, 6)
        ReDim refsave(NELT, 6)
        '------------------------------------------------------------------------------------------
        '材料特性（Em,AA,AI）
        '------------------------------------------------------------------------------------------
        For i = 1 To MATEL
            dat = sr.ReadLine : sbuf = dat.Split(delim)
            For j = 1 To 3
                wm(i, j) = CDbl(sbuf(j - 1))
            Next j
        Next i
        '------------------------------------------------------------------------------------------
        '要素構成節点と材料種別（node-1,node-2,matno）
        '------------------------------------------------------------------------------------------
        For ne = 1 To NELT
            dat = sr.ReadLine : sbuf = dat.Split(delim)
            kakom(ne, 1) = CInt(sbuf(0))
            kakom(ne, 2) = CInt(sbuf(1))
            matno(ne) = CInt(sbuf(2))
            For i = 1 To MATEL
                If i = matno(ne) Then
                    Em(ne) = wm(i, 1)
                    AA(ne) = wm(i, 2)
                    AI(ne) = wm(i, 3)
                End If
            Next i
        Next ne
        '--------------------------------------------------------------------------------------
        '節点座標（x,y）
        '--------------------------------------------------------------------------------------
        For i = 1 To NODT
            dat = sr.ReadLine : sbuf = dat.Split(delim)
            x(i) = CDbl(sbuf(0))
            y(i) = CDbl(sbuf(1))
        Next i
        '--------------------------------------------------------------------------------------
        '変位既知節点番号，既知節点変位量
        '--------------------------------------------------------------------------------------
        If 0 < KOX Then
            For i = 1 To KOX
                dat = sr.ReadLine : sbuf = dat.Split(delim)
                nokx(i) = CInt(sbuf(0))
            Next i
        End If
        If 0 < KOY Then
            For i = 1 To KOY
                dat = sr.ReadLine : sbuf = dat.Split(delim)
                noky(i) = CInt(sbuf(0))
            Next i
        End If
        If 0 < KOZ Then
            For i = 1 To KOZ
                dat = sr.ReadLine : sbuf = dat.Split(delim)
                nokz(i) = CInt(sbuf(0))
            Next i
        End If
        '--------------------------------------------------------------------------------------
        '節点番号,x方向荷重,y方向荷重,モーメント荷重
        '--------------------------------------------------------------------------------------
        For i = 1 To nt : df(i) = 0.0 : Next i
        If 0 < NF Then
            For i = 1 To NF
                dat = sr.ReadLine : sbuf = dat.Split(delim)
                n = CInt(sbuf(0))
                df(3 * n - 2) = CDbl(sbuf(1))
                df(3 * n - 1) = CDbl(sbuf(2))
                df(3 * n) = CDbl(sbuf(3))
            Next i
        End If
        sr.Close()
        '*****************************************************
        '入力確認出力
        '*****************************************************
        Call INPDAT(fnameW, strcom, NODT, NELT, MATEL, KOX, KOY, KOZ, NF, deltaS, _
                    x, y, df, nokx, noky, nokz, kakom, Em, AA, AI, matno)
        '*****************************************************
        'メイン処理計算
        '*****************************************************
        '-----------------------------------------------------
        '変位既知節点自由度記憶
        '-----------------------------------------------------
        For i = 1 To nt : index(i) = i : Next i
        If 0 < KOX Then
            For i = 1 To KOX : n = 3 * nokx(i) - 2 : index(n) = 0 : Next i
        End If
        If 0 < KOY Then
            For i = 1 To KOY : n = 3 * noky(i) - 1 : index(n) = 0 : Next i
        End If
        If 0 < KOZ Then
            For i = 1 To KOZ : n = 3 * nokz(i) : index(n) = 0 : Next i
        End If
        mm = 0
        For i = 1 To nt
            If index(i) > 0 Then
                mm = mm + 1 : index(mm) = i
            End If
        Next i
        '-----------------------------------------------------
        '全体剛性行列・全体荷重ベクトルクリア
        '-----------------------------------------------------
        For i = 1 To nt
            ff(i) = 0.0
            dist(i) = 0.0
            dsav(i) = 0.0
            dis0(i) = 0.0
            disr(i) = 0.0
            reac(i) = 0.0
            ftvec(i) = 0.0
        Next i
        For ne = 1 To NELT
            For i = 1 To 6
                reforce(ne, i) = 0.0
                refsave(ne, i) = 0.0
            Next i
        Next ne
        '-----------------------------------------------------
        '方法1スケーリングパラメータ設定
        '-----------------------------------------------------
        ftest = 0.0
        For i = 1 To nt
            ftest = ftest + df(i) * df(i)
        Next i
        Spara = Math.Sqrt(coefAB / ftest)
        ftest = 0.0
        '*****************************************************
        '荷重増分過程
        '*****************************************************
        jflagc = 0 '正常終了フラグ初期化
        iload = 0
        deltaS = deltaS0
        ib = 0
        fend0 = 0.0
        jdginv = 0
        Do
            iload = iload + 1
            jdg = 0
again:
            itaration = 0
            nocon = 0
            lam = 0.0
            dlam = 0.0
            dtest0 = 0.0
            ftest0 = 0.0
            For ne = 1 To NELT : For i = 1 To 6 : reforce(ne, i) = refsave(ne, i) : Next i, ne
            For i = 1 To nt : dist(i) = dsav(i) : Next i
            '初期処理 {df}=(KT){dis0}
            For i = 1 To nt : For j = 1 To nt : tsm(i, j) = 0.0 : Next j, i
            For ne = 1 To NELT
                FRAMESM(ne, kakom, x, y, dist, Em, AA, AI, reforce, sm)
                ASMAT(ne, 2, 3, kakom, sm, tsm)
            Next ne
            For i = 1 To nt : reac(i) = 0.0 : dis0(i) = 0.0 : work(i) = 0.0 : Next i
            For i = 1 To mm
                ia = index(i) : work(i) = df(ia) : reac(i) = work(i)
                For j = 1 To mm : ja = index(j) : tsm(i, j) = tsm(ia, ja) : Next j
            Next i
            Select Case jdg
                Case 0 'Cholesky
                    ib = IBAND(mm, tsm, fact0)
                    BAND(mm, ib, tsm)
                    jdg = CHBAND10(mm, ib, tsm) '三角行列
                    CHBAND11(mm, ib, tsm, reac) '前進消去・後退代入
                Case 1 'Gauss-Jordan
                    jdginv = MATINV(mm, tsm)
                    If jdginv = 1 Then GoTo epart
                    For i = 1 To mm
                        reac(i) = 0.0
                        For j = 1 To mm
                            reac(i) = reac(i) + tsm(i, j) * work(j)
                        Next j
                    Next i
            End Select
            For i = 1 To mm : ia = index(i) : dis0(ia) = reac(i) : Next i
            Select Case kprog
                Case 1 : dlam = CAL_LAM01(nt, deltaS, Spara, dis0, df)
                Case 2 : dlam = CAL_LAM02(nt, deltaS, coefAB, Spara, dis0, df)
            End Select
            For i = 1 To nt : dis(i) = dlam * dis0(i) : Next i
            '====================================
            '収束過程
            '====================================
            Do
                itaration = itaration + 1
                lam = lam + dlam
                For i = 1 To nt : dist(i) = dist(i) + dis(i) : ftvec(i) = 0.0 : Next i
                '内力計算
                For ne = 1 To NELT
                    CALST(ne, kakom, x, y, dis, dist, Em, AA, AI, fevec, reforce)
                    ASVEC(ne, 2, 3, kakom, fevec, ftvec)
                Next ne
                '不平衡力計算
                For i = 1 To nt : reac(i) = ff(i) + lam * df(i) - ftvec(i) : ftvec(i) = reac(i) : Next i
                For i = 1 To mm : ia = index(i) : reac(i) = reac(ia) : dis(i) = dis(ia) : Next i
                '収束判定
                dtest = 0.0 : ftest = 0.0
                For i = 1 To mm
                    dtest = dtest + Math.Abs(dis(i))
                    ftest = ftest + Math.Abs(reac(i))
                Next i
                dtest0 = dtest0 + dtest
                ftest0 = ftest0 + ftest
                '=====================================================
                If dtest < 0.0000000001 Then Exit Do
                If dtest < dtest0 * 0.0001 And ftest < ftest0 * 0.0001 Then Exit Do
                '=====================================================
                For i = 1 To nt
                    dis0(i) = 0.0 : disr(i) = 0.0 : reac(i) = 0.0
                    For j = 1 To nt : tsm(i, j) = 0.0 : Next j
                Next i
                For ne = 1 To NELT
                    FRAMESM(ne, kakom, x, y, dist, Em, AA, AI, reforce, sm)
                    ASMAT(ne, 2, 3, kakom, sm, tsm)
                Next ne
                For i = 1 To mm

                    ia = index(i) : reac(i) = df(ia) : work(i) = reac(i)
                    For j = 1 To mm : ja = index(j) : tsm(i, j) = tsm(ia, ja) : Next j
                Next i
                Select Case jdg
                    Case 0 'Cholesky
                        ib = IBAND(mm, tsm, fact0)
                        BAND(mm, ib, tsm)
                        jdg = CHBAND10(mm, ib, tsm)
                        If jdg = 1 Then GoTo again
                        CHBAND11(mm, ib, tsm, reac) '前進消去・後退代入
                        For i = 1 To mm : ia = index(i) : dis0(ia) = reac(i) : reac(i) = ftvec(ia) : Next i
                        CHBAND11(mm, ib, tsm, reac) '前進消去・後退代入
                    Case 1 'Gauss-Jordan
                        jdginv = MATINV(mm, tsm)
                        If jdginv = 1 Then GoTo epart
                        For i = 1 To mm
                            reac(i) = 0.0
                            For j = 1 To mm
                                reac(i) = reac(i) + tsm(i, j) * work(j)
                            Next j
                        Next i
                        For i = 1 To mm : ia = index(i) : dis0(ia) = reac(i) : work(i) = ftvec(ia) : Next i
                        For i = 1 To mm
                            reac(i) = 0.0
                            For j = 1 To mm
                                reac(i) = reac(i) + tsm(i, j) * work(j)
                            Next j
                        Next i
                End Select
                For i = 1 To mm : ia = index(i) : disr(ia) = reac(i) : Next i
                dlam = CAL_LAMNEW(nt, Spara, dis0, disr, df)
                For i = 1 To nt : dis(i) = dlam * dis0(i) + disr(i) : Next i
            Loop While itaration <= maxita
            '====================================
            If maxita < itaration Then nocon = 1
            '釣合点確定
            For i = 1 To nt : ff(i) = ff(i) + lam * df(i) : dsav(i) = dist(i) : Next i  '荷重変位確定
            For ne = 1 To NELT : For i = 1 To 6 : refsave(ne, i) = reforce(ne, i) : Next i, ne  '要素内力確定
            Call DATAOUT(fnameW, NODT, NELT, iload, itaration, nocon, deltaS, lam, dtest, ftest, x, y, dist, ff, ftvec, reforce)
            '荷重・変位が動かなくなったら計算終了
            fend1 = 0.0 : dend1 = 0.0
            For i = 1 To nt
                fend1 = fend1 + Math.Abs(ff(i))
                dend1 = dend1 + Math.Abs(dist(i))
            Next i
            If 2 <= iload Then
                If ((1.0 - 0.000001) * fend0 < fend1 And fend1 < (1.0 + 0.000001) * fend0) And ((1.0 - 0.000001) * fend0 < fend1 And fend1 < (1.0 + 0.000001) * fend0) Then
                    jflagc = 1 '荷重・変位変化停止
                    GoTo epart
                End If
            End If
            fend0 = fend1 : dend0 = dend1
            '弧長変更判定
            '初回の増分荷重に対する増分荷重の比率が1/100未満となったら弧長を coefS 倍にして次回増分荷重を算定
            If iload = 1 Then
                fjdg0 = 0.0
                For i = 1 To nt
                    If fjdg0 < Math.Abs(lam * df(i)) Then fjdg0 = Math.Abs(lam * df(i))
                Next i
            End If
            fjdg = 0.0
            For i = 1 To nt
                If fjdg < Math.Abs(lam * df(i)) Then fjdg = Math.Abs(lam * df(i))
            Next i
            If fjdg < fjdg0 * 0.01 Then deltaS = coefS * deltaS
            If deltaS0 * 20.0 < deltaS Then deltaS = deltaS0 * 20.0
        Loop While nocon = 0 And iload < maxload
        '*****************************************************
epart:
        If nocon = 1 Then jflagc = 2 '未収束
        If jdginv = 1 Then jflagc = 3 '剛性行列特異
        '-----------------------------------------------------
        '時間計測完了
        '-----------------------------------------------------
        entime = DateTime.Now
        tspan = entime - sttime

        sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default)
        Select Case jflagc
            Case 0 : dat = "jflagc=" & jflagc & " : Normal completion"
            Case 1 : dat = "jflagc=" & jflagc & " : Stand still"
            Case 2 : dat = "jflagc=" & jflagc & " : No convergence"
            Case 3 : dat = "jflagc=" & jflagc & " : Singular matrix"
        End Select
        sw.WriteLine(dat)
        dat = "NODT=" & NODT.ToString() & "  nt=" & nt.ToString() & "  mm=" & mm.ToString() & "  ib=" & ib.ToString()
        sw.WriteLine(dat)
        dat = "deltaS=" & deltaS.ToString("g") & "  jdg=" & jdg.ToString() & "  iload=" & iload.ToString() & "  itaration=" & itaration.ToString()
        sw.WriteLine(dat)
        dat = "Calculation time=" & tspan.TotalSeconds.ToString("0.00sec")
        sw.WriteLine(dat)
        dat = entime.ToString()
        sw.WriteLine(dat)
        sw.Close()

    End Sub


    Private Sub INPDAT(ByRef fnameW As String, ByVal strcom As String, _
                   ByVal NODT As Integer, ByVal NELT As Integer, ByVal MATEL As Integer, _
                   ByVal KOX As Integer, ByVal KOY As Integer, ByVal KOZ As Integer, _
                   ByVal NF As Integer, ByVal deltaS As Double, _
                   ByVal x() As Double, ByRef y() As Double, ByRef f() As Double, _
                   ByRef nokx() As Integer, ByRef noky() As Integer, ByRef nokz() As Integer, _
                   ByRef kakom(,) As Integer, ByRef Em() As Double, ByRef AA() As Double, ByRef AI() As Double, _
                   ByRef matno() As Integer)

        '入力確認出力
        Dim sw As System.IO.StreamWriter
        Dim dat As String
        Dim i As Integer
        Dim k As Integer
        Dim l As Integer
        Dim m As Integer
        Dim ne As Integer
        sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default)
        dat = strcom : sw.WriteLine(dat)
        dat = "NODT,NELT,MATEL,KOX,KOY,KOZ,NF" : sw.WriteLine(dat)
        dat = NODT.ToString("0")
        dat = dat & "," & NELT.ToString("0")
        dat = dat & "," & MATEL.ToString("0")
        dat = dat & "," & KOX.ToString("0")
        dat = dat & "," & KOY.ToString("0")
        dat = dat & "," & KOZ.ToString("0")
        dat = dat & "," & NF.ToString("0")
        sw.WriteLine(dat)
        dat = "*input-data" : sw.WriteLine(dat)
        dat = "node,x,y,fx,fy,fz,x-fix,y-fix,z-fix" : sw.WriteLine(dat)
        For i = 1 To NODT
            k = 0 : l = 0 : m = 0
            If 0 < KOX Then
                For j = 1 To KOX
                    If i = nokx(j) Then k = 1
                Next j
            End If
            If 0 < KOY Then
                For j = 1 To KOY
                    If i = noky(j) Then l = 1
                Next j
            End If
            If 0 < KOZ Then
                For j = 1 To KOZ
                    If i = nokz(j) Then m = 1
                Next j
            End If
            dat = i.ToString("0")
            dat = dat & "," & x(i).ToString("E")
            dat = dat & "," & y(i).ToString("E")
            dat = dat & "," & f(3 * i - 2).ToString("E")
            dat = dat & "," & f(3 * i - 1).ToString("E")
            dat = dat & "," & f(3 * i).ToString("E")
            dat = dat & "," & k.ToString("0")
            dat = dat & "," & l.ToString("0")
            dat = dat & "," & m.ToString("0")
            sw.WriteLine(dat)
        Next i
        dat = "element,node-1,node-2,E,A,I,matno" : sw.WriteLine(dat)
        For ne = 1 To NELT
            dat = ne.ToString("0")
            dat = dat & "," & kakom(ne, 1).ToString("0")
            dat = dat & "," & kakom(ne, 2).ToString("0")
            dat = dat & "," & Em(ne).ToString("E")
            dat = dat & "," & AA(ne).ToString("E")
            dat = dat & "," & AI(ne).ToString("E")
            dat = dat & "," & matno(ne).ToString("0")
            sw.WriteLine(dat)
        Next ne
        sw.Close()

    End Sub


    Private Sub DATAOUT(ByVal fnameW As String, ByVal NODT As Integer, ByVal NELT As Integer, _
                     ByVal iload As Integer, ByVal itaration As Integer, ByVal nocon As Integer, _
                     ByVal deltaS As Double, ByVal lam As Double, ByVal dtest As Double, ByVal ftest As Double, _
                     ByRef x() As Double, ByRef y() As Double, ByRef dist() As Double, _
                     ByRef ff() As Double, ByRef ftvec() As Double, ByRef reforce(,) As Double)

        '結果出力（Appendモード）
        Dim i As Integer
        Dim ne As Integer
        Dim sw As System.IO.StreamWriter
        Dim dat As String

        sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default)
        Select Case nocon
            Case 0 : dat = "***** Output *****" : sw.WriteLine(dat)
            Case 1 : dat = "##### Output #####" : sw.WriteLine(dat)
        End Select

        dat = "iload,itaration,nocon,deltaS,lam,dtest,ftest" : sw.WriteLine(dat)
        dat = iload.ToString("0")
        dat = dat & "," & itaration.ToString("0")
        dat = dat & "," & nocon.ToString("0")
        dat = dat & "," & deltaS.ToString("E")
        dat = dat & "," & lam.ToString("E")
        dat = dat & "," & dtest.ToString("E")
        dat = dat & "," & ftest.ToString("E")
        sw.WriteLine(dat)
        dat = "*Node value" : sw.WriteLine(dat)
        dat = "node,x-coord,y-coord,x-dir,y-dir,z-dir,ff-x,ff-y,ff-z,ftvec-x,ftvec-y,ftvec-z" : sw.WriteLine(dat)
        For i = 1 To NODT
            dat = i.ToString("0")
            dat = dat & "," & x(i).ToString("E")
            dat = dat & "," & y(i).ToString("E")
            dat = dat & "," & dist(3 * i - 2).ToString("E")
            dat = dat & "," & dist(3 * i - 1).ToString("E")
            dat = dat & "," & dist(3 * i).ToString("E")
            dat = dat & "," & ff(3 * i - 2).ToString("E")
            dat = dat & "," & ff(3 * i - 1).ToString("E")
            dat = dat & "," & ff(3 * i).ToString("E")
            dat = dat & "," & ftvec(3 * i - 2).ToString("E")
            dat = dat & "," & ftvec(3 * i - 1).ToString("E")
            dat = dat & "," & ftvec(3 * i).ToString("E")
            sw.WriteLine(dat)
        Next i
        dat = "*Stress-resultant" : sw.WriteLine(dat)
        dat = "element,Ni,Si,Mi,Nj,Sj,Mj" : sw.WriteLine(dat)
        For ne = 1 To NELT
            dat = ne.ToString("0")
            For i = 1 To 6
                dat = dat & "," & reforce(ne, i).ToString("E")
            Next i
            sw.WriteLine(dat)
        Next ne
        sw.Close()

    End Sub


    Private Function CAL_LAM01(ByVal nt As Integer, ByVal deltaS As Double, ByVal Spara As Double, ByRef dis0() As Double, ByRef df() As Double) As Double
        Dim i As Integer
        Dim sum1 As Double
        Dim sum2 As Double

        sum1 = 0.0 : sum2 = 0.0
        For i = 1 To nt
            sum1 = sum1 + dis0(i) * dis0(i)
            sum2 = sum2 + df(i) * df(i)
        Next i
        Return Math.Sqrt(deltaS * deltaS / (sum1 + Spara * Spara * sum2))
    End Function


    Private Function CAL_LAM02(ByVal nt As Integer, ByVal deltaS As Double, ByVal coefAB As Double, ByRef Spara As Double, ByRef dis0() As Double, ByRef df() As Double) As Double
        Dim i As Integer
        Dim sum1 As Double
        Dim sum2 As Double
        Dim sum3 As Double

        sum1 = 0.0 : sum2 = 0.0 : sum3 = 0.0
        For i = 1 To nt
            sum1 = sum1 + dis0(i) * dis0(i)
            sum2 = sum2 + df(i) * df(i)
            sum3 = sum3 + dis0(i) * df(i)
        Next i
        Spara = 1.0 / coefAB * sum3 / sum2
        Return coefAB * Math.Sqrt(deltaS * deltaS / (sum1 + Spara * Spara * sum2))
    End Function


    Private Function CAL_LAMNEW(ByVal nt As Integer, ByVal Spara As Double, ByRef dis0() As Double, ByRef disr() As Double, ByRef df() As Double) As Double
        Dim i As Integer
        Dim sum1 As Double
        Dim sum2 As Double
        Dim sum3 As Double

        sum1 = 0.0
        sum2 = 0.0
        sum3 = 0.0
        For i = 1 To nt
            sum1 = sum1 + dis0(i) * disr(i)
            sum2 = sum2 + dis0(i) * dis0(i)
            sum3 = sum3 + df(i) * df(i)
        Next i

        Return -1.0 * sum1 / (sum2 + Spara * Spara * sum3)
    End Function


    Private Sub FRAMESM(ByVal ne As Integer, ByRef kakom(,) As Integer,
                        ByRef x() As Double, ByRef y() As Double, ByRef dist() As Double, _
                        ByRef Em() As Double, ByRef AA() As Double, ByRef AI() As Double, _
                        ByRef reforce(,) As Double, ByRef sm(,) As Double)

        Dim i As Integer : Dim j As Integer : Dim k As Integer

        Dim hen(6, 6) As Double
        Dim sme(6, 6) As Double
        Dim smv(6, 6) As Double
        Dim w(6, 6) As Double
        Dim AL As Double : Dim P As Double
        Dim xx1 As Double : Dim yy1 As Double : Dim xx2 As Double : Dim yy2 As Double

        For i = 1 To 6
            For j = 1 To 6
                smv(i, j) = 0.0
                sm(i, j) = 0.0
        Next j, i
        '軸力計算
        i = kakom(ne, 1)
        j = kakom(ne, 2)
        xx1 = x(i) + dist(3 * i - 2)
        yy1 = y(i) + dist(3 * i - 1)
        xx2 = x(j) + dist(3 * j - 2)
        yy2 = y(j) + dist(3 * j - 1)
        AL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1))
        P = 0.5 * (reforce(ne, 4) - reforce(ne, 1))

        '幾何剛性行列
        smv(1, 1) = 0.0
        smv(2, 1) = 0.0
        smv(2, 2) = 6.0 / 5 / AL
        smv(3, 1) = 0.0
        smv(3, 2) = 1.0 / 10.0
        smv(3, 3) = 2.0 * AL / 15.0
        smv(4, 1) = 0.0
        smv(4, 2) = 0.0
        smv(4, 3) = 0.0
        smv(4, 4) = 0.0
        smv(5, 1) = 0.0
        smv(5, 2) = -6.0 / 5.0 / AL
        smv(5, 3) = -1.0 / 10.0
        smv(5, 4) = 0.0
        smv(5, 5) = 6.0 / 5.0 / AL
        smv(6, 1) = 0.0
        smv(6, 2) = 1.0 / 10.0
        smv(6, 3) = -AL / 30.0
        smv(6, 4) = 0.0
        smv(6, 5) = -1.0 / 10.0
        smv(6, 6) = 2.0 * AL / 15.0
        For i = 1 To 6
            For j = i To 6
                smv(i, j) = smv(j, i)
        Next j, i

        '線形弾性剛性行列
        Call ELMSM(ne, xx1, yy1, xx2, yy2, Em, AA, AI, sme)
        For i = 1 To 6
            For j = 1 To 6
                sm(i, j) = sme(i, j) + P * smv(i, j)
        Next j, i

        Call HENMX(ne, xx1, yy1, xx2, yy2, hen)
        '(w)=(T)'*(k)
        For i = 1 To 6
            For j = 1 To 6
                w(i, j) = 0.0
                For k = 1 To 6
                    w(i, j) = w(i, j) + hen(k, i) * sm(k, j)
        Next k, j, i
        '(w)*(T)=(T)'(k,T)：(sm)は全体座標系での要素剛性行列
        For i = 1 To 6
            For j = 1 To 6
                sm(i, j) = 0.0
                For k = 1 To 6
                    sm(i, j) = sm(i, j) + w(i, k) * hen(k, j)
        Next k, j, i

    End Sub


    Private Sub ELMSM(ByVal ne As Integer, ByVal xx1 As Double, ByVal yy1 As Double, ByVal xx2 As Double, ByVal yy2 As Double, _
                      ByRef Em() As Double, ByRef AA() As Double, ByRef AI() As Double, ByRef sme(,) As Double)

        '要素剛性行列
        Dim i As Integer : Dim j As Integer
        Dim AL As Double : Dim EA As Double : Dim EI As Double

        AL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1))
        EA = Em(ne) * AA(ne)
        EI = Em(ne) * AI(ne)
        For i = 1 To 6
            For j = 1 To 6
                sme(i, j) = 0.0
        Next j, i
        sme(1, 1) = EA / AL
        sme(2, 2) = 12.0 * EI / (AL * AL * AL)
        sme(3, 2) = 6.0 * EI / (AL * AL)
        sme(3, 3) = 4.0 * EI / AL
        sme(4, 1) = -EA / AL
        sme(4, 4) = sme(1, 1)
        sme(5, 2) = -sme(2, 2)
        sme(5, 3) = -sme(3, 2)
        sme(5, 5) = sme(2, 2)
        sme(6, 2) = sme(3, 2)
        sme(6, 3) = 2.0 * EI / AL
        sme(6, 5) = -sme(3, 2)
        sme(6, 6) = sme(3, 3)
        '対称項セット
        For i = 1 To 6
            For j = 1 To i
                sme(j, i) = sme(i, j)
        Next j, i

    End Sub


    Private Sub CALST(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef x() As Double, ByRef y() As Double, _
                      ByRef dis() As Double, ByRef dist() As Double, _
                      ByRef Em() As Double, ByRef AA() As Double, ByRef AI() As Double, _
                      ByRef fevec() As Double, ByRef reforce(,) As Double)

        '断面力算定
        Dim i As Integer : Dim k As Integer
        Dim ii As Integer : Dim jj As Integer : Dim ia As Integer
        Dim hen(6, 6) As Double
        Dim sme(6, 6) As Double
        Dim AL As Double : Dim BL As Double
        Dim deltai As Double : Dim deltaj As Double
        Dim wd(6) As Double : Dim wde(6) As Double : Dim ws(6) As Double
        Dim xx1 As Double : Dim yy1 As Double : Dim xx2 As Double : Dim yy2 As Double

        Call ROTATE(ne, kakom, x, y, dis, dist, deltai, deltaj)

        '要素節点変位増分（全体座標系）
        For i = 1 To 2
            ia = kakom(ne, i)
            wd(3 * i - 2) = dis(3 * ia - 2)
            wd(3 * i - 1) = dis(3 * ia - 1)
            wd(3 * i) = dis(3 * ia)
        Next i
        ii = kakom(ne, 1)
        jj = kakom(ne, 2)
        '前回要素座標系からの変位増分
        xx1 = x(ii) + dist(3 * ii - 2) - dis(3 * ii - 2)
        yy1 = y(ii) + dist(3 * ii - 1) - dis(3 * ii - 1)
        xx2 = x(jj) + dist(3 * jj - 2) - dis(3 * jj - 2)
        yy2 = y(jj) + dist(3 * jj - 1) - dis(3 * jj - 1)
        HENMX(ne, xx1, yy1, xx2, yy2, hen)
        For i = 1 To 6
            wde(i) = 0.0
            For k = 1 To 6
                wde(i) = wde(i) + hen(i, k) * wd(k)
        Next k, i
        '前回部材長
        xx1 = x(ii) + dist(3 * ii - 2) - dis(3 * ii - 2)
        yy1 = y(ii) + dist(3 * ii - 1) - dis(3 * ii - 1)
        xx2 = x(jj) + dist(3 * jj - 2) - dis(3 * jj - 2)
        yy2 = y(jj) + dist(3 * jj - 1) - dis(3 * jj - 1)
        BL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1))
        '今回部材長
        xx1 = x(ii) + dist(3 * ii - 2)
        yy1 = y(ii) + dist(3 * ii - 1)
        xx2 = x(jj) + dist(3 * jj - 2)
        yy2 = y(jj) + dist(3 * jj - 1)
        AL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1))
        wd(1) = 0.0
        wd(2) = 0.0
        wd(3) = deltai
        wd(4) = AL - BL
        wd(5) = 0.0
        wd(6) = deltaj
        '要素内力増分(前回微小変形要素剛性使用)
        xx1 = x(ii) + dist(3 * ii - 2) - dis(3 * ii - 2)
        yy1 = y(ii) + dist(3 * ii - 1) - dis(3 * ii - 1)
        xx2 = x(jj) + dist(3 * jj - 2) - dis(3 * jj - 2)
        yy2 = y(jj) + dist(3 * jj - 1) - dis(3 * jj - 1)
        Call ELMSM(ne, xx1, yy1, xx2, yy2, Em, AA, AI, sme)
        For i = 1 To 6
            ws(i) = 0.0
            For k = 1 To 6
                ws(i) = ws(i) + sme(i, k) * wd(k)
        Next k, i
        '断面力保存
        For i = 1 To 6
            reforce(ne, i) = reforce(ne, i) + ws(i)
        Next i
        '全体座標系での変形後要素内力
        xx1 = x(ii) + dist(3 * ii - 2)
        yy1 = y(ii) + dist(3 * ii - 1)
        xx2 = x(jj) + dist(3 * jj - 2)
        yy2 = y(jj) + dist(3 * jj - 1)
        Call HENMX(ne, xx1, yy1, xx2, yy2, hen)
        For i = 1 To 6
            fevec(i) = 0.0
            For k = 1 To 6
                fevec(i) = fevec(i) + hen(k, i) * reforce(ne, k)
        Next k, i

    End Sub


    Private Sub ROTATE(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef x() As Double, ByRef y() As Double, _
        ByRef dis() As Double, ByRef dist() As Double, ByRef deltai As Double, ByRef deltaj As Double)

        '剛体回転除去
        Dim i As Integer : Dim k As Integer
        Dim ii As Integer : Dim jj As Integer : Dim ia As Integer

        Dim hen(6, 6) As Double
        Dim wd(6) As Double : Dim wde(6) As Double
        Dim xx1 As Double : Dim yy1 As Double
        Dim xx2 As Double : Dim yy2 As Double

        Dim BL As Double : Dim Ti As Double : Dim Tj As Double
        Dim theta1i As Double : Dim theta1j As Double
        Dim theta2i As Double : Dim theta2j As Double

        ii = kakom(ne, 1)
        jj = kakom(ne, 2)
        '前回における初期部材座標系からの変位量
        For i = 1 To 2
            ia = kakom(ne, i)
            wd(3 * i - 2) = dist(3 * ia - 2) - dis(3 * ia - 2)
            wd(3 * i - 1) = dist(3 * ia - 1) - dis(3 * ia - 1)
            wd(3 * i) = dist(3 * ia) - dis(3 * ia)
        Next i
        '初期部材長さ
        xx1 = x(ii)
        yy1 = y(ii)
        xx2 = x(jj)
        yy2 = y(jj)
        BL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1))
        Call HENMX(ne, xx1, yy1, xx2, yy2, hen)
        For i = 1 To 6
            wde(i) = 0.0
            For k = 1 To 6
                wde(i) = wde(i) + hen(i, k) * wd(k)
        Next k, i
        '前回の内力計算用回転角
        Ti = Math.Tan(wde(3))
        Tj = Math.Tan(wde(6))
        theta1i = ((BL + wde(4) - wde(1)) * Ti - (wde(5) - wde(2))) / ((BL + wde(4) - wde(1)) + (wde(5) - wde(2)) * Ti)
        theta1j = ((BL + wde(4) - wde(1)) * Tj - (wde(5) - wde(2))) / ((BL + wde(4) - wde(1)) + (wde(5) - wde(2)) * Tj)

        '今回における初期部材座標系からの変位量
        For i = 1 To 2
            ia = kakom(ne, i)
            wd(3 * i - 2) = dist(3 * ia - 2)
            wd(3 * i - 1) = dist(3 * ia - 1)
            wd(3 * i) = dist(3 * ia)
        Next i
        For i = 1 To 6
            wde(i) = 0.0
            For k = 1 To 6
                wde(i) = wde(i) + hen(i, k) * wd(k)
        Next k, i
        '前回の内力計算用回転角
        Ti = Math.Tan(wde(3))
        Tj = Math.Tan(wde(6))
        theta2i = ((BL + wde(4) - wde(1)) * Ti - (wde(5) - wde(2))) / ((BL + wde(4) - wde(1)) + (wde(5) - wde(2)) * Ti)
        theta2j = ((BL + wde(4) - wde(1)) * Tj - (wde(5) - wde(2))) / ((BL + wde(4) - wde(1)) + (wde(5) - wde(2)) * Tj)

        '内力計算用回転角増分
        deltai = theta2i - theta1i
        deltaj = theta2j - theta1j

    End Sub


    Private Sub HENMX(ByVal ne As Integer, ByRef xx1 As Double, ByVal yy1 As Double, ByVal xx2 As Double, ByVal yy2 As Double, ByRef hen(,) As Double)
        '座標変換行列
        Dim i As Integer
        Dim j As Integer
        Dim AL As Double
        Dim cs As Double
        Dim sn As Double

        AL = Math.Sqrt((xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1))
        cs = (xx2 - xx1) / AL
        sn = (yy2 - yy1) / AL
        For i = 1 To 6
            For j = 1 To 6
                hen(i, j) = 0.0
        Next j, i
        hen(1, 1) = cs
        hen(1, 2) = sn
        hen(2, 1) = -sn
        hen(2, 2) = cs
        hen(3, 3) = 1.0
        For i = 1 To 3
            For j = 1 To 3
                hen(i + 3, j + 3) = hen(i, j)
        Next j, i
    End Sub


    Private Sub ASMAT(ByVal ne As Integer, ByVal nod As Integer, ByVal nhen As Integer, _
                      ByRef kakom(,) As Integer, ByRef sm(,) As Double, ByRef tsm(,) As Double)
        '剛性行列組立(nod：1要素節点数，nhen：1節点の自由度)
        Dim i As Integer : Dim j As Integer : Dim k As Integer : Dim l As Integer
        Dim ki As Integer : Dim kj As Integer : Dim ks As Integer : Dim js As Integer
        Dim kik As Integer : Dim kjl As Integer : Dim ksk As Integer : Dim jsl As Integer

        For i = 1 To nod
            For j = 1 To nod
                ki = (kakom(ne, i) - 1) * nhen
                kj = (kakom(ne, j) - 1) * nhen
                ks = (i - 1) * nhen
                js = (j - 1) * nhen
                For k = 1 To nhen
                    For l = 1 To nhen
                        kik = ki + k : kjl = kj + l : ksk = ks + k : jsl = js + l
                        tsm(kik, kjl) = tsm(kik, kjl) + sm(ksk, jsl)
                Next l, k
        Next j, i

    End Sub


    Private Sub ASVEC(ByVal ne As Integer, ByVal nod As Integer, ByVal nhen As Integer, _
                   ByRef kakom(,) As Integer, ByRef fevec() As Double, ByRef ftvec() As Double)
        '荷重ベクトル組立(nod：1要素節点数，nhen：1節点の自由度)
        Dim i As Integer : Dim k As Integer
        Dim ki As Integer : Dim ks As Integer : Dim kik As Integer : Dim ksk As Integer

        For i = 1 To nod
            ki = (kakom(ne, i) - 1) * nhen
            ks = (i - 1) * nhen
            For k = 1 To nhen
                kik = ki + k : ksk = ks + k
                ftvec(kik) = ftvec(kik) + fevec(ksk)
            Next k
        Next i

    End Sub


    Private Function CHBAND10(ByVal n1 As Integer, ByVal m1 As Integer, ByRef array(,) As Double) As Integer
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim mm As Integer
        Dim mn As Integer
        Dim z As Double

        array(1, 1) = Math.Sqrt(array(1, 1))
        For j = 2 To m1
            array(1, j) = array(1, j) / array(1, 1)
        Next j
        For i = 2 To n1
            If n1 - i + 1 > m1 Then
                mm = m1
            Else
                mm = n1 - i + 1
            End If
            For j = 1 To mm
                z = array(i, j)
                If j <> m1 Then
                    If m1 - j + 1 < i Then
                        mn = m1 - j + 1
                    Else
                        mn = i
                    End If
                    For k = 2 To mn
                        z = z - array(i - k + 1, k) * array(i - k + 1, j + k - 1)
                    Next k
                    If j > 1 Then
                        array(i, j) = z / array(i, 1)
                    Else
                        If z < 0.0 Then Return 1
                        array(i, 1) = Math.Sqrt(z)
                    End If
                Else
                    array(i, j) = z / array(i, 1)
                End If
        Next j, i
        Return 0

    End Function


    Private Sub CHBAND11(ByVal n1 As Integer, ByVal m1 As Integer, ByRef array(,) As Double, ByRef vector() As Double)
        Dim i As Integer
        Dim j As Integer
        Dim mm As Integer
        Dim z As Double

        vector(1) = vector(1) / array(1, 1)
        For i = 2 To n1
            z = vector(i)
            If i > m1 Then
                mm = m1
            Else
                mm = i
            End If
            For j = 2 To mm
                z = z - array(i - j + 1, j) * vector(i - j + 1)
            Next j
            vector(i) = z / array(i, 1)
        Next i

        vector(n1) = vector(n1) / array(n1, 1)
        For i = n1 - 1 To 1 Step -1
            z = vector(i)
            For j = 2 To m1
                If i + j - 1 > n1 Then Exit For
                z = z - array(i, j) * vector(i + j - 1)
            Next j
            vector(i) = z / array(i, 1)
        Next i

    End Sub


    Private Function IBAND(ByVal mm As Integer, ByRef array(,) As Double, ByVal fact0 As Double) As Integer
        Dim i As Integer
        Dim j As Integer
        '-----------------------------------------------------
        'バンド幅計算
        '-----------------------------------------------------
        IBAND = 1
        For i = 1 To mm - 1
            For j = mm To i + 1 Step -1
                If fact0 <= Math.Abs(array(i, j)) Then Exit For
            Next j
            If IBAND <= j - i + 1 Then IBAND = j - i + 1
        Next i
        If mm < IBAND Then IBAND = mm

    End Function


    Private Sub BAND(ByVal mm As Integer, ByVal ib As Integer, ByRef array(,) As Double)
        Dim i As Integer
        Dim j As Integer
        Dim mn As Integer
        Dim work() As Double
        ReDim work(mm)
        '-----------------------------------------------------
        'フルマトリックスのバンド化記憶
        '-----------------------------------------------------
        For i = 2 To mm
            If i + ib - 1 < mm Then
                mn = ib
            Else
                mn = mm - i + 1
            End If
            For j = 1 To ib
                work(j) = 0.0
            Next j
            For k = 1 To mn
                work(k) = array(i, i + k - 1)
            Next k
            For j = 1 To ib
                array(i, j) = work(j)
            Next j
        Next i

    End Sub


    Private Sub IERROR1(ByVal fnameW As String)
        '結果出力（Appendモード）
        Dim sw As System.IO.StreamWriter
        Dim dat As String
        sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default)
        dat = "area of element is less or equal to 0" : sw.WriteLine(dat)
        sw.Close()
        MessageBox.Show("要素面積<=0", "異常終了通知")
        Me.Close()

    End Sub


    Private Sub IERROR2(ByVal fnameW As String)
        '結果出力（Appendモード）
        Dim sw As System.IO.StreamWriter
        Dim dat As String
        sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default)
        dat = "diagonal value is less or equal to 0" : sw.WriteLine(dat)
        sw.Close()
        MessageBox.Show("対角項<=0", "異常終了通知")
        Me.Close()

    End Sub

    Private Function MATINV(ByVal n As Integer, ByRef ary(,) As Double) 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 Return 0
                        If ipivot(k) - 1 < 0 Then
                            If Math.Abs(amax) - Math.Abs(ary(j, k)) < 0 Then
                                irow = j
                                icol = k
                                amax = ary(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 = ary(irow, l)
                    ary(irow, l) = ary(icol, l)
                    ary(icol, l) = s
                Next l
            End If
            index(i, 1) = irow
            index(i, 2) = icol
            pivot(i) = ary(icol, icol)
            If Math.Abs(pivot(i)) < factor Then Return 1
            If determ <= 1.0E+30 Then
                determ = determ * pivot(i)
            End If
            ary(icol, icol) = 1.0
            For l = 1 To n
                ary(icol, l) = ary(icol, l) / pivot(i)
            Next l

            For l1 = 1 To n
                If l1 - icol <> 0 Then
                    t = ary(l1, icol)
                    ary(l1, icol) = 0.0
                    For l = 1 To n
                        ary(l1, l) = ary(l1, l) - ary(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 = ary(k, jrow)
                    ary(k, jrow) = ary(k, jcol)
                    ary(k, jcol) = s
                Next k
            End If
        Next i

        Return 0

    End Function

End Class