forked from MinMa1122/GaussianElimination
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgauss_solve.py
executable file
·143 lines (113 loc) · 4.2 KB
/
gauss_solve.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#----------------------------------------------------------------
# File: gauss_solve.py
#----------------------------------------------------------------
#
# Author: Marek Rychlik ([email protected])
# 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
import numpy as np
# Load the shared library
lib = ctypes.CDLL('./libgauss.so')
def lu(A, use_c=False):
""" Accepts a list of lists A of floats and
it returns (L, U) - the LU-decomposition as a tuple.
"""
n = len(A)
if use_c:
# Create a 2D array in Python and flatten it
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
L = [
[modified_array_2d[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) ] + [ modified_array_2d[i][j] for j in range(i, n) ]
for i in range(n)
]
else:
# Initialize matrices
L = np.zeros((n,n))
U = np.zeros((n,n))
for k in range(n):
for i in range(k, n):
U[k, i] = A[k][i] - np.dot(L[k,:k],U[:k,i])
for i in range(k+1, n):
L[i, k] = (A[i][k] - np.dot(L[i,:k],U[:k,k]))/U[k, k]
L[k, k] = 1
L = L.tolist()
U = U.tolist()
return L, U
def plu(A, use_c=False):
"""
Accepts list of lists A of floats and returns
P: list of integers (permutation)
L: lower triangular list of lists
U: upper triangular list of lists
"""
n = len(A)
if use_c:
flat_array_2d = [item for row in A for item in row]
A_c_version = (ctypes.c_double * len(flat_array_2d))(*flat_array_2d)
P_c_version = (ctypes.c_int * n)()
# Define function signature
lib.plu.argtypes = (ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_int))
# Do PA=LU in C to c_array_2d
lib.plu(n, A_c_version, P_c_version)
# Convert A back to a 2D Python list of lists
modified_array_2d = [
[A_c_version[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
L = [
[modified_array_2d[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) ] + [ modified_array_2d[i][j] for j in range(i, n) ]
for i in range(n)
]
P = list(P_c_version)
# Use Python version
else:
# Initialize matrices
P_np = np.linspace(0, n-1, n)
L = np.zeros(shape=(n,n))
U = np.array(A)
for k in range(n-1):
# Find the row index r of the pivot element
r = np.argmax(np.abs(U[k:, k])) + k
# Make row/element swaps
U[[k, r]] = U[[r, k]]
P_np[r], P_np[k] = P_np[k], P_np[r]
L[[k, r], 0:k] = L[[r, k], 0:k]
for i in range(k+1, n):
L[i, k] = (U[i, k] / U[k, k])
U[i] = U[i] - L[i, k] * U[k]
# Add ones to diagonal of L and convert matrices to correct return form
L = (L + np.eye(n)).tolist()
P = P_np.astype(int).tolist()
U = U.tolist()
return P, L, U
if __name__ == "__main__":
A = [[4.0, 9.0, 10.0],
[14.0, 30.0, 34.0],
[2.0, 3.0, 3.0]];
L1, U1 = lu(A, use_c=False)
P, L, U = plu(A, use_c=True)