gauss_solve.py

#----------------------------------------------------------------
# File:     gauss_solve.py
#----------------------------------------------------------------
#
# Author:   Marek Rychlik (rychlik@arizona.edu)
# Date:     Thu Sep 26 10:38:32 2024
# Copying:  (C) Marek Rychlik, 2020. All rights reserved.
# 
#----------------------------------------------------------------
# A Python wrapper module around the C library libgauss.so

import ctypes

gauss_library_path = './libgauss.so'

def unpack(A):
    """ Extract L and U parts from A, fill with 0's and 1's """
    n = len(A)
    L = [[A[i][j] for j in range(i)] + [1] + [0 for j in range(i+1, n)]
         for i in range(n)]

    U = [[0 for j in range(i)] + [A[i][j] for j in range(i, n)]
         for i in range(n)]

    return L, U

def lu_c(A):
    """ Accepts a list of lists A of floats and
    it returns (L, U) - the LU-decomposition as a tuple.
    """
    # Load the shared library
    lib = ctypes.CDLL(gauss_library_path)

    # Create a 2D array in Python and flatten it
    n = len(A)
    flat_array_2d = [item for row in A for item in row]

    # Convert to a ctypes array
    c_array_2d = (ctypes.c_double * len(flat_array_2d))(*flat_array_2d)

    # Define the function signature
    lib.lu_in_place.argtypes = (ctypes.c_int, ctypes.POINTER(ctypes.c_double))

    # Modify the array in C (e.g., add 10 to each element)
    lib.lu_in_place(n, c_array_2d)

    # Convert back to a 2D Python list of lists
    modified_array_2d = [
        [c_array_2d[i * n + j] for j in range(n)]
        for i in range(n)
    ]

    # Extract L and U parts from A, fill with 0's and 1's
    return unpack(modified_array_2d)

def lu_python(A):
    n = len(A)
    for k in range(n):
        for i in range(k,n):
            for j in range(k):
                A[k][i] -= A[k][j] * A[j][i]
        for i in range(k+1, n):
            for j in range(k):
                A[i][k] -= A[i][j] * A[j][k]
            A[i][k] /= A[k][k]

    return unpack(A)


def lu(A, use_c=False):
    if use_c:
        return lu_c(A)
    else:
        return lu_python(A)


def plu_c(A):
    """ Accepts a list of lists A of floats and
    it returns (P, L, U) - the LU-decomposition as a tuple.
    """
    # Load the shared library
    lib = ctypes.CDLL(gauss_library_path)

    # Create a 2D array in Python and flatten it
    n = len(A)
    flat_array_2d = [item for row in A for item in row]
    p_array = list(range(n))

    # Convert to a ctypes array
    c_array_2d = (ctypes.c_double * len(flat_array_2d))(*flat_array_2d)
    c_p_array = (ctypes.c_int * len(p_array))(*p_array)

    # Define the function signature
    lib.plu_in_place.argtypes = (ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_int))

    # Modify the array in C (e.g., add 10 to each element)
    lib.plu_in_place(n, c_array_2d, c_p_array)

    # Convert back to a 2D Python list of lists
    modified_array_2d = [
        [c_array_2d[i * n + j] for j in range(n)]
        for i in range(n)
    ]
    modified_p_array = [c_p_array[i] for i in range(n)]

    L, U = unpack(modified_array_2d)

    # Extract L and U parts from A, fill with 0's and 1's
    return modified_p_array, L, U


def plu_python(A):
    n = len(A)
    P = list(range(n))

    for k in range(n):
        # Find the pivot element
        pivot = A[k][k]
        r = k

        # Find pivot
        for i in range(k + 1, n):
            if abs(pivot) < abs(A[i][k]):
                pivot = A[i][k]
                r = i

        # row exchange
        if r != k:
            P[k], P[r] = P[r], P[k]
            # Swap rows k and r in L, only cols 1 to k-1
            for j in range(k):
                A[k][j], A[r][j] = A[r][j], A[k][j]

            # In this version, we swap the rest of the row, too
            for j in range(k, n):
                A[k][j], A[r][j] = A[r][j], A[k][j]

        for i in range(k+1, n):
            A[i][k] /= A[k][k]
            for j in range(k+1, n):
                A[i][j] -= A[i][k] * A[k][j]

    L, U = unpack(A)
    return P, L, U


def plu(A, use_c=False):
    if use_c:
        return plu_c(A)
    else:
        return plu_python(A)


if __name__ == "__main__":

    def get_A():
        """ Make a test matrix """
        A = [[2.0, 3.0, -1.0],
             [4.0, 1.0,  4.0],
             [-2.0, 7.0, 2.0]]
        return A

    # Must re-initialize A as it was destroyed
    A = get_A()

    L, U = lu(A,use_c=False)
    print("LU - Python version")
    print(L)
    print(U)

    # Must re-initialize A as it was destroyed
    A = get_A()


    L, U = lu(A, use_c=True)
    print("LU - C version")
    print(L)
    print(U)

    # Must re-initialize A as it was destroyed
    A = get_A()
    P, L, U = plu(A)
    print("PLU - Python  version")
    print(P)
    print(L)
    print(U)

    # Must re-initialize A as it was destroyed
    A = get_A()
    P, L, U = plu(A, use_c=True)
    print("PLU - C  version")
    print(P)
    print(L)
    print(U)