lab3 oksana

using namespace std;

#include <iostream>
#include <cstdio>
#include <vector>
#include <cmath>

void PrintMatrix(vector<vector<double>> matrix) {
    for (int i = 0; i < matrix.size(); i++) {
        for (int j = 0; j < matrix[i].size(); j++) {
            printf("%12f", matrix[i][j]);
        }
        cout << endl;
    }
    cout << endl;
}

void DeleteRow(vector<vector<double>> &matrix, int line) {
    if (line >= matrix.size()) return;
    matrix.erase(matrix.begin() + line);
}

void DeleteColumn(vector<vector<double>> &matrix, int column) {
    if (column >= matrix[0].size()) return;
    for (int i = 0; i < matrix.size(); i++) {
        matrix[i].erase(matrix[i].begin() + column);
    }
}

void ChooseMainElement(vector<vector<double>> matrix) {
    cout << "Choose primary element method" << endl;
    PrintMatrix(matrix);

    vector<vector<double>> matrix1;
    vector<int> pos;

    int N = matrix.size();
    int M = matrix[0].size();

    for (int k = 0; k < M - 2; k++) {
        int p = 0, q = 0;

        double a = matrix[0][0];
        for (int i = 0; i < N - k; i++) {
            for (int j = 0; j < M - k - 1; j++) {
                if (fabs(matrix[i][j]) > fabs(a)) {
                    a = matrix[i][j];
                    p = i;
                    q = j;
                }
            }
        }

        vector<double> m;
        for (int i = 0; i < (N - k); i++) {
            m.push_back(-matrix[i][q] / a);
        }


        for (int i = 0; i < N - k; i++) {
            if (i == p) continue;

            for (int j = 0; j < M - k; j++) {
                matrix[i][j] += matrix[p][j] * m[i];
            }
        }
        pos.insert(pos.begin(), q);
        matrix1.insert(matrix1.begin(), matrix[p]);

        DeleteRow(matrix, p);
        DeleteColumn(matrix, q);
    }

    pos.insert(pos.begin(), 0);
    matrix1.insert(matrix1.begin(), matrix[0]);

    vector<double> X;
    for (int i = 0; i < matrix1.size(); i++) {
        double s = 0;
        for (int j = 0, index = 0; j < matrix1[i].size() - 1; j++) {
            if (j == pos[i]) continue;
            s += X[index++] * matrix1[i][j];
        }
        X.push_back((matrix1[i][matrix1[i].size() - 1] - s) / matrix1[i][pos[i]]);
    }



    printf("%f\t%f\t%f\t%f\n", X[0], X[1], X[3], X[2]);
}

void SimpleIteration(vector<vector<double>> matrix) {
    cout << endl << "Simple iteration method" << endl;
    PrintMatrix(matrix);

    int N = matrix.size();
    int M = matrix[0].size();

    double q = 0, result[N], x0[N];

    for (int i = 0; i < N; i++) {
        double a = matrix[i][i];
        matrix[i][M - 1] /= a;
        for (int j = 0; j < N; j++) {
            if (j != i) {
                matrix[i][j] /= a * -1;
            } else {
                matrix[i][j] = 0;
            }
        }
    }
    for (int i = 0; i < N; i++) {
        double a = 0;
        for (int j = 0; j < N; j++) {
            a += abs(matrix[i][j]);
        }
        if (q < a) {
            q = a;
        }
    }
    for (int i = 0; i < N; i++) {
        x0[i] = matrix[i][M - 1];
    }

    printf("-------------------------------------\n");
    printf("|%7s |%7s |%7s |%7s |\n", "x1", "x2", "x3", "x4");
    printf("-------------------------------------\n");

    float t = abs((1 - q) * 1e-2 / q);
    while (true) {
        printf("|");
        for (int i = 0; i < N; i++) {
            result[i] = matrix[i][M - 1];
            for (int j = 0; j < N; j++) {
                result[i] += matrix[i][j] * x0[j];
            }
            printf("%7.3f |", result[i]);
        }
        printf("\n-------------------------------------\n");
        double a = abs(result[0] - x0[0]);
        for (int i = 1; i < N; i++) {
            if (abs(result[i] - x0[i]) > a) {
                a = abs(result[i] - x0[i]);
            }
        }

        if (a <= t) break;

        for (int i = 0; i < N; i++) {
            x0[i] = result[i];
        }
    }
}

vector<double> operator-(const vector<double> &l, const vector<double> &r) {
    vector<double> result(l.size());
    for (int i = 0; i < l.size(); i++) {
        result[i] = l[i] - r[i];
    }
    return result;
}

vector<double> operator+(const vector<double> &l, const vector<double> &r) {
    vector<double> result(l.size());
    for (int i = 0; i < l.size(); i++) {
        result[i] = l[i] + r[i];
    }
    return result;
}

vector<double> operator*(const vector<double> &l, const double &r) {
    vector<double> result(l.size());
    for (int i = 0; i < l.size(); i++) {
        result[i] = l[i] * r;
    }
    return result;
}

int main() {
    vector<vector<double>> AB = {
        {  3, 19, 11,  8, 149 },
        {  9, 31,  3, 18, 257 },
        { 11,  7, 32, 13, 143 },
        { 12, 19, 12,  5, 144 }
    };
    ChooseMainElement(AB);

    vector<vector<double>> matrix = {
        AB[3] - AB[0],
        AB[0],
        AB[3] * -1 + AB[2],
        AB[1] + AB[2] - AB[3] - AB[0]
    };
    /*vector<vector<double>> matrix = {
        {  9,   0,  1, -3,  -5 }, //[3] - [0]
        {  3,  19, 11,  8, 149 }, //[0]
        { -1, -12, 20,  8,  -1 }, //[3] * -1 + [2]
        {  5,   0, 11, 18, 107 }  //[1] + [2] - [3] - [0] (мало бути 12 а не 11)
    };*/

    SimpleIteration(matrix);
}