Gaussian-Seidel method to solve a system of linear equations. [Innopolis University] Linear Algebra 2023

main.cpp

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <stdlib.h>

using namespace std;

typedef long long ll;

struct Matrix;
struct Vector;


/* Square matrix */
struct Matrix
{
    ll size;
    vector<vector<double>> matrix;
};

struct Matrix* matrix_init();
struct Matrix* matrix_init(ll size);
struct Matrix* matrix_init(ll size, const vector<vector<double>> &);
struct Matrix* matrix_init_ident(ll size);
struct Matrix* matrix_copy(const Matrix *);
struct Matrix* matrix_make_diag(const Matrix *);
struct Matrix* matrix_make_ltr(const Matrix *); /* lower triangular */
struct Matrix* matrix_make_utr(const Matrix *); /* upper triangular */
struct Matrix* matrix_inv_diag(const Matrix *);
struct Matrix* matrix_inv_ltr(const Matrix *);
struct Matrix* matrix_sum_matrix(const Matrix *, const Matrix *);
struct Matrix* matrix_dif_matrix(const Matrix *, const Matrix *);
struct Matrix* matrix_mul_matrix(const Matrix *, const Matrix *);
struct Vector* matrix_mul_vector(const Matrix *, const Vector *);
struct Matrix* matrix_mul_scalar(Matrix *, double);
bool matrix_is_diagonally_dominant(Matrix* );
void matrix_print(const Matrix *);


/* 1-dimensional vector */
struct Vector
{
    ll dimension;
    vector<double> values;
};

struct Vector* vector_init();
struct Vector* vector_init(ll dim);
struct Vector* vector_init(ll dim, const vector<double> &);
struct Vector* vector_copy(const Vector *);
struct Vector* vector_dif_vector(const Vector *, const Vector *);
double vector_norm(const Vector*);
void vector_print(const Vector *);
void vector_vprint(const Vector *);

/****************/
/* MATRIX START */
/****************/
struct Matrix* matrix_init()
{
    struct Matrix* matrix = (struct Matrix*)malloc(sizeof(struct Matrix));
    return matrix;
}

struct Matrix* matrix_init(ll size)
{
    struct Matrix* matrix = (struct Matrix*)malloc(sizeof(struct Matrix));
    matrix->matrix = vector<vector<double>>(size);
    for (ll i = 0; i < size; ++i)
    {
        matrix->matrix[i] = vector<double>(size);
        for (ll j = 0; j < size; ++j)
        {
            matrix->matrix[i][j] = 0;
        }
    }
    matrix->size = size;

    return matrix;
}

struct Matrix* matrix_init(ll size, const vector<vector<double>> & vec2d)
{
    struct Matrix* matrix = (struct Matrix*)malloc(sizeof(struct Matrix));
    matrix->matrix = vector<vector<double>>(size);
    for (ll i = 0; i < size; ++i)
    {
        matrix->matrix[i] = vector<double>(size);
        for (ll j = 0; j < size; ++j)
        {
            matrix->matrix[i][j] = vec2d[i][j];
        }
    }
    matrix->size = size;

    return matrix;
}

struct Matrix* matrix_init_ident(ll size)
{
    struct Matrix* matrix = matrix_init(size);
    for (ll i = 0; i < matrix->size; ++i)
    {
        matrix->matrix[i][i] = 1.;
    }

    return matrix;
}

struct Matrix* matrix_copy(const Matrix* origin)
{
    struct Matrix* matrix = matrix_init();
    matrix->size = origin->size;
    matrix->matrix = origin->matrix;
    return matrix;
}

struct Matrix* matrix_make_diag(const Matrix* origin)
{
    ll size = origin->size;
    struct Matrix* matrix = matrix_init(size);
    for (ll i = 0; i < size; ++i)
    {
        matrix->matrix[i][i] = origin->matrix[i][i];
    }

    return matrix;
}

struct Matrix* matrix_make_ltr(const Matrix* origin)
{
    struct Matrix* ltr = matrix_init(origin->size);

    for (ll i = 0; i < origin->size; ++i)
    {
        for (ll j = 0; j < origin->size; ++j)
        {
            if (i >= j)
            {
                ltr->matrix[i][j] = origin->matrix[i][j];
            }
        }
    }

    return ltr;
} 

struct Matrix* matrix_make_utr(const Matrix* origin)
{
    struct Matrix* utr = matrix_init(origin->size);

    for (ll i = 0; i < origin->size; ++i)
    {
        for (ll j = 0; j < origin->size; ++j)
        {
            if (i < j)
            {
                utr->matrix[i][j] = origin->matrix[i][j];
            }
        }
    }

    return utr;
}

struct Matrix* matrix_inv_diag(const Matrix* diag)
{
    struct Matrix* matrix = matrix_init(diag->size);
    matrix->matrix = diag->matrix; /* copy values */

    for (ll i = 0; i < matrix->size; ++i)
    {
        for (ll j = 0; j < matrix->size; ++j)
        {
            if (i == j)
            {
                double el = matrix->matrix[i][j];
                matrix->matrix[i][j] = 1. / el;
            }
        }
    }

    return matrix;
}

/* Using Gaussian elimination */
struct Matrix* matrix_inv_ltr(const Matrix* ltr)
{
    ll size = ltr->size;
    struct Matrix* res = matrix_init_ident(ltr->size);

    for (ll i = 0; i < size; ++i)
    {
        res->matrix[i][i] = 1. / ltr->matrix[i][i];
        for (ll j = i + 1; j < size; ++j)
        {
            double s = 0;
            for (ll k = i; k < j; ++k)
            {
                s -= ltr->matrix[j][k] * res->matrix[k][i];
            }
            res->matrix[j][i] = s / ltr->matrix[j][j];
        }
    }

    return res;
}

struct Matrix* matrix_dif_matrix(const Matrix* lhs, const Matrix* rhs)
{
    struct Matrix* res = matrix_init(lhs->size);

    for (ll i = 0; i < res->size; ++i)
    {
        for (ll j = 0; j < res->size; ++j)
        {
            double c = lhs->matrix[i][j] - rhs->matrix[i][j];
            res->matrix[i][j] = c;
        }
    }

    return res;
}

struct Matrix* matrix_mul_matrix(const Matrix* lhs, const Matrix* rhs)
{
    struct Matrix* res = matrix_init(lhs->size);
    for (ll r1 = 0; r1 < res->size; ++r1)
    {
        for (ll c2 = 0; c2 < res->size; ++c2)
        {
            for (ll c1 = 0; c1 < res->size; ++c1)
            {
                double c = lhs->matrix[r1][c1] * rhs->matrix[c1][c2];
                res->matrix[r1][c2] += c;
            }
        }
    }

    return res;
}

struct Vector* matrix_mul_vector(const Matrix* matrix, const Vector* vec)
{
    struct Vector* res = vector_init(vec->dimension); /* NxN x Nx1 => Nx1, cause squared */
    for (ll i = 0; i < matrix->size; ++i)
    {
        for (ll j = 0; j < vec->dimension; ++j)
        {
            double c = matrix->matrix[i][j] * vec->values[j];
            res->values[i] += c;
        }
    }

    return res;
}

bool matrix_is_diagonally_dominant(Matrix* matrix)
{
    for (ll i = 0; i < matrix->size; ++i)
    {
        for (ll j = 0; j < matrix->size; ++j)
        {
            double d = matrix->matrix[i][j];
            if (i == j)
            {
                double s = 0.;
                for (ll k = 0; k < matrix->size; ++k)
                {
                    if (k != j)
                    {
                        s += matrix->matrix[i][k];
                    }
                }

                if (d < s)
                {
                    return false;
                }
            }
        }
    }

    return true;
}

void matrix_print(const Matrix* matrix)
{
    for (ll i = 0; i < matrix->size; ++i)
    {
        for (ll j = 0; j < matrix->size; ++j)
        {
            cout << matrix->matrix[i][j] << " ";
        }
        cout << endl;
    }
}
/**************/
/* MATRIX END */
/**************/

/****************/
/* VECTOR START */
/****************/
struct Vector* vector_init()
{
    struct Vector* vec = (struct Vector*)malloc(sizeof(struct Vector));
    return vec;
}

struct Vector* vector_init(ll dim)
{
    struct Vector* vec = vector_init();
    vec->values = vector<double>(dim);
    vec->dimension = dim;
    return vec;
}

struct Vector* vector_init(ll dim, const vector<double>& vec)
{
    struct Vector* res = vector_init();
    res->dimension = dim;
    res->values = vector<double>(dim);
    for (ll i = 0; i < dim; ++i)
    {
        res->values[i] = vec[i];
    }

    return res;
}

struct Vector* vector_copy(const Vector* vec)
{
    struct Vector* res = vector_init();
    res->dimension = vec->dimension;
    res->values = vec->values;
    return res;
}

struct Vector* vector_dif_vector(const Vector* lhs, const Vector* rhs)
{
    struct Vector* vec = vector_init(lhs->dimension);

    for (ll i = 0; i < lhs->dimension; ++i)
    {
        vec->values[i] = lhs->values[i] - rhs->values[i];
    }

    return vec;
}

double vector_norm(const Vector* vec)
{
    double res = 0.;
    for (ll i = 0; i < vec->dimension; ++i)
    {
        res += vec->values[i] * vec->values[i];
    }

    return sqrt(res);
}

void vector_print(const Vector* vec)
{
    for (ll i = 0; i < vec->dimension; ++i)
    {
        cout << vec->values[i] << " ";
    }
    cout << endl;
}

void vector_vprint(const Vector* vec)
{
    for (ll i = 0; i < vec->dimension; ++i)
    {
        cout << vec->values[i] << endl;
    }
}
/**************/
/* VECTOR END */
/**************/

void print_vec(const vector<double>& vec)
{
    for (ll i = 0; i < vec.size(); ++i)
    {
        cout << vec[i] << " ";
    }
}

void print_vec2d(const vector<vector<double>>& vec)
{
    for (ll i = 0; i < vec.size(); ++i)
    {
        for (ll j = 0; j < vec.size(); ++j)
        {
            cout << vec[i][j] << " ";
        }
        cout << endl;
    }
}

// α = I - (D_1 * A)
// β = D_1 * b
// X0 = β
// X(i+1) = B_1 * (b - C*X(i) )
// A, b : they are given


int main(int argc, char **argv)
{
    /* Gauss–Seidel method
     * 
     * also known as the Liebmann method or 
     * the method of successive displacement, 
     * is an iterative method used to solve 
     * a system of linear equations.
     * https://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method
     */

    ll size; cin >> size;
    vector<vector<double>> m_values(size);
    for (ll i = 0; i < size; ++i)
    {
        m_values[i] = vector<double>(size);
        for (ll j = 0; j < size; ++j)
        {
            cin >> m_values[i][j];
        }
    }
    
    ll dim; cin >> dim;
    vector<double> v_values(dim);
    for (ll i = 0; i < dim; ++i)
    {
        cin >> v_values[i];
    }
    double accur; cin >> accur;


    cout << setprecision(4) << fixed;

    struct Matrix* A = matrix_init(size, m_values);
    struct Matrix* I = matrix_init_ident(size);
    struct Matrix* D = matrix_make_diag(A);
    struct Matrix* D_1 = matrix_inv_diag(D);
    struct Vector* b = vector_init(dim, v_values);

    struct Matrix* D_1_mul_A = matrix_mul_matrix(D_1, A);
    struct Matrix* alpha = matrix_dif_matrix(I, D_1_mul_A);

    struct Matrix* B = matrix_make_ltr(A);
    struct Matrix* B_1 = matrix_inv_ltr(B);
    struct Matrix* C = matrix_make_utr(A);

    struct Matrix* aB = matrix_make_ltr(alpha);
    struct Matrix* aC = matrix_make_utr(alpha);
    struct Matrix* I_aB = matrix_dif_matrix(I, aB);
    struct Matrix* I_aB_1 = matrix_inv_ltr(I_aB);
    struct Vector* beta = matrix_mul_vector(D_1, b);
    
    if (!matrix_is_diagonally_dominant(A))
    {
        cout << "The method is not applicable!" << endl;
        return 0;
    }

    cout << "beta:" << endl;
    vector_vprint(beta);
    cout << "alpha:" << endl;
    matrix_print(alpha);
    cout << "B:" << endl;
    matrix_print(aB);
    cout << "C:" << endl;
    matrix_print(aC);
    cout << "I-B:" << endl;
    matrix_print(I_aB);
    cout << "(I-B)_-1:" << endl;
    matrix_print(I_aB_1);

    double eps = 1.;
    struct Vector* X_prev = vector_init(dim);
    struct Vector* X_curr = vector_copy(beta);
    struct Vector* X_dif = vector_dif_vector(X_curr, X_prev);
    struct Vector* C_mul_X_curr = matrix_mul_vector(C, X_curr);
    struct Vector* b_dif_CX_curr = vector_dif_vector(b, C_mul_X_curr);
    ll idx = 0;
    for (idx = 0; eps > accur; ++idx)
    {
        if (eps <= accur)
        {
            break;
        }

        cout << "x(" << idx << "):" << endl;
        vector_vprint(X_curr);
    
        X_prev = X_curr;

        C_mul_X_curr = matrix_mul_vector(C, X_curr);
        b_dif_CX_curr = vector_dif_vector(b, C_mul_X_curr);
        X_curr = matrix_mul_vector(B_1, b_dif_CX_curr);

        X_dif = vector_dif_vector(X_curr, X_prev); 
        eps = vector_norm(X_dif);
        cout << "e: " << eps << endl;
    }
    cout << "x(" << idx << "):" << endl;
    vector_vprint(X_curr);

    return 0;
}

Makefile

run: build
    ./main.o

build:
    g++ main.cpp -o main.o -std=c++11 -O0