Iteracii_the_same

#define _CRT_SECURE_NO_WARNINGS
#include <conio.h>
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <math.h>
using namespace std;
std::ifstream f("in.txt");
std::ofstream out("1.txt");
double eps = 1.0e-12;
#define PI 3.14159265

class Matrix
{
private:
    int m;
    int n;
    double **arr; // двумерный массив
public:
    void Create() {
        arr = new double*[n];
        for (int i = 0; i < n; i++)
            arr[i] = new double[m];
    }
    Matrix(int n, int m)
    {
        this->n = n;
        this->m = m;
        Create();
    }
    void del() {
        for (int i = 0; i < m; i++)
            delete[] arr[i];
        delete[] arr;
    }
    int Get_n() {
        return n;
    }
    int Get_m() {
        return m;
    }
    double Get_arr(int i, int j) {
        return arr[i][j];
    }
    void Copy(Matrix &right) {
        del();
        Create();
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                arr[i][j] = right.arr[i][j];
                cout << arr[i][j] << " ";
            }
            cout << "\n";
        }
        cout << "\n";
    };
    friend void fout(Matrix A);
    Matrix operator= (Matrix &right);
    Matrix operator= (Matrix right);
    Matrix operator+= (Matrix right);
    Matrix operator-= (Matrix right);
    friend Matrix operator+ (Matrix left, Matrix right);
    friend Matrix operator- (Matrix left, Matrix right);
    friend void Mfile(int a, Matrix &A, Matrix &b);
    friend Matrix operator*(Matrix A, Matrix B);
    friend Matrix operator*(double a, Matrix A);
    friend double Gershgorin(Matrix A, double *m, double *M);
    friend double Norm2(Matrix A);
    friend void Cheb_ord(Matrix &A);
    friend void Cheb(double *t, Matrix &P, double m, double M);
    friend Matrix I(int n);
};

void fout(Matrix A) {
    for (int i = 0; i < A.n; i++) {
        for (int j = 0; j < A.m; j++) {
            out << A.arr[i][j] << " ";
        }
        out << "\n";
    }
    out << "\n";
}

Matrix Matrix::operator= (Matrix &right) {
    if (m != right.m || n != right.n) {
        cout << "error \n";
        exit(1);
    }
    if (this != &right)
    Copy(right);
    return *this;
}

Matrix Matrix::operator= (Matrix right) {
    if (m != right.m || n != right.n) {
        cout << "error \n";
        exit(1);
    }
    for (int i = 0; i < n; i++)
        for (int j = 0; j < m; j++)
            arr[i][j] = right.arr[i][j];
    return *this;
}

Matrix Matrix::operator+= (Matrix right)
{
    if (m != right.m || n != right.n) {
        cout << "error \n";
        exit(1);
    }
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            arr[i][j] += right.arr[i][j];
            cout << arr[i][j];
            cout << " ";
        }
        cout << "\n";
    }
    cout << "\n";
    return *this;
}

Matrix Matrix::operator-= (Matrix right)
{
    if (m != right.m || n != right.n) {
        cout << "error \n";
        exit(1);
    }
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            arr[i][j] -= right.arr[i][j];
            cout << arr[i][j] << " ";
        }
        cout << "\n";
    }
    cout << "\n";
    return *this;
}

Matrix operator+(Matrix left, Matrix right) {
    if (left.n != right.n || right.m != left.m) {
        cout << "error \n";
        exit(1);
    }
    Matrix add(left.n, left.m);
    for (int i = 0; i < left.n; i++)
        for (int j = 0; j < left.m; j++)
            add.arr[i][j] = left.arr[i][j] + right.arr[i][j];
    return add;
}

Matrix operator*(Matrix A, Matrix B) {
    if (A.m != B.n) {
        cout << "error";
        exit(1);
    }
    Matrix mult(A.n, B.m);
    for (int i = 0; i < A.n; i++){
        for (int k = 0; k < B.m; k++) {
            double t = 0;
            for (int j = 0; j < A.n; j++) t += A.arr[i][j] * B.arr[j][k];
            mult.arr[i][k] = t;
            cout << mult.arr[i][k] << " ";
        }
        cout << "\n";
    }
    cout << "\n";
    return mult;
}

Matrix operator*(double a, Matrix A) {
    Matrix mult(A.n, A.m);
    for (int i = 0; i < A.n; i++) {
        for (int j = 0; j < A.m; j++) {
            mult.arr[i][j] = a * A.arr[i][j];
            cout << mult.arr[i][j] << " ";
        }
        cout << "\n";
    }
    return mult;
}

Matrix operator-(Matrix left, Matrix right) {
    if (left.n != right.n || left.m != right.m) {
        cout << "error";
        exit(1);
    }
    Matrix sub = (-1)*right;
    sub += left;
    return sub;
}

void Mfile(int N, Matrix &A, Matrix &b) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            f >> A.arr[i][j];
            cout << A.arr[i][j] << " ";
        }
        cout << "\n";
    }
    for (int i = 0; i < N; i++) {
        f >> b.arr[i][0];
        cout << b.arr[i][0] << " ";
    };
    cout << "\n";
}

Matrix I(int n) {
    Matrix I(n, n);
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < i; j++) {
            I.arr[i][j] = 0;
            cout << I.arr[i][j] << " ";
        }
        I.arr[i][i] = 1;
        cout << I.arr[i][i] << " ";
        for (int j = i+1; j < n; j++) {
            I.arr[i][j] = 0;
            cout << I.arr[i][j] << " ";
        }
        cout << "\n";
    }
    return I;
}

double Gershgorin(Matrix A, double *m, double *M) {
    if (A.n != A.m) {
        cout << "error";
        exit(1);
    }
    double t = 0;
    *m = A.arr[0][0], *M = A.arr[0][0];
    for (int i = 0; i < A.n; i++) {
        t = 0;
        for (int j=0; j<i; j++) t += fabs(A.arr[i][j]);
        for (int j = i+1; j < A.n; j++) t += fabs(A.arr[i][j]);
        if (*m > A.arr[i][i] - t) *m = A.arr[i][i] - t;
        if (*M < A.arr[i][i] + t) *M = A.arr[i][i] + t;
    }
    if (*m + *M == 0) {
        cout << "error";
        exit(1);
    }
    return 2 / (*m + *M);
}

double Norm2(Matrix A) {
    double t = 0;
    for (int i = 0; i < A.n; i++) {
        t += A.arr[i][0]*A.arr[i][0];
    }
    double M = sqrt(t);
    return M;
}

void Cheb_ord(Matrix &P) { //0 строка - параметры P[1], 1 - P[2], 2- P[4] и тд
    int p = 1;
    P.arr[0][0] = 1;
    cout << p << P.arr[0][0] << "\n";
    for (int i = 1; i < P.n; i++) {
        p *= 2;
        for (int j = 0; j < p; j+=2) {
            P.arr[i][j] = P.arr[i-1][j/2];
            P.arr[i][j + 1] = 2*p - P.arr[i][j];
            cout << P.arr[i][j] << " " << P.arr[i][j + 1] << " ";
        }
        cout << p << "\n";
    }
}

void Cheb(double *t, Matrix &P, double m, double M) {
    Cheb_ord(P);
    int p = P.n - 1, p2 = P.m;
    double k, a = m + M, b = M - m;
    for (int i = 0; i < p2; i++) {
        k = P.arr[p][i] * PI / (2*p2);
        t[i] = 2 / (a - b * cos(k));
        cout << t[i] << " ";
    }
    cout << "\n";
}

int main()
{
    int N, p;
    f >> N;
    f >> p;
    Matrix E = I(N);
    Matrix A(N, N), b(N, 1);
    Mfile(N, A, b);
    double m, M;
    double tau=Gershgorin(A, &m, &M);
    cout << tau;
    Matrix X0 = 0 * b;
    out << " Linear equation A*x = b. We initialize X[0] to the solution of the equation.\n";
    out << "Matrix A = \n";
    fout(A);
    out << "Vector b = \n";
    fout(b);
    out << "Initial guess X[0] = \n";
    fout(X0);
    Matrix A1 = E + (-tau)*A;
    Matrix b1 = tau * b;
    Matrix X1 = A1 * X0 + b1;
    Matrix Res = A * X1 - b;
    double x = Norm2(Res);
    int k = 1;
    out << "1) Norm of Res[" << k << "] = " << x << "\n";
    fout(X1);
    while (x >eps) {
        if (k % 2) {// то есть мы по X1 строим X2, но чтобы обойтись двумя векторами, стираем и заполняем заново X0;
            X0 = A1 * X1 + b1; //X[n+1]=(E-tau*A)*X[n]+tau*b; n+1 - количество итераций - одну мы выполнили до while;
            Res = A * X0 - b;
        }
        else {
            X1 = A1 * X0 + b1;
            Res = A * X1 - b;
        }
        x = Norm2(Res);
        k++;
        out << k << ") Norm of Res[" << k << "] = " << x << "\n";
    }
    out << "Fixed-point iteration method with " << k << " iterations gives the solution:\n";
    if (k % 2) fout(X0); //решение методом простых итераций
    else fout(X1); // или это
    cout << "\n\n\n";

    int p2 = (int)pow(2, p);
    Matrix P(p+1, p2);
    double *t= new double[p2];
    Cheb(t, P, m, M);
    X0 = 0 * b;
    X1 = X0 + t[0] * b - t[0] * A * X0;
    Res = A * X1 - b;
    x = Norm2(Res);
    k = 1;
    int k1 = 1;
    out << "1) Norm of Res[" << k << "] = " << x << "\n";
    fout(X1);
    while (x > eps) {
        k1 = (k - 1) % p2;
        if (k % 2) {// то есть мы по X1 строим X2, но чтобы обойтись двумя векторами, стираем и заполняем заново X0;
            X0 = X1 + t[k1] * b - t[k1] * A * X1;
            Res = A * X0 - b;
        }
        else {
            X1 = X0 + t[k1] * b - t[k1] * A * X0;
            Res = A * X1 - b;
        }
        x = Norm2(Res);
        k++;
        out << "Parameter " << t[k1] << "\n";
        out << k << ") Norm of Res[" << k << "] = " << x << "\n";
    }
    out << "Chebyshev-point iteration method with " << k << " iterations gives the solution:\n";
    if (k % 2) fout(X0); //решение методом чебышёвских параметров
    else fout(X1); // или это
    X0.del();
    X1.del();

    out.close();
    return 0;
}