Some examples of matrices operations:

A is: [[1, 2], [3, 4]]
B is: [[5, 6], [7, 8]]

type(A) is: <class 'matrix.Matrix'>
A.__class__ is: <class 'matrix.Matrix'>

A[0, 1] is: 2
B[1, 1] is: 8

A[0, 1] becomes 0.
A[0, 1] is now: 0
A is now: [[1, 0], [3, 4]]

A == A is: True
A == B is: False

A + B is: [[6, 6], [10, 12]]
B + A is: [[6, 6], [10, 12]]

A + 3 is: [[4, 3], [6, 7]]
1 + B is: [[6, 7], [8, 9]]

A += 1 can be done
A is now: [[2, 1], [4, 5]]

A - B is: [[-3, -5], [-3, -3]]
B - A is: [[3, 5], [3, 3]]

-A is: [[-2, -1], [-4, -5]]
B + (-A) is: [[3, 5], [3, 3]]

+A is: [[2, 1], [4, 5]]

A * B is: [[17, 20], [55, 64]]
B * A is: [[34, 35], [46, 47]]

A * 3 is: [[6, 3], [12, 15]]
2 * B is: [[10, 12], [14, 16]]

A is: [[2, 1], [4, 5]]

A ** 3 is: [[44, 43], [172, 173]]
B ** 91 is: [[260886190336827603639069235910825462911444049857238777869408910350902285458764592605977965777341722785, 303823113391902795973729664639518907663320331121475544074064063760904586750217533165849296350745720242], [354460298957219928636017942079438725607207052975054801419741407721055351208587122026824179075870006949, 412797747032779001625934068230584916743104215417976549906440942231354578833873359188902613952714582906]]

THETA is: [[3, 0], [0, -1]]

THETA.exp() is: [[20.0855369232, 0.0], [0.0, 0.367879441171]]
Conversely diag([exp(3), exp(-1)]) is: [[20.0855369232, 0], [0, 0.367879441171]]

exp(diag([3, -1])) ~= diag([exp(3), exp(-1)]): True

With epsilon = 1e-16:
exp(diag([3, -1])) ~= diag([exp(3), exp(-1)]): False
exp(diag([3, -1])) == diag([exp(3), exp(-1)]): True

A is: [[2, 1], [4, 5]]

exp(A) is: [[82.8603841611, 80.1421023326], [320.568409331, 323.286691159]]

A - 3 is: [[-1, -2], [1, 2]]
1 - B is: [[-4, -5], [-6, -7]]

A is: [[2, 1], [4, 5]] and B is: [[5, 6], [7, 8]]
A.multiply_elementwise(B) is: [[10, 6], [28, 40]]
B.multiply_elementwise(A) is: [[10, 6], [28, 40]]

A.count(1) is: 1
eye(4).count(1) is: 4

1 in A is: True
0 not in B is: True

I = eye(3) is: [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
I.n is: 3
I.m is: 3
len(I) is: 9
len(I) is: (3, 3)
I.count(0) is: 6

diag([1, 2, 3, 4]) is: [[1, 0, 0, 0], [0, 2, 0, 0], [0, 0, 3, 0], [0, 0, 0, 4]]

O = ones(3, 4) is: [[1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]]
I * O is: [[1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]]

diag([1+1j, 2+2j, 3+3j, 4+4j]) is: [[(1+1j), 0, 0, 0], [0, (2+2j), 0, 0], [0, 0, (3+3j), 0], [0, 0, 0, (4+4j)]]
D.real is: [[1.0, 0, 0, 0], [0, 2.0, 0, 0], [0, 0, 3.0, 0], [0, 0, 0, 4.0]]
D.imag is: [[1.0, 0, 0, 0], [0, 2.0, 0, 0], [0, 0, 3.0, 0], [0, 0, 0, 4.0]]
D.conjugate() is: [[(1-1j), 0, 0, 0], [0, (2-2j), 0, 0], [0, 0, (3-3j), 0], [0, 0, 0, (4-4j)]]

A is: [[2, 1], [4, 5]]
A.T is: [[2, 4], [1, 5]]

A is: [[2, 1], [4, 5]]
A.col(0) is: [[2], [4]]
A.col(1) is: [[1], [5]]
A.row(0) is: [[2, 1]]
A.row(1) is: [[4, 5]]

A is: [[2, 1], [4, 5]]
For f: x -> x² + x + 1, A.map(f) is: [[7, 3], [21, 31]]

A is: [[2, 1], [4, 5]]
Iterating over values of a matrix: 
for v in A: print v
2
1
4
5

B is: [[5, 6], [7, 8]]
list(B) is: [5, 6, 7, 8]
tuple(B) is: (5, 6, 7, 8)

U is: [[2, 3], [7, -1]]
2×U² + 4×U + I_n = [[59, 18], [42, 41]]

A random matrix of size (3, 2) and integer coefficients between -5 and 5 is: [[4, 5], [0, 1], [2, -2]]

A random matrix of size (2, 4) and float coefficients between -2 and 2 is: [[1.37269238284, 1.29686098896, -1.99040629531, 0.780250833011], [0.919050216176, -1.05031431628, 1.74817738252, 0.554349698283]]

A == -(-A) is: True
A is -(-A) is: False

A1 is: [[1, 1], [1, 0]]
B1 is: [[0, 0], [1, 0]]
A1 and B1 is: [[0, 0], [1, 0]]
A1 or B1 is: [[1, 1], [1, 0]]
not A1 is: False  ==> WARNING is not the negation element-wise!

WARNING: 'or' and 'and' work weirdly with numbers

A2 is: [[7, 7], [7, 0]]
B2 is: [[0, 0], [-6, 0]]
A2 and B2 is: [[0, 0], [-6, 0]]
A2 or B2 is: [[7, 7], [7, 0]]

A is: [[2, 1], [4, 5]]
A % 3 is: [[2, 1], [1, 2]]
A % 2 is: [[0, 1], [0, 1]]

B is: [[5, 6], [7, 8]]
B % 4 is: [[1, 2], [3, 0]]
B % 2 is: [[1, 0], [1, 0]]

ones(4) is: [[1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1], [1, 1, 1, 1]]
ones(2, 3) is: [[1, 1, 1], [1, 1, 1]]

zeros(4) is: [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]
zeros(2, 3) is: [[0, 0, 0], [0, 0, 0]]

Z is: [[-2, 7], [-10, 6], [-10, -4], [-9, -10]]
Vector v is: [1, 1]
Z.dot(v) is: [[5], [-4], [-14], [-19]]

A is: [[2, 1], [4, 5]]
A.norm() is: 6.78232998313
A.norm(p=1) is: 12 which is like sum(A): 12
A.norm(3) is: 5.82847668325
norm(A) is: 6.78232998313
norm(A, 42) is: 5.00001012703

Random U is: [[0, 14, 7, -13, 20], [12, -16, 15, 8, -3], [11, 15, -5, -14, 11], [7, 9, 6, -3, 10], [4, -9, -12, -14, -12]]
abs(U) is: [[0, 14, 7, 13, 20], [12, 16, 15, 8, 3], [11, 15, 5, 14, 11], [7, 9, 6, 3, 10], [4, 9, 12, 14, 12]]

eye(12) is: [[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]]
eye(12).trace() is: 12
trace(ones(5)) is: 5

A is: [[2, 1], [4, 5]]
A.is_symetric is: False
A.is_anti_symetric is: False
A.is_diagonal is: False
A.is_hermitian is: False

V is: [[1j, 0j, 0j, 0j], [0j, 1j, 0j, 0j], [0j, 0j, 1j, 0j], [0j, 0j, 0j, 1j]]
V.is_square is: True
V.is_symetric is: True
V.is_anti_symetric is: False
V.is_diagonal is: True
V.is_hermitian is: False

X is: [[3, 0], [1, -9]]
X.trace() is: -6
X.is_upper is: False
X.is_lower is: True
X.T.is_upper is: True
X.T.is_lower is: False
X.is_zero is: False

B = mat_from_f(lambda i, j: '{},{}'.format(i, j), 13) is: [[0,0, 0,1, 0,2, 0,3, 0,4, 0,5, 0,6, 0,7, 0,8, 0,9, 0,10, 0,11, 0,12], [1,0, 1,1, 1,2, 1,3, 1,4, 1,5, 1,6, 1,7, 1,8, 1,9, 1,10, 1,11, 1,12], [2,0, 2,1, 2,2, 2,3, 2,4, 2,5, 2,6, 2,7, 2,8, 2,9, 2,10, 2,11, 2,12], [3,0, 3,1, 3,2, 3,3, 3,4, 3,5, 3,6, 3,7, 3,8, 3,9, 3,10, 3,11, 3,12], [4,0, 4,1, 4,2, 4,3, 4,4, 4,5, 4,6, 4,7, 4,8, 4,9, 4,10, 4,11, 4,12], [5,0, 5,1, 5,2, 5,3, 5,4, 5,5, 5,6, 5,7, 5,8, 5,9, 5,10, 5,11, 5,12], [6,0, 6,1, 6,2, 6,3, 6,4, 6,5, 6,6, 6,7, 6,8, 6,9, 6,10, 6,11, 6,12], [7,0, 7,1, 7,2, 7,3, 7,4, 7,5, 7,6, 7,7, 7,8, 7,9, 7,10, 7,11, 7,12], [8,0, 8,1, 8,2, 8,3, 8,4, 8,5, 8,6, 8,7, 8,8, 8,9, 8,10, 8,11, 8,12], [9,0, 9,1, 9,2, 9,3, 9,4, 9,5, 9,6, 9,7, 9,8, 9,9, 9,10, 9,11, 9,12], [10,0, 10,1, 10,2, 10,3, 10,4, 10,5, 10,6, 10,7, 10,8, 10,9, 10,10, 10,11, 10,12], [11,0, 11,1, 11,2, 11,3, 11,4, 11,5, 11,6, 11,7, 11,8, 11,9, 11,10, 11,11, 11,12], [12,0, 12,1, 12,2, 12,3, 12,4, 12,5, 12,6, 12,7, 12,8, 12,9, 12,10, 12,11, 12,12]]
B[0, 0] is: 0,0

Slicing with the first index: row sub-vectors.
B[1:8:2, 0] is: [[1,0], [3,0], [5,0], [7,0]]
B[0::2, 0] is: [[0,0], [2,0], [4,0], [6,0], [8,0], [10,0], [12,0]]
B[:2, 0] is: [[0,0], [1,0]]
B[:10:3, 0] is: [[0,0], [3,0], [6,0], [9,0]]
B[::3, 0] is: [[0,0], [3,0], [6,0], [9,0], [12,0]]
B[:, 0] is: [[0,0], [1,0], [2,0], [3,0], [4,0], [5,0], [6,0], [7,0], [8,0], [9,0], [10,0], [11,0], [12,0]]
B[3:, 0] is: [[3,0], [4,0], [5,0], [6,0], [7,0], [8,0], [9,0], [10,0], [11,0], [12,0]]

Slicing with the second index: column sub-vectors.
B[0, 1:8:2] is: [[0,1, 0,3, 0,5, 0,7]]
B[0, 0::2] is: [[0,0, 0,2, 0,4, 0,6, 0,8, 0,10, 0,12]]
B[0, :2] is: [[0,0, 0,1]]
B[0, :10:3] is: [[0,0, 0,3, 0,6, 0,9]]
B[0, ::3] is: [[0,0, 0,3, 0,6, 0,9, 0,12]]
B[0, :] is: [[0,0, 0,1, 0,2, 0,3, 0,4, 0,5, 0,6, 0,7, 0,8, 0,9, 0,10, 0,11, 0,12]]
B[0, 3:] is: [[0,3, 0,4, 0,5, 0,6, 0,7, 0,8, 0,9, 0,10, 0,11, 0,12]]

Slicing with the two indeces: sub-matrices.
B[:4, 1:8:2] is: [[0,1, 0,3, 0,5, 0,7], [1,1, 1,3, 1,5, 1,7], [2,1, 2,3, 2,5, 2,7], [3,1, 3,3, 3,5, 3,7]]
B[::4, 0::2] is: [[0,0, 0,2, 0,4, 0,6, 0,8, 0,10, 0,12], [4,0, 4,2, 4,4, 4,6, 4,8, 4,10, 4,12], [8,0, 8,2, 8,4, 8,6, 8,8, 8,10, 8,12], [12,0, 12,2, 12,4, 12,6, 12,8, 12,10, 12,12]]
B[5:, :2] is: [[5,0, 5,1], [6,0, 6,1], [7,0, 7,1], [8,0, 8,1], [9,0, 9,1], [10,0, 10,1], [11,0, 11,1], [12,0, 12,1]]
B[1::2, :10:3] is: [[1,0, 1,3, 1,6, 1,9], [3,0, 3,3, 3,6, 3,9], [5,0, 5,3, 5,6, 5,9], [7,0, 7,3, 7,6, 7,9], [9,0, 9,3, 9,6, 9,9], [11,0, 11,3, 11,6, 11,9]]
B[1:5, ::3] is: [[1,0, 1,3, 1,6, 1,9, 1,12], [2,0, 2,3, 2,6, 2,9, 2,12], [3,0, 3,3, 3,6, 3,9, 3,12], [4,0, 4,3, 4,6, 4,9, 4,12]]
B[:, :] is: [[0,0, 0,1, 0,2, 0,3, 0,4, 0,5, 0,6, 0,7, 0,8, 0,9, 0,10, 0,11, 0,12], [1,0, 1,1, 1,2, 1,3, 1,4, 1,5, 1,6, 1,7, 1,8, 1,9, 1,10, 1,11, 1,12], [2,0, 2,1, 2,2, 2,3, 2,4, 2,5, 2,6, 2,7, 2,8, 2,9, 2,10, 2,11, 2,12], [3,0, 3,1, 3,2, 3,3, 3,4, 3,5, 3,6, 3,7, 3,8, 3,9, 3,10, 3,11, 3,12], [4,0, 4,1, 4,2, 4,3, 4,4, 4,5, 4,6, 4,7, 4,8, 4,9, 4,10, 4,11, 4,12], [5,0, 5,1, 5,2, 5,3, 5,4, 5,5, 5,6, 5,7, 5,8, 5,9, 5,10, 5,11, 5,12], [6,0, 6,1, 6,2, 6,3, 6,4, 6,5, 6,6, 6,7, 6,8, 6,9, 6,10, 6,11, 6,12], [7,0, 7,1, 7,2, 7,3, 7,4, 7,5, 7,6, 7,7, 7,8, 7,9, 7,10, 7,11, 7,12], [8,0, 8,1, 8,2, 8,3, 8,4, 8,5, 8,6, 8,7, 8,8, 8,9, 8,10, 8,11, 8,12], [9,0, 9,1, 9,2, 9,3, 9,4, 9,5, 9,6, 9,7, 9,8, 9,9, 9,10, 9,11, 9,12], [10,0, 10,1, 10,2, 10,3, 10,4, 10,5, 10,6, 10,7, 10,8, 10,9, 10,10, 10,11, 10,12], [11,0, 11,1, 11,2, 11,3, 11,4, 11,5, 11,6, 11,7, 11,8, 11,9, 11,10, 11,11, 11,12], [12,0, 12,1, 12,2, 12,3, 12,4, 12,5, 12,6, 12,7, 12,8, 12,9, 12,10, 12,11, 12,12]]
B[6:, 3:] is: [[6,3, 6,4, 6,5, 6,6, 6,7, 6,8, 6,9, 6,10, 6,11, 6,12], [7,3, 7,4, 7,5, 7,6, 7,7, 7,8, 7,9, 7,10, 7,11, 7,12], [8,3, 8,4, 8,5, 8,6, 8,7, 8,8, 8,9, 8,10, 8,11, 8,12], [9,3, 9,4, 9,5, 9,6, 9,7, 9,8, 9,9, 9,10, 9,11, 9,12], [10,3, 10,4, 10,5, 10,6, 10,7, 10,8, 10,9, 10,10, 10,11, 10,12], [11,3, 11,4, 11,5, 11,6, 11,7, 11,8, 11,9, 11,10, 11,11, 11,12], [12,3, 12,4, 12,5, 12,6, 12,7, 12,8, 12,9, 12,10, 12,11, 12,12]]

Modifying a slice with the first index: row sub-vectors.

Modifying a slice (with the first index) with a constant value: 0
B[1:8:2, 0] is: [[1,0], [3,0], [5,0], [7,0]]
B[1:8:2, 0] is: [[0], [0], [0], [0]]

Modifying a slice (with the first index) with a list: [1, 1, 1, 1, 1, 1, 1]
B[0::2, 0] is: [[0,0], [2,0], [4,0], [6,0], [8,0], [10,0], [12,0]]
B[0::2, 0] is: [[1], [1], [1], [1], [1], [1], [1]]

Modifying a slice (with the first index) with a row vector: [[6, 6, 6, 6]]
B[:10:3, 0] is: [[1], [0], [1], [9,0]]
B[:10:3, 0] is: [[6], [0], [[[6, 6, 6, 6]]], [9,0]]
WARNING: Modifying a slice with a row vector (list of 1 list) fails.
B[:10:3, 0] is fixed: [[6], [0], [6], [9,0]]

Modifying a slice (with the first index) with a column vector: [[-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j]]
B[:, 2] is: [[0,2], [1,2], [2,2], [3,2], [4,2], [5,2], [6,2], [7,2], [8,2], [9,2], [10,2], [11,2], [12,2]]
B[:, 2] is: [[-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j]]

Modifying a slice (with the first index) with a row vector (as a Matrix object): [[-9, -9, -9, -9]]
B[:10:3, 0] is: [[6], [0], [6], [9,0]]
B[:10:3, 0] is: [[-9], [0], [[[-9, -9, -9, -9]]], [9,0]]
WARNING: Modifying a slice with a row vector (list of 1 list) fails.
B[:10:3, 0] is fixed: [[-9], [0], [-9], [9,0]]

Modifying a slice (with the first index) with a column vector (as a Matrix object): [[4], [4], [4], [4], [4], [4], [4], [4], [4], [4], [4], [4], [4]]
B[:, 2] is: [[-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j], [-3j]]
B[:, 2] is: [[4], [4], [4], [4], [4], [4], [4], [4], [4], [4], [4], [4], [4]]

Gauss elimination examples:

A is: [[2, 1], [4, 5]]
A.gauss() row echelon form of A is: [[2, 1], [0, 3]]
det_A is: 6
A.det is: 6
A.rank is: 2
det(A) is: 6
rank(A) is: 2

Changing A with: A[1, :] = A[0, :]
Now A is: [[2, 1], [2, 1]]
A.det is: 0
A.rank is: 1

Z = zeros(3) is: [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
Z.det is: 0
Z.rank is: 0


Z2 is: [[0, 0, 0], [0, 0, 1], [0, 0, 0]]

Trying to find the 0-th pivot:
With these indeces: [0, 1, 2]
And that array: [0, 0, 0]
_argmax has given i_max = 0, and c[k, i_max] = 0 (with k = 0).
WARNING: A.gauss() might have been called on a singular matrix. FIXME remove these warnings

Trying to find the 1-th pivot:
With these indeces: [1, 2]
And that array: [0, 0, 1]
_argmax has given i_max = 2, and c[k, i_max] = 1 (with k = 1).
c.col(i_max) is: [[0], [1], [0]]
c.col(k) is: [[0], [0], [0]]
i_max = 2, k = 1.
We swapped two different lines (2 and 1), the determinant will be multiplied by -1.
Gauss Elimination: using the 1th line (L_1 = [[0, 1, 0]]).
  We use 1 as a pivot.
 Operation L_2' <-- L_2 - gamma * L_1
  with gamma = 0.0
  with old L_2 = [[0, 0, 0]]
  with new L_2' = [[0, 0, 0]]

Trying to find the 2-th pivot:
With these indeces: [2]
And that array: [0, 0, 0]
_argmax has given i_max = 2, and c[k, i_max] = 0 (with k = 2).
WARNING: A.gauss() might have been called on a singular matrix. FIXME remove these warnings
Z2.det is: 0
Z2.rank is: 1


With this matrix A0: [[2, -1, 0], [-1, 2, -1], [0, -1, 2]] we try the extended Gauss-Jordan algorithm.
A0.map(Fraction) is: [[2, -1, 0], [-1, 2, -1], [0, -1, 2]]
A0.gauss_jordan(inv=True): row echelon form of A0 is: [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
A0.gauss_jordan(inv=True): inverse of A0 is: [[3/4, 1/2, 1/4], [1/2, 1, 1/2], [1/4, 1/2, 3/4]]
det(A0) is: 4

A0 * inv_A0 is: [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
inv_A0 * A0 is: [[1, 0, 0], [0, 1, 0], [0, 0, 1]]


Looking for a non-singular random matrix R of size (4, 4):
Trying [[-2, 6, 4, 8], [6, 1, 9, 2], [-7, -1, -5, 8], [-10, -4, 2, -4]]
R.map(Fraction) is: [[-2, 6, 4, 8], [6, 1, 9, 2], [-7, -1, -5, 8], [-10, -4, 2, -4]]
Note: we use Fraction in order to be exact and not numerically approximative !
R.gauss_jordan(inv=True): row echelon form of R is: [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]
R.gauss_jordan(inv=True): inverse of R is: [[-91/1752, 37/876, 1/146, -121/1752], [269/1752, -119/876, -19/146, -37/1752], [3/146, 6/73, -1/73, 4/73], [-47/3504, 125/1752, 31/292, -101/3504]]
det(R) is: 7008

R * inv_R is: [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]
inv_R * R is: [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]

For A: [[1, 2], [3, 4]]

For i = 0 and j = 0
  - The (i,j) minor is 4
  - The (i,j) cofactor is 4

For i = 0 and j = 1
  - The (i,j) minor is 3
  - The (i,j) cofactor is -3

For i = 1 and j = 0
  - The (i,j) minor is 2
  - The (i,j) cofactor is -2

For i = 1 and j = 1
  - The (i,j) minor is 1
  - The (i,j) cofactor is 1
So the co-matrix of A is [[4, -3], [-2, 1]]

For A: [[-1, 4], [4, 1]]

For i = 0 and j = 0
  - The (i,j) minor is 1
  - The (i,j) cofactor is 1

For i = 0 and j = 1
  - The (i,j) minor is 4
  - The (i,j) cofactor is -4

For i = 1 and j = 0
  - The (i,j) minor is 4
  - The (i,j) cofactor is -4

For i = 1 and j = 1
  - The (i,j) minor is -1
  - The (i,j) cofactor is -1
So the co-matrix of A is [[1, -4], [-4, -1]]

For A: [[1, 2, -1, 0], [2, 4, -2, -1], [-3, -5, 6, 1], [-1, 2, 8, -2]]
[["<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>"], ["<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>"], ["<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>"], ["<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>"]]
We compute the permuted LU decomposition...

We have computed the permutation matrix P: [[1, 0, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0], [0, 1, 0, 0]]
[["<type int>", "<type int>", "<type int>", "<type int>"], ["<type int>", "<type int>", "<type int>", "<type int>"], ["<type int>", "<type int>", "<type int>", "<type int>"], ["<type int>", "<type int>", "<type int>", "<type int>"]]

And the lower triangular matrix L: [[1, 0, 0, 0], [2, 1, 0, 1], [-3, 1/4, 1, 0], [-1, 0, 0, 0]]
L.type():
 [["<type int>", "<type int>", "<type int>", "<type int>"], ["<class matrix.Fraction>", "<type int>", "<type int>", "<type int>"], ["<class matrix.Fraction>", "<class matrix.Fraction>", "<type int>", "<type int>"], ["<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>", "<type int>"]]
L.is_lower: False
L.is_upper: False

And the upper triangular matrix  U: [[1, 2, -1, 0], [0, 4, 7, -2], [0, 0, 5/4, 3/2], [0, 0, 0, -1]]
U.type():
 [["<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>"], ["<class matrix.Fraction>", "<class fractions.Fraction>", "<class fractions.Fraction>", "<class fractions.Fraction>"], ["<class matrix.Fraction>", "<class matrix.Fraction>", "<class fractions.Fraction>", "<class fractions.Fraction>"], ["<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>", "<class fractions.Fraction>"]]
U.is_lower: False
U.is_upper: True

So P * A is: [[1, 2, -1, 0], [-1, 2, 8, -2], [-3, -5, 6, 1], [2, 4, -2, -1]]
So L * U is: [[1, 2, -1, 0], [2, 8, 5, -3], [-3, -5, 6, 1], [-1, -2, 1, 0]]
P * A == L * U: False


Examples are done. Are you satisfied?