using System;

namespace ConsoleApplication5
{
    class Program
    {
        static void Main(string[] args)
        {
            int size;

            Console.WriteLine("Gauss'es method of LU.\nEnter system dimension: ");
            size = int.Parse(Console.ReadLine());

            double[,] matrixU = new double[size, size];
            double[,] matrixL = new double[size, size];

            for (int i = 0; i < size; i++)
            {
                for (int j = 0; j < size; j++)
                {
                    matrixL[i,j] = (i == j ? 1 : 0);
                }
            }
            for (int i = 0; i < size; i++)
            {
                Console.WriteLine("Enter " + (i + 1) + " row: ");

                for (int j = 0; j < size; j++)
                {
                    matrixU[i,j] = double.Parse(Console.ReadLine());
                }
            }
            if (!getLUDecomposition(ref matrixU, ref matrixL, 0, size))
                Console.WriteLine("LU-decomposition for this matrix not found");
            else
            {
                Console.WriteLine("L-matrix is:\n");
                for (int i = 0; i < size; i++)
                {
                    for (int j = 0; j < size; j++)
                    {
                        Console.WriteLine(matrixL[i, j] +  " ");
                    }
                    Console.WriteLine();
                }
                Console.WriteLine("\nU-matrix is:\n");
                for (int i = 0; i < size; i++)
                {
                    for (int j = 0; j < size; j++)
                    {
                        Console.WriteLine(matrixU[i,j] + " ");
                    }
                    Console.WriteLine();
                } 
            }

        }

        static bool getLUDecomposition(ref double[,] matrixU, ref double[,] matrixL, int currRowAndColumn, int numberOfEquation)
        {
                if (matrixU[currRowAndColumn, currRowAndColumn] == 0
                    && currRowAndColumn < (numberOfEquation - 1))
                    return false;
            double temp;

            for (int i = currRowAndColumn + 1; i < numberOfEquation; i++)
            {
                matrixL[i, currRowAndColumn] = matrixU[i, currRowAndColumn] / matrixU[currRowAndColumn, currRowAndColumn];
                temp = -matrixL[i, currRowAndColumn];
                for (int k = currRowAndColumn; k < numberOfEquation; k++)
                {
                    matrixU[i, k] = matrixU[i, k] + matrixU[currRowAndColumn, k] * temp;
                }
            }
            return getLUDecomposition(ref matrixU, ref  matrixL, currRowAndColumn + 1, numberOfEquation);
        }
    }
}