-
Notifications
You must be signed in to change notification settings - Fork 0
/
gauss.py
90 lines (70 loc) · 2.22 KB
/
gauss.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
import numpy as np
from mpi4py import MPI
from utils import rank_from_row
def gauss(data, comm, n):
"""
Inverse rows using Gauss-Jordan elimination
:param data: list of rows, with row index as last element of row
:param comm: communicator
:param n: matrix dimension
:return: inversed rows
"""
def _eliminate(l, r, cur):
"""
Do gauss elimination step
:param l: left-side row
:param r: right-side row
:param cur: current processed index
:return: None
"""
for indx in indexes:
if indx == cur:
continue
rhs[indx] -= r * lhs[indx, cur]
lhs[indx] -= l * lhs[indx, cur]
def _send(indx, rnk):
"""
Send changed row to other process
:param indx: index of changed row
:param rnk: rank of process that change row
:return: None
"""
comm.Bcast([lhs[indx], MPI.DOUBLE], root=rnk)
comm.Bcast([rhs[indx], MPI.DOUBLE], root=rnk)
_eliminate(lhs[indx], rhs[indx], indx)
def _receive(indx, rnk):
"""
Receive changed row from current active process
:param indx: index of changed row
:param rnk: rank of process that change row
:return: None
"""
comm.Bcast([l_row, MPI.DOUBLE], root=rnk)
comm.Bcast([r_row, MPI.DOUBLE], root=rnk)
_eliminate(l_row, r_row, indx)
indexes = []
lhs = np.zeros((n, n))
rhs = np.identity(n)
for line in data:
ind = int(line[-1])
indexes.append(ind)
lhs[ind] = line[:-1]
l_row = np.zeros(n, dtype=np.float64)
r_row = np.zeros(n, dtype=np.float64)
for i in range(n): # forward elimination
rank = rank_from_row(n, comm.size, i)
if i in indexes:
rhs[i] /= lhs[i, i]
lhs[i] /= lhs[i, i]
_send(i, rank)
else:
_receive(i, rank)
comm.Barrier()
for i in range(n - 1, -1, -1): # back substitution
rank = rank_from_row(n, comm.size, i)
if i in indexes:
_send(i, rank)
else:
_receive(i, rank)
comm.Barrier()
return [rhs.tolist()[ind] [ind] for ind in indexes]