Gaussian_elimination/Gauss.py at master · Megaco/Gaussian_elimination · GitHub

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
__author__ = 'Megaco'

def pprint(A):
    n = len(A)
    for i in range(0, n):
        line = ""
        for j in range(0, n + 1):
            line += str(A[i][j]) + "\t"
            if j == n - 1:
                line += "| "
        print(line)
    print("")


def gauss(A):
    n = len(A)

    for i in range(0, n):
        # Search for maximum in this column
        maxEl = abs(A[i][i])
        maxRow = i
        for k in range(i + 1, n):
            if abs(A[k][i]) > maxEl:
                maxEl = abs(A[k][i])
                maxRow = k

        # Swap maximum row with current row (column by column)
        for k in range(i, n + 1):
            tmp = A[maxRow][k]
            A[maxRow][k] = A[i][k]
            A[i][k] = tmp

        # Make all rows below this one 0 in current column
        for k in range(i + 1, n):
            try:
                c = -A[k][i] / A[i][i]
            except ZeroDivisionError as err:
                print("NO")
                exit(0)
            for j in range(i, n + 1):
                if i == j:
                    A[k][j] = 0
                else:
                    A[k][j] += c * A[i][j]

    # Solve equation Ax=b for an upper triangular matrix A
    x = [0 for i in range(n)]
    for i in range(n - 1, -1, -1):
        try:
            x[i] = A[i][n] / A[i][i]
        except ZeroDivisionError as err:
            print("NO")
            exit(0)
        for k in range(i - 1, -1, -1):
            A[k][n] -= A[k][i] * x[i]
    return x


if __name__ == "__main__":
    from fractions import Fraction

    n, neiz = input().strip().split(' ')
    n = int(n)
    neiz = int(neiz)
    if n < neiz:
        print("INF")
        exit(0)

    # n = int(input())
    # n=int(n)

    A = [[0 for j in range(neiz+1)] for i in range(n)]
    # Read input data
    for i in range(0, n):
        line = map(Fraction, input().strip().split(" "))

        for j, el in enumerate(line):
            A[i][j] = el

# input()
# line = input().split(" ")
# lastLine = map(Fraction, line)
# for i in range(0, n):
#     A[i][n] = lastLine[i]

# Print input
# pprint(A)

# Calculate solution
x = gauss(A)

# Print result
# line = "Result:\t"


print("YES")
line = ""
x = [float(x) for x in x]
for i in range(0, n):
    line += str(x[i]) + " "
print(line)