python --- 实现LU=PA分解（部分主元的高斯消去）/带行交换的LU分解

题设： Pascal矩阵 进行 LU=PA分解

import numpy as np
import random
from scipy.linalg import lu

#Вычисление значения определителя матрицы.
def det(A):
    if len(A) <= 0:
        return None
    elif len(A) == 1:
        return A[0][0]
    else:
        s = 0
        for i in range(len(A)):
            n = [[row[a] for a in range(len(A)) if a != i] for row in A[1:]]  # 这里生成余子式
            s += A[0][i] * det(n) * (-1) ** (i )
        return s

#Judge whether the k-th SequentialPrincipal Minor of A is non-zero,
# and If it is non-zero, the LU decomposition condition is satisfied
def Master_Sequential(A,n):
    for i in range(0,n):
        Master = np.zeros([i+1,i+1])
        for row in range(0,i+1):
            for a in range(0,i+1):
                Master[row][a]=A[row][a]
        if det(Master)==0:
            done=False
            return done
#-LU разложение с выбором главного элемента по строке (схема частичного выбора)
def LU_with_partial_pivoting(A):
    A = np.array(A, dtype=np.float64)
    L = np.zeros_like(A)
    m = A.shape[0]

    P = np.identity(m)

    for i in range(m-1):
        #Get the maximum number of rows of absolute value elements in this column
        idx = i + np.argmax(abs(A[i:, i]))
        #Swap rows of a, l, P
        if idx != i & idx != m-1:
            A[[i, idx], :] = A[[idx, i], :]
            L[[i, idx], :] = L[[idx, i], :]
            P[[i, idx], :] = P[[idx, i], :]
        #
        L[i, i] = 1
        for j in range(i+1, m):
            #Calculate the values of the L and u elements
            L[j, i] = A[j, i] / A[i, i]
            A[j, :] = A[j, :] - L[j, i] * A[i, :]
    #Finally, deal with the main diagonal element
    L[i+1, i+1] = 1
    #Get u
    U = np.array(A)
    print('L = \n', L)
    print('U = \n', U)
    print('P = \n', P)

#LU- разложением без выбора главного элемента.
def LU_without_P(A):
    n=len(A[0])
    L = np.zeros([n,n])
    U = np.zeros([n, n])
    for i in range(n):
        L[i][i]=1
        if i==0:
            U[0][0] = A[0][0]
            for j in range(1,n):
                U[0][j]=A[0][j]
                L[j][0]=A[j][0]/U[0][0]
        else:
                for j in range(i, n):#U
                    temp=0
                    for k in range(0, i):
                        temp = temp+L[i][k] * U[k][j]
                    U[i][j]=A[i][j]-temp
                for j in range(i+1, n):#L
                    temp = 0
                    for k in range(0, i ):
                        temp = temp + L[j][k] * U[k][i]
                    L[j][i] = (A[j][i] - temp)/U[i][i]
    return L,U


#Создайте матрицу Pascal произвольного размера.
def pascalmatrix(n):
    C = [[0 for x in range(2 * n + 1)]

            for y in range(2 * n + 1)]

    for i in range(2 * n + 1):
        for j in range(min(i, 2 * n) + 1):
            if (j == 0 or j == i):
                C[i][j] = 1
            else:
                C[i][j] = (C[i - 1][j - 1] + C[i - 1][j])

    P = np.zeros(shape=(n,n))
    for i in range(n):
        for j in range(n):
            P[i,j] = C[i + j][i]

    return P

#Проверьте выполнение det() = 1 и зависимость определителя от .
def check_det():
    TF = True
    for i in range(1, 10):
        if det(pascalmatrix(i)) != 1:
            TF = False
    return TF



#Test program
n= random.randint(2,6)
print('n = ',n)
A = pascalmatrix(n)
print("Pascal Matrix = \n",A)
LU_with_partial_pivoting(A)

if(check_det()):
    print("\033[31mdet(A) == 1 and не зависимость определителя от n\033[0m")


print("\033[5;31;46mПовторите эксперименты с 𝐿𝑈- разложением без выбора главного элемента =>\033[0m")


if Master_Sequential(A,n) != False:
        L,U=LU_without_P(A)
        print('L =：\n',L, '\nU =：\n',U)
else:
        print('LU decomposition conditions are not satisfied。')

运行结果如下：