def lup(A):
    # функция для LUP-разложения матрицы A
    # возвращает L, U и столбец престановок pi

    # копирование матрицы (для предотвращения изменения внешней матрицы)
    A = [list(row) for row in A]
    n = len(A) # число строк в мтр A
    pi = [i for i in range(n)] # создание столбца перестановок

    for k in range(n):
        p = 0.0
        k1 = k

        for i in range(k, n): # для i от k до n-1
            if abs(A[i][k]) > p:
                p = abs(A[i][k])
                k1 = i

        if p == 0.0:
            print("error") # выдает ошибку при обнаружении вырожденной матрицы

        pi[k], pi[k1] = pi[k1], pi[k]  # обмен значений pi[k] и pi[k1]
        A[k], A[k1] = A[k1], A[k]      # обмен строк A[k] и A[k1]

        for i in range(k + 1, n):
            A[i][k] = A[i][k] / A[k][k] # деление на ведущий элемент
            for j in range(k + 1, n):
                A[i][j] = A[i][j] - A[i][k] * A[k][j] # вычитание из подматрицы

    L = [[0] * n for i in range(n)]# создание пустых матриц L, U
    U = [[0] * n for i in range(n)]

    # далее формирование матриц
    for i in range(n):
        for j in range(n):
            if i > j:
                L[i][j] = A[i][j]
            else:
                U[i][j] = A[i][j]

    for i in range(n):
        L[i][i] = 1.0

    return pi, L, U

def solve(A, b):
    # функция для нахождения корней х с помощью lup

    pi, L, U = lup(A) # вычисление мтр с помощью функции lup() для мтр а

    n = len(L) # число строк в матрице L
    x = [0] * n  # создание пустых векторов x, y для записи решения
    y = [0] * n

    for i in range(n):
        s = sum(L[i][j] * y[j] for j in range(i)) # суммирование по j до i-1
        y[i] = b[pi[i]] - s                       # прямая подстановка

    for i in reversed(range(n)): #разворот для прохода от n-1 до 0
        s = sum(U[i][j] * x[j] for j in range(i + 1, n)) # суммирование от j = i+1 до n
        x[i] = (y[i] - s) / U[i][i] #  обратная подстановка

    return x


A = [[2.0, 1.0, 3.0], [1.0, 1.0, 1.0], [3.0, 1.0, 1.0]]
b = [13.0, 6.0, 8.0]
print(solve(A, b)) # выведет [1.0, 2.0, 3.0]