Inverse matrix

using CSparse.Double;
using CSparse.Double.Factorization;
using CSparse.Storage;
using System;

/// <summary>
/// Compute the inverse of given matrix.
/// </summary>
public static CompressedColumnStorage<double> Inverse(CompressedColumnStorage<double> A)
{
    int size = A.RowCount;

    if (A.ColumnCount != size)
    {
        throw new ArgumentException("Matrix must be square.");
    }

    var lu = SparseLU.Create(A, ColumnOrdering.MinimumDegreeAtPlusA, 1.0);

    var b = new double[size];
    var x = new double[size];

    // Result might be dense. We use a sparse storage with memory fully
    // allocated for simplicity.
    var inverse = new CoordinateStorage<double>(size, size, size * size);

    for (int i = 0; i < size; i++)
    {
        b[i] = 1.0;

        lu.Solve(b, x);

        for (int j = 0; j < size; j++)
        {
            inverse.At(j, i, x[j]);
        }

        b[i] = 0.0;
    }

    return SparseMatrix.OfIndexed(inverse);
}