Обратная матрица

/// <summary>
        /// Нахождение единичной матрицы размера n x n
        /// </summary>
        /// <param name="n">Размер матрицы</param>
        /// <returns>Единичная матрица данного размера</returns>
        public Matrix IdentityMatrix(int n)
        {
            Matrix result = new Matrix(n, n);
            for (int i = 0; i < n; i++)
                for (int j = 0; j < n; j++)
                    result[i, j] = i == j ? 1 : 0; // если диагональный элемент, то 1, иначе 0
            return result;
        }
 
        /// <summary>
        /// Перестановка строк по заданным индексам
        /// </summary>
        /// <param name="rowIndex1">Первый индекс сткроки</param>
        /// <param name="rowIndex2">Второй индекс строки</param>
        public void SwapRows(int rowIndex1, int rowIndex2)
        {
            for (int i = 0; i < Columns; i++) // перебор всех столбцов и замена значений в них
            {
                double temp = matrix[rowIndex1, i];
                matrix[rowIndex1, i] = matrix[rowIndex2, i];
                matrix[rowIndex2, i] = temp;
            }
        }
 
        /// <summary>
        /// Приведение матрицы к LUP виду
        /// </summary>
        /// <param name="result">Результирующая матрица</param>
        /// <param name="p">Матрица перестановок</param>
        public void LUP(ref Matrix result, ref Matrix p)
        {
            if (result == null) throw new ArgumentNullException("result");
            if (p == null) throw new ArgumentNullException("p");
            int n = Rows; // n - размерность исходной матрицы
 
            result = this;
 
            //загружаем в матрицу p единичную матрицу
            p = IdentityMatrix(n);
 
            for (int i = 0; i < n; i++)
            {
                // поиск опорного элемента (максимального по модулю)
                double pivotValue = 0; // опорный элемент
                int pivot = -1;         // индекс опорного элемента
                for (int row = i; row < n; row++)
                {
                    if (Math.Abs(result[row, i]) > pivotValue)
                    {
                        pivotValue = Math.Abs(result[row, i]);
                        pivot = row;
                    }
                }
                if (pivotValue == 0)
                    throw new Exception("Матрица вырождена. " + "Обратная ей матрица не может быть найдена. "
                        + "Попробуйте ввести другие значения.");
 
                // меняем местами i-ю строку и строку с опорным элементом
                p.SwapRows(pivot, i);
                result.SwapRows(pivot, i);
                for (int j = i + 1; j < n; j++)
                {
                    result[j, i] /= result[i, i];
                    for (int k = i + 1; k < n; k++)
                        result[j, k] -= result[j, i] * result[i, k];
                }
            }
            // теперь матрица result = L + U - E
        }
 
        /// <summary>
        /// Нахождение обратной матрицы
        /// </summary>
        /// <returns>Обратная матрица</returns>
        public Matrix Inverse()
        {
            if (Rows != Columns) throw new Exception("Матрица должна быть квадратной!");
            // n - размерность квадратной матрицы Matrix1
            int n = Rows;//  
            // при инициализации задается размерность nxn
            Matrix result = new Matrix(n, n);
            Matrix p = new Matrix(n, n); // матрица перестановок
            Matrix c = new Matrix(n, n); // матрица для LUP разложения
            // предполагается что в результате следующего вызова матрица C = L + U - E
            LUP(ref c, ref p);
            for (int k = n - 1; k >= 0; k--)
            {
                result[k, k] = 1;
                for (int j = n - 1; j > k; j--)
                    result[k, k] -= c[k, j] * result[j, k];
                result[k, k] /= c[k, k];
                for (int i = k - 1; i >= 0; i--)
                {
                    for (int j = n - 1; j > i; j--)
                    {
                        result[i, k] -= c[i, j] * result[j, k];
                        result[k, i] -= c[j, i] * result[k, j];
                    }
                    result[i, k] /= c[i, i];
                }
            }
            return result * p; // обратная матрица
        }

private void FillResult(Matrix matrix)
        {
            dataResult.ColumnCount = matrix.Columns;
            dataResult.RowCount = matrix.Rows;
 
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    dataResult[j, i].Value = matrix[i, j];
                }
            }
        }