VB.Net Tutorial/Data Type/Matrix — различия между версиями

Материал из VB Эксперт
Перейти к: навигация, поиск
м (1 версия)
 
(нет различий)

Текущая версия на 15:54, 26 мая 2010

Calculate the determinant of an array

<source lang="vbnet">" 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</source>

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

Display a Matrix

<source lang="vbnet">" 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</source>

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

Invert a matrix

<source lang="vbnet">" 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</source>

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

<source lang="vbnet">" 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</source>

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