using System;

namespace Science.Mathematics.LinearAlgebra
{
    public class FactorizationAeqCCt
    {
        private Matrix C;
        private double[,] a;
        public FactorizationAeqCCt(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];
            C = new Matrix(A.SizeOfRow, A.SizeOfColumn);
        }
        public void Compute()
        {
            Cholesky obj = new Cholesky(a);

            for (int i = 0; i < a.GetLength(0); i++)
            {
                for (int j = 0; j < a.GetLength(1); j++)
                {
                    C[i, j] = obj.LowerTriangularMatrix[i, j];
                }
            }   
        }
        public Matrix LowerTriangular
        {
            get { return C; }
        }
    }


    class Cholesky
    {
        private int n;
        private double[,] el;
        public double[,] LowerTriangularMatrix
        {
            get { return el; }
        }
        public Cholesky(double[,] a)
        {
            n = a.GetLength(0);
            el = a;
            int i, j, k;
            double sum;

            for (i = 0; i < n; i++)
            {
                for (j = i; j < n; j++)
                {
                    for (sum = el[i, j], k = i - 1; k >= 0; k--) sum -= el[i, k] * el[j, k];
                    if (i == j)
                    {
                        if (sum <= 0.0) try { throw new Science.Error(); }
                            catch (Science.Error e)
                            {
                                e.Write("Cholesky failed");
                            }
                        el[i, i] = Math.Sqrt(sum);
                    }
                    else el[j, i] = sum / el[i, i];
                }
            }
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < i; j++)
                {
                    el[j, i] = 0.0;
                }
            }
        }
        public void solve(double[] b, double[] x)
        {
            int i, k;
            double sum;

            for (i = 0; i < n; i++)
            {
                for (sum = b[i], k = i - 1; k >= 0; k--) sum -= el[i, k] * x[k];
                x[i] = sum / el[i, i];
            }
            for (i = n - 1; i >= 0; i--)
            {
                for (sum = x[i], k = i + 1; k < n; k++) sum -= el[k, i] * x[k];
                x[i] = sum / el[i, i];
            }
        }
        public void elmult(double[] y, double[] b)
        {
            int i, j;
            for (i = 0; i < n; i++)
            {
                b[i] = 0.0;
                for (j = 0; j <= i; j++) b[i] += el[i, j] * y[j];
            }
        }
        public void elsolve(double[] b, double[] y)
        {
            int i, j;
            double sum;
            for (i = 0; i < n; i++)
            {
                for (sum = b[i], j = 0; j < i; j++) sum -= el[i, j] * y[j];
                y[i] = sum / el[i, i];
            }
        }
        public void inverse(double[,] ainv)
        {
            int i, j, k;
            double sum;
            for (i = 0; i < n; i++)
            {
                for (j = 0; j <= i; j++)
                {
                    sum = (i == j ? 1.0 : 0.0);
                    for (k = i - 1; k >= j; k--) sum -= el[i, k] * ainv[j, k];
                    ainv[j, i] = sum / el[i, i];
                }
            }
            for (i = n - 1; i >= 0; i--)
            {
                for (j = 0; j <= i; j++)
                {
                    sum = (i < j ? 0.0 : ainv[j, i]);
                    for (k = i + 1; k < n; k++) sum -= el[k, i] * ainv[j, k];
                    ainv[i, j] = ainv[j, i] = sum / el[i, i];
                }
            }
        }
        public double logdet()
        {
            double sum = 0.0;
            for (int i = 0; i < n; i++) sum += Math.Log(el[i, i]);
            return 2.0 * sum;
        }
    }
}