Attribute VB_Name = "LINEAR2"

' Copyright (c) 1997 C. David Eagle Jr.

DefInt I-N
DefDbl A-H, O-Z

Static Sub LINEAR2(N, A(), B(), IER)

    ' Solution of a system of linear equations subroutine

    ' Solves [ A ] * { X } = { B } using the Gauss-Jordan method

    ' Input

    '  N   = number of equations
    '  A() = matrix of coefficients   (N rows by N columns)
    '  B() = right hand column vector (N rows)

    ' Output

    '  A() = inverse of incoming matrix [ A ] (N rows by N columns)
    '  B() = solution vector of linear system (N rows)
    '  IER = error flag
    '    0 = no error
    '    1 = singular matrix

    ReDim IPIVOT(N), INDEXR(N), INDEXC(N)

    IER = 0

    For j = 1 To N
        IPIVOT(j) = 0
    Next j

    For i = 1 To N
        PMAX = 0#

        For j = 1 To N
            If (IPIVOT(j) <> 1) Then
               For K = 1 To N
                   If (IPIVOT(K) = 0) Then
                      If (Abs(A(j, K)) >= PMAX) Then
                         PMAX = Abs(A(j, K))
                         IR = j
                         IC = K
                      End If
                   ElseIf (IPIVOT(K) > 1) Then
                      IER = 1
                      GoTo EXITSUB
                   End If
               Next K
            End If
        Next j

        IPIVOT(IC) = IPIVOT(IC) + 1
       
        If (IR <> IC) Then
           For L = 1 To N
               TMP = A(IR, L)
               A(IR, L) = A(IC, L)
               A(IC, L) = TMP
           Next L

           TMP = B(IR)
           B(IR) = B(IC)
           B(IC) = TMP
        End If

        INDEXR(i) = IR
        INDEXC(i) = IC
       
        If (A(IC, IC) = 0#) Then
           IER = 1
           GoTo EXITSUB
        End If
       
        PIVINV = 1# / A(IC, IC)
        A(IC, IC) = 1#

        For L = 1 To N
            A(IC, L) = A(IC, L) * PIVINV
        Next L

        B(IC) = B(IC) * PIVINV

        For LL = 1 To N
            If (LL <> IC) Then
               TMP = A(LL, IC)
               A(LL, IC) = 0#
               For L = 1 To N
                   A(LL, L) = A(LL, L) - A(IC, L) * TMP
               Next L
               B(LL) = B(LL) - B(IC) * TMP
            End If
        Next LL
    Next i

    For L = N To 1 Step -1
        If (INDEXR(L) <> INDEXC(L)) Then
           For K = 1 To N
               TMP = A(K, INDEXR(L))
               A(K, INDEXR(L)) = A(K, INDEXC(L))
               A(K, INDEXC(L)) = TMP
           Next K
        End If
    Next L

EXITSUB:
    Erase IPIVOT, INDEXC, INDEXR

End Sub