A sub-routine for the in-situ inversion of a square matrix

When a programmer deals with mathematics, one of the most common routines added to his/her math library is matrix inversion. The following code implements a very efficient in-situ matrix inversion (meaning that no memory is allocated for another matrix to hold the result; the inverse of the matrix is built gradually in the matrix itself!). The code is based on an algorithm published in Numerical Recipes.

Classic VB (vb6) code

The code is written in classic vb (vb6). I am afraid that I do not have the time to convert it to .NET. Have fun:

Public Function TLMatrixInvertGaussJordan(mMatrix1() As Double, gDim As Long) As Long


               On Error GoTo Er


    'Matrix inversion with Gauss-Jordan elimination and full pivoting.

    'mMatrix1 is transformed gradually to its inverse for reasons of memory usage :-)

    'if you need to keep the original matrix, copy it to another matrix before inverting!

    'Code taken from Numerical Recipes and modified for Visual Basic

    'For this version we removed the solution vectors, because we only need the inverse of the matrix



    Dim N As Long

    Dim nP As Long

    Dim i As Long

    Dim iCol As Long

    Dim iRow As Long

    Dim J As Long

    Dim K As Long

    Dim L As Long

    Dim lL As Long

    Dim Indxc() As Long

    Dim Indxr() As Long

    Dim Ipiv() As Long


    Dim Big As Double

    Dim Dum As Double

    Dim piVinv As Double



        N = gDim


        ReDim Indxc(N) As Long

        ReDim Indxr(N) As Long

        ReDim Ipiv(N) As Long


        For i = 1 To N 'do 22 i=1,n


            Big = 0


            For J = 1 To N 'do 13 j=1,n


                If Ipiv(J) <> 1 Then 'if(ipiv(j).ne.1)then


                    For K = 1 To N 'do 12 k=1,n


                        If Ipiv(K) = 0 Then 'if (ipiv(k).eq.0) then


                            If Abs(mMatrix1(J, K)) >= Big Then ' if (abs(a(j,k)).ge.big)then


                                Big = Abs(mMatrix1(J, K)) 'Big = Abs(a(J, K))

                                iRow = J

                                iCol = K


                            End If


                        End If


                    Next K '12


                End If


            Next J '13


            Ipiv(iCol) = Ipiv(iCol) + 1


            If iRow <> iCol Then 'If (iRow.ne.iCol) Then


                For L = 1 To N 'do 14 l=1,n


                    Dum = mMatrix1(iRow, L) '  Dum = a(iRow, L)

                    mMatrix1(iRow, L) = mMatrix1(iCol, L) 'a(irow,l)=a(icol,l)

                    mMatrix1(iCol, L) = Dum


                Next L '14


                'do 15 l=1,m

                'Dum = b(iRow, L)

                'b(iRow, L) = b(iCol, L)

                'b(iCol, L) = Dum

                'enddo 15


            End If


            Indxr(i) = iRow

            Indxc(i) = iCol


            If Abs(mMatrix1(iRow, iCol)) <= 0.000000000000001 Then 'if (a(icol,icol).eq.0.) pause

                Debug.Print "ERROR : TLMatrixInvertGaussJordan : Singular Matrix."

                Exit Function

            End If


            piVinv = 1 / mMatrix1(iCol, iCol) 'pivinv=1./a(icol,icol)


            mMatrix1(iCol, iCol) = 1 'a(iCol, iCol) = 1


            For L = 1 To N 'do 16 l=1,n

                mMatrix1(iCol, L) = mMatrix1(iCol, L) * piVinv 'a(icol,l)=a(icol,l)*pivinv

            Next L '16


            'do 17 l=1,m

            'b(iCol, L) = b(iCol, L) * piVinv

            'enddo 17


            For lL = 1 To N 'do 21 ll=1,n


                If lL <> iCol Then 'if(ll.ne.icol)then


                    Dum = mMatrix1(lL, iCol) 'dum=a(ll,icol)


                    mMatrix1(lL, iCol) = 0 'a(ll,icol)=0.


                    For L = 1 To N 'do 18 l=1,n

                        mMatrix1(lL, L) = mMatrix1(lL, L) - mMatrix1(iCol, L) * Dum 'a(ll,l)=a(ll,l)-a(icol,l)*dum

                    Next L ' 18


                    'do 19 l=1,m

                    'b(lL, L) = b(lL, L) - b(iCol, L) * Dum

                    'enddo 19


                End If


            Next lL '21


        Next i '22


        For L = N To 1 Step -1 'do 24 l=n,1,-1


            If Indxr(L) <> Indxc(L) Then 'if(indxr(l).ne.indxc(l))then


                For K = 1 To N 'do 23 k=1,n


                    Dum = mMatrix1(K, Indxr(L)) 'dum=a(k,indxr(l))

                    mMatrix1(K, Indxr(L)) = mMatrix1(K, Indxc(L)) 'a(k,indxr(l))=a(k,indxc(l))

                    mMatrix1(K, Indxc(L)) = Dum 'a(k,indxc(l))=dum


                Next K '23


            End If


        Next L '24


    TLMatrixInvertGaussJordan = 1

    Exit Function


    Debug.Print "ERROR : TLMatrixInvertGaussJordan : Unspecified error."

End Sub