lab3.cpp

#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <vector>

using namespace std;

typedef vector<double> Vector;
typedef vector<Vector> Martix;

void print(Martix matrix, int n, int m) {
    int i, j;
    for (i = 0; i < n; i++){
        for (j = 0; j < m; j++){
            printf("%8.2f", matrix [i][j]);
        }
        printf("\n");
    }
    return;
}

void CompletePivoting(Martix matrix, int n, int m) {
    printf("Complete pivoting method\n");
    print(matrix, n, m);

    int p, q, line_new = 0;
    double M[n], t[n][m], x[n];

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

        for (int j = 0; j < m; j++) {
            t[line_new][j] = matrix[p][j];
        }
        line_new++;

        for (int i = 0; i < n; i++) {
            M[i] = -matrix[i][q] / a;
        }

        for (int i = 0; i < n; i++) {
            for (int j = 0; j < m; j++) {
                if (i != p) {
                    matrix[i][j] += matrix[p][j] * M[i];
                }
            }
        }

        for (int j = 0; j < m; j++) {
            matrix[p][j] = 0;
        }

        for (int i = 0; i < n; i++) {
            matrix[i][q] = 0;
        }

    }

    for (int i = 0; i < n; i++) {
        if (matrix[i][m - 1] != 0) {
            p = i;
        }
    }

    for (int j = 0; j < m; j++) {
        t[line_new][j] = matrix[p][j];
    }

    for (int j = 0; j < m-1; j++) {
        x[j] = 0;
    }

    for (int i = n - 1; i >= 0; i--){
        double R = t[i][m - 1];
        for (int j = 0; j < m - 1; j++) {
            if ((t[i][j] != 0) && (x[j] != 0)) {
                R -= t[i][j] * x[j];
            }
        }
        for (int j = 0; j < m - 1; j++) {
            if ((t[i][j] != 0) && (x[j] == 0)) {
                x[j] = R / t[i][j];
            }
        }
    }

    for (int i = 0; i < n; i++) {
        printf("x%i = %f\n", i, x[i]);
    }
}

double* SimpleIteration(Martix matrix, int n, int m) {
    printf("Simple iteration method\n");
    print(matrix, n, m);

    double x[n], x0[n];
    for (int i = 0; i < n; i++) {
        double a = 1.0 / matrix[i][i];
        for (int j = 0; j < m; j++) {
            matrix[i][j] = matrix[i][j] * a;
        }
    }

    double q = 0;
    for (int i = 0; i < n; i++) {
        double s = 0;
        for (int j = 0; j < m - 1; j++) {
            if (i != j) {
                s += fabs(matrix[i][j]);
            }
        }
        if (s > q) q = s;
    }

    for (int i = 0; i < n; i++) {
        x0[i] = matrix[i][m-1];
    }

    for (int i = 0; i < n; i++) {
        x[i] = matrix[i][m-1];
        for (int j = 0; j < m - 1; j++) {
            if (i != j) {
                x[i] -= matrix[i][j] * x0[j];
            }
        }
    }

    for (int i = 0; i < n; i++) {
        x0[i] = x[i];
    }

    int k = 2;
    double m_norma_x = 0;
    double error = 1e-8 * fabs(((1 - q) / q));
    do{
        for (int i = 0; i < n; i++) {
            x[i] = matrix[i][m - 1];
            for (int j = 0; j < m - 1; j++) {
                if (i != j) {
                    x[i] -= matrix[i][j] * x0[j];
                }
            }
        }

        for (int i = 0; i < n; i++)
            x0[i] = x[i] - x0[i];

        m_norma_x = fabs(x0[0]);
        for (int i = 1; i < n; i++)
            if (fabs(x0[i]) > m_norma_x)
                m_norma_x = fabs(x0[i]);

        for (int i = 0; i < n; i++) {
            x0[i] = x[i];
        }

        k++;
    } while (m_norma_x > error);

    for (int i = 0; i < n; i++) {
        printf("x%i = %f\n", i, x[i]);
    }
}

Vector operator+(const Vector &a, const Vector &b) {
    Vector r(a.size());
    for (int i = 0; i < a.size(); i++) {
        r[i] = a[i] + b[i];
    }
    return r;
}

Vector operator-(const Vector &a, const Vector &b) {
    Vector r(a.size());
    for (int i = 0; i < a.size(); i++) {
        r[i] = a[i] - b[i];
    }
    return r;
}

Vector operator*(const Vector &a, const double &b) {
    Vector r(a.size());
    for (int i = 0; i < a.size(); i++) {
        r[i] = a[i] * b;
    }
    return r;
}

int main() {
    Martix matrix1 = {
        {13, 14, 17, 14, 146},
        { 9, 26,  7,  9, 135},
        { 8,  4, 25, 11, 119},
        {15, 11,  5, 16, 109}
    };

    Martix matrix2 = {
        matrix1[0] + matrix1[3] - matrix1[2] - matrix1[1],
        matrix1[1],
        matrix1[0],
        matrix1[3],
    };

    CompletePivoting(matrix1, 4, 5);
    printf("\n");
    SimpleIteration(matrix2, 4, 5);
}