Recipe 6.35. Listing of the MatrixHelper Class

Sample code folder: Chapter 06\Matrix

Following is the full code for the MatrixHelper class described in Recipes 6.31, 6.32, 6.33 through 6.34:

 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