//Факторизация матрицы или LU-разложение

        void factorization(int n)

        {

            void Start () 

            decimal macheps = 1 / 1000000;

            int tmp;

            decimal sum = 0;

            for (int k = 0; k < n; k++)

            {

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

                {

                    sum = 0;

                    for (int p = 0; p < k; p++)

                    {

                        sum += matrix[k, q[p]] * matrix[p, q[j]];

                    }

                    matrix[k, q[j]] = matrix[k, q[j]] - sum;

                }

                decimal max = matrix[k, q[k]];

                int max_j = k;

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

                {

                    if (Math.Abs(matrix[k, q[j]]) > Math.Abs(max))

                    {

                        max = matrix[k, q[j]];

                        max_j = j;

                    }

                    if (max_j != k)

                    {

                        tmp = q[max_j];

                        q[max_j] = q[k];

                        q[k] = tmp;

                        sign = -sign;

                    }

                }

                if (matrix[k, q[k]] == 0)

                {

                    MessageBox.Show("Матрица вырожденная!");

                    break;

                }

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

                {

                    sum = 0;

                    for (int p = 0; p < k; p++)

                    {

                        sum += matrix[i, q[p]] * matrix[p, q[k]];

                    }

                    matrix[i, q[k]] = (matrix[i, q[k]] - sum) / matrix[k, q[k]];

                }

            }