using System;

namespace Science.Mathematics.LinearAlgebra
{
	/// <summary>
	/// Matrix
	/// </summary>
	public class Matrix
	{
        // Norm, Frobenius Norm
		protected int sr, sc;
		public int SizeOfRow
		{		    
			get{return sr;}
		}
		public int SizeOfColumn
		{		    
			get{return sc;}
		}

        public Matrix()
        {
        }
     
        protected double[,] el;
    	public Matrix(int numberOfRows, int numberOfColumns)
        {
            sr = numberOfRows;
            sc = numberOfColumns;
            el = new double[numberOfRows, numberOfColumns];
		}
		public Matrix(double[,] x)
		{
			sr = x.GetLength(0);
			sc = x.GetLength(1);
			el = new double[sr,sc];
			for(int k = 0; k < sr; k++)
				for(int kk = 0; kk < sc; kk++)
					el[k,kk] = x[k,kk];
		}
        public Matrix(Vector[] columnVector)
        {
            sr = columnVector[0].Size;
            sc = columnVector.Length;
            el = new double[sr, sc];
            for (int k = 0; k < sr; k++)
                for (int kk = 0; kk < sc; kk++)
                    el[k, kk] = columnVector[kk][k];
        }
		public double this[int whereRow, int whereColumn]
		{
			get{return el[whereRow, whereColumn];}
			set{el[whereRow, whereColumn]=value;}
		}

        public int Rank
        {
            get
            {
                FactorizationEAeqR obj = new FactorizationEAeqR(this);
                obj.Compute();
                int nofpivot = 0;
                for (int i = 0; i < this.SizeOfRow; i++)
                    for (int j = 0; j < this.SizeOfColumn; j++)
                        if (obj.ReducedRowEchelonForm[i, j] > 1.0E-12)
                        {
                            nofpivot++;
                            break;
                        }
                return nofpivot;
            }
        }

		public Matrix Transpose
		{
			get
			{
				Matrix m = new Matrix(sc,sr);
				for(int k = 0; k < sr; k++)
					for(int kk = 0; kk < sc; kk++)
						m[kk,k] = el[k,kk];
				return m;
			}
		}

        public double Determinant
        {
            get
            {
                if (sr != sc) try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("Square matrix should be considered.");
                    }
                double[,] ell = new double[sr, sc];
                for (int j = 0; j < sr; j++)
                    for (int i = 0; i < sc; i++) ell[j, i] = el[j, i];
                LUdcmp obj = new LUdcmp(ell);
                return obj.det();
            }
        }

        public double Trace
        {
            get
            {
                if (sr != sc) try { throw new Science.Error(); }
                    catch (Science.Error e)
                    {
                        e.Write("Square matrix should be considered.");
                    } 
                double sum = 0.0;
                for (int k = 0; k < this.SizeOfRow; k++)
                    sum += this.el[k, k];
                return sum;

            }
        }

        public Matrix 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[sr, sc];
                for (int j = 0; j < sr; j++)
                    for (int i = 0; i < sc; i++) ell[j, i] = el[j, i];
                LUdcmp obj = new LUdcmp(ell);
                double[,] a = obj.inverse();
                Matrix im = new Matrix(sr, sc);
                for (int j = 0; j < sr; j++)
                    for (int i = 0; i < sc; i++) im[j, i] = a[j, i];
                return im;
            }
        }

        public Matrix Similar(Matrix M)
        {
            return M.Inverse * this * M;            
        }

        public Matrix Pseudoinverse
        {
            get
            {
                FactorizationAeqUSVt obj = new FactorizationAeqUSVt(this);
                obj.Compute();
                Matrix mi = new Matrix(sc, sc);
                for (int i = 0; i < sc; i++)
                {
                    if (Math.Abs(obj.Diagonal[i, i]) < 1.0E-13) mi[i, i] = 0.0;
                    else mi[i, i] = 1.0 / obj.Diagonal[i, i];
                }
                Matrix pim = obj.Orthogonal * mi * obj.ColumnOrthogonal.Transpose;
                return pim;
            }
        }

        private Complex[] e;
        private VectorComplex[] v;
        public virtual void Diagonalize()
        {
            if (sr != sc) try { throw new Science.Error(); }
                catch (Science.Error ee)
                {
                    ee.Write("Square matrix should be considered.");
                }
            FactorizationAeqSLSi obj = new FactorizationAeqSLSi(this);
            obj.Compute();
            e = new Complex[sr];
            v = new VectorComplex[sr];
            for (int i = 0; i < sr; i++)
            {
                e[i] = obj.Eigenvalues[i, i];
                v[i] = new VectorComplex(sr);
                for (int j = 0; j < sr; j++) v[i][j] = obj.Eigenvectors[j, i];
            }
        }
        public Complex[] EigenvalueComplex
        {
            get { return e; }
        }
        public VectorComplex[] EigenvectorComplex
        {
            get { return v; }
        }

		public static Matrix operator *(double a, Matrix c)
        {
            Matrix b = new Matrix(c.SizeOfRow, c.SizeOfColumn);
            for (int k = 0; k < c.SizeOfRow; k++)
                for (int j = 0; j < c.SizeOfColumn; j++)
                    b[k, j] = a * c[k, j];
            return b;
        }
		public static Matrix operator *(Matrix c, double a)
		{
			return a*c;
		}
		public static Matrix operator +(Matrix a, Matrix c)
        {
            Matrix b = new Matrix(a.SizeOfRow, a.SizeOfColumn);
            for (int k = 0; k < a.SizeOfRow; k++)
                for (int j = 0; j < a.SizeOfColumn; j++)
                    b[k, j] = a[k, j] + c[k, j];
            return b;
        }
		public static Matrix operator -(Matrix a, Matrix c)
		{
            Matrix b = new Matrix(a.SizeOfRow, a.SizeOfColumn);
            for (int k = 0; k < a.SizeOfRow; k++)
                for (int j = 0; j < a.SizeOfColumn; j++)
                    b[k, j] = a[k, j] - c[k, j];
            return b;
		}
        public static Matrix operator -(Matrix a)
        {
            Matrix b = new Matrix(a.SizeOfRow, a.SizeOfColumn);
            for (int k = 0; k < a.SizeOfRow; k++)
                for (int j = 0; j < a.SizeOfColumn; j++)
                    b[k, j] = - a[k, j];
            return b;
        }
	
		public static Vector operator *(Matrix a, Vector c)
		{
			Vector b = new Vector(a.SizeOfRow);
			for(int k = 0; k < b.Size; k++)
			{
				double sum = 0.0;
				for(int s = 0; s < c.Size; s++)
					sum += a[k,s]*c[s];
				b[k] = sum;
	    	}	
			return b;
		}
		public static Vector operator *(Vector c, Matrix a)
        {
            Vector b = new Vector(a.SizeOfColumn);
            for (int k = 0; k < b.Size; k++)
            {
                double sum = 0.0;
                for (int s = 0; s < c.Size; s++)
                    sum += c[s] * a[s, k];
                b[k] = sum;
            }
            return b;
        }
		public static Matrix operator *(Matrix a, Matrix c)
        {
            Matrix b = new Matrix(a.SizeOfRow, c.SizeOfColumn);
            double sum;
            for (int k = 0; k < b.SizeOfRow; k++)
                for (int j = 0; j < b.SizeOfColumn; j++)
                {
                    sum = 0.0;
                    for (int s = 0; s < a.SizeOfColumn; s++)
                        sum += a[k, s] * c[s, j];
                    b[k, j] = sum;
                }
            return b;
        }

        public void Chop()
        {
            for (int k = 0; k < sr; k++)
            {
                for (int kk = 0; kk < sc; kk++)
                {
                    if (el[k, kk] < 0.00000001 && el[k, kk] > -0.00000001) el[k, kk] = 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;
        }        
	}
}