#include <iostream>
#include <math.h>
#include <cmath>
#include <stdlib.h>
#include <fstream>
using namespace std;

//Нормирует вектор y размерности N
void normalize(double* y, int N)
{
    int i;
    double norm, result = 0;
    for (i = 0; i < N; i++)
        result += y[i] * y[i];
    norm = sqrt(result);

    for (i = 0; i < N; i++)
        y[i] = y[i] / norm;
}

//Записывает в res результат умножения матрицы А на столбец х
void Multiply_Ax(double** A, double* x, double* res, int N)
{
    int i, j;

    for (i = 0; i < N; i++) {
        res[i] = 0;
        for (j = 0; j < N; j++) {
            res[i] += A[i][j] * x[j];
        }
    }
}

void Gauss_Seidel_method(double** A, double* B, double* x, int N, double deviation)
{
    int i, j, k;
    double tmp, sum;

    //Перестановка строк в случае, если есть нулевой диаг элемент
    for (i = 0; i < N; i++) {
        if (A[i][i] == 0) {
            for (j = i + 1; j < N; j++) {
                if (A[j][i] != 0) {
                    for (k = 0; k < N; k++) {
                        tmp = A[j][k];
                        A[j][k] = A[i][k];
                        A[i][k] = tmp;
                    }
                    tmp = x[i];
                    x[i] = x[j];
                    x[j] = tmp;
                    tmp = B[i];
                    B[i] = B[j];
                    B[j] = tmp;
                }
            }
        }
    }

    //Реализация алгоритма
    k = 0;
    int flag = 1;
    while (flag > 0 && k < 15) {
        flag = 0;
        for (i = 0; i < N; i++) {
            for (j = 0; j < N; j++) {
                if (j != i)
                    sum += A[i][j] * x[j];
            }
            tmp = B[i];
            x[i] = (tmp - sum) / A[i][i];
            sum = 0;
        }
        k++;

        //Проверка невязки (B - A*x')
        for (i = 0; i < N; i++) {
            tmp = B[i];
            for (j = 0; j < N; j++)
                tmp -= A[i][j] * x[j];

            if (abs(tmp) > deviation) {
                flag = 1;
            }
        }
    }
}

void Gradient_descent_method(double** A, double* B, double* x, int N, double deviation)
{
    int i, j, k, flag = 1;
    double tau = 0, sum = 0;
    double denominator = 0, numerator = 0, tmp;
    double* r = new double[N];
    double* res = new double[N];

    k = 0;
    while (flag > 0 && k < 50) {
        flag = 0;

        for (i = 0; i < N; i++) {
            sum = 0;
            for (j = 0; j < N; j++)
                sum += A[i][j] * x[j];
            r[i] = sum - B[i];
            cout << "r[" << i << "] = " << r[i] << endl;
        }

        Multiply_Ax(A, r, res, N);
        for (j = 0; j < N; j++) {
            cout << "res " << res[j] << endl;
            denominator += res[j] * r[j];
            numerator += r[j] * r[j];
        }
        cout << "denominator = " << denominator << endl;
        cout << "numerator = " << numerator << endl;
        tau = numerator / denominator;
        denominator = 0;
        numerator = 0;
        cout << "tau " << tau << endl;

        for (i = 0; i < N; i++) {
            x[i] -= tau * r[i];
            cout << "x[" << i << "] = " << x[i] << endl;
        }
        //Проверка невязки (B - A*x' = 0)
        for (i = 0; i < N; i++) {
            tmp = B[i];
            for (j = 0; j < N; j++)
                tmp -= A[i][j] * x[j];
           // cout << "difference : " << tmp << endl;
            if (abs(tmp) > deviation)
                flag = 1;
            break;
        }
        k++;
    }
}

//Степенной метод нахождения макс собств числа
double Power_iteration(double** A, int N, double deviation)
{
    int flag = 1, i, j, k = 0;
    double lambda = 1, result;

    double* x = new double[N];
    double* y = new double[N];
    for (i = 0; i < N; i++)
        x[i] = 0;
    x[0] = 1;

    k = 0;
    while (flag > 0 && k < 10) {
        flag = 0;

        Multiply_Ax(A, x, y, N); //Строим новое приближение собств вектора у: y = A*x

        //Проверка на погрешность Ax - x*lambda = 0
        for (i = 0; i < N; i++)
            if (abs(y[i] - lambda * x[i]) > deviation) //Хотя бы одна компонента > deviation то продолжаем искать
                flag = 1;

        normalize(y, N); //нормируем найденный вектор

        for (i = 0; i < N; i++) {
            lambda += x[i] * y[i]; //Строим новое приблежение макс. собств числа: lambda = <x,y>/norma
            x[i] = y[i]; //Перезаписываем x
        }

        k++; //Кол-во итераций, ограничим их, чтобы не ушло в бесконечность если не работает алгоритм
    }
}

int main()
{
    setlocale(LC_ALL, "RUS");

    int i, j;
    int N = 3;
    double deviation = 1e-5;

    double* B = new double[N];
    double* x = new double[N];
    double** A;
    A = new double*[N];
    for (i = 0; i < N; i++) {
        A[i] = new double[N];
        x[i] = 1; //Первое приближение
    }

    A[2][0] = 3;
    A[2][1] = 6;
    A[2][2] = 4;
    A[0][0] = 0;
    A[0][1] = 5;
    A[0][2] = 7;
    A[1][0] = 0;
    A[1][1] = 0;
    A[1][2] = 2;
    B[2] = 5;
    B[0] = 2;
    B[1] = 3;

    //Gauss_Seidel_method(A, B, x, N, deviation); //Работает

    //Вывод
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++)
            cout << A[i][j] << "\t";
        cout << " = " << B[i];
        cout << "; x[" << i << "] = " << x[i] << endl;
    } cout << endl << endl;

    Gradient_descent_method(A, B, x, N, deviation);

    //Вывод
    cout << endl << endl;
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++)
            cout << A[i][j] << "\t";
        cout << " = " << B[i];
        cout << "; x[" << i << "] = " << x[i] << endl;
    }

    cout << endl;

    return 0;
}