lab3.cpp

using namespace std;

#include <vector>
#include <iostream>
#include <iomanip>
#include <cmath>

template <int N>
struct Array {
    double data[N];

    double &operator[](int index) {
        return data[index];
    }

    const double &operator[](int index) const {
        return data[index];
    }

    int size() {
        return N;
    }
};

template<int N>
Array<N> operator+(const Array<N> &a, const Array<N> &b) {
    Array<N> r;
    for (int i = 0; i < N; i++) {
        r[i] = a[i] + b[i];
    }
    return r;
}

template<int N>
Array<N> operator-(const Array<N> &a, const Array<N> &b) {
    Array<N> r;
    for (int i = 0; i < N; i++) {
        r[i] = a[i] - b[i];
    }
    return r;
}

template<int N>
Array<N> operator*(const Array<N> &a, const double &b) {
    Array<N> r;
    for (int i = 0; i < N; i++) {
        r[i] = a[i] * b;
    }
    return r;
}

template <int N, int M>
struct Matrix {
    Array<M> data[N];

    Array<M> &operator[](int index) {
        return data[index];
    }
};


template <int N, int M>
void print(Matrix<N, M> A) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < M; j++) {
            cout << setw(5) << A[i][j] << ' ';
        }
        cout << endl;
    }
}

template <int N, int M>
void single_division(Matrix<N, M> A) {
    cout << "Single Division method" << endl;
    print(A);

    Array<N> x;

    for (int j = 0; j < N; j++) {
        double a = A[j][j];
        for (int i = j; i < N + 1; i++) {
            A[j][i] /= a;
        }
        for (int k = j + 1; k < N; k++) {
            double a = A[k][j];
            for (int l = j; l < N + 1; l++) {
                A[k][l] -= A[j][l] * a;
            }
        }
    }
    for (int j = N - 1; j >= 0; j--) {
        double s = 0;
        for (int i = N - 1; i > j; i--) {
            s += A[j][i] * x[i];
        }
        x[j] = A[j][N] - s;
    }

    for (int i = 0; i < x.size(); i++) {
        cout << "x[" << (i + 1) << "] = " << x[i] << endl;
    }
}

template <int N, int M>
void gauss_seidel(Matrix<N, M> matrix, double eps) {
    cout << "Gauss-Seidel method" << endl;
    print(matrix);

    Array<N> x, x0, beta;

    for (int i = 0; i < N; i++) {
        x[i] = beta[i] = matrix[i][N] / matrix[i][i];
    }

    Matrix<N, N> alfa;
    for(int i = 0; i < N; i++) {
        for(int j = 0; j < N; j++) {
            alfa[i][j] = (i != j) ? (-matrix[i][j] / matrix[i][i]) : 0;
        }
    }

    double q = 0;
    for(int i = 0; i < N; i++) {
        double s = 0;
        for(int j = 0; j < N; j++) {
            s += fabs(alfa[i][j]);
        }
        if(s > fabs(q)) q = s;
    }

    double epsilon = fabs((1 - q) * eps / q);
    while (true) {
        for (int i = 0; i < N; i++) {
            x0[i] = x[i];

            double s = 0;
            for (int j = 0; j < N; j++) {
                s += alfa[i][j] * x[j];
            }
            x[i] = beta[i] + s;
        }


        double error = 0;
        for(int i = 0; i < N; i++) {
            double delta = fabs(x0[i] - x[i]);
            if (delta > error) error = delta;
        }
        if (error <= epsilon) break;
    }
    for (int i = 0; i < x.size(); i++) {
        cout << "x[" << (i + 1) << "] = " << x[i] << endl;
    }
}

int main() {
    Matrix<4, 5> AB = {{
        {14, 19, 15, 4, 180},
        {17, 33, 5, 10, 208},
        {11, 6, 28, 10, 230},
        { 6, 19, 3, 13, 149}
    }};

    Matrix<4, 5> T = {
        (AB[0] - AB[3]) * 2 - AB[2],
        AB[1],
        AB[2],
        AB[3] * 3 - AB[0] * 4,
    };


    single_division(AB);

    cout << endl;

    gauss_seidel(T, 1e-20);

}