using System;

namespace Science.Mathematics.LinearAlgebra
{
    public class MatrixComplex
    {
        protected int sr, sc;
		public int SizeOfRow
		{		    
			get{return sr;}
		}
		public int SizeOfColumn
		{		    
			get{return sc;}
		}
		protected Complex[,] cel;
        public MatrixComplex(int numberOfRows, int numberOfColumns)
        {
            sr = numberOfRows;
            sc = numberOfColumns;
            cel = new Complex[sr, sc];
            for (int k = 0; k < sr; k++)
                for (int kk = 0; kk < sc; kk++)
                    cel[k, kk] = new Complex(0.0, 0.0);
        }
		public MatrixComplex(Complex[,] x)
		{
			sr = x.GetLength(0);
			sc = x.GetLength(1);
			cel = new Complex[sr,sc];
			for(int k = 0; k < sr; k++)
				for(int kk = 0; kk < sc; kk++)
					cel[k,kk] = new Complex(x[k,kk].Real,x[k,kk].Imaginary);
	    }
		public Complex this[int whereRow, int whereColumn]
		{
			get{return cel[whereRow, whereColumn];}
			set{cel[whereRow, whereColumn]=value;}
		}	

        public MatrixComplex Transpose
		{
            get
            {
                MatrixComplex m = new MatrixComplex(sc, sr);
                for (int k = 0; k < sr; k++)
                    for (int kk = 0; kk < sc; kk++)
                    {
                        m[kk, k].Real = cel[k, kk].Real;
                        m[kk, k].Imaginary = cel[k, kk].Imaginary;
                    }
                return m;
            }
		}

        public MatrixComplex ConjugateTranspose
        {
            get
            {
                MatrixComplex m = new MatrixComplex(sc, sr);
                for (int k = 0; k < sr; k++)
                    for (int kk = 0; kk < sc; kk++)
                    {
                        m[kk, k].Real = cel[k, kk].Real;
                        m[kk, k].Imaginary = - cel[k, kk].Imaginary;
                    }
                return m;
            }
        }

        public MatrixComplex Inverse
        {
            get
            {
                if (sr != sc) try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("Square matrix should be considered.");
                    }
                double[,] ell = new double[2 * sr, 2 * sc];
                for (int i = 0; i < sr; i++)
                    for (int j = 0; j < sc; j++)
                    {
                        ell[i, j] = cel[i, j].Real;
                        ell[sr + i, j] = cel[i, j].Imaginary;
                        ell[i, sc + j] = -cel[i, j].Imaginary;
                        ell[sr + i, sc + j] = cel[i, j].Real;
                    }
                LUdcmp obj = new LUdcmp(ell);
                double[,] a = obj.inverse();
                MatrixComplex im = new MatrixComplex(sr, sc);
                for (int j = 0; j < sr; j++)
                    for (int i = 0; i < sc; i++) im[j, i] = new Complex(a[j, i],a[j+sr, i]);
                return im;
            }
        }

        public static MatrixComplex operator *(double a, MatrixComplex c)
        {
            MatrixComplex b = new MatrixComplex(c.SizeOfRow, c.SizeOfColumn);
            for (int k = 0; k < b.SizeOfRow; k++)
                for (int j = 0; j < b.SizeOfColumn; j++)
                {
                    b[k, j].Real = a * c[k, j].Real;
                    b[k, j].Imaginary = a * c[k, j].Imaginary;
                }
            return b;
        }
        public static MatrixComplex operator *(MatrixComplex c, double a)
		{
			return a*c;
		}
        public static MatrixComplex operator *(Complex a, MatrixComplex c)
        {
            MatrixComplex b = new MatrixComplex(c.SizeOfRow, c.SizeOfColumn);
            for (int k = 0; k < b.SizeOfRow; k++)
                for (int j = 0; j < b.SizeOfColumn; j++)
                {
                    Complex t = a * c[k, j];
                    b[k, j].Real = t.Real;
                    b[k, j].Imaginary = t.Imaginary;
                }
            return b;
        }
        public static MatrixComplex operator *(MatrixComplex c, Complex a)
		{
			return a*c;
		}
        public static MatrixComplex operator +(MatrixComplex a, MatrixComplex c)
        {
            MatrixComplex b = new MatrixComplex(c.SizeOfRow, c.SizeOfColumn);
            for (int k = 0; k < b.SizeOfRow; k++)
                for (int j = 0; j < b.SizeOfColumn; j++)
                {
                    Complex t = a[k, j] + c[k, j];
                    b[k, j].Real = t.Real;
                    b[k, j].Imaginary = t.Imaginary;
                }
            return b;
        }
        public static MatrixComplex operator -(MatrixComplex a, MatrixComplex c)
		{
            MatrixComplex b = new MatrixComplex(c.SizeOfRow, c.SizeOfColumn);
            for (int k = 0; k < b.SizeOfRow; k++)
                for (int j = 0; j < b.SizeOfColumn; j++)
                {
                    Complex t = a[k, j] - c[k, j];
                    b[k, j].Real = t.Real;
                    b[k, j].Imaginary = t.Imaginary;
                }
            return b;
		}

        public static MatrixComplex operator -(MatrixComplex a)
        {
            MatrixComplex b = new MatrixComplex(a.SizeOfRow, a.SizeOfColumn);
            for (int k = 0; k < a.SizeOfRow; k++)
                for (int j = 0; j < a.SizeOfColumn; j++)
                {
                    b[k, j].Real =  - a[k, j].Real;
                    b[k, j].Imaginary = - a[k, j].Imaginary;
                }
            return b;
        }

        public static VectorComplex operator *(MatrixComplex a, VectorComplex c)
        {
            VectorComplex b = new VectorComplex(a.SizeOfRow);
            for (int k = 0; k < b.Size; k++)
            {
                Complex sum = new Complex(0.0, 0.0);
                for (int s = 0; s < c.Size; s++)
                    sum = sum + a[k, s] * c[s];
                b[k].Real = sum.Real;
                b[k].Imaginary = sum.Imaginary;
            }
            return b;
        }
        public static VectorComplex operator *(VectorComplex c, MatrixComplex a)
		{
            VectorComplex b = new VectorComplex(a.SizeOfRow);
            for (int k = 0; k < b.Size; k++)
            {
                Complex sum = new Complex(0.0, 0.0);
                for (int s = 0; s < c.Size; s++)
                    sum = sum + c[s] * a[s, k] ;
                b[k].Real = sum.Real;
                b[k].Imaginary = sum.Imaginary;
            }
            return b;
		}
        public static MatrixComplex operator *(MatrixComplex a, MatrixComplex c)
        {
            MatrixComplex b = new MatrixComplex(a.SizeOfRow, c.SizeOfColumn);
            for (int k = 0; k < b.SizeOfRow; k++)
                for (int j = 0; j < b.SizeOfColumn; j++)
                {
                    Complex sum = new Complex(0.0, 0.0);
                    for (int s = 0; s < a.SizeOfColumn; s++)
                    {
                        sum = sum + a[k, s] * c[s, j];
                    }
                    b[k, j].Real = sum.Real;
                    b[k, j].Imaginary = sum.Imaginary;
                }
            return b;
        }

        public void Chop()
        {
            for (int k = 0; k < sr; k++)
                for (int kk = 0; kk < sc; kk++)
                {
                    if (cel[k, kk].Real < 0.00000001 && cel[k, kk].Real > -0.00000001)
                        cel[k, kk].Real = 0.0;
                    if (cel[k, kk].Imaginary < 0.00000001 && cel[k, kk].Imaginary > -0.00000001)
                        cel[k, kk].Imaginary = 0.0;
                }
        }

        public override string ToString()
        {
            string res = "";
            for (int i = 0; i < this.SizeOfRow; i++)
            {
                for (int j = 0; j < this.SizeOfColumn; j++)
                    res += this[i, j].ToString() + "     ";
                res += "\r\n";
            }
            res += "\r\n";
            return res;
        } 
	}
}