Option Explicit On
Option Strict On
Public Class Form1

    Private Sub Form1_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load
        Label1.Dock = DockStyle.Fill : Label1.TextAlign = ContentAlignment.MiddleCenter
        Label3.Dock = DockStyle.Fill : Label3.TextAlign = ContentAlignment.MiddleCenter
        Label5.Dock = DockStyle.Fill : Label5.TextAlign = ContentAlignment.MiddleCenter
        Label7.Dock = DockStyle.Fill : Label7.TextAlign = ContentAlignment.MiddleCenter
        Label9.Dock = DockStyle.Fill : Label9.TextAlign = ContentAlignment.MiddleCenter
        Label2.Dock = DockStyle.Fill : Label2.TextAlign = ContentAlignment.MiddleLeft
        Label4.Dock = DockStyle.Fill : Label4.TextAlign = ContentAlignment.MiddleLeft
        Label6.Dock = DockStyle.Fill : Label6.TextAlign = ContentAlignment.MiddleLeft
        Label8.Dock = DockStyle.Fill : Label8.TextAlign = ContentAlignment.MiddleLeft
        Label10.Dock = DockStyle.Fill : Label10.TextAlign = ContentAlignment.MiddleLeft

        Label1.Text = "入力ファイル"
        Label3.Text = "出力ファイル"
        Label5.Text = "開始時刻"
        Label7.Text = "完了時刻"
        Label9.Text = "計算時間"
        Label2.Text = ""
        Label4.Text = ""
        Label6.Text = ""
        Label8.Text = ""
        Label10.Text = ""
        Label11.Text = "nt,mm,ib"
    End Sub


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


    Private Sub MAIN()
        '２次元飽和定常浸透流解析
        Dim i As Integer : Dim j As Integer : Dim k As Integer
        Dim n As Integer
        Dim ia As Integer
        Dim ja As Integer
        Dim ne As Integer
        Dim sr As System.IO.StreamReader
        Dim dat As String
        Dim sbuf() As String
        Dim delim() As Char = {","c}
        Dim fnameR As String = ""
        Dim fnameW As String = ""

        Dim strcom As String          '書出用コメント
        Dim idan As Integer           '0：水平断面，1：鉛直断面
        Dim NODT As Integer           '節点総数
        Dim NELT As Integer           '要素総数
        Dim MATEL As Integer          '材料種類数
        Dim KOH As Integer            '全水頭指定節点数
        Dim KOQ As Integer            '流量指定節点数
        Dim nt As Integer             '全自由度（総節点数×１[１節点自由度]）
        Dim mm As Integer             '温度固定点処理後自由度
        Dim ib As Integer             'バンド幅
        Dim kakom(,) As Integer       '要素構成接点番号
        Dim Akx() As Double           '要素x方向透水係数
        Dim Aky() As Double           '要素y方向透水係数
        Dim matno() As Integer        '材料種別No
        Dim wm(,) As Double           '作業用材料物性記憶
        Dim x() As Double             '節点x座標
        Dim y() As Double             '節点y座標
        Dim nokh() As Integer         '全水頭指定節点番号
        Dim nokq() As Integer         '流量指定節点番号
        Dim Hinp() As Double          '指定全水頭入力
        Dim Qinp() As Double          '指定流量入力
        Dim tak(,) As Double          '全体透水行列
        Dim eak(3, 3) As Double       '要素透水行列
        Dim wak(,) As Double          '水頭指定節点対応列記憶
        Dim hvec() As Double          '全体全水頭ベクトル
        Dim qvec() As Double          '全体節点流量ベクトル
        Dim reac() As Double          '作業用ベクトル
        Dim vx() As Double            'x方向流速ベクトル
        Dim vy() As Double            'y方向流速ベクトル
        Dim index() As Integer        '温度固定点処理用
        Dim sttime As Date            '計算開始時刻
        Dim entime As Date            '計算終了時刻
        Dim tspan As TimeSpan         '計算時間
        Dim fact0 As Double = 1.0E-30 '0判定
        Dim elarea As Double

        '-----------------------------------------------------
        '入出力ファイル指定
        '-----------------------------------------------------
        OpenFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory()
        If OpenFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fnameR = OpenFileDialog1.FileName
        If SaveFileDialog1.ShowDialog() = Windows.Forms.DialogResult.OK Then fnameW = SaveFileDialog1.FileName
        Label2.Text = fnameR
        Label4.Text = fnameW
        '-----------------------------------------------------
        '時間計測開始
        '-----------------------------------------------------
        sttime = DateTime.Now
        Label6.Text = sttime.ToString
        '*****************************************************
        'データ入力
        '*****************************************************
        sr = New System.IO.StreamReader(fnameR, System.Text.Encoding.Default)
        '--------------------------------------------------------------------------------------
        '書出用コメント入力
        '--------------------------------------------------------------------------------------
        dat = sr.ReadLine : sbuf = dat.Split(delim) : strcom = sbuf(0)
        '------------------------------------------------------------------------------------------
        '解析断面(0：水平，1：鉛直)，節点数，要素数，材料種類数,全水頭固定節点数，流量指定節点数
        '------------------------------------------------------------------------------------------
        dat = sr.ReadLine : sbuf = dat.Split(delim)
        idan = CInt(sbuf(0))
        NODT = CInt(sbuf(1))
        NELT = CInt(sbuf(2))
        MATEL = CInt(sbuf(3))
        KOH = CInt(sbuf(4))
        KOQ = CInt(sbuf(5))
        nt = NODT * 1
        '------------------------------------------------------------------------------------------
        '配列寸法宣言
        '------------------------------------------------------------------------------------------
        ReDim kakom(NELT, 3)
        ReDim Akx(NELT)
        ReDim Aky(NELT)
        ReDim matno(NELT)
        ReDim wm(MATEL, 2)
        ReDim x(NODT)
        ReDim y(NODT)
        ReDim nokh(KOH)
        ReDim nokq(KOQ)
        ReDim Hinp(KOH)
        ReDim Qinp(KOQ)
        ReDim vx(NELT)
        ReDim vy(NELT)
        ReDim index(nt)
        ReDim tak(nt, nt)
        ReDim hvec(nt)
        ReDim qvec(nt)
        ReDim reac(nt)
        ReDim wak(nt, KOH)
        '------------------------------------------------------------------------------------------
        '材料特性（Akx,Aky）
        '------------------------------------------------------------------------------------------
        For i = 1 To MATEL
            dat = sr.ReadLine : sbuf = dat.Split(delim)
            For j = 1 To 2
                wm(i, j) = CDbl(sbuf(j - 1))
            Next j
        Next i
        '------------------------------------------------------------------------------------------
        '要素構成節点と材料種別（node-1,node-2,node-3,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))
            kakom(ne, 3) = CInt(sbuf(2))
            matno(ne) = CInt(sbuf(3))
            For i = 1 To MATEL
                If i = matno(ne) Then
                    Akx(ne) = wm(i, 1)
                    Aky(ne) = wm(i, 2)
                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 < KOH Then
            For i = 1 To KOH
                dat = sr.ReadLine : sbuf = dat.Split(delim)
                nokh(i) = CInt(sbuf(0)) : Hinp(i) = CDbl(sbuf(1))
            Next i
        End If
        '------------------------------------------------------------------------------------------
        '流量指定節点番号，指定流量
        '------------------------------------------------------------------------------------------
        For i = 1 To nt : qvec(i) = 0.0 : Next i
        If 0 < KOQ Then
            For i = 1 To KOQ
                dat = sr.ReadLine : sbuf = dat.Split(delim)
                nokq(i) = CInt(sbuf(0)) : Qinp(i) = CDbl(sbuf(1))
            Next i
        End If
        sr.Close()
        '-----------------------------------------------------
        '入力確認出力
        '-----------------------------------------------------
        Call INPDAT(fnameW, strcom, idan, NODT, NELT, MATEL, KOH, KOQ, x, y, _
                    nokh, nokq, kakom, Akx, Aky, matno, Hinp, Qinp)
        '要素面積確認と異常終了通知
        For ne = 1 To NELT
            i = kakom(ne, 1) : j = kakom(ne, 2) : k = kakom(ne, 3)
            elarea = 0.5 * ((x(k) - x(j)) * y(i) + (x(i) - x(k)) * y(j) + (x(j) - x(i)) * y(k))
            If elarea < fact0 Then Call IERROR1(fnameW)
        Next ne
        '*****************************************************
        'メイン処理計算
        '*****************************************************
        '-----------------------------------------------------
        '全体行列クリア
        '-----------------------------------------------------
        For i = 1 To nt
            reac(i) = 0.0
            hvec(i) = 0.0
            qvec(i) = 0.0
            For j = 1 To nt
                tak(i, j) = 0.0
        Next j, i
        '-----------------------------------------------------
        '全体行列組立
        '-----------------------------------------------------
        For ne = 1 To NELT
            Call STIFF_S3(ne, kakom, x, y, Akx, Aky, eak)
            Call ASMAT(ne, 3, 1, kakom, eak, tak)
        Next ne
        '-----------------------------------------------------
        '全水頭指定節点対応列記憶，前処理
        '-----------------------------------------------------
        If 0 < KOH Then
            For k = 1 To KOH
                For i = 1 To nt
                    wak(i, k) = tak(i, nokh(k))
            Next i, k
            For k = 1 To KOH
                For i = 1 To nt
                    hvec(i) = hvec(i) + Hinp(k) * wak(i, k)
            Next i, k
        End If
        If 0 < KOQ Then
            For k = 1 To KOQ
                qvec(nokq(k)) = Qinp(k)
            Next k
        End If
        For i = 1 To nt
            qvec(i) = qvec(i) - hvec(i)
        Next i
        '-----------------------------------------------------
        '全水頭指定節点処理
        '-----------------------------------------------------
        For i = 1 To nt : index(i) = i : Next i
        If 0 < KOH Then
            For i = 1 To KOH : n = nokh(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 mm
            ia = index(i)
            reac(i) = qvec(ia)
            For j = 1 To mm
                ja = index(j)
                tak(i, j) = tak(ia, ja)
        Next j, i
        '-----------------------------------------------------
        'バンド幅計算とバンド行列記憶
        '-----------------------------------------------------
        ib = IBAND(mm, tak, fact0)
        Call BAND(mm, ib, tak)
        '-----------------------------------------------------
        '対角項確認と異常終了通知
        '-----------------------------------------------------
        For i = 1 To mm
            If Math.Abs(tak(i, 1)) < fact0 Then Call IERROR2(fnameW)
        Next i
        '-----------------------------------------------------
        'コレスキー法
        '-----------------------------------------------------
        Call CHBAND10(mm, ib, tak)       '三角行列
        Call CHBAND11(mm, ib, tak, reac)
        '-----------------------------------------------------
        '全水頭算定
        '-----------------------------------------------------
        For i = 1 To nt : hvec(i) = 0.0 : Next i
        For i = 1 To mm : ia = index(i) : hvec(ia) = reac(i) : Next i
        For k = 1 To KOH
            hvec(nokh(k)) = Hinp(k)
        Next k
        '-----------------------------------------------------
        '節点流量算定
        '-----------------------------------------------------
        For i = 1 To nt : qvec(i) = 0.0 : Next i
        If 0 < KOQ Then
            For k = 1 To KOQ
                qvec(nokq(k)) = Qinp(k)
            Next k
        End If
        If 0 < KOH Then
            For k = 1 To KOH
                For i = 1 To nt
                    qvec(nokh(k)) = qvec(nokh(k)) + wak(i, k) * hvec(i)
            Next i, k
        End If
        '-----------------------------------------------------
        '要素平均流速算定
        '-----------------------------------------------------
        For ne = 1 To NELT
            Call CALV_S3(ne, kakom, x, y, hvec, Akx, Aky, vx, vy)
        Next ne
        '-----------------------------------------------------
        '時間計測完了
        '-----------------------------------------------------
        entime = DateTime.Now
        tspan = entime - sttime
        Label8.Text = entime.ToString
        Label10.Text = tspan.TotalSeconds.ToString("0.00sec")
        '*****************************************************
        '結果出力
        '*****************************************************
        Call OUTDAT(fnameW, idan, NODT, NELT, x, y, hvec, qvec, vx, vy)

        Label11.Text = "nt=" & nt.ToString() & "  mm=" & mm.ToString() & "  ib=" & ib.ToString()
        MessageBox.Show("計算完了", "通知")

    End Sub


    Private Sub INPDAT(ByRef fnameW As String, ByVal strcom As String, ByVal idan As Integer, _
                       ByVal NODT As Integer, ByVal NELT As Integer, ByVal MATEL As Integer, _
                       ByVal KOH As Integer, ByVal KOQ As Integer, _
                       ByVal x() As Double, ByRef y() As Double, _
                       ByRef nokh() As Integer, ByRef nokq() As Integer, _
                       ByRef kakom(,) As Integer, _
                       ByRef Akx() As Double, ByRef Aky() As Double, ByRef matno() As Integer, _
                       ByRef Hinp() As Double, ByRef Qinp() As Double)

        '入力確認出力
        Dim sw As System.IO.StreamWriter
        Dim dat As String
        Dim i As Integer
        Dim ne As Integer
        Dim w1() As String
        Dim w2() As String

        ReDim w1(NODT)
        ReDim w2(NODT)

        For i = 1 To NODT
            w1(i) = "--"
            w2(i) = "--"
        Next i
        If 0 < KOH Then
            For i = 1 To KOH
                w1(nokh(i)) = Hinp(i).ToString("E")
            Next i
        End If
        If 0 < KOQ Then
            For i = 1 To KOQ
                w2(nokq(i)) = Qinp(i).ToString("E")
            Next i
        End If

        sw = New System.IO.StreamWriter(fnameW, False, System.Text.Encoding.Default)
        dat = strcom : sw.WriteLine(dat)
        dat = "idan,NODT,NELT,MATEL,KOH,KOQ" : sw.WriteLine(dat)
        dat = idan.ToString("0")
        dat = dat & "," & NODT.ToString("0")
        dat = dat & "," & NELT.ToString("0")
        dat = dat & "," & MATEL.ToString("0")
        dat = dat & "," & KOH.ToString("0")
        dat = dat & "," & KOQ.ToString("0")
        sw.WriteLine(dat)
        dat = "*input-data" : sw.WriteLine(dat)
        dat = "node,x,y,Hinp,Qinp" : 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 & "," & w1(i)
            dat = dat & "," & w2(i)
            sw.WriteLine(dat)
        Next i
        dat = "element,node-1,node-2,node-3,Akx,Aky,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 & "," & kakom(ne, 3).ToString("0")
            dat = dat & "," & Akx(ne).ToString("E")
            dat = dat & "," & Aky(ne).ToString("E")
            dat = dat & "," & matno(ne).ToString
            sw.WriteLine(dat)
        Next ne
        sw.Close()

    End Sub


    Private Sub OUTDAT(ByRef fnameW As String, ByVal idan As Integer, ByVal NODT As Integer, ByVal NELT As Integer, _
                       ByVal x() As Double, ByRef y() As Double, _
                       ByRef hvec() As Double, ByRef qvec() As Double, _
                       ByRef vx() As Double, ByRef vy() As Double)

        '入力確認出力
        Dim sw As System.IO.StreamWriter
        Dim dat As String
        Dim i As Integer
        Dim ne As Integer

        sw = New System.IO.StreamWriter(fnameW, True, System.Text.Encoding.Default)
        dat = "*output-data" : sw.WriteLine(dat)
        dat = "node,x,y,hvec,qvec,P-head" : sw.WriteLine(dat)
        '節点番号,x座標,y座標,全水頭,節点流量,圧力水頭
        '接点流量は，流入が正・流出が負の値
        For i = 1 To NODT
            dat = i.ToString("0")
            dat = dat & "," & x(i).ToString("E")
            dat = dat & "," & y(i).ToString("E")
            dat = dat & "," & hvec(i).ToString("E")
            dat = dat & "," & qvec(i).ToString("E")
            Select Case idan
                Case 0 : dat = dat & "," & "--"
                Case 1 : dat = dat & "," & (hvec(i) - y(i)).ToString("E")
            End Select
            sw.WriteLine(dat)
        Next i
        dat = "element,vx,vy" : sw.WriteLine(dat)
        For ne = 1 To NELT
            dat = ne.ToString("0")
            dat = dat & "," & vx(ne).ToString("E")
            dat = dat & "," & vy(ne).ToString("E")
            sw.WriteLine(dat)
        Next ne
        dat = "Calculation time=" & Label10.Text : sw.WriteLine(dat)
        sw.Close()

    End Sub


    Private Sub STIFF_S3(ByVal ne As Integer, ByRef kakom(,) As Integer, _
                      ByRef x() As Double, ByRef y() As Double, _
                      ByRef Akx() As Double, ByRef Aky() As Double, ByRef eak(,) As Double)

        '要素透水行列
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim b1 As Double
        Dim b2 As Double
        Dim b3 As Double
        Dim c1 As Double
        Dim c2 As Double
        Dim c3 As Double
        Dim Delta As Double
        Dim wx(3, 3) As Double
        Dim wy(3, 3) As Double

        For i = 1 To 3
            For j = 1 To 3
                eak(i, j) = 0.0
                wx(i, j) = 0.0
                wy(i, j) = 0.0
        Next j, i
        i = kakom(ne, 1)
        j = kakom(ne, 2)
        k = kakom(ne, 3)
        Delta = 0.5 * ((x(k) - x(j)) * y(i) + (x(i) - x(k)) * y(j) + (x(j) - x(i)) * y(k))
        b1 = y(j) - y(k) : c1 = x(k) - x(j)
        b2 = y(k) - y(i) : c2 = x(i) - x(k)
        b3 = y(i) - y(j) : c3 = x(j) - x(i)
        wx(1, 1) = b1 * b1
        wx(1, 2) = b1 * b2
        wx(1, 3) = b1 * b3
        wx(2, 1) = b2 * b1
        wx(2, 2) = b2 * b2
        wx(2, 3) = b2 * b3
        wx(3, 1) = b3 * b1
        wx(3, 2) = b3 * b2
        wx(3, 3) = b3 * b3
        wy(1, 1) = c1 * c1
        wy(1, 2) = c1 * c2
        wy(1, 3) = c1 * c3
        wy(2, 1) = c2 * c1
        wy(2, 2) = c2 * c2
        wy(2, 3) = c2 * c3
        wy(3, 1) = c3 * c1
        wy(3, 2) = c3 * c2
        wy(3, 3) = c3 * c3
        For i = 1 To 3
            For j = 1 To 3
                eak(i, j) = Akx(ne) / 4.0 / Delta * wx(i, j) + Aky(ne) / 4.0 / Delta * wy(i, j)
        Next j, i


    End Sub


    Private Sub CALV_S3(ByVal ne As Integer, ByRef kakom(,) As Integer, ByRef x() As Double, ByRef y() As Double, _
                     ByRef hvec() As Double, ByRef Akx() As Double, ByRef Aky() As Double, _
                     ByRef vx() As Double, ByRef vy() As Double)
        '流速ベクトル算定
        Dim i As Integer
        Dim j As Integer
        Dim k As Integer
        Dim b1 As Double
        Dim b2 As Double
        Dim b3 As Double
        Dim c1 As Double
        Dim c2 As Double
        Dim c3 As Double
        Dim Delta As Double
        Dim h1 As Double
        Dim h2 As Double
        Dim h3 As Double

        i = kakom(ne, 1)
        j = kakom(ne, 2)
        k = kakom(ne, 3)
        Delta = 0.5 * ((x(k) - x(j)) * y(i) + (x(i) - x(k)) * y(j) + (x(j) - x(i)) * y(k))
        b1 = y(j) - y(k) : c1 = x(k) - x(j)
        b2 = y(k) - y(i) : c2 = x(i) - x(k)
        b3 = y(i) - y(j) : c3 = x(j) - x(i)
        h1 = hvec(i)
        h2 = hvec(j)
        h3 = hvec(k)
        vx(ne) = -Akx(ne) / 2.0 / Delta * (b1 * h1 + b2 * h2 + b3 * h3)
        vy(ne) = -Aky(ne) / 2.0 / Delta * (c1 * h1 + c2 * h2 + c3 * h3)

    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 CHBAND10(ByVal n1 As Integer, ByVal m1 As Integer, ByRef array(,) As Double)
        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
                        array(i, 1) = Math.Sqrt(z)
                    End If
                Else
                    array(i, j) = z / array(i, 1)
                End If
        Next j, i

    End Sub


    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

End Class