Ускорение/оптимизация метода. Решение СЛАУ методом Гаусса

/// <summary>
        /// Решение слау методом прямого и обратного хода Гаусса
        /// </summary>
        /// <param name="matrix">augment(Глобальная матрица жёсткости, Глобальный грузовой вектор)</param>
        /// <param name="n">кол-во неизвестных</param>
        /// <param name="x">Вектор неизвестных</param>
        /// <param name="solutions">Пустая строка, покажет результат решения(0,1 или бесконечноть)</param>
        /// <returns></returns>
        public static int SolveSLAE5(double[,] matrix, int n, out double[] x, out string solutions)
        {
            //Объявляем переменные:
            int i;//Переменная, номер строки.
            int j;//Переменная, номер столбца.
            int m;//Размерность системы (учитываем вектор-расширение матрицы).
 
            m = n + 1;
 
            int k;
            int nl;//Переменная номер строки с ведущим элементом.
            double o = 0;//Перемнная для поиска самого большого по модулю элемента в матрице. 
            double p = 0;//Перемнная для поиска самого большого по модулю элемента в матрице.
            double tmp = 0;//Временная перемнная для прямого хода метода Гаусса.
 
            for (k = 0; k < n; k++)
            {
                //Поиск ведущего элемента:
                p = 0;
                nl = k;
                for (i = k; i < n; i++)
                {
                    o = Math.Abs(matrix[i,k]);
                    if (Math.Abs(o) > Math.Abs(p))
                    {
                        p = Math.Abs(o);//Нашли значение самого большого элемента.
                        nl = i;//Нашли номер строки с самым большим по модулю элементом (ведущим элементом).
                    }
                }
 
                //Перестановка строк:
                if (nl != k)
                {
                    for (i = k; i < n + 1; i++)
                    {
                        double tmp_line = 0;//Временная переменная для перестановки строк местами.
                        tmp_line = matrix[nl,i];
                        matrix[nl,i] = matrix[k,i];
                        matrix[k,i] = tmp_line;
                    }
                }
 
                //Прямой ход, приведение к верхнетреугольному виду:
 
                tmp = matrix[k,k];//Первый элемент в текущей матрице.
 
                if (tmp != 0)
                {
                    for (j = k; j < n + 1; j++)//Нормализация.
                    {
                        matrix[k,j] = matrix[k,j] / tmp;
                    }
                    for (i = k + 1; i < n; i++)
                    {
                        tmp = matrix[i,k];//Главный элемент.
                        for (j = k; j < n + 1; j++)//Из следующей строки вычитаем первую, умноженную на главный элемент.
                        {
                            matrix[i,j] = matrix[i,j] - tmp * matrix[k,j];
                        }
                    }
                }
            }
            //Обратный ход:
            for (k = n - 1; k >= 0; k--)
            {
                for (i = 0; i < k; i++)
                {
                    tmp = matrix[i,k];//Множитель.
                    for (j = k; j < n + 1; j++)//Из следующей строки вычитаем первую, умноженную на главный элемент.
                    {
                        matrix[i,j] = matrix[i,j] - tmp * matrix[k,j];
                    }
                }
            }
 
            bool nenulinline = false;//Переменная, которая обозначает, что в строке есть ненулевой элемент.
            bool netresh = false;//Переменная, которая обозначает, что система не имеет решений.
            bool beskresh = false;//Переменная, которая обозначает, что система имеет бесконечное множество решений.
 
            for (k = 0; k < n; k++)
            {
                if (matrix[k,k] == 0)//Если на диагонали 0, значит проверяем строку.
                {
                    nenulinline = false;//Обнуляем переменную.
                    for (j = k + 1; j < n; j++)
                    {
                        if (matrix[k,j] != 0)//Если в строке есть ненулевой элемент значит обозначим это через флаг.
                        {
                            nenulinline = true;
                        }
                    }
                    if (!nenulinline)//Если в строке все нули, то проверям свободный член.
                    {
                        if (matrix[k,n] != 0)//Если не равно нулю, то решений нет.
                        {
                            netresh = true;
                        }
                        else//А иначе данное уравнение решается всегда.
                        {
                            beskresh = true;
                        }
                    }
                }
            }
            x = new double[n];
            solutions = "one sulution";
            //Выводим решения
            if (netresh == true)
            {
                solutions = "Null solution";
            }
            else
            {
                if (beskresh == true)
                {
                    solutions = "Infinit solution.";
                }
                else
                {
                    for (i = 0; i < n; i++)
                    {
                        x[i] = matrix[i,n];
                    }
                }
            }
            return 0;
        }

Parallel.For(0, Length, i =>
            {
                X[i] = new double[Length]; //Вектор неизвестных
                MyMath.SolveSLAE(Array, Array.GetLength(0), out X[i], out sol); //Решение СЛАУ 
                                                                           //Array*X = последний столбец Array
                                                                           //Array - массив размерностью (n x n+1)
            });