using System using System Collections Generic using System Linq using

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

namespace QR {
    class Program {
        public static Matrix QR(Matrix A) {
            int n = A.GetLength(0),
                m = A.GetLength(1);
            Matrix y,
                   Q = new Matrix(n, m).MakeIdentity();
            for (int i = 0; i < n - 1; ++i) {
                y = A.GetColumn(i, i);
                Matrix curV = new Matrix(n - i, m - i).MakeIdentity(),
                       curQ = new Matrix(n, m).MakeIdentity(),
                       omega = new Matrix(n - i, 1),
                       e = new Matrix(n - i, 1); e[0, 0] = 1;
                double yNorm = Math.Sqrt((y.GetTransposed() * y)[0, 0]),
                       omegaNorm = 0.0;
                int sign = Math.Sign(y[0, 0]);

                omega = y + sign * (yNorm * e);
                omegaNorm = Math.Sqrt((omega.GetTransposed() * omega)[0, 0]);
                omega *= (1.0 / omegaNorm);
                curV = curV - (omega * omega.GetTransposed()) * 2.0;
                curQ.SetMinor2(curV, i);

                A = curQ * A;
                Q = curQ * Q;
                
            }
            return Q.GetTransposed();
        }

        public static Matrix GetSolve(Matrix A, Matrix b, double eps) {
            int n = A.GetLength(0);
            Matrix Q = QR(A),
                   R = Q.GetTransposed() * A,
                   x = new Matrix(n, 1),
                   nb = new Matrix(n, 1);

            Console.WriteLine("Матрица Q:");
            Console.WriteLine(Q);
            Console.WriteLine("Матрица R:");
            Console.WriteLine(R);

            for (int i = 0; i < n; i++) {
                nb[i, 0] = 0;
                for (int j = 0; j < n; j++)
                    nb[i, 0] += b[j, 0] * Q[j, i];
            }

            for (int i = n - 1; i >= 0; --i) {
                double s = 0.0;
                for (int j = i + 1; j < n; ++j)
                    s += R[i, j] * x[j, 0];
                    x[i, 0] = (nb[i, 0] - s) / R[i, i];
            }

            
            Console.WriteLine("Решение СЛАУ:");
            Console.WriteLine(x);
            return x;
        }

        public static Matrix TryNewSolve(Matrix A, Matrix b, double eps) {
            int n = A.GetLength(0);
            Matrix Q = QR(A),
                   R = Q.GetTransposed() * A,
                   x = new Matrix(n, 1),
                   nb = new Matrix(n, 1);

            for (int i = 0; i < n; i++) {
                nb[i, 0] = 0;
                for (int j = 0; j < n; j++)
                    nb[i, 0] += b[j, 0] * Q[j, i];
            }

            x[A.GetLength(0) - 1, 0] = 1.0;

            for (int i = n - 2; i >= 0; --i) {
                double s = 0.0;
                for (int j = i + 1; j < n; ++j)
                    s += R[i, j] * x[j, 0];
                x[i, 0] = (nb[i, 0] - s) / R[i, i];
            }

            for (int i = 0; i < n; ++i)
                x[i, 0] /= Math.Sqrt(x.GetSummColumn(0, 0));
            return x;
        }

        public static List<double> GetEigenvalues(Matrix A, double eps) {
            Matrix reserve = new Matrix(A);
            List<double> retValue = new List<double>();
            int countIter = 0;
            int indColumn = A.GetLength(0) - 2;
            while (true) {
                bool stop = true;
                for (int i = indColumn; i >= 0; --i)
                    if (A.GetColumn(i, i + 1).GetNorm() >= eps)
                        stop = false;
                if (stop)
                    break;
                if (countIter > 100) {
                    countIter = 0;
                    indColumn--;
                    A = reserve;
                }
                Matrix Q = QR(A),
                       R = Q.GetTransposed() * A;
                A = R * Q;
                countIter++;
            }
            for (int i = 0; i < A.GetLength(0) - 1; ++i) {
                if (A.GetColumn(i, i + 1).GetNorm() < eps) {
                    retValue.Add(A[i, i]);
                }
                if (retValue.Count == A.GetLength(0) - 1)
                    retValue.Add(A[A.GetLength(0) - 1, A.GetLength(1) - 1]);
            }
            
            return retValue;
        }

        public static List<Matrix> GetEigenvectors(Matrix A, List<double> eigenvalues) {
            List<Matrix> retValue = new List<Matrix>();
            Matrix b = new Matrix(A.GetLength(0), 1);
            for (int i = 0; i < A.GetLength(0); ++i)
                b[i, 0] = 0;
            Matrix reserve = new Matrix(A);
            for (int i = 0; i < eigenvalues.Count; ++i) {
                A = new Matrix(reserve);
                
                A = A - A.MakeIdentity() * eigenvalues[i];
                retValue.Add(TryNewSolve(A, b, 0.001));
            }
            return retValue;
        }

        static void Main(string[] args) {
            Random rand = new Random();
            Matrix A = new Matrix(new double[,] {{21, 18, 6},
                                                 {18, 17, 5},
                                                 {6, 5, 2}});
            Matrix b = new Matrix(new double[,] {{1}, 
                                                 {2}, 
                                                 {3}});

            

            Console.WriteLine("Исходная матрица А:");
            Console.WriteLine(A);
            Console.WriteLine("Вектор свободных членов:");
            Console.WriteLine(b);
            
            GetSolve(A, b, 0.00001);
            List<double> eigenValues = GetEigenvalues(A, 0.0000001);
            List<Matrix> eigenVectors = GetEigenvectors(A, eigenValues);

            for (int i = 0; i < eigenVectors.Count; ++i) {
                Console.WriteLine("Собственное значение : {0:f2}", eigenValues[i]);
                Console.WriteLine("Собственный вектор :");
                Console.WriteLine(eigenVectors[i]);
                Console.WriteLine();
            }
           
            
            
            
        }
    }
}