I finally found, that it can be done using **LU decomposition**. Here the **U** matrix represents the reduced form of the linear system.

```
from numpy import array
from scipy.linalg import lu
a = array([[2.,4.,4.,4.],[1.,2.,3.,3.],[1.,2.,2.,2.],[1.,4.,3.,4.]])
pl, u = lu(a, permute_l=True)
```

Then `u`

reads

```
array([[ 2., 4., 4., 4.],
[ 0., 2., 1., 2.],
[ 0., 0., 1., 1.],
[ 0., 0., 0., 0.]])
```

Depending on the solvability of the system this matrix has an upper triangular or trapezoidal structure. In the above case a line of zeros arises, as the matrix has only rank `3`

.