#include <iostream>
#include <vector>
#include <string>
#include<time.h>
#include<cmath>

using namespace std;
//////////////////////////////////////////////////////////////////////////

    class Matrix{

    public:
        int m;
        int n;
        vector<vector<double> > vvi;
        Matrix() {}
        Matrix(int m, int n) {
            this->m = m;
            this->n = n;
            vvi.resize(m);
            for (int i = 0; i < m; i++) {
                vvi[i].resize(n);
            }
        }
        void getPrint() {
            for (int i = 0; i < this->m; i++)
            {
                cout << "| ";
                for (int j = 0; j < this->n; j++)
                    cout << "" << this->vvi[i][j] << " ";
                cout << "|";
                //cout << "" << endl;
                cout << "" << endl;
            }
            cout << endl;
        }
        void setValue() {
            for (int i = 0; i < this->m; i++)
                for (int j = 0; j < this->n; j++)
                    cin >> this->vvi[i][j];

        }
        void fillMatrix() {
            for (int i = 0; i < this->m; i++)
                for (int j = 0; j < this->n; j++)
                    this->vvi[i][j] = rand() % 2;

        }

    };
    Matrix E(Matrix& a) {
        Matrix e(a.m,a.m);
        for (int i = 0; i < e.m; i++)
            for (int j = 0; j < e.n; j++)
                e.vvi[i][j]=0;
        for (int i = 0; i < e.m; i++)
            e.vvi[i][i] = 1;
        return e;
    }

    Matrix E(int m) {
        Matrix e(m, m);
        for (int i = 0; i < m; i++)
            for (int j = 0; j < m; j++)
                e.vvi[i][j] = 0;
        for (int i = 0; i < m; i++)
            e.vvi[i][i] = 1;
        return e;
    }
    template <typename T>
    Matrix Sum(T& a, T& b) {
        if ((a.n != b.n) or (a.m != b.m))
        {
        }   //cout << "Некорректное выражение!" << endl;
        else
        {
            int n = a.n;
            int m = a.m;
            Matrix c(m, n);
            for (int i = 0; i < m; i++)
            {
                for (int j = 0; j < n; j++)
                    c.vvi[i][j] = a.vvi[i][j] + b.vvi[i][j];
            }
            return c;
        }
    }
    template <typename T>
    Matrix Mult(T& a, T& b) {
        int m = a.m;
        int n = b.n;
        int f = b.m;
        Matrix c(m, n);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++)
            {
                c.vvi[i][j] = 0;
                for (int k = 0; k < f; k++)
                    c.vvi[i][j] += a.vvi[i][k] * b.vvi[k][j];
            }
        }
        return c;
    }

    Matrix Transp(Matrix& a) {
        Matrix c(a.n, a.m);
        for (int i = 0; i < a.m; i++)
            for (int j = 0; j < a.n; j++)
                c.vvi[j][i] = a.vvi[i][j];
        return c;
    }



    Matrix GetMatr(Matrix& a, int i, int j) {
        int ki, kj, di, dj;
        Matrix p(a.m - 1, a.n - 1);
        di = 0;
        for (ki = 0; ki < a.m - 1; ki++) { // проверка индекса строки
            if (ki == i) di = 1;
            dj = 0;
            for (kj = 0; kj < a.m - 1; kj++) { // проверка индекса столбца
                if (kj == j) dj = 1;
                p.vvi[ki][kj] = a.vvi[ki + di][kj + dj];
            }
        }
        return p;
    }



    Matrix operator * (Matrix a, Matrix b) {
        int m = a.m;
        int n = b.n;
        int f = b.m;
        Matrix c(m, n);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++)
            {
                c.vvi[i][j] = 0;
                for (int k = 0; k < f; k++)
                    c.vvi[i][j] += a.vvi[i][k] * b.vvi[k][j];
            }
        }
        return c;
    }



    double det(Matrix a) {
        int i, j, k, n;
        double d;
        Matrix p(a.m, a.n);
        j = 0; d = 0;
        k = 1; //(-1) в степени i
        n = a.m - 1;
        if (a.m < 1) cout << "Определитель вычислить невозможно!";
        if (a.m == 1) {
            d = a.vvi[0][0];
            return d;
        }
        if (a.m == 2) {
            d = a.vvi[0][0] * a.vvi[1][1] - (a.vvi[1][0] * a.vvi[0][1]);
            return d;
        }
        if (a.m > 2) {
            for (i = 0; i < a.m; i++) {
                p = GetMatr(a, i, 0);
                //cout << a.vvi[i][j] << endl;
                if (a.vvi[i][0] == 0) continue;
                if(a.vvi[i][0] != 0) d += k * a.vvi[i][0] * det(p);
                k = -k;
            }
        }
        return d;
    }

    ////////////////////////////////////////////////////////////////////////////

    class Vector :public Matrix {
    public:
        Vector(int m) {
            this->m = m;
            this->n = 1;
            vvi.resize(m);
            for (int i = 0; i < m; i++) {
                vvi[i].resize(n);
            }
        }
    };
    Matrix transV(Vector a) {
        Matrix x(a.m, 1);
        for (int i = 0; i < x.m; i++)
            x.vvi[i][0] = a.vvi[i][0];
        return x;
    }
    Vector transM(Matrix a) {
        Vector x(a.m);
        for (int i = 0; i < x.m; i++)
            x.vvi[i][0] = a.vvi[i][0];
        return x;
    }
    template <typename T>
    T operator + (T a, T b) {
        int n = a.n;
        int m = a.m;
        for (int i = 0; i < m; i++)
        {
            for (int j = 0; j < n; j++)
                a.vvi[i][j] += b.vvi[i][j];
        }
        return a;

    }
    template <typename T>
    T operator - (T a, T b) {
        int n = a.n;
        int m = a.m;
        for (int i = 0; i < m; i++)
        {
            for (int j = 0; j < n; j++)
            a.vvi[i][j] -= b.vvi[i][j];
            }
            return a;
        }

    template <typename T>
    T operator * (double a, T x) {
        for (int i = 0; i < x.m; i++)
            for (int j = 0; j < x.n; j++)
                x.vvi[i][j] = x.vvi[i][j] * a;
        return x;

    }
    Vector operator * (Matrix a, Vector b) {
        Matrix f = a * transV(b);
        Vector x = transM(f);
        return x;
    }


    Vector getVector(Matrix& a, int n) {
        Vector z(a.m);
        for (int i = 0; i < a.m; i++) {
            z.vvi[i][0] = a.vvi[i][n];
            //cout << "Типа скопировали " << z.vvi[i][0] << endl;
        }
        //getPrint(z);
        return z;
    }

    double dotProduct(Vector a, Vector b) {// Скалярное произведение
        double d = 0;
        for (int i = 0; i < a.m; i++)
            d += a.vvi[i][0] * b.vvi[i][0];
        return d;
    }

    Matrix addVector(Matrix& a, Vector& b) {
        a.n = a.n + 1;
        for (int i = 0; i < a.m; i++)
            a.vvi[i].push_back(b.vvi[i][0]);
        return a;
    }

    void setColumns(Matrix& a, Vector& b, int k) {
        for (int i = 0; i < a.m; i++)
            a.vvi[i][k] = b.vvi[i][0];

    }
    ////////////////////////////////////////////////////////////////////////////
    //Метод Гаусса

    Matrix makeDiag(Matrix& a) {
        int i, j, m = 0, n = 0;
        double k;
        double c;

        for (int n = 0; n < a.n; n++)
        {
            if ((m < a.m) and (a.n > n)) {
                k = a.vvi[m][n];
                if (k == 0) goto end1;
                for (i = 0; i < a.n; i++)
                    a.vvi[m][i] = a.vvi[m][i] / k;//поделили на элемент и получили единицу
                for (i = 0; i < a.m; i++)//вычитаем строчку из строчек
                    if (i != m)
                    {
                        c = a.vvi[i][n];
                        for (j = 0; j < a.n; j++)
                            a.vvi[i][j] = a.vvi[i][j] - c * a.vvi[m][j];
                    }
            end1:
                m++;
                //getPrint(a);
            }
        }
        for (i = 0; i < a.m; i++)
            for (j = 0; j < a.n; j++)
                if (a.vvi[i][j] == -0)
                    a.vvi[i][j] = 0;
        return a;

    }//

    Matrix Invert(Matrix& a) {
        Matrix b(a.m, 2 * a.n);
        for (int i = 0; i < a.m; i++)
            for (int j = 0; j < a.n; j++)
                b.vvi[i][j] = a.vvi[i][j];
        for (int i = 0; i < a.m; i++)
            b.vvi[i][i + a.n] = 1;
        b = makeDiag(b);
        Matrix c(a.m, a.n);
        for (int i = 0; i < a.m; i++)
            for (int j = 0; j < a.n; j++)
                c.vvi[i][j] = b.vvi[i][j + a.n];
        return c;
    }
    template <typename T>
    bool operator == (T a, T b) {
        int k = 1;
        for (int i = 0; i < a.m; i++)
            for (int j = 0; j < a.n; j++)
                if (a.vvi[i][j] != b.vvi[i][j])
                    k = 0;
        if ((a.m != b.m) or (a.n != b.n))
            k = 0;
        return k;
    }
    template <typename T>
    bool operator != (T a, T b) {
        int k = 1;
        for (int i = 0; i < a.m; i++)
            for (int j = 0; j < a.n; j++)
                if (a.vvi[i][j] == b.vvi[i][j])
                    k = 0;
        return k;
    }
    bool isOrth(Matrix& a) {
        Matrix b = Transp(a);
        if (((a * b) == E(a)) == 1 )
            return true;
        else
            return false;
    }

    Vector solveGauss(Matrix& a, Vector& f) {
        a = addVector(a, f);
        a = makeDiag(a);
        Vector x = getVector(a, a.n - 1);
        cout << "Solution:" << endl;
        return x;
    }

    bool isSymetric(Matrix& a) {
        int k = 1;
        for (int i = 0; i < a.m - 1; i++)
            for (int j = 1; j < a.n; j++)
                if (a.vvi[i][j] != a.vvi[j][i])
                    k = 0;
        if (a.m != a.n)
            k = 0;
        return k;
    }

    bool isPosDetermined(Matrix& a) {// Реализация критерия Сильвестра
        int k = 1;
        for (int i = 0; i < a.m; i++)
        {
            Matrix tmp(i + 1, i + 1);
            for (int k = 0; k < i + 1; k++)
                for (int q = 0; q < i + 1; q++)
                    tmp.vvi[k][q] = a.vvi[k][q];
            if ((det(tmp) > 0) == 0)
                k = 0;


        }
        return k;
    }

    ////////////////////////////////////////////////////////////////////////////
    //Правило Крамера
    Vector solveCramer(Matrix& a, Vector& b) {
        Vector res(a.m);
        Vector r(a.m);
        double Delta = det(a);
        if (Delta != 0)
            for (int i = 0; i < a.n; i++)
            {
                r = getVector(a, i);
                setColumns(a, b, i);
                //cout << "a" << endl;
                double delta = det(a);
                cout << "delta" << delta << endl;
                res.vvi[i][0] = delta / Delta;
                for (int j = 0; j < a.m; j++)
                    a.vvi[j][i] = r.vvi[j][0];
                //cout << res.vvi[i][0] << endl;
            }
        cout << "Solution is" << endl;
        return res;

    }
    ////////////////////////////////////////////////////////////////////////////
    // Метод сопряженных градиентов
    Vector solveCG(Matrix a, Vector b) {

            Vector x(a.m);
            int i = 0;
            double alpha = 0;
            double betha = 0;
            for (int i = 0; i < x.m; i++)
                x.vvi[i][0] = 2;
            Vector r = b - a * x;
            Vector tmp(r.m);
            Vector z = r;
            while (i<100){
                 //начало i-ой итерации

            alpha = dotProduct(r, r) / dotProduct(a*z, z);
            x = x + alpha * z;
            tmp = r;
            r = r - alpha * (a*z);
            betha = alpha = dotProduct(r, r) / dotProduct(tmp, tmp);
            z = r + betha * z;
            i++;
            x.getPrint();
            }
            for (int j = 0; j < x.m; j++)
                if (abs(x.vvi[j][0]) < 0.001)// с точностью до 8 знака
                    x.vvi[j][0] = 0;
            return x;

    }


    ////////////////////////////////////////////////////////////////////////////

    int main()
    {
        unsigned int start_time = clock();
        setlocale(LC_ALL, "ru");
        Matrix A(3, 3);
        A.setValue();
        Vector b(A.m);
        b.setValue();
        Vector x = solveCG(A, b);
        x.getPrint();
        unsigned int end_time = clock(); // конечное время
        unsigned int search_time = end_time - start_time;
        system("pause");
        return 0;
    }