linear-system-solver

About

This is a CPython library for solving linear systems of equations with the Jacobi iteration method and the Gauss-Seidel method. The library consists of two member functions that take in two lists A and b and finds the solution x to Ax = b. Given the nature of the algorithms, a solution can only be found if the matrix is diagonally dominant() with no diagonal zeroes. There seems to be a way to determine if a matrix is possibly diagonalizably dominant and transform it from one that is not initially, but for this project all inputs must be diagonally dominant. The algorithms are implemented in C and use CPython to wrap the functions into a python library.

Code

Build system

The build system is called with python build. pyproject.toml requires setuptools, which then calls setup.py, which instantiates the module as an Extension module with the provided source files.

Module definition

The module itself is defined with the PyMODINIT_FUNC definition, which returns PyModule_Create on the PyModuleDef. The methods are listed in the PyMethodDef, which wrap the source functions by parsing the arguments with PyArg_ParseTuple.

CPython and PyObjects

To call the functions that contain the algorithms, we can create PyObject arguments to use as parameters. First, the C code calls Py_Initialize to start the python interpreter. Then, the values can be built with PyList_New and Py_BuildValue. After, PyList_GetItem and PyFloat_AsDouble are called to retrieve the returned values. The PyObjects will be destructed dusing Py_Finalize.

This is then compiled with

gcc $(python3-config --cflags --embed --ldflags)

Algorithms

I implemented the algorithms as is on wikipedia, with the new estimates of x formulated using the iterative formulas. The A, b, and x used are heap allocated and are freed at the end of the function. Similar to the wikipedia implementation, I determined a convergence was reached when all values of x_i differ from the previous value by less than a given value.

To use

Build

make build

Test

make test

Call

source .venv/bin/activate; python3 source-file.py

# source-file.py
from solverlib import jacobi, gauss_seidel

A = [
    [ 15., 2., 0.0, 6.5, 4. ],
    [ 2., 21.5, -5., 3., 0.0 ],
    [ -9., 4., 18., 1., -2.4, ],
    [ 0.0, -4., -1., 16., 3. ],
    [ 5., 3.1, 6., 1., 17., ]
]

b = [
    -53.,
    103.,
    23.1,
    70.,
    -16.1
]

x = jacobi(A, b)
y = gauss_seidel(A, b)

print([round(a, 2) for a in x])
print([round(b, 2) for b in y])

# output -
# [ -6.5, 4.0, -3.0, 5.0, 1.0 ]
# [ -6.5, 4.0, -3.0, 5.0, 1.0 ]

Sources

Jacobi Method- Wikipedia

Gauss-Seidel Method- Wikipedia