using System;

namespace Science.Mathematics.LinearAlgebra
{
    public class FactorizationPAeqLU
    {
        private Matrix U;
        private Matrix L;
        private Matrix P;
        private double d;
        private double[,] a;
        public FactorizationPAeqLU(Matrix A)
        {
            a = new double[A.SizeOfRow, A.SizeOfColumn];
            for (int i = 0; i < a.GetLength(0); i++)
                for (int j = 0; j < a.GetLength(1); j++)
                    a[i, j] = A[i, j];
            U = new Matrix(A.SizeOfRow, A.SizeOfColumn);
            L = new Matrix(A.SizeOfRow, A.SizeOfColumn);
            P = new Matrix(A.SizeOfRow, A.SizeOfColumn);
            for (int i = 0; i < P.SizeOfRow; i++) P[i, i] = 1.0;
        }
        public void Compute()
        {
            LUdcmp obj = new LUdcmp(a);

            for (int i = 0; i < a.GetLength(0); i++)
            {
                for (int j = 0; j < a.GetLength(1); j++)
                {
                    if (j > i) U[i, j] = obj.LU[i, j];
                    else if (j == i) { U[i, j] = obj.LU[i, j]; L[i, j] = 1.0; }
                    else L[i, j] = obj.LU[i, j];
                }
                if (obj.RowPermutation[i] > i)
                {
                    for (int j = 0; j < a.GetLength(1); j++)
                    {
                        double dum = P[i, j];
                        P[i, j] = P[obj.RowPermutation[i], j];
                        P[obj.RowPermutation[i], j] = dum;
                    }
                }
            }
            d = obj.Sign;
        }
        public Matrix UpperTriangular
        {
            get { return U; }
        }
        public Matrix LowerTriangular
        {
            get { return L; }
        }
        public Matrix RowPermutation
        {
            get { return P; }
        }
        public double Sign
        {
            get { return d; }
        }
    }

    class LUdcmp
    {
        int n;
        private double[,] lu;
        private int[] indx;
        private double d;
        private double[,] aref;

        public double[,] LU  // output
        {
            get { return lu; }  
        }
        public int[] RowPermutation // output
        {
            get { return indx; }
        }
        public double Sign // output
        {
            get { return d; }
        }
        public LUdcmp(double[,] a)   // a input
        {
            n = a.GetLength(0);
            lu = a;
            aref = a;
            indx = new int[n];
            double TINY = 1.0e-40;
            int i, imax = 0, j, k;
            double big, temp;
            double[] vv = new double[n];
            d = 1.0;
            for (i = 0; i < n; i++)
            {
                big = 0.0;
                for (j = 0; j < n; j++)
                    if ((temp = Math.Abs(lu[i, j])) > big) big = temp;
                if (big == 0.0) try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("Singular matrix in LUdcmp");
                    }
                vv[i] = 1.0 / big;
            }
            for (k = 0; k < n; k++)
            {
                big = 0.0;
                for (i = k; i < n; i++)
                {
                    temp = vv[i] * Math.Abs(lu[i, k]);
                    if (temp > big)
                    {
                        big = temp;
                        imax = i;
                    }
                }
                if (k != imax)
                {
                    for (j = 0; j < n; j++)
                    {
                        temp = lu[imax, j];
                        lu[imax, j] = lu[k, j];
                        lu[k, j] = temp;
                    }
                    d = -d;
                    vv[imax] = vv[k];
                }
                indx[k] = imax;
                if (lu[k, k] == 0.0) lu[k, k] = TINY;
                for (i = k + 1; i < n; i++)
                {
                    temp = lu[i, k] /= lu[k, k];
                    for (j = k + 1; j < n; j++) lu[i, j] -= temp * lu[k, j];
                }
            }
        }

        public void solve(double[] b, double[] x) // b input and x output
        {
            int i, ii = 0, ip, j;
            double sum;
            if (b.Length != n || x.Length != n) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Bad sizes");
                }
            for (i = 0; i < n; i++) x[i] = b[i];
            for (i = 0; i < n; i++)
            {
                ip = indx[i];
                sum = b[ip];
                x[ip] = x[i];
                if (ii != 0) for (j = ii - 1; j < i; j++) sum -= lu[i, j] * x[j];
                else if (sum != 0.0) ii = i + 1;
                x[i] = sum;
            }
            for (i = n - 1; i >= 0; i--)
            {
                sum = x[i];
                for (j = i + 1; j < n; j++) sum -= lu[i, j] * x[j];
                x[i] = sum / lu[i, i];
            }
        }

        public void solve(double[,] b, double[,] x) // b input and x output
        {
            int i, j, m = b.GetLength(1);
            if (b.GetLength(0) != n || x.GetLength(0) != n
                || b.GetLength(1) != x.GetLength(1)) try { throw new Science.Error(); }
                catch (Science.Error e)
                {
                    e.Write("Bad sizes");
                }
            double[] xx = new double[n];
            for (j = 0; j < m; j++)
            {
                for (i = 0; i < n; i++) xx[i] = b[i, j];
                solve(xx, xx);
                for (i = 0; i < n; i++) x[i, j] = xx[i];
            }
        }

        public double[,] inverse()
        {
            double[,] ainv = new double[n, n];
            for (int i = 0; i < n; i++) ainv[i, i] = 1.0;
            solve(ainv, ainv);
            return ainv;
        }

        public double det()
        {
            double dd = d;
            for (int i = 0; i < n; i++) dd *= lu[i, i];
            return dd;
        }

        public void mprove(double[] b, double[] x)
        {
            int j, i;
            double sdp;
            double[] r = new Double[n];

            for (i = 0; i < n; i++)
            {
                sdp = -b[i];
                for (j = 0; j < n; j++) sdp += aref[i, j] * x[j];
                r[i] = sdp;
            }
            solve(r, r);
            for (i = 0; i < n; i++) x[i] -= r[i];
        }
    }
}