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;
}

3. Seidel Method

Time limit 1 second
Memory limit 64Mb
Input standard input or input.txt
Output standard output or output.txt
Write a computer program in C++ programming language to solve the given system of linear algebraic equations with the use of iterative Jacobi method.

Input format
The input contains:

A square matrix A (in element-wise manner with the dimension firstly) as in the previous exercises.
A vector of free coefficients b (in element-wise manner with the dimension firstly).
The approximation accuracy .

Output format
The output contains:

The string "The method is not applicable!"

or

Matrix , entitled "alpha:"
Vector , entitled "beta:"
Matrix B, entitled "B:"
Matrix C, entitled "C:"
Matrix I−B, entitled "I-B:"
Matrix (I−B)-1, entitled "(I-B)_-1:"

Sample 1
Input

3
10 1 1
2 10 1
2 2 10
3
12 13 14
0.01

Output

beta:
1.2000
1.3000
1.4000
alpha:
0.0000 -0.1000 -0.1000
-0.2000 0.0000 -0.1000
-0.2000 -0.2000 0.0000
B:
0.0000 0.0000 0.0000
-0.2000 0.0000 0.0000
-0.2000 -0.2000 0.0000
C:
0.0000 -0.1000 -0.1000
0.0000 0.0000 -0.1000
0.0000 0.0000 0.0000
I-B:
1.0000 0.0000 0.0000
0.2000 1.0000 0.0000
0.2000 0.2000 1.0000
(I-B)_-1:
1.0000 0.0000 0.0000
-0.2000 1.0000 0.0000
-0.1600 -0.2000 1.0000
x(0):
1.2000
1.3000
1.4000
e: 0.5694
x(1):
0.9300
0.9740
1.0192
e: 0.0770
x(2):
1.0007
0.9979
1.0003
e: 0.0021
x(3):
1.0002
0.9999
1.0000

Sample 2
Input

3
2 10 3
10 1 1
2 2 10
3
12 13 14
0.01

Output
The method is not applicable!

Notes
The type of the elements of the matrix is double. Print the result using the formatting to 4 digits after the floating point.

Hints

how to solve jacobi and seidel tasks in agla/assignment2 using matrix language in less than an hour so you can enjoy your time solving other assignments? (and no, this is not a research question)

Jacobi:
α = I - (D_1 * A)
β = D_1 * b
X0 = β
X(i+1) = X(i) + D_1 * (b - A*X(i) )

such that:
D : the diagonal matrix which contains the diagonal elements from A
D_1 : the inverse matrix of D
I : the identity matrix
A, b : they are given

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

such that:
D : the diagonal matrix which contains the diagonal elements from A
D_1 : the inverse matrix of D
B : the lower triangular part of A (with the diagonal)
C : the upper triangular part of A (without the diagonal)
B_1 : the inverse matrix of B

note : for the required B and C to be printed at first, they have to be taken from α instead of A (i really have no idea why)

finally, epsilon could be calculated by finding the norm of the vector X(i+1) - X(i)
and don't forget to check whether the method is applicable -> A have to be diagonally dominant (each element belongs to the diagonal in the absolute value must be greater than or equal to the absolute value of the sum of the other elements in the same row)

for better understanding for the above formulas and how they could be derived:
https://www.youtube.com/watch?v=VH0TZlkZPRo
https://www.youtube.com/watch?v=vN7fDiNxjss

Makefile

run: build
    ./main.o

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