VB.Net Tutorial/Data Type/Matrix

Calculate the determinant of an array

" Quote from
"Visual Basic 2005 Cookbook Solutions for VB 2005 Programmers
"by Tim Patrick (Author), John Craig (Author)
"# Publisher: O"Reilly Media, Inc. (September 21, 2006)
"# Language: English
"# ISBN-10: 0596101775
"# ISBN-13: 978-0596101770
Public Class Tester
    Public Shared Sub Main
        Dim matrixA(,) As Double = { _
            {1, 2, 3}, _
            {5, 4, 6}, _
            {9, 7, 8}}
        Dim determinant As Double = MatrixHelper.Determinant(matrixA)
        Console.WriteLine(MatrixHelper.MakeDisplayable(matrixA) & _
            vbNewLine & vbNewLine & "Determinant: " & _
            determinant.ToString)
     End Sub
End Class

Module MatrixHelper
    Public Function MakeDisplayable(ByVal sourceMatrix(,) As Double) As String
        " ----- Prepare a multi-line string that shows the contents
        "       of a matrix, a 2D array.
        Dim rows As Integer
        Dim cols As Integer
        Dim eachRow As Integer
        Dim eachCol As Integer
        Dim result As New System.Text.StringBuilder
        " ----- Process all rows of the matrix, generating one
        "       output line per row.
        rows = UBound(sourceMatrix, 1) + 1
        cols = UBound(sourceMatrix, 2) + 1
        For eachRow = 0 To rows - 1
            " ----- Process each column of the matrix on a single
            "       row, separating values by commas.
            If (eachRow > 0) Then result.AppendLine()
            For eachCol = 0 To cols - 1
                " ----- Add a single matrix element to the output.
                If (eachCol > 0) Then result.Append(",")
                result.Append(sourceMatrix(eachRow, eachCol).ToString)
            Next eachCol
        Next eachRow
        " ----- Finished.
        Return result.ToString
    End Function
    Public Function MakeDisplayable(ByVal sourceArray() As Double) As String
        " ----- Present an array as multiple lines of output.
        Dim result As New System.Text.StringBuilder
        Dim scanValue As Double
        For Each scanValue In sourceArray
            result.AppendLine(scanValue.ToString)
        Next scanValue
        Return result.ToString
    End Function
    Public Function Inverse(ByVal sourceMatrix(,) As Double) As Double(,)
        " ----- Build a new matrix that is the mathematical inverse
        "       of the supplied matrix. Multiplying a matrix and its
        "       inverse together will give the identity matrix.
        Dim eachCol As Integer
        Dim eachRow As Integer
        Dim rowsAndCols As Integer
        " ----- Determine the size of each dimension of the matrix.
        "       Only square matrices can be inverted.
        If (UBound(sourceMatrix, 1) <> UBound(sourceMatrix, 2)) Then
            Throw New Exception("Matrix must be square.")
        End If
        Dim rank As Integer = UBound(sourceMatrix, 1)
        " ----- Clone a copy of the matrix (not just a new reference).
        Dim workMatrix(,) As Double = _
            CType(sourceMatrix.Clone, Double(,))
        " ----- Variables used for backsolving.
        Dim destMatrix(rank, rank) As Double
        Dim rightHandSide(rank) As Double
        Dim solutions(rank) As Double
        Dim rowPivots(rank) As Integer
        Dim colPivots(rank) As Integer
        " ----- Use LU decomposition to form a triangular matrix.
        workMatrix = FormLU(workMatrix, rowPivots, colPivots, rowsAndCols)
        " ----- Backsolve the triangular matrix to get the inverted
        "       value for each position in the final matrix.
        For eachCol = 0 To rank
            rightHandSide(eachCol) = 1
            BackSolve(workMatrix, rightHandSide, solutions, rowPivots, colPivots)
            For eachRow = 0 To rank
                destMatrix(eachRow, eachCol) = solutions(eachRow)
                rightHandSide(eachRow) = 0
            Next eachRow
        Next eachCol
        " ----- Return the inverted matrix result.
        Return destMatrix
    End Function
    Public Function Determinant(ByVal sourceMatrix(,) As Double) As Double
        " ----- Calculate the determinant of a matrix.
        Dim result As Double
        Dim pivots As Integer
        Dim count As Integer
        " ----- Only calculate the determinants of square matrices.
        If (UBound(sourceMatrix, 1) <> UBound(sourceMatrix, 2)) Then
            Throw New Exception( _
                "Determinant only calculated for square matrices.")
        End If
        Dim rank As Integer = UBound(sourceMatrix, 1)
        " ----- Make a copy of the matrix so we can work inside of it.
        Dim workMatrix(rank, rank) As Double
        Array.Copy(sourceMatrix, workMatrix, sourceMatrix.Length)
        " ----- Use LU decomposition to form a triangular matrix.
        Dim rowPivots(rank) As Integer
        Dim colPivots(rank) As Integer
        workMatrix = FormLU(workMatrix, rowPivots, colPivots, count)
        " ----- Get the product at each of the pivot points.
        result = 1
        For pivots = 0 To rank
            result *= workMatrix(rowPivots(pivots), colPivots(pivots))
        Next pivots
        " ----- Determine the sign of the result using LaPlace"s formula.
        result = (-1) ^ count * result
        Return result
    End Function
    Public Function SimultEq(ByVal sourceEquations(,) As Double, _
            ByVal sourceRHS() As Double) As Double()
        " ----- Use matrices to solve simultaneous equations.
        Dim rowsAndCols As Integer
        " ----- The matrix must be square and the array size must match.
        Dim rank As Integer = UBound(sourceEquations, 1)
        If (UBound(sourceEquations, 2) <> rank) Or _
                (UBound(sourceRHS, 1) <> rank) Then
            Throw New Exception("Size problem for simultaneous equations.")
        End If
        " ----- Create some arrays for doing all of the work.
        Dim coefficientMatrix(rank, rank) As Double
        Dim rightHandSide(rank) As Double
        Dim solutions(rank) As Double
        Dim rowPivots(rank) As Integer
        Dim colPivots(rank) As Integer
        " ----- Make copies of the original matrices so we don"t
        "       mess them up.
        Array.Copy(sourceEquations, coefficientMatrix, sourceEquations.Length)
        Array.Copy(sourceRHS, rightHandSide, sourceRHS.Length)
        " ----- Use LU decomposition to form a triangular matrix.
        coefficientMatrix = FormLU(coefficientMatrix, rowPivots, _
            colPivots, rowsAndCols)
        " ----- Find the unique solution for the upper-triangle.
        BackSolve(coefficientMatrix, rightHandSide, solutions, _
            rowPivots, colPivots)
        " ----- Return the simultaneous equations result in an array.
        Return solutions
    End Function
    Private Function FormLU(ByVal sourceMatrix(,) As Double, _
            ByRef rowPivots() As Integer, ByRef colPivots() As Integer, _
            ByRef rowsAndCols As Integer) As Double(,)
        " ----- Perform an LU (lower and upper) decomposition of a matrix,
        "       a modified form of Gaussian elimination.
        Dim eachRow As Integer
        Dim eachCol As Integer
        Dim pivot As Integer
        Dim rowIndex As Integer
        Dim colIndex As Integer
        Dim bestRow As Integer
        Dim bestCol As Integer
        Dim rowToPivot As Integer
        Dim colToPivot As Integer
        Dim maxValue As Double
        Dim testValue As Double
        Dim oldMax As Double
        Const Deps As Double = 0.0000000000000001
        " ----- Determine the size of the array.
        Dim rank As Integer = UBound(sourceMatrix, 1)
        Dim destMatrix(rank, rank) As Double
        Dim rowNorm(rank) As Double
        ReDim rowPivots(rank)
        ReDim colPivots(rank)
        " ----- Make a copy of the array so we don"t mess it up.
        Array.Copy(sourceMatrix, destMatrix, sourceMatrix.Length)
        " ----- Initialize row and column pivot arrays.
        For eachRow = 0 To rank
            rowPivots(eachRow) = eachRow
            colPivots(eachRow) = eachRow
            For eachCol = 0 To rank
                rowNorm(eachRow) += Math.Abs(destMatrix(eachRow, eachCol))
            Next eachCol
            If (rowNorm(eachRow) = 0) Then
                Throw New Exception("Cannot invert a singular matrix.")
            End If
        Next eachRow
        " ----- Use Gauss-Jordan elimination on the matrix rows.
        For pivot = 0 To rank - 1
            maxValue = 0
            For eachRow = pivot To rank
                rowIndex = rowPivots(eachRow)
                For eachCol = pivot To rank
                    colIndex = colPivots(eachCol)
                    testValue = Math.Abs(destMatrix(rowIndex, colIndex)) _
                        / rowNorm(rowIndex)
                    If (testValue > maxValue) Then
                        maxValue = testValue
                        bestRow = eachRow
                        bestCol = eachCol
                    End If
                Next eachCol
            Next eachRow
            " ----- Detect a singular, or very nearly singular, matrix.
            If (maxValue = 0) Then
                Throw New Exception("Singular matrix used for LU.")
            ElseIf (pivot > 1) Then
                If (maxValue < (Deps * oldMax)) Then
                    Throw New Exception("Non-invertible matrix used for LU.")
                End If
            End If
            oldMax = maxValue
            " ----- Swap row pivot values for the best row.
            If (rowPivots(pivot) <> rowPivots(bestRow)) Then
                rowsAndCols += 1
                Swap(rowPivots(pivot), rowPivots(bestRow))
            End If
            " ----- Swap column pivot values for the best column.
            If (colPivots(pivot) <> colPivots(bestCol)) Then
                rowsAndCols += 1
                Swap(colPivots(pivot), colPivots(bestCol))
            End If
            " ----- Work with the current pivot points.
            rowToPivot = rowPivots(pivot)
            colToPivot = colPivots(pivot)
            " ----- Modify the remaining rows from the pivot points.
            For eachRow = (pivot + 1) To rank
                rowIndex = rowPivots(eachRow)
                destMatrix(rowIndex, colToPivot) = _
                    -destMatrix(rowIndex, colToPivot) / _
                    destMatrix(rowToPivot, colToPivot)
                For eachCol = (pivot + 1) To rank
                    colIndex = colPivots(eachCol)
                    destMatrix(rowIndex, colIndex) += _
                        destMatrix(rowIndex, colToPivot) * _
                        destMatrix(rowToPivot, colIndex)
                Next eachCol
            Next eachRow
        Next pivot
        " ----- Detect a non-invertible matrix.
        If (destMatrix(rowPivots(rank), colPivots(rank)) = 0) Then
            Throw New Exception("Non-invertible matrix used for LU.")
        ElseIf (Math.Abs(destMatrix(rowPivots(rank), colPivots(rank))) / _
                rowNorm(rowPivots(rank))) < (Deps * oldMax) Then
            Throw New Exception("Non-invertible matrix used for LU.")
        End If
        " ----- Success. Return the LU triangular matrix.
        Return destMatrix
    End Function
    Private Sub Swap(ByRef firstValue As Integer, ByRef secondValue As Integer)
        " ----- Reverse the values of two reference integers.
        Dim holdValue As Integer
        holdValue = firstValue
        firstValue = secondValue
        secondValue = holdValue
    End Sub
    Private Sub BackSolve(ByVal sourceMatrix(,) As Double, _
            ByVal rightHandSide() As Double, ByVal solutions() As Double, _
            ByRef rowPivots() As Integer, ByRef colPivots() As Integer)
        " ----- Solve an upper-right-triangle matrix.
        Dim pivot As Integer
        Dim rowToPivot As Integer
        Dim colToPivot As Integer
        Dim eachRow As Integer
        Dim eachCol As Integer
        Dim rank As Integer = UBound(sourceMatrix, 1)
        " ----- Work through all pivot points. This section builds
        "       the "B" in the AX=B formula.
        For pivot = 0 To (rank - 1)
            colToPivot = colPivots(pivot)
            For eachRow = (pivot + 1) To rank
                rowToPivot = rowPivots(eachRow)
                rightHandSide(rowToPivot) += _
                    sourceMatrix(rowToPivot, colToPivot) _
                    * rightHandSide(rowPivots(pivot))
            Next eachRow
        Next pivot
        " ----- Now solve for each X using the general formula
        "       x(i) = (b(i) - summation(a(i,j)x(j)))/a(i,i)
        For eachRow = rank To 0 Step -1
            colToPivot = colPivots(eachRow)
            rowToPivot = rowPivots(eachRow)
            solutions(colToPivot) = rightHandSide(rowToPivot)
            For eachCol = (eachRow + 1) To rank
                solutions(colToPivot) -= _
                    sourceMatrix(rowToPivot, colPivots(eachCol)) _
                    * solutions(colPivots(eachCol))
            Next eachCol
            solutions(colToPivot) /= sourceMatrix(rowToPivot, colToPivot)
        Next eachRow
    End Sub
End Module

1,2,3
5,4,6
9,7,8
Determinant: 15

Display a Matrix

" Quote from
"Visual Basic 2005 Cookbook Solutions for VB 2005 Programmers
"by Tim Patrick (Author), John Craig (Author)
"# Publisher: O"Reilly Media, Inc. (September 21, 2006)
"# Language: English
"# ISBN-10: 0596101775
"# ISBN-13: 978-0596101770
Public Class Tester
    Public Shared Sub Main
        Dim matrixA(,) As Double = { _
            {4, 5, 6}, _
            {7, 8, 9}, _
            {3, 2, 1}}
        Console.WriteLine(MatrixHelper.MakeDisplayable(matrixA))
     End Sub
End Class

Module MatrixHelper
    Public Function MakeDisplayable(ByVal sourceMatrix(,) As Double) As String
        " ----- Prepare a multi-line string that shows the contents
        "       of a matrix, a 2D array.
        Dim rows As Integer
        Dim cols As Integer
        Dim eachRow As Integer
        Dim eachCol As Integer
        Dim result As New System.Text.StringBuilder
        " ----- Process all rows of the matrix, generating one
        "       output line per row.
        rows = UBound(sourceMatrix, 1) + 1
        cols = UBound(sourceMatrix, 2) + 1
        For eachRow = 0 To rows - 1
            " ----- Process each column of the matrix on a single
            "       row, separating values by commas.
            If (eachRow > 0) Then result.AppendLine()
            For eachCol = 0 To cols - 1
                " ----- Add a single matrix element to the output.
                If (eachCol > 0) Then result.Append(",")
                result.Append(sourceMatrix(eachRow, eachCol).ToString)
            Next eachCol
        Next eachRow
        " ----- Finished.
        Return result.ToString
    End Function
    Public Function MakeDisplayable(ByVal sourceArray() As Double) As String
        " ----- Present an array as multiple lines of output.
        Dim result As New System.Text.StringBuilder
        Dim scanValue As Double
        For Each scanValue In sourceArray
            result.AppendLine(scanValue.ToString)
        Next scanValue
        Return result.ToString
    End Function
    Public Function Inverse(ByVal sourceMatrix(,) As Double) As Double(,)
        " ----- Build a new matrix that is the mathematical inverse
        "       of the supplied matrix. Multiplying a matrix and its
        "       inverse together will give the identity matrix.
        Dim eachCol As Integer
        Dim eachRow As Integer
        Dim rowsAndCols As Integer
        " ----- Determine the size of each dimension of the matrix.
        "       Only square matrices can be inverted.
        If (UBound(sourceMatrix, 1) <> UBound(sourceMatrix, 2)) Then
            Throw New Exception("Matrix must be square.")
        End If
        Dim rank As Integer = UBound(sourceMatrix, 1)
        " ----- Clone a copy of the matrix (not just a new reference).
        Dim workMatrix(,) As Double = _
            CType(sourceMatrix.Clone, Double(,))
        " ----- Variables used for backsolving.
        Dim destMatrix(rank, rank) As Double
        Dim rightHandSide(rank) As Double
        Dim solutions(rank) As Double
        Dim rowPivots(rank) As Integer
        Dim colPivots(rank) As Integer
        " ----- Use LU decomposition to form a triangular matrix.
        workMatrix = FormLU(workMatrix, rowPivots, colPivots, rowsAndCols)
        " ----- Backsolve the triangular matrix to get the inverted
        "       value for each position in the final matrix.
        For eachCol = 0 To rank
            rightHandSide(eachCol) = 1
            BackSolve(workMatrix, rightHandSide, solutions, rowPivots, colPivots)
            For eachRow = 0 To rank
                destMatrix(eachRow, eachCol) = solutions(eachRow)
                rightHandSide(eachRow) = 0
            Next eachRow
        Next eachCol
        " ----- Return the inverted matrix result.
        Return destMatrix
    End Function
    Public Function Determinant(ByVal sourceMatrix(,) As Double) As Double
        " ----- Calculate the determinant of a matrix.
        Dim result As Double
        Dim pivots As Integer
        Dim count As Integer
        " ----- Only calculate the determinants of square matrices.
        If (UBound(sourceMatrix, 1) <> UBound(sourceMatrix, 2)) Then
            Throw New Exception( _
                "Determinant only calculated for square matrices.")
        End If
        Dim rank As Integer = UBound(sourceMatrix, 1)
        " ----- Make a copy of the matrix so we can work inside of it.
        Dim workMatrix(rank, rank) As Double
        Array.Copy(sourceMatrix, workMatrix, sourceMatrix.Length)
        " ----- Use LU decomposition to form a triangular matrix.
        Dim rowPivots(rank) As Integer
        Dim colPivots(rank) As Integer
        workMatrix = FormLU(workMatrix, rowPivots, colPivots, count)
        " ----- Get the product at each of the pivot points.
        result = 1
        For pivots = 0 To rank
            result *= workMatrix(rowPivots(pivots), colPivots(pivots))
        Next pivots
        " ----- Determine the sign of the result using LaPlace"s formula.
        result = (-1) ^ count * result
        Return result
    End Function
    Public Function SimultEq(ByVal sourceEquations(,) As Double, _
            ByVal sourceRHS() As Double) As Double()
        " ----- Use matrices to solve simultaneous equations.
        Dim rowsAndCols As Integer
        " ----- The matrix must be square and the array size must match.
        Dim rank As Integer = UBound(sourceEquations, 1)
        If (UBound(sourceEquations, 2) <> rank) Or _
                (UBound(sourceRHS, 1) <> rank) Then
            Throw New Exception("Size problem for simultaneous equations.")
        End If
        " ----- Create some arrays for doing all of the work.
        Dim coefficientMatrix(rank, rank) As Double
        Dim rightHandSide(rank) As Double
        Dim solutions(rank) As Double
        Dim rowPivots(rank) As Integer
        Dim colPivots(rank) As Integer
        " ----- Make copies of the original matrices so we don"t
        "       mess them up.
        Array.Copy(sourceEquations, coefficientMatrix, sourceEquations.Length)
        Array.Copy(sourceRHS, rightHandSide, sourceRHS.Length)
        " ----- Use LU decomposition to form a triangular matrix.
        coefficientMatrix = FormLU(coefficientMatrix, rowPivots, _
            colPivots, rowsAndCols)
        " ----- Find the unique solution for the upper-triangle.
        BackSolve(coefficientMatrix, rightHandSide, solutions, _
            rowPivots, colPivots)
        " ----- Return the simultaneous equations result in an array.
        Return solutions
    End Function
    Private Function FormLU(ByVal sourceMatrix(,) As Double, _
            ByRef rowPivots() As Integer, ByRef colPivots() As Integer, _
            ByRef rowsAndCols As Integer) As Double(,)
        " ----- Perform an LU (lower and upper) decomposition of a matrix,
        "       a modified form of Gaussian elimination.
        Dim eachRow As Integer
        Dim eachCol As Integer
        Dim pivot As Integer
        Dim rowIndex As Integer
        Dim colIndex As Integer
        Dim bestRow As Integer
        Dim bestCol As Integer
        Dim rowToPivot As Integer
        Dim colToPivot As Integer
        Dim maxValue As Double
        Dim testValue As Double
        Dim oldMax As Double
        Const Deps As Double = 0.0000000000000001
        " ----- Determine the size of the array.
        Dim rank As Integer = UBound(sourceMatrix, 1)
        Dim destMatrix(rank, rank) As Double
        Dim rowNorm(rank) As Double
        ReDim rowPivots(rank)
        ReDim colPivots(rank)
        " ----- Make a copy of the array so we don"t mess it up.
        Array.Copy(sourceMatrix, destMatrix, sourceMatrix.Length)
        " ----- Initialize row and column pivot arrays.
        For eachRow = 0 To rank
            rowPivots(eachRow) = eachRow
            colPivots(eachRow) = eachRow
            For eachCol = 0 To rank
                rowNorm(eachRow) += Math.Abs(destMatrix(eachRow, eachCol))
            Next eachCol
            If (rowNorm(eachRow) = 0) Then
                Throw New Exception("Cannot invert a singular matrix.")
            End If
        Next eachRow
        " ----- Use Gauss-Jordan elimination on the matrix rows.
        For pivot = 0 To rank - 1
            maxValue = 0
            For eachRow = pivot To rank
                rowIndex = rowPivots(eachRow)
                For eachCol = pivot To rank
                    colIndex = colPivots(eachCol)
                    testValue = Math.Abs(destMatrix(rowIndex, colIndex)) _
                        / rowNorm(rowIndex)
                    If (testValue > maxValue) Then
                        maxValue = testValue
                        bestRow = eachRow
                        bestCol = eachCol
                    End If
                Next eachCol
            Next eachRow
            " ----- Detect a singular, or very nearly singular, matrix.
            If (maxValue = 0) Then
                Throw New Exception("Singular matrix used for LU.")
            ElseIf (pivot > 1) Then
                If (maxValue < (Deps * oldMax)) Then
                    Throw New Exception("Non-invertible matrix used for LU.")
                End If
            End If
            oldMax = maxValue
            " ----- Swap row pivot values for the best row.
            If (rowPivots(pivot) <> rowPivots(bestRow)) Then
                rowsAndCols += 1
                Swap(rowPivots(pivot), rowPivots(bestRow))
            End If
            " ----- Swap column pivot values for the best column.
            If (colPivots(pivot) <> colPivots(bestCol)) Then
                rowsAndCols += 1
                Swap(colPivots(pivot), colPivots(bestCol))
            End If
            " ----- Work with the current pivot points.
            rowToPivot = rowPivots(pivot)
            colToPivot = colPivots(pivot)
            " ----- Modify the remaining rows from the pivot points.
            For eachRow = (pivot + 1) To rank
                rowIndex = rowPivots(eachRow)
                destMatrix(rowIndex, colToPivot) = _
                    -destMatrix(rowIndex, colToPivot) / _
                    destMatrix(rowToPivot, colToPivot)
                For eachCol = (pivot + 1) To rank
                    colIndex = colPivots(eachCol)
                    destMatrix(rowIndex, colIndex) += _
                        destMatrix(rowIndex, colToPivot) * _
                        destMatrix(rowToPivot, colIndex)
                Next eachCol
            Next eachRow
        Next pivot
        " ----- Detect a non-invertible matrix.
        If (destMatrix(rowPivots(rank), colPivots(rank)) = 0) Then
            Throw New Exception("Non-invertible matrix used for LU.")
        ElseIf (Math.Abs(destMatrix(rowPivots(rank), colPivots(rank))) / _
                rowNorm(rowPivots(rank))) < (Deps * oldMax) Then
            Throw New Exception("Non-invertible matrix used for LU.")
        End If
        " ----- Success. Return the LU triangular matrix.
        Return destMatrix
    End Function
    Private Sub Swap(ByRef firstValue As Integer, ByRef secondValue As Integer)
        " ----- Reverse the values of two reference integers.
        Dim holdValue As Integer
        holdValue = firstValue
        firstValue = secondValue
        secondValue = holdValue
    End Sub
    Private Sub BackSolve(ByVal sourceMatrix(,) As Double, _
            ByVal rightHandSide() As Double, ByVal solutions() As Double, _
            ByRef rowPivots() As Integer, ByRef colPivots() As Integer)
        " ----- Solve an upper-right-triangle matrix.
        Dim pivot As Integer
        Dim rowToPivot As Integer
        Dim colToPivot As Integer
        Dim eachRow As Integer
        Dim eachCol As Integer
        Dim rank As Integer = UBound(sourceMatrix, 1)
        " ----- Work through all pivot points. This section builds
        "       the "B" in the AX=B formula.
        For pivot = 0 To (rank - 1)
            colToPivot = colPivots(pivot)
            For eachRow = (pivot + 1) To rank
                rowToPivot = rowPivots(eachRow)
                rightHandSide(rowToPivot) += _
                    sourceMatrix(rowToPivot, colToPivot) _
                    * rightHandSide(rowPivots(pivot))
            Next eachRow
        Next pivot
        " ----- Now solve for each X using the general formula
        "       x(i) = (b(i) - summation(a(i,j)x(j)))/a(i,i)
        For eachRow = rank To 0 Step -1
            colToPivot = colPivots(eachRow)
            rowToPivot = rowPivots(eachRow)
            solutions(colToPivot) = rightHandSide(rowToPivot)
            For eachCol = (eachRow + 1) To rank
                solutions(colToPivot) -= _
                    sourceMatrix(rowToPivot, colPivots(eachCol)) _
                    * solutions(colPivots(eachCol))
            Next eachCol
            solutions(colToPivot) /= sourceMatrix(rowToPivot, colToPivot)
        Next eachRow
    End Sub
End Module

4,5,6
7,8,9
3,2,1

Invert a matrix

" Quote from
"Visual Basic 2005 Cookbook Solutions for VB 2005 Programmers
"by Tim Patrick (Author), John Craig (Author)
"# Publisher: O"Reilly Media, Inc. (September 21, 2006)
"# Language: English
"# ISBN-10: 0596101775
"# ISBN-13: 978-0596101770

Public Class Tester
    Public Shared Sub Main
        Dim matrixA(,) As Double = { _
            {1, 3, 3}, _
            {2, 4, 3}, _
            {1, 3, 4}}
        Dim matrixB(,) As Double = MatrixHelper.Inverse(matrixA)
        Console.WriteLine(MatrixHelper.MakeDisplayable(matrixA) & _
            vbNewLine & vbNewLine & "Inverse: " & _
            vbNewLine & MatrixHelper.MakeDisplayable(matrixB))
     End Sub
End Class

Module MatrixHelper
    Public Function MakeDisplayable(ByVal sourceMatrix(,) As Double) As String
        " ----- Prepare a multi-line string that shows the contents
        "       of a matrix, a 2D array.
        Dim rows As Integer
        Dim cols As Integer
        Dim eachRow As Integer
        Dim eachCol As Integer
        Dim result As New System.Text.StringBuilder
        " ----- Process all rows of the matrix, generating one
        "       output line per row.
        rows = UBound(sourceMatrix, 1) + 1
        cols = UBound(sourceMatrix, 2) + 1
        For eachRow = 0 To rows - 1
            " ----- Process each column of the matrix on a single
            "       row, separating values by commas.
            If (eachRow > 0) Then result.AppendLine()
            For eachCol = 0 To cols - 1
                " ----- Add a single matrix element to the output.
                If (eachCol > 0) Then result.Append(",")
                result.Append(sourceMatrix(eachRow, eachCol).ToString)
            Next eachCol
        Next eachRow
        " ----- Finished.
        Return result.ToString
    End Function
    Public Function MakeDisplayable(ByVal sourceArray() As Double) As String
        " ----- Present an array as multiple lines of output.
        Dim result As New System.Text.StringBuilder
        Dim scanValue As Double
        For Each scanValue In sourceArray
            result.AppendLine(scanValue.ToString)
        Next scanValue
        Return result.ToString
    End Function
    Public Function Inverse(ByVal sourceMatrix(,) As Double) As Double(,)
        " ----- Build a new matrix that is the mathematical inverse
        "       of the supplied matrix. Multiplying a matrix and its
        "       inverse together will give the identity matrix.
        Dim eachCol As Integer
        Dim eachRow As Integer
        Dim rowsAndCols As Integer
        " ----- Determine the size of each dimension of the matrix.
        "       Only square matrices can be inverted.
        If (UBound(sourceMatrix, 1) <> UBound(sourceMatrix, 2)) Then
            Throw New Exception("Matrix must be square.")
        End If
        Dim rank As Integer = UBound(sourceMatrix, 1)
        " ----- Clone a copy of the matrix (not just a new reference).
        Dim workMatrix(,) As Double = _
            CType(sourceMatrix.Clone, Double(,))
        " ----- Variables used for backsolving.
        Dim destMatrix(rank, rank) As Double
        Dim rightHandSide(rank) As Double
        Dim solutions(rank) As Double
        Dim rowPivots(rank) As Integer
        Dim colPivots(rank) As Integer
        " ----- Use LU decomposition to form a triangular matrix.
        workMatrix = FormLU(workMatrix, rowPivots, colPivots, rowsAndCols)
        " ----- Backsolve the triangular matrix to get the inverted
        "       value for each position in the final matrix.
        For eachCol = 0 To rank
            rightHandSide(eachCol) = 1
            BackSolve(workMatrix, rightHandSide, solutions, rowPivots, colPivots)
            For eachRow = 0 To rank
                destMatrix(eachRow, eachCol) = solutions(eachRow)
                rightHandSide(eachRow) = 0
            Next eachRow
        Next eachCol
        " ----- Return the inverted matrix result.
        Return destMatrix
    End Function
    Public Function Determinant(ByVal sourceMatrix(,) As Double) As Double
        " ----- Calculate the determinant of a matrix.
        Dim result As Double
        Dim pivots As Integer
        Dim count As Integer
        " ----- Only calculate the determinants of square matrices.
        If (UBound(sourceMatrix, 1) <> UBound(sourceMatrix, 2)) Then
            Throw New Exception( _
                "Determinant only calculated for square matrices.")
        End If
        Dim rank As Integer = UBound(sourceMatrix, 1)
        " ----- Make a copy of the matrix so we can work inside of it.
        Dim workMatrix(rank, rank) As Double
        Array.Copy(sourceMatrix, workMatrix, sourceMatrix.Length)
        " ----- Use LU decomposition to form a triangular matrix.
        Dim rowPivots(rank) As Integer
        Dim colPivots(rank) As Integer
        workMatrix = FormLU(workMatrix, rowPivots, colPivots, count)
        " ----- Get the product at each of the pivot points.
        result = 1
        For pivots = 0 To rank
            result *= workMatrix(rowPivots(pivots), colPivots(pivots))
        Next pivots
        " ----- Determine the sign of the result using LaPlace"s formula.
        result = (-1) ^ count * result
        Return result
    End Function
    Public Function SimultEq(ByVal sourceEquations(,) As Double, _
            ByVal sourceRHS() As Double) As Double()
        " ----- Use matrices to solve simultaneous equations.
        Dim rowsAndCols As Integer
        " ----- The matrix must be square and the array size must match.
        Dim rank As Integer = UBound(sourceEquations, 1)
        If (UBound(sourceEquations, 2) <> rank) Or _
                (UBound(sourceRHS, 1) <> rank) Then
            Throw New Exception("Size problem for simultaneous equations.")
        End If
        " ----- Create some arrays for doing all of the work.
        Dim coefficientMatrix(rank, rank) As Double
        Dim rightHandSide(rank) As Double
        Dim solutions(rank) As Double
        Dim rowPivots(rank) As Integer
        Dim colPivots(rank) As Integer
        " ----- Make copies of the original matrices so we don"t
        "       mess them up.
        Array.Copy(sourceEquations, coefficientMatrix, sourceEquations.Length)
        Array.Copy(sourceRHS, rightHandSide, sourceRHS.Length)
        " ----- Use LU decomposition to form a triangular matrix.
        coefficientMatrix = FormLU(coefficientMatrix, rowPivots, _
            colPivots, rowsAndCols)
        " ----- Find the unique solution for the upper-triangle.
        BackSolve(coefficientMatrix, rightHandSide, solutions, _
            rowPivots, colPivots)
        " ----- Return the simultaneous equations result in an array.
        Return solutions
    End Function
    Private Function FormLU(ByVal sourceMatrix(,) As Double, _
            ByRef rowPivots() As Integer, ByRef colPivots() As Integer, _
            ByRef rowsAndCols As Integer) As Double(,)
        " ----- Perform an LU (lower and upper) decomposition of a matrix,
        "       a modified form of Gaussian elimination.
        Dim eachRow As Integer
        Dim eachCol As Integer
        Dim pivot As Integer
        Dim rowIndex As Integer
        Dim colIndex As Integer
        Dim bestRow As Integer
        Dim bestCol As Integer
        Dim rowToPivot As Integer
        Dim colToPivot As Integer
        Dim maxValue As Double
        Dim testValue As Double
        Dim oldMax As Double
        Const Deps As Double = 0.0000000000000001
        " ----- Determine the size of the array.
        Dim rank As Integer = UBound(sourceMatrix, 1)
        Dim destMatrix(rank, rank) As Double
        Dim rowNorm(rank) As Double
        ReDim rowPivots(rank)
        ReDim colPivots(rank)
        " ----- Make a copy of the array so we don"t mess it up.
        Array.Copy(sourceMatrix, destMatrix, sourceMatrix.Length)
        " ----- Initialize row and column pivot arrays.
        For eachRow = 0 To rank
            rowPivots(eachRow) = eachRow
            colPivots(eachRow) = eachRow
            For eachCol = 0 To rank
                rowNorm(eachRow) += Math.Abs(destMatrix(eachRow, eachCol))
            Next eachCol
            If (rowNorm(eachRow) = 0) Then
                Throw New Exception("Cannot invert a singular matrix.")
            End If
        Next eachRow
        " ----- Use Gauss-Jordan elimination on the matrix rows.
        For pivot = 0 To rank - 1
            maxValue = 0
            For eachRow = pivot To rank
                rowIndex = rowPivots(eachRow)
                For eachCol = pivot To rank
                    colIndex = colPivots(eachCol)
                    testValue = Math.Abs(destMatrix(rowIndex, colIndex)) _
                        / rowNorm(rowIndex)
                    If (testValue > maxValue) Then
                        maxValue = testValue
                        bestRow = eachRow
                        bestCol = eachCol
                    End If
                Next eachCol
            Next eachRow
            " ----- Detect a singular, or very nearly singular, matrix.
            If (maxValue = 0) Then
                Throw New Exception("Singular matrix used for LU.")
            ElseIf (pivot > 1) Then
                If (maxValue < (Deps * oldMax)) Then
                    Throw New Exception("Non-invertible matrix used for LU.")
                End If
            End If
            oldMax = maxValue
            " ----- Swap row pivot values for the best row.
            If (rowPivots(pivot) <> rowPivots(bestRow)) Then
                rowsAndCols += 1
                Swap(rowPivots(pivot), rowPivots(bestRow))
            End If
            " ----- Swap column pivot values for the best column.
            If (colPivots(pivot) <> colPivots(bestCol)) Then
                rowsAndCols += 1
                Swap(colPivots(pivot), colPivots(bestCol))
            End If
            " ----- Work with the current pivot points.
            rowToPivot = rowPivots(pivot)
            colToPivot = colPivots(pivot)
            " ----- Modify the remaining rows from the pivot points.
            For eachRow = (pivot + 1) To rank
                rowIndex = rowPivots(eachRow)
                destMatrix(rowIndex, colToPivot) = _
                    -destMatrix(rowIndex, colToPivot) / _
                    destMatrix(rowToPivot, colToPivot)
                For eachCol = (pivot + 1) To rank
                    colIndex = colPivots(eachCol)
                    destMatrix(rowIndex, colIndex) += _
                        destMatrix(rowIndex, colToPivot) * _
                        destMatrix(rowToPivot, colIndex)
                Next eachCol
            Next eachRow
        Next pivot
        " ----- Detect a non-invertible matrix.
        If (destMatrix(rowPivots(rank), colPivots(rank)) = 0) Then
            Throw New Exception("Non-invertible matrix used for LU.")
        ElseIf (Math.Abs(destMatrix(rowPivots(rank), colPivots(rank))) / _
                rowNorm(rowPivots(rank))) < (Deps * oldMax) Then
            Throw New Exception("Non-invertible matrix used for LU.")
        End If
        " ----- Success. Return the LU triangular matrix.
        Return destMatrix
    End Function
    Private Sub Swap(ByRef firstValue As Integer, ByRef secondValue As Integer)
        " ----- Reverse the values of two reference integers.
        Dim holdValue As Integer
        holdValue = firstValue
        firstValue = secondValue
        secondValue = holdValue
    End Sub
    Private Sub BackSolve(ByVal sourceMatrix(,) As Double, _
            ByVal rightHandSide() As Double, ByVal solutions() As Double, _
            ByRef rowPivots() As Integer, ByRef colPivots() As Integer)
        " ----- Solve an upper-right-triangle matrix.
        Dim pivot As Integer
        Dim rowToPivot As Integer
        Dim colToPivot As Integer
        Dim eachRow As Integer
        Dim eachCol As Integer
        Dim rank As Integer = UBound(sourceMatrix, 1)
        " ----- Work through all pivot points. This section builds
        "       the "B" in the AX=B formula.
        For pivot = 0 To (rank - 1)
            colToPivot = colPivots(pivot)
            For eachRow = (pivot + 1) To rank
                rowToPivot = rowPivots(eachRow)
                rightHandSide(rowToPivot) += _
                    sourceMatrix(rowToPivot, colToPivot) _
                    * rightHandSide(rowPivots(pivot))
            Next eachRow
        Next pivot
        " ----- Now solve for each X using the general formula
        "       x(i) = (b(i) - summation(a(i,j)x(j)))/a(i,i)
        For eachRow = rank To 0 Step -1
            colToPivot = colPivots(eachRow)
            rowToPivot = rowPivots(eachRow)
            solutions(colToPivot) = rightHandSide(rowToPivot)
            For eachCol = (eachRow + 1) To rank
                solutions(colToPivot) -= _
                    sourceMatrix(rowToPivot, colPivots(eachCol)) _
                    * solutions(colPivots(eachCol))
            Next eachCol
            solutions(colToPivot) /= sourceMatrix(rowToPivot, colToPivot)
        Next eachRow
    End Sub
End Module

1,3,3
2,4,3
1,3,4
Inverse:
-3.5,1.5,1.5
2.5,-0.5,-1.5
-1,0,1

Solve equations using matrices

" Quote from
"Visual Basic 2005 Cookbook Solutions for VB 2005 Programmers
"by Tim Patrick (Author), John Craig (Author)
"# Publisher: O"Reilly Media, Inc. (September 21, 2006)
"# Language: English
"# ISBN-10: 0596101775
"# ISBN-13: 978-0596101770
Public Class Tester
    Public Shared Sub Main
        Dim matrixA(,) As Double = { _
            {1, 1, 1, 1}, _
            {1, 5, 10, 25}, _
            {0, 5, 10, 0}, _
            {0, 0, 10, 25}}
        Dim arrayB() As Double = {18, 223, 70, 200}
        Dim arrayC() As Double = MatrixHelper.SimultEq(matrixA, arrayB)
        Console.WriteLine(MatrixHelper.MakeDisplayable(matrixA) & _
           vbNewLine & vbNewLine & MatrixHelper.MakeDisplayable(arrayB) & _
           vbNewLine & vbNewLine & "Simultaneous Equations Solution:" & _
           vbNewLine & MatrixHelper.MakeDisplayable(arrayC))
     End Sub
End Class

Module MatrixHelper
    Public Function MakeDisplayable(ByVal sourceMatrix(,) As Double) As String
        " ----- Prepare a multi-line string that shows the contents
        "       of a matrix, a 2D array.
        Dim rows As Integer
        Dim cols As Integer
        Dim eachRow As Integer
        Dim eachCol As Integer
        Dim result As New System.Text.StringBuilder
        " ----- Process all rows of the matrix, generating one
        "       output line per row.
        rows = UBound(sourceMatrix, 1) + 1
        cols = UBound(sourceMatrix, 2) + 1
        For eachRow = 0 To rows - 1
            " ----- Process each column of the matrix on a single
            "       row, separating values by commas.
            If (eachRow > 0) Then result.AppendLine()
            For eachCol = 0 To cols - 1
                " ----- Add a single matrix element to the output.
                If (eachCol > 0) Then result.Append(",")
                result.Append(sourceMatrix(eachRow, eachCol).ToString)
            Next eachCol
        Next eachRow
        " ----- Finished.
        Return result.ToString
    End Function
    Public Function MakeDisplayable(ByVal sourceArray() As Double) As String
        " ----- Present an array as multiple lines of output.
        Dim result As New System.Text.StringBuilder
        Dim scanValue As Double
        For Each scanValue In sourceArray
            result.AppendLine(scanValue.ToString)
        Next scanValue
        Return result.ToString
    End Function
    Public Function Inverse(ByVal sourceMatrix(,) As Double) As Double(,)
        " ----- Build a new matrix that is the mathematical inverse
        "       of the supplied matrix. Multiplying a matrix and its
        "       inverse together will give the identity matrix.
        Dim eachCol As Integer
        Dim eachRow As Integer
        Dim rowsAndCols As Integer
        " ----- Determine the size of each dimension of the matrix.
        "       Only square matrices can be inverted.
        If (UBound(sourceMatrix, 1) <> UBound(sourceMatrix, 2)) Then
            Throw New Exception("Matrix must be square.")
        End If
        Dim rank As Integer = UBound(sourceMatrix, 1)
        " ----- Clone a copy of the matrix (not just a new reference).
        Dim workMatrix(,) As Double = _
            CType(sourceMatrix.Clone, Double(,))
        " ----- Variables used for backsolving.
        Dim destMatrix(rank, rank) As Double
        Dim rightHandSide(rank) As Double
        Dim solutions(rank) As Double
        Dim rowPivots(rank) As Integer
        Dim colPivots(rank) As Integer
        " ----- Use LU decomposition to form a triangular matrix.
        workMatrix = FormLU(workMatrix, rowPivots, colPivots, rowsAndCols)
        " ----- Backsolve the triangular matrix to get the inverted
        "       value for each position in the final matrix.
        For eachCol = 0 To rank
            rightHandSide(eachCol) = 1
            BackSolve(workMatrix, rightHandSide, solutions, rowPivots, colPivots)
            For eachRow = 0 To rank
                destMatrix(eachRow, eachCol) = solutions(eachRow)
                rightHandSide(eachRow) = 0
            Next eachRow
        Next eachCol
        " ----- Return the inverted matrix result.
        Return destMatrix
    End Function
    Public Function Determinant(ByVal sourceMatrix(,) As Double) As Double
        " ----- Calculate the determinant of a matrix.
        Dim result As Double
        Dim pivots As Integer
        Dim count As Integer
        " ----- Only calculate the determinants of square matrices.
        If (UBound(sourceMatrix, 1) <> UBound(sourceMatrix, 2)) Then
            Throw New Exception( _
                "Determinant only calculated for square matrices.")
        End If
        Dim rank As Integer = UBound(sourceMatrix, 1)
        " ----- Make a copy of the matrix so we can work inside of it.
        Dim workMatrix(rank, rank) As Double
        Array.Copy(sourceMatrix, workMatrix, sourceMatrix.Length)
        " ----- Use LU decomposition to form a triangular matrix.
        Dim rowPivots(rank) As Integer
        Dim colPivots(rank) As Integer
        workMatrix = FormLU(workMatrix, rowPivots, colPivots, count)
        " ----- Get the product at each of the pivot points.
        result = 1
        For pivots = 0 To rank
            result *= workMatrix(rowPivots(pivots), colPivots(pivots))
        Next pivots
        " ----- Determine the sign of the result using LaPlace"s formula.
        result = (-1) ^ count * result
        Return result
    End Function
    Public Function SimultEq(ByVal sourceEquations(,) As Double, _
            ByVal sourceRHS() As Double) As Double()
        " ----- Use matrices to solve simultaneous equations.
        Dim rowsAndCols As Integer
        " ----- The matrix must be square and the array size must match.
        Dim rank As Integer = UBound(sourceEquations, 1)
        If (UBound(sourceEquations, 2) <> rank) Or _
                (UBound(sourceRHS, 1) <> rank) Then
            Throw New Exception("Size problem for simultaneous equations.")
        End If
        " ----- Create some arrays for doing all of the work.
        Dim coefficientMatrix(rank, rank) As Double
        Dim rightHandSide(rank) As Double
        Dim solutions(rank) As Double
        Dim rowPivots(rank) As Integer
        Dim colPivots(rank) As Integer
        " ----- Make copies of the original matrices so we don"t
        "       mess them up.
        Array.Copy(sourceEquations, coefficientMatrix, sourceEquations.Length)
        Array.Copy(sourceRHS, rightHandSide, sourceRHS.Length)
        " ----- Use LU decomposition to form a triangular matrix.
        coefficientMatrix = FormLU(coefficientMatrix, rowPivots, _
            colPivots, rowsAndCols)
        " ----- Find the unique solution for the upper-triangle.
        BackSolve(coefficientMatrix, rightHandSide, solutions, _
            rowPivots, colPivots)
        " ----- Return the simultaneous equations result in an array.
        Return solutions
    End Function
    Private Function FormLU(ByVal sourceMatrix(,) As Double, _
            ByRef rowPivots() As Integer, ByRef colPivots() As Integer, _
            ByRef rowsAndCols As Integer) As Double(,)
        " ----- Perform an LU (lower and upper) decomposition of a matrix,
        "       a modified form of Gaussian elimination.
        Dim eachRow As Integer
        Dim eachCol As Integer
        Dim pivot As Integer
        Dim rowIndex As Integer
        Dim colIndex As Integer
        Dim bestRow As Integer
        Dim bestCol As Integer
        Dim rowToPivot As Integer
        Dim colToPivot As Integer
        Dim maxValue As Double
        Dim testValue As Double
        Dim oldMax As Double
        Const Deps As Double = 0.0000000000000001
        " ----- Determine the size of the array.
        Dim rank As Integer = UBound(sourceMatrix, 1)
        Dim destMatrix(rank, rank) As Double
        Dim rowNorm(rank) As Double
        ReDim rowPivots(rank)
        ReDim colPivots(rank)
        " ----- Make a copy of the array so we don"t mess it up.
        Array.Copy(sourceMatrix, destMatrix, sourceMatrix.Length)
        " ----- Initialize row and column pivot arrays.
        For eachRow = 0 To rank
            rowPivots(eachRow) = eachRow
            colPivots(eachRow) = eachRow
            For eachCol = 0 To rank
                rowNorm(eachRow) += Math.Abs(destMatrix(eachRow, eachCol))
            Next eachCol
            If (rowNorm(eachRow) = 0) Then
                Throw New Exception("Cannot invert a singular matrix.")
            End If
        Next eachRow
        " ----- Use Gauss-Jordan elimination on the matrix rows.
        For pivot = 0 To rank - 1
            maxValue = 0
            For eachRow = pivot To rank
                rowIndex = rowPivots(eachRow)
                For eachCol = pivot To rank
                    colIndex = colPivots(eachCol)
                    testValue = Math.Abs(destMatrix(rowIndex, colIndex)) _
                        / rowNorm(rowIndex)
                    If (testValue > maxValue) Then
                        maxValue = testValue
                        bestRow = eachRow
                        bestCol = eachCol
                    End If
                Next eachCol
            Next eachRow
            " ----- Detect a singular, or very nearly singular, matrix.
            If (maxValue = 0) Then
                Throw New Exception("Singular matrix used for LU.")
            ElseIf (pivot > 1) Then
                If (maxValue < (Deps * oldMax)) Then
                    Throw New Exception("Non-invertible matrix used for LU.")
                End If
            End If
            oldMax = maxValue
            " ----- Swap row pivot values for the best row.
            If (rowPivots(pivot) <> rowPivots(bestRow)) Then
                rowsAndCols += 1
                Swap(rowPivots(pivot), rowPivots(bestRow))
            End If
            " ----- Swap column pivot values for the best column.
            If (colPivots(pivot) <> colPivots(bestCol)) Then
                rowsAndCols += 1
                Swap(colPivots(pivot), colPivots(bestCol))
            End If
            " ----- Work with the current pivot points.
            rowToPivot = rowPivots(pivot)
            colToPivot = colPivots(pivot)
            " ----- Modify the remaining rows from the pivot points.
            For eachRow = (pivot + 1) To rank
                rowIndex = rowPivots(eachRow)
                destMatrix(rowIndex, colToPivot) = _
                    -destMatrix(rowIndex, colToPivot) / _
                    destMatrix(rowToPivot, colToPivot)
                For eachCol = (pivot + 1) To rank
                    colIndex = colPivots(eachCol)
                    destMatrix(rowIndex, colIndex) += _
                        destMatrix(rowIndex, colToPivot) * _
                        destMatrix(rowToPivot, colIndex)
                Next eachCol
            Next eachRow
        Next pivot
        " ----- Detect a non-invertible matrix.
        If (destMatrix(rowPivots(rank), colPivots(rank)) = 0) Then
            Throw New Exception("Non-invertible matrix used for LU.")
        ElseIf (Math.Abs(destMatrix(rowPivots(rank), colPivots(rank))) / _
                rowNorm(rowPivots(rank))) < (Deps * oldMax) Then
            Throw New Exception("Non-invertible matrix used for LU.")
        End If
        " ----- Success. Return the LU triangular matrix.
        Return destMatrix
    End Function
    Private Sub Swap(ByRef firstValue As Integer, ByRef secondValue As Integer)
        " ----- Reverse the values of two reference integers.
        Dim holdValue As Integer
        holdValue = firstValue
        firstValue = secondValue
        secondValue = holdValue
    End Sub
    Private Sub BackSolve(ByVal sourceMatrix(,) As Double, _
            ByVal rightHandSide() As Double, ByVal solutions() As Double, _
            ByRef rowPivots() As Integer, ByRef colPivots() As Integer)
        " ----- Solve an upper-right-triangle matrix.
        Dim pivot As Integer
        Dim rowToPivot As Integer
        Dim colToPivot As Integer
        Dim eachRow As Integer
        Dim eachCol As Integer
        Dim rank As Integer = UBound(sourceMatrix, 1)
        " ----- Work through all pivot points. This section builds
        "       the "B" in the AX=B formula.
        For pivot = 0 To (rank - 1)
            colToPivot = colPivots(pivot)
            For eachRow = (pivot + 1) To rank
                rowToPivot = rowPivots(eachRow)
                rightHandSide(rowToPivot) += _
                    sourceMatrix(rowToPivot, colToPivot) _
                    * rightHandSide(rowPivots(pivot))
            Next eachRow
        Next pivot
        " ----- Now solve for each X using the general formula
        "       x(i) = (b(i) - summation(a(i,j)x(j)))/a(i,i)
        For eachRow = rank To 0 Step -1
            colToPivot = colPivots(eachRow)
            rowToPivot = rowPivots(eachRow)
            solutions(colToPivot) = rightHandSide(rowToPivot)
            For eachCol = (eachRow + 1) To rank
                solutions(colToPivot) -= _
                    sourceMatrix(rowToPivot, colPivots(eachCol)) _
                    * solutions(colPivots(eachCol))
            Next eachCol
            solutions(colToPivot) /= sourceMatrix(rowToPivot, colToPivot)
        Next eachRow
    End Sub
End Module

1,1,1,1
1,5,10,25
0,5,10,0
0,0,10,25
18
223
70
200

Simultaneous Equations Solution:
3
4
5
6