LU Decomposition and Linear Equation Solving in c#

LU decomposition is used to write a matrix in Upper and Lower triangle forms which is used to solve linear equations.Here is my code for LU decomposition and Linear Equation Solving, it is based on some very old fortran code, i think lapack may have the original sources but they are hard to understand.I added a few comments but to really understand what is going on i can recommend using  pen and paper and doing with hand.

static void Main(string[] args)

{

    double[,] lu = new double[3, 3] { { 2, 1, -1 }, { -3, -1, 2 }, { -2, 1, 2 } };

    double[] b = new double[] { 8, -11, -3 };

    double[] x = new double[3];

    int[] indx = new int[b.Length];

    Program p = new Program();

    p.LUDecompose(ref lu, ref indx);

    p.Solve(ref b, ref x, indx, lu);

}

private void LUDecompose(ref double[,] lu, ref int[] indx)

{

    int i, imax = 0, j, k, n = lu.GetLength(0);

    double big, temp;

    double[] vv = new double[n];

    //preChecks

    if (lu.GetLength(0) != lu.GetLength(1) || lu.GetLength(0) != indx.Length)

        throw new Exception("matrix dimension problem only use square matrices");

    //for each row find the absolute value of the greatest cell and store in vv

    for (i = 0; i < n; i++)

    {

        big = 0.0;

        for (j = 0; j < n; j++)

            if ((temp = Math.Abs(lu[i, j])) > big) big = temp;

        if (big == 0.0)

        throw new Exception("singular matrix");

        vv[i] = 1.0 / big;//calculate scaling and save

    }

    //k is for colums start with the left look for the columns under the diagonal for the biggest value want to move the largest over diagonal

    for (k = 0; k < n; k++)//find the largest pivot element 

    {

        big = 0.0;

        for (i = k; i < n; i++)

        {

            temp = vv[i] * Math.Abs(lu[i, k]);

            if (temp > big)

            {

                big = temp;

                imax = i;

            }

        }

        if (k != imax)//do we need a row change 

        {

            for (j = 0; j < n; j++)// counter for the colums

            {

                temp = lu[imax, j];// change the rows

                lu[imax, j] = lu[k, j];

                lu[k, j] = temp;

            }

            vv[imax] = vv[k];

        }

        indx[k] = imax;

        for (i = k + 1; i < n; i++)

        {

            temp = lu[i, k] /= lu[k, k];//divide pilot element

            for (j = k + 1; j < n; j++)

                lu[i, j] -= temp * lu[k, j];

        }

    }

}

private void Solve(ref double[] b, ref double[] x, int[] indx, double[,] lu)

{

    if (b.Length != lu.GetLength(0) || x.Length != lu.GetLength(0))

        throw new Exception("vector dimension problem");

    int n = lu.GetLength(0);

    int i, ii = 0, ip, j;

    double sum = 0;

    for (i = 0; i < n; i++) x[i] = b[i];

    for (i = 0; i < n; i++)

    {

        ip = indx[i];

        sum = x[ip];

        x[ip] = x[i];

        if (ii != 0)

            for (j = ii - 1; j < i; j++) sum -= lu[i, j] * x[j];

        else if (sum != 0.0)

            ii = i + 1;

        x[i] = sum;

    }

    for (i = n - 1; i >= 0; i--)

    {

        sum = x[i];

        for (j = i + 1; j < n; j++) sum -= lu[i, j] * x[j];

        x[i] = sum / lu[i, i];

    }

}

The code also has an example how to use b is the vector right hand side of the equation x is the solution vector, lu is the matrix.