SALE

dec.py

def LU_decomposition(A):
    """
    input:
        A - матрица
        
    output:
        (LU, P)
        LU - матрица разложения
        P - перестановка строк
        Q - перестановка столбцов
    """
    
    size = A.shape # размер матрицы
    P = Permutation(size[0]) # перестановка строк
    Q = Permutation(size[1]) # перестановка столбцов
    
    LU = A.copy()
    
    for i in range(min(*size)):
        max_idx = np.argmax(LU[P.P[i:]][:, Q.P[i:]])
        if eq(LU[P.P[i:]][:, Q.P[i:]].flat[max_idx], 0):
            break
            
        P.swap(i, i + max_idx // (size[1]-i))
        Q.swap(i, i + max_idx % (size[1]-i))   
    
        #j = Q[i] - столбец ведущего элемента
        #t = P[i] - строка ведущего элемента

        t = P[i]
        j = Q[i]
    
        if i + 1 < P.size():
            for r in P[i+1:]:
                LU[r][j] = LU[r][j] / LU[t][j]
                LU[r][Q[i+1:]] = LU[r][Q[i+1:]] - LU[t][Q[i+1:]] * LU[r][j]
    LU = LU[P.P][:, Q.P]

        
    return LU, P, Q

P.py

class Permutation(object):
    """
     Класс, описывающий перестановку строк в матрице
     
     P - массив перестановки
     even - чётность перстановки
         1  - чётна
         -1 - нечётна
    """
    def __init__(self, n):
        self.P = np.arange(n)
        self.even = 1
        
    def __getitem__(self, i):
        return self.P[i]
    
    def swap(self, i, j):
        """
        Поменять строки i и j местами
        """
        if i != j:
            self.P[i], self.P[j] = self.P[j], self.P[i]
            self.even *= -1
            
    def size(self):
        return self.P.size
    
    def __str__(self):
        return self.P.__str__()

SLAE.py

def SLAE_triangle(A, B, diagonal_is_ones = False):
    """
    Решение СЛАУ, у которых матрица А имеет вид верхней правой треугольной матрицы
    
    diagonal_is_ones = True, если матрица унитреугольная
    """
    b = B.copy()
    n = A.shape[1]
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = b[i] / (A[i][i] if not diagonal_is_ones else 1)
        if i > 0:
            b = b - b[i]*(A[:, i] / (A[i][i] if not diagonal_is_ones else 1))
    return x
  
  
  def SLAE(LU, b, P: Permutation, Q: Permutation):
    """
    Решение СЛАУ по LU разложению для невырожденных матриц
    """
    b = b[P.P]
    
    # если размерность столбца b не соотносится с числом строк
    if LU.shape[0] != b.size:
        b = b[:LU.shape[0]]

    y = SLAE_triangle(LU[::-1, ::-1], b[::-1], True)[::-1] # Ly = Pb
    
    x = SLAE_triangle(LU, y) # Ux = y
    X = np.zeros(Q.size())
    X[:x.size] = x
    return X[Q.inverse()]
  
  
  def SLAE2(A, b):
    """Решение СЛАУ для любых матриц"""
    def rank(LU):
        for i in range(min(*LU.shape)):
            if eq(LU[i,i], 0):
                return i
        return min(*LU.shape)
    
    LU, P, Q = LU_decomposition(A)
    print("LU=", LU)
    if eq(find_det(LU, P, Q), 0):
        print("det(A) = 0")
        LUb, _, _ = LU_decomposition(np.concatenate((A, b.reshape((-1, 1))), axis=1))
        rank_A = rank(LU)
        print(f"rank(A) = {rank_A}")
        if rank_A != rank(LUb):
            print("Система несовместна")
            return False
        print("Система совместна")
        return SLAE(LU[:rank_A, :rank_A], b, P, Q)
        
        
    else:
        print("det(A) != 0")
        return SLAE(LU, b, P, Q)