using System;

namespace Science.Mathematics.LinearAlgebra
{
    public class FactorizationEAeqR
    {
        private Matrix R;
        private Matrix E;
        private double[,] a;
        private double[,] b;
        private double[,] ao;
        public FactorizationEAeqR(Matrix A)
        {
            E = new Matrix(A.SizeOfRow, A.SizeOfRow);
            R = new Matrix(A.SizeOfRow, A.SizeOfColumn);
            a = new double[A.SizeOfRow, A.SizeOfColumn];
            ao = new double[A.SizeOfRow, A.SizeOfColumn];
            b = new double[A.SizeOfRow, A.SizeOfRow];
            for (int i = 0; i < a.GetLength(0); i++)
                for (int j = 0; j < a.GetLength(1); j++)
                {
                    a[i, j] = A[i, j];
                    ao[i, j] = a[i, j];
                }
            for (int i = 0; i < b.GetLength(0); i++) b[i, i] = 1.0;
        }
        public void Compute()
        {
            Gaussj obj = new Gaussj();
            obj.gaussj(a, b);
            for (int i = 0; i < b.GetLength(0); i++)
                for (int j = 0; j < b.GetLength(1); j++)
                    E[i, j] = b[i, j];
            for (int i = 0; i < a.GetLength(0); i++)
            {
                double sum;
                for (int j = 0; j < a.GetLength(1); j++)
                {
                    sum = 0.0;
                    for (int k = 0; k < b.GetLength(1); k++) sum += b[i, k] * ao[k, j];
                    R[i, j] = sum;
                }
            }
        }
        public Matrix ReducedRowEchelonForm
        {
            get { return R; }
        }
        public Matrix Elimination
        {
            get { return E; }
        }
    }

    class Gaussj
    {
        public Gaussj()
        {
        }
        private void SWAP(ref double a, ref double b)
        {
            double temp = a;
            a = b;
            b = temp;
        }
        public void gaussj(double[,] a, double[,] b)  // a (m by n) and b (m by m) are input, also output
        {
            int m = a.GetLength(0);
            int n = a.GetLength(1);
            int i, j, l, ll;
            int irow = 0;
            double big, dum, pivinv;
            int np = 0;
            for (i = 0; i < n; i++)
            {
                big = 0.0;
                for (j = np; j < m; j++)
                    if (Math.Abs(a[j, i]) >= big)
                    {
                        big = Math.Abs(a[j, i]);
                        irow = j;                               
                    }       
     
                if (irow != np)
                {
                    for (l = 0; l < n; l++)
                    {
                        SWAP(ref a[irow, l], ref a[np, l]);
                    }
                    for (l = 0; l < m; l++)
                    {
                        SWAP(ref b[irow, l], ref b[np, l]);
                    }
                }
                if (Math.Abs(a[np, i]) > 1.0E-13)
                {
                    pivinv = 1.0 / a[np, i];
                    a[np, i] = 1.0;
                    for (l = 0; l < n; l++) a[np, l] *= pivinv;
                    for (l = 0; l < m; l++) b[np, l] *= pivinv;
                    for (ll = 0; ll < m; ll++)
                        if (ll != np)
                        {
                            dum = a[ll, i];
                            a[ll, i] = 0.0;
                            for (l = 0; l < n; l++) a[ll, l] -= a[np, l] * dum;
                            for (l = 0; l < m; l++) b[ll, l] -= b[np, l] * dum;
                        }
                }
                else
                {
                    --np;
                }
                ++np;
                if (np >= m) break;
            }
        }
    }      
}