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.08553692318766, 0.0], [0.0, 0.36787944117144245]]
Conversely diag([exp(3), exp(-1)]) is: [[20.085536923187668, 0], [0, 0.36787944117144233]]

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.86038416110793, 80.14210233264886], [320.56840933059544, 323.2866911590545]]

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: [[-2, 0], [-1, -3], [3, -2]]

A random matrix of size (2, 4) and float coefficients between -2 and 2 is: [[-1.1218240398994217, 1.2308776567513258, -0.4184343804448649, -0.09548031576957028], [1.7877586493153328, -0.7592711065155795, -0.3564851818977819, -1.4761697508978053]]

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: [[6, 9], [7, 9], [0, 0], [-6, -2]]
Vector v is: [1, 1]
Z.dot(v) is: [[15], [16], [0], [-8]]

A is: [[2, 1], [4, 5]]
A.norm() is: 6.782329983125268
A.norm(p=1) is: 12 which is like sum(A): 12
A.norm(3) is: 5.828476683251456
norm(A) is: 6.782329983125268
norm(A, 42) is: 5.000010127030902

Random U is: [[0, -12, -13, -3, -18], [-1, 9, 11, 11, 9], [7, -17, 4, -17, -11], [-16, 1, -7, -5, 10], [17, 5, -12, -11, 3]]
abs(U) is: [[0, 12, 13, 3, 18], [1, 9, 11, 11, 9], [7, 17, 4, 17, 11], [16, 1, 7, 5, 10], [17, 5, 12, 11, 3]]

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: [[(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-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: [[(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-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: [[(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-3j)], [(-0-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 [[8, -5, -4, 7], [-9, -7, -2, -5], [-2, 9, -10, -10], [-2, 10, 3, 8]]
R.map(Fraction) is: [[8, -5, -4, 7], [-9, -7, -2, -5], [-2, 9, -10, -10], [-2, 10, 3, 8]]
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: [[59/19419, -609/6473, -26/19419, -1226/19419], [-508/19419, -242/6473, 553/19419, 682/19419], [-1790/19419, -394/6473, -1186/19419, -655/19419], [1321/19419, 298/6473, -253/19419, 1514/19419]]
det(R) is: 19419

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]]
[["<class int>", "<class int>", "<class int>", "<class int>"], ["<class int>", "<class int>", "<class int>", "<class int>"], ["<class int>", "<class int>", "<class int>", "<class int>"], ["<class int>", "<class int>", "<class int>", "<class 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():
 [["<class int>", "<class int>", "<class int>", "<class int>"], ["<class matrix.Fraction>", "<class int>", "<class int>", "<class int>"], ["<class matrix.Fraction>", "<class matrix.Fraction>", "<class int>", "<class int>"], ["<class matrix.Fraction>", "<class matrix.Fraction>", "<class matrix.Fraction>", "<class 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?