Sub main()
'-----------------------------------------------------------------------
' Use GAUSSJ of "Numerical Recipes" to solve a set of linear algebraic equations
' Ax=b
' Instructor: Nam Sun Wang
'-----------------------------------------------------------------------
' Either declare the subroutines or make sure the paremeters types
' are passed correctly.  Otherwise, results are highly unpredictable.

      Const NP = 100
      Dim A(NP, NP), B(NP, 1)      'B must be two-dimensional to match the notation in GAUSSJ.

' Dimension of A -------------------------------------------------------
      N = 60
      
' Input matrix/vector coefficients -------------------------------------
      For i = 1 To N
        For j = 1 To N
'         A(1, 1) = Cells(5, 2)
          A(i, j) = Cells(i + 4, j + 1)
        Next
        B(i, 1) = 0
      Next
      
' Solve the linear set of equations ------------------------------------
      Call GAUSSJ(A(), N, NP, B(), 1, 1)

' Output results (Upon returning from GAUSSJ, inverse of A is stored in A, and solution x in B)
      For i = 1 To N
        For j = 1 To N
'         Cells(127,2) = A(i, j)
          Cells(i + 126, j + 1) = A(i, j)
        Next
      Next

End Sub


'-----------------------------------------------------------------------
' GAUSSJ of "Numerical Recipes"
'-----------------------------------------------------------------------
Sub GAUSSJ(A(), N, NP, B(), M, MP)
ReDim IPIV(N), INDXR(N), INDXC(N)    'Used to be Dim
For j = 1 To N
  IPIV(j) = 0
Next j
For i = 1 To N
  BIG = 0!
  For j = 1 To N
    If IPIV(j) <> 1 Then
      For K = 1 To N
        If IPIV(K) = 0 Then
          If Abs(A(j, K)) >= BIG Then
            BIG = Abs(A(j, K))
            irow = j
            ICOL = K
          End If
        ElseIf IPIV(K) > 1 Then
'          Print "Singular matrix"      'no printing in Excel
           Call MsgBox("Singular matrix")
          Exit Sub
        End If
      Next K
    End If
  Next j
  IPIV(ICOL) = IPIV(ICOL) + 1
  If irow <> ICOL Then
    For L = 1 To N
      DUM = A(irow, L)
      A(irow, L) = A(ICOL, L)
      A(ICOL, L) = DUM
    Next L
    For L = 1 To M
      DUM = B(irow, L)
      B(irow, L) = B(ICOL, L)
      B(ICOL, L) = DUM
    Next L
  End If
  INDXR(i) = irow
  INDXC(i) = ICOL
'  If A(ICOL, ICOL) = 0! Then Print "Singular matrix.": Exit Sub
   If A(ICOL, ICOL) = 0! Then Call MsgBox("Singular matrix"): Exit Sub
  PIVINV = 1! / A(ICOL, ICOL)
  A(ICOL, ICOL) = 1!
  For L = 1 To N
    A(ICOL, L) = A(ICOL, L) * PIVINV
  Next L
  For L = 1 To M
    B(ICOL, L) = B(ICOL, L) * PIVINV
  Next L
  For LL = 1 To N
    If LL <> ICOL Then
      DUM = A(LL, ICOL)
      A(LL, ICOL) = 0!
      For L = 1 To N
        A(LL, L) = A(LL, L) - A(ICOL, L) * DUM
      Next L
      For L = 1 To M
        B(LL, L) = B(LL, L) - B(ICOL, L) * DUM
      Next L
    End If
  Next LL
Next i
For L = N To 1 Step -1
  If INDXR(L) <> INDXC(L) Then
    For K = 1 To N
      DUM = A(K, INDXR(L))
      A(K, INDXR(L)) = A(K, INDXC(L))
      A(K, INDXC(L)) = DUM
    Next K
  End If
Next L
Erase INDXC, INDXR, IPIV
End Sub