-
Notifications
You must be signed in to change notification settings - Fork 1
/
gauss_jordan.py
88 lines (71 loc) · 2.24 KB
/
gauss_jordan.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
import fractions
import numpy as np
def show_matrix(B):
matrix_print = []
for i in range(B.shape[0]):
row = "["
for j in range(B.shape[1]):
if B[i, j].denominator == 1:
row += "%2d " % (B[i, j].numerator,)
else:
row += "%2d/%2d " % (B[i, j].numerator, B[i, j].denominator)
row += "]"
matrix_print.append(row)
return matrix_print
def show(A, B):
print("")
print_A = show_matrix(A)
print_B = show_matrix(B)
for row_A, row_B in zip(print_A, print_B):
print(row_A + " " + row_B)
def map_map(f, xss):
return [list(map(f, xs)) for xs in xss]
A = np.array(map_map(fractions.Fraction,
[[6,3,4,5],
[1,2,2,1],
[2,4,3,2],
[3,3,4,2]]), dtype=fractions.Fraction)
def gauss_jordan(A):
A = A.copy()
num_rows, num_columns = A.shape
# identity matrix with the same dimension as A
E = np.array([[fractions.Fraction(1 if i == j else 0) for i in range(num_rows)] for j in range(num_columns)])
# a "DSL" to invert the matrix
def at(i, j):
return A[i, j]
def sub(i, j):
A[i, :] -= A[j, :]
E[i, :] -= E[j, :]
def div(i, x):
A[i, :] /= x
E[i, :] /= x
def swap_row(i, j):
A[i, :], A[j, :] = A[j, :], A[i, :]
E[i, :], E[j, :] = E[j, :], E[i, :]
def get_non_zero_at_throu_row_swapping(index):
for i in range(index, num_rows):
if at(i, index) != 0:
break
else:
raise ValueError("SINGULAR MATRIX")
swap_row(index, i)
# convert the matrix in a upper triangular matrix
for i in range(0, num_rows):
show(A, E)
get_non_zero_at_throu_row_swapping(i)
div(i, at(i, i))
for j in range(i + 1, num_rows):
if at(j, i) != 0:
div(j, at(j, i))
sub(j, i)
# convert the upper triangular matrix into a identity matrix
for i in range(num_rows - 1, -1, -1):
show(A, E)
div(i, at(i, i))
for j in range(i - 1, -1, -1):
if at(j, i) != 0:
div(j, at(j, i))
sub(j, i)
show(A, E)
return E
gauss_jordan(A)