Matrix inversion in C#

I was talking to someone about my demos the other day, and he mentioned that he’d tried to find a simple matrix inversion routine in C# on the web, and not found one he could easily understand. So, for what it’s worth, I thought I’d post the one I use. It uses the System.Numerics namespace for the Complex structure for complex numbers, and a straight-forward Gaussian elimination technique with pivoting. It returns “true” if the inversion was successful, and “false” if the matrix is singular. Somewhat controversially, it does carry on with the inversion when the matrix is singular, but in these cases the result will normally be rubbish.

The theory behind this technique is described in the book chapter on Gaussian Elimination, now available in the Linear Algebra section of the book.

Just in case anyone finds it interesting / useful, the code is here:

// The complex matrix inversion routine:
internal Boolean MyInvertMatrix(Complex[,] A, out Complex[,] invA)
{
    // There are faster ways to do this, but for simplicity
    // and to get something working quickly, I'll just write a 
    // simple Gaussian-elimination matrix inverter here, and look 
    // at speeding things up later.

    // Keep a record to see if there is a sensible inverse:
    Boolean bNotIllConditioned = true;

    // If the matrix is not square, there is no inverse:
    if (A.GetLength(0) != A.GetLength(1))
    {
        invA = A; // Have to assign something to invA before returning
        return false;
    }

    // This routine destroys the matrix it's working on, so I'll first 
    // make a copy of it, as well as setting up the output matrix invA 
    // as a unit matrix of the appropriate size:
    int dimension = A.GetLength(0);
    Complex[,] working = (Complex[,])A.Clone();
    Complex[,] inverse = new Complex[dimension, dimension];
    // C# will set the initial values to zero, so to create a unit
    // matrix, I just need to fill in the diagonal elements:
    for (int loop = 0; loop < dimension; loop++) inverse[loop, loop] = 1.0;

    // OK, first convert working to upper triangular form:
    for (int loop = 0; loop < dimension; loop++) // for each row
    {
        int currentRow = loop;

        // First step is pivoting: make sure the biggest element
        // remaining in any column is on the next row.  First, find
        // the biggest element remaining in the current column:
        double biggestSoFar = 0.0; int biggestRow = currentRow;
        for (int x = currentRow; x < dimension; x++)
        {
            double sizeOfThis = working[x, currentRow].Magnitude;
            if (sizeOfThis > biggestSoFar)
            {
                biggestSoFar = sizeOfThis;
                biggestRow = x;
            }
        }

        // and if this is not at the top, swop the rows of working
        // and inverse around until it is:
        if (biggestRow != currentRow)
        {
            Complex temp;
            for (int lop = currentRow; lop < dimension; lop++)
            {
                temp = working[currentRow, lop];
                working[currentRow, lop] = working[biggestRow, lop];
                working[biggestRow, lop] = temp;
            }
            for (int lop = 0; lop < dimension; lop++)
            {
                temp = inverse[currentRow, lop];
                inverse[currentRow, lop] = inverse[biggestRow, lop];
                inverse[biggestRow, lop] = temp;
            }
        }

        // Then, go down the matrix subtracting as necessary
        // to get rid of the lower-triangular elements:
        for (int lop = currentRow + 1; lop < dimension; lop++)
        {
            // Matrix might be ill-conditioned.  I should check:
            if (working[currentRow, currentRow] == 0)
            {
                bNotIllConditioned = false;
                working[currentRow, currentRow] = Globals.TINYDOUBLE;
            }
            Complex factor = working[lop, currentRow] / working[currentRow, currentRow];

            // If the matrix is fairly sparse (quite common for this
            // application), it might make sense to check that the 
            // lower elements are not already zero before doing all
            // the scaling and replacing:
            if (factor != 0.0)
            {
                // Only have to do from current row on in working, but due
                // to pivoting, might have to do the entire row in inverse:
                for (int lp = currentRow; lp < dimension; lp++)
                    working[lop, lp] -= factor * working[currentRow, lp];
                for (int lp = 0; lp < dimension; lp++)
                    inverse[lop, lp] -= factor * inverse[currentRow, lp];
            }
        }
        // That's it for this row, now on to the next one...
    }

    // Now with the working matrix in upper-triangular form, continue the same
    // process amongst the upper-triangular elements to convert working into
    // diagonal form:
    for (int loop = dimension - 1; loop >= 0; loop--) // for each row
    {
        int currentRow = loop;

        // Matrix might be ill-conditioned.  I should check:
        if (working[currentRow, currentRow] == 0)
        {
            bNotIllConditioned = false;
            working[currentRow, currentRow] = Globals.TINYDOUBLE;
        }

        // Then, go up the matrix subtracting as necessary to get 
        // rid of the remaining upper-triangular elements:
        for (int lop = currentRow - 1; lop >= 0; lop--)
        {
            Complex factor = working[lop, currentRow] / working[currentRow, currentRow];

            // There's only one element in working to change (the other elements
            // in the row of working are all zero), and that will always be set
            // to zero; but you might have to do the entire row in inverse:
            working[lop, currentRow] = 0.0;

            if (factor != 0.0)
            {
                for (int lp = 0; lp < dimension; lp++)
                {
                    inverse[lop, lp] -= factor * inverse[currentRow, lp];
                }
            }
        }
        // That's it for this row, now on to the next one...
    }

    // Should now have working as a diagonal matrix.  Final thing is 
    // to scale all the rows:
    for (int loop = 0; loop < dimension; loop++)
    {
        Complex scale = working[loop, loop];
        for (int lop = 0; lop < dimension; lop++) inverse[loop, lop] /= scale;
    }

    // That's it.  inverse should now be the inverse of the original matrix.
    invA = inverse;
    return bNotIllConditioned;
}