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];
}
}