using System;

namespace Science.Mathematics.LinearAlgebra
{
	/// <summary>
	/// HermitianMatrix
	/// </summary>
	public class HermitianMatrix : MatrixComplex
	{
		public HermitianMatrix(int nr, int nc) : base(nr,nc)
		{
            Check();
		} 
		public HermitianMatrix(Complex[,] x) : base(x)
	    {
            Check();
        }
        private void Check()
        {
			for(int i = 0; i < this.SizeOfRow; i++)
				for(int j = i; j < this.SizeOfColumn; j++)
				{
					if ( this.cel[j,i].Real != this.cel[i,j].Real)
					{
                        try { throw new Science.Error(); }
                        catch (Science.Error e)
                        {
                            e.Write("Something Wrong");
                        }
					}
		            if (this.cel[j,i].Imaginary != - this.cel[i,j].Imaginary)
					{
                        try { throw new Science.Error(); }
                        catch (Science.Error e)
                        {
                            e.Write("Something Wrong");
                        } 
					}
				}
		}

        private double[] e;
        private VectorComplex[] v;
        public void Diagonalize()
        {
            SymmetricMatrix a = new SymmetricMatrix(2 * sr, 2 * sc);
            for (int i = 0; i < this.SizeOfRow; i++)
                for (int j = 0; j < this.SizeOfColumn; j++)
                {
                    a[i, j] = this.cel[i, j].Real;
                    a[this.SizeOfRow + i, j] = this.cel[i, j].Imaginary;
                    a[i, this.SizeOfColumn + j] = -this.cel[i, j].Imaginary;
                    a[this.SizeOfRow + i, this.SizeOfColumn + j] = this.cel[i, j].Real;
                }

            FactorizationAeqQLQt obj = new FactorizationAeqQLQt(a);
            obj.Compute();
            e = new double[sr];
            v = new VectorComplex[sr];
            for (int i = 0; i < sr; i++)
            {
                e[i] = obj.Eigenvalues[2*i, 2*i];
                v[i] = new VectorComplex(sr);
                for (int j = 0; j < sr; j++)
                {
                    v[i][j].Real = obj.Eigenvectors[j, 2*i];
                    v[i][j].Imaginary = obj.Eigenvectors[j + sr, 2*i];
                }
            }
        }
        public double[] Eigenvalue
        {
            get { return e; }
        }
        public VectorComplex[] EigenvectorComplex
        {
            get { return v; }
        }
	}
}