using System using System Collections Generic using System Linq using

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Methods_First_Task
{

    #region Matrix_Class
    class Matrix
    {
        double[,] data;
        int n, k;

        public Matrix(int n, int k)
        {
            this.n = n;
            this.k = k;
            data = new double[n, k];
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < k; ++j)
                    data[i, j] = 0;
        }

        public Matrix(Matrix c)
        {
            this.n = c.n;
            this.k = c.k;
            data = new double[n, k];
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < k; ++j)
                    data[i, j] = c.data[i, j];
        }

        public int getColumns()
        {
            return k;
        }
        public int getRows()
        {
            return n;
        }

        public void makeIdentity()
        {
            if (n != k)
            {
                Console.WriteLine("MAKE IDENTITY ERROR: N != K");
                return;
            }
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < k; ++j)
                {
                    data[i, j] = ((i == j) ? 1 : 0);
                }
        }


        public Matrix(double[,] d)
        {
            n = d.GetLength(0);
            k = d.GetLength(1);
            data = new double[n, k];
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < n; ++j)
                    data[i, j] = d[i, j];
        }

        public void setValue(double v, int i, int j)
        {
            data[i, j] = v;
        }


        public double getValue(int i, int j)
        {
            return data[i, j];
        }

        static public Matrix operator *(Matrix inp1, Matrix inp2)
        {
            if (inp1.k != inp2.n)
            {
                Console.WriteLine("MATRIX MULT ERROR");
                return new Matrix(1, 1);
            }
            Matrix rs = new Matrix(inp1.n, inp2.k);
            for (int i = 0; i < inp1.n; ++i)
            {
                for (int j = 0; j < inp2.k; ++j)
                {
                    double tmp = 0;
                    for (int c = 0; c < inp1.k; ++c)
                        tmp += inp1.data[i, c] * inp2.data[c, j];
                    rs.setValue(tmp, i, j);
                }
            }

            return rs;
        }

        static public Matrix operator *(Matrix inp1, double alpha)
        {

            Matrix rs = new Matrix(inp1.n, inp1.k);
            for (int i = 0; i < inp1.n; ++i)
            {
                for (int j = 0; j < inp1.k; ++j)
                {
                    rs.setValue(inp1.data[i, j] * alpha, i, j);
                }
            }

            return rs;
        }

        static public Matrix operator -(Matrix inp1, double alpha)
        {

            Matrix rs = new Matrix(inp1.n, inp1.k);
            for (int i = 0; i < inp1.n; ++i)
            {
                for (int j = 0; j < inp1.k; ++j)
                {
                    rs.setValue(inp1.data[i, j] - alpha, i, j);
                }
            }

            return rs;
        }

        static public Matrix operator -(Matrix inp1, Matrix inp2)
        {
            if (inp1.n != inp2.n || inp1.k != inp2.k)
            {
                Console.WriteLine("Operator- ERROR: Matr size not equal");
                return new Matrix(1, 1);
            }

            Matrix rs = new Matrix(inp1.n, inp1.k);
            for (int i = 0; i < inp1.n; ++i)
            {
                for (int j = 0; j < inp1.k; ++j)
                {
                    rs.setValue(inp1.data[i, j] - inp2.data[i, j], i, j);
                }
            }

            return rs;
        }

        static public Matrix operator +(Matrix inp1, Matrix inp2)
        {
            if (inp1.n != inp2.n || inp1.k != inp2.k)
            {
                Console.WriteLine("Operator+ ERROR: Matr size not equal");
                return new Matrix(1, 1);
            }

            Matrix rs = new Matrix(inp1.n, inp1.k);
            for (int i = 0; i < inp1.n; ++i)
            {
                for (int j = 0; j < inp1.k; ++j)
                {
                    rs.setValue(inp1.data[i, j] + inp2.data[i, j], i, j);
                }
            }

            return rs;
        }

        /// Tensor multiplication of vectors
        static public Matrix operator ^(Matrix inp1, Matrix inp2)
        {
            if (inp1.k != inp2.n || inp1.n != inp2.k)
            {
                Console.WriteLine("Operator^ ERROR: Matr size not equal");
                return new Matrix(1, 1);
            }

            Matrix rs = new Matrix(inp1.n, inp2.k);
            for (int i = 0; i < inp1.n; ++i)
            {
                for (int j = 0; j < inp2.k; ++j)
                {
                    rs.setValue(inp1.data[i, 0] * inp2.data[0, j], i, j);
                }
            }

            return rs;
        }

        //транспонированная матрица
        public Matrix getTranspozed()
        {
            Matrix ret = new Matrix(k, n);
            for (int i = 0; i < n; ++i)
            {
                for (int j = 0; j < k; ++j)
                {
                    ret.setValue(data[i, j], j, i);
                }
            }
            return ret;
        }

        public void transpoze()
        {
            double[,] ret = new double[k, n];
            for (int i = 0; i < n; ++i)
            {
                for (int j = 0; j < k; ++j)
                {
                    ret[j, i] = data[i, j];
                }
            }
            data = ret;
        }

        public Matrix getSubMatrix(int si, int sj, int ei, int ej)
        {
            if (si >= n || sj >= k || ei >= n || ej >= k ||
                si < 0 || sj < 0 || ei < 0 || ej < 0)
            {
                Console.WriteLine("getSubMatrix RANGE ERROR");
                return new Matrix(1, 1);
            }

            Matrix retMtr = new Matrix(ei - si + 1, ej - sj + 1);
            int fi = 0;
            int fj = 0;
            for (int i = si; i <= ei; ++i)
            {
                fj = 0;
                for (int j = sj; j <= ej; ++j)
                {
                    retMtr.setValue(data[i, j], fi, fj);
                    ++fj;
                }
                ++fi;
            }
            return retMtr;
        }

        public bool isZeroMatrix()
        {
            for (int i = 0; i < n; ++i)
                for (int j = 0; j < k; ++j)
                    if (data[i, j] != 0)
                        return false;
            return true;
        }

        public void printMatrix()
        {
            for (int i = 0; i < n; ++i)
            {
                for (int j = 0; j < k; ++j)
                {

                    Console.Write("{0, 3:f} ", data[i, j]);
                }
                Console.WriteLine();

            }
            Console.WriteLine(); Console.WriteLine(); Console.WriteLine();

        }

        public double getMaximumAboveDiag(ref int ii, ref int jj)
        {
            double val = -1000 * 1000 * 1000;

            for (int i = 0; i < n; ++i)
            {
                for (int j = i + 1; j < k; ++j)
                {
                    if (val < Math.Abs(data[i, j]))
                    {
                        val = Math.Abs(data[i, j]);
                        ii = i;
                        jj = j;
                    }
                }
            }
            return val;
        }
    } 
    #endregion
    class Methods
    {
        public double[,] rightSide;
        int str = 0;
        int stl = 0;
        public Methods(int str, int stl)
        {
            rightSide = Methods.rightSideMatrix(str, stl);
            this.str = str;
        }
        // единичная матрица
        public static double[,] rightSideMatrix(int str, int stl)
        {
            double[,] rightSide = new double[str, stl];
            for (int i = 0; i < str; i++)
            {
                for (int j = 0; j < stl; j++)
                {
                    if (i == j)
                        rightSide[i, j] = 1;
                    else
                        rightSide[i, j] = 0;
                }
            }
            return rightSide;
        }
        //public double[,] equationMatrix = new double[3, 3] { { 3, 45, 1 }, { 2, 15, 3 }, { 9, 8, 12 } };
        // матрица A
        public double[,] equationMatrix = new double[3, 3] { { 5, 1, 2 }, { 1, 4, 1 }, { 2, 1, 3 } };

        // Метод Гаусса
        public void gauss_solver()
        {
            int n = str;//this.rightSide.Length;
            for (int i = 0; i < n; ++i)
                for (int j = i + 1; j < n; ++j)
                {
                    double value = equationMatrix[i, i];
                    double koeff = equationMatrix[j, i] / value;
                    for (int c = 0; c < n; ++c)
                        equationMatrix[j, c] -= koeff * equationMatrix[i, c];
                    // для всех элементов строки
                    for (int k = 0; k < n; k++)
                    {
                        rightSide[j, k] -= koeff * rightSide[i, k];
                    }
                }
            for (int i = n - 1; i >= 0; --i)
                for (int j = i - 1; j >= 0; --j)
                {
                    double value = equationMatrix[i, i];
                    double koeff = equationMatrix[j, i] / value;
                    for (int c = 0; c < n; ++c)
                        equationMatrix[j, c] -= koeff * equationMatrix[i, c];
                    for (int k = 0; k < n; k++)
                    {
                        rightSide[j, k] -= koeff * rightSide[i, k];
                    }
                }
            for (int i = 0; i < n; ++i)
            {

                for (int k = 0; k < n; k++)
                {
                    rightSide[i, k] = rightSide[i, k] / equationMatrix[i, i];
                }
            }

        }
        static double getNorm(Matrix a)
        {
            return Math.Sqrt((a.getTranspozed() * a).getValue(0, 0));
        }

        //public double 
        //Matrix[] triada = new Matrix[3];

        // метод скалярных произведений
        public void smAlgorithm(Matrix A, Matrix y, ref Matrix x, ref double v, bool secondIter = false)
        {
            int n = A.getColumns();
            double EPS = 0.000001;

            double s = (y.getTranspozed() * y).getValue(0, 0);
            double yNorm = Math.Sqrt(s);
            Matrix lx = y * (1.0f / yNorm);
            double lamdda = 0;
            double lastLamda = 0;
            for (int iter = 0; iter < 100000; ++iter)
            {
                y = A * lx;
                if (secondIter)
                    y = y - x * ((x.getTranspozed() * y).getValue(0, 0));
                s = (y.getTranspozed() * y).getValue(0, 0);
                double t = (lx.getTranspozed() * y).getValue(0, 0);
                yNorm = Math.Sqrt(s);
                lx = y * (1.0 / yNorm);
                lastLamda = lamdda;
                lamdda = s / t;

                if (Math.Abs(lastLamda - lamdda) < EPS)
                    break;
            }
            v = lamdda;
            x = lx;
        }

        public void findSecondBySm(Matrix A)
        {
            Matrix y = new Matrix(A.getColumns(), 1);
            Matrix x1 = new Matrix(A.getColumns(), 1);
            double vl = 0;

            for (int i = 0; i < A.getColumns(); ++i) y.setValue(1, i, 0);
            smAlgorithm(A, y, ref x1, ref vl);
            Console.WriteLine("Максимально собственное значение: {0}: ", vl);
            Console.WriteLine("Собственный вектор максимального числа:");
            x1.printMatrix();
            for (int i = 0; i < A.getColumns(); ++i) y.setValue(123, i, 0);

            y = y - x1 * ((x1.getTranspozed() * y).getValue(0, 0));
            smAlgorithm(A, y, ref x1, ref vl, true);

            //Console.WriteLine("adfasdf");
            Console.WriteLine("Второе по величине собственное число: {0}", vl);
            //Console.WriteLine(vl);
            Console.WriteLine("Его собственный вектор:");
            x1.printMatrix();

        }

        // Метод Якоби
        public static void Jacobi(int n, double[][] A, double[] F, double[] X)
        {
            double EPS = 0.000001;
            double[] TempX = new double[n];
            double norm;
            for (int i = 0; i < n; i++)
            {
                TempX[i] = X[i];
            }
            do
            {
                for (int i = 0; i < n; i++)
                {
                    TempX[i] = F[i];
                    for (int g = 0; g < n; g++)
                    {
                        if (i != g)
                            TempX[i] -= A[i][g] * X[g];
                    }
                    TempX[i] /= A[i][i];
                }
                norm = Math.Abs(X[0] - TempX[0]);
                for (int h = 0; h < n; h++)
                {
                    if (Math.Abs(X[h] - TempX[h]) > norm)
                        norm = Math.Abs(X[h] - TempX[h]);
                    X[h] = TempX[h];
                }
            }
            while (norm > EPS);
        }
        // Метод Зейделя
        static bool ZeidelEnd(double[] xk, double[] xp, int n, double EPS)
        {
            double norm = 0;
            for (int i = 0; i < n; i++)
            {
                norm += (xk[i] - xp[i]) * (xk[i] - xp[i]);
            }
            if (Math.Sqrt(norm) >= EPS)
                return false;
            return true;
        }
        public static void Zeidel(int n, double[,] a, double[] x, double[] p, double[] rightSide)
        {
            do
            {
                for (int i = 0; i < n; i++)
                {
                    p[i] = x[i];
                }
                for (int i = 0; i < n; i++)
                {
                    double var = 0;
                    for (int j = 0; j < i; j++)
                    {
                        var += (a[i, j] * x[j]);
                    }
                    for (int j = i + 1; j < n; j++)
                    {
                        var += (a[i, j] * p[j]);
                    }
                    x[i] = (rightSide[i] - var) / a[i,i];
                }
            }
            while(!ZeidelEnd(x, p, n, 0.0001));
        }

    }
    class Program
    {
        static void Main(string[] args)
        {
            double[,] test = Methods.rightSideMatrix(3, 3);
            Methods a = new Methods(3, 3);
            double[,] test2 = a.rightSide;
            a.gauss_solver();
            Console.WriteLine("Gauss:");
            for (int i = 0; i < 3; i++)
            {
                for (int j = 0; j < 3; j++)
                {
                    Console.Write(a.rightSide[i, j].ToString("0.00"));
                    Console.Write(" ");
                }
                Console.WriteLine();
            }
            double[,] res = new double[3, 3] { { 5, 1, 2 }, { 1, 4, 1 }, { 2, 1, 3 } };
            Matrix eqMatrix = new Matrix(res);
            a.findSecondBySm(eqMatrix);

            // Метод Якоби
            //double[][] JA = new double[3][] { new double[3]{ 1, 1, 1 }, new double[3]{4, 2, 1 }, new double[3]{9,3,1 } };
            //double[] JFree = new double[] { 0, 1, 3 };
            //double[] JRes = new double[3] {1,1,1};
            //Methods.Jacobi(3, JA, JFree, JRes);
            //Console.WriteLine("{0} {1} {2}", JRes[0], JRes[1], JRes[2]);

            //Метод Зейделя 
            double[,] ZA = new double[3, 3] { { 1, 1, 1 }, { 4, 2, 1 }, { 9, 3, 1 } };
            double[] ZFree = new double[] { 0, 1, 3 };
            double[] ZRes = new double[3] { 0, 0, 0 };
            double[] Zpred = new double[3];
            Methods.Zeidel(3, ZA, ZRes, Zpred, ZFree);
            
            Console.ReadLine();
        }
    }
}