e2

#https://scicomp.stackexchange.com/questions/10815/sort-of-problems-where-sor-is-faster-than-gauss-seidel
import numpy as np

# Declaring the A and b 
A = np.matrix([[0.0,3.0,-1.0,8.0],[-1.0,11.0,-1.0,3.0],[2.0,-1.0,10.0,-1.0],[10.0,-1.0,2.0,0.0]])
b = np.array([15.0,25.0,-11.0,6.0])

#debug = False
debug = True


"""
  Parameters
  ----------
  A  : Matrix to be be checked  
  
  Returns
  -------  
  True incase it is strictly dominant  
  using definition in https://mathworld.wolfram.com/DiagonallyDominantMatrix.html#:~:text=If%20a%20matrix%20is%20strictly,of%20its%20eigenvalues%20are%20negative.
  np.all(diagonal > others.sum(axis=1)) for strictly diagonally dominant 
  np.all(diagonal >= others.sum(axis=1)) for diagonally dominant 
"""
def is_strictly_diagonally_dominant_np(A):
    A = np.abs(np.asarray(A))  # converts lists to numpy arrays
    diagonal = np.diag(A)  # 
    others = A - np.diag(diagonal)  # matrix 'A' but with zeros along the diagonal
    return np.all(diagonal > others.sum(axis=1)) # tests if each element from diagonal is greater than the sum of the remaining elements in the row. 


"""
  Parameters
  ----------
  A  : Matrix to be pivoted 
  
  Returns
  -------  
  A  : Pivoted matrix 
  piv: Pivote 

"""
def pivot_matrix(A):
    n = A.shape[0] # Get shape of matrix in this case 4 
    piv = np.arange(0,n) # piv array index [0 1 2 3]  
    if(debug): print("Pivoting matrix: A")
    if(debug): print(A)
    if(debug): print("Setting intial piv array :")
    if(debug): print(piv)
    if(debug): print("----------------------------")
    for k in range(n-1):
        max_row_index = np.argmax(abs(A[k:n,k])) + k 
        if(debug): print("Detected max row index at row #", max_row_index , "-> switching row #" , max_row_index , " with row #" , k  )
        piv[[k,max_row_index]] = piv[[max_row_index,k]] 
        A[[k,max_row_index]] = A[[max_row_index,k]]
        if(debug): print(A)
        if(debug): print(piv)

    return [A,piv]


def forward_substitution(L,b):
    y = np.full_like(b,0) 
    for i in range(L.shape[0]):
        y[i] = b[i]
        for j in range(i):
            y[i] -= L[i,j]*y[j]
        y[i] = y[i] / L[i,i]
    return y

def backward_substitution(U,y):
    for i in range(U.shape[0]-1,-1,-1): 
        for j in range(i+1, U.shape[1]):
            y[i] -= U[i,j]*y[j]
        y[i] = y[i]/U[i,i]
    return y

def doolittle(A):
    
  n = len(A)
  L = np.zeros([n,n])
  U = np.zeros([n,n])
        
  for z in range(n):
    L[z, z] = 1 # Diagonal elements are all ones for Lower trangular matrix
    U[z, z] = (A[z, z] - np.dot(L[z, :z], U[:z, z])) # For taking care of first eliment in diagonal order 

    # Since first element is already taken care of start with z+1 row and z+1 col.
    for i in range(z+1, n):
      # Algebraically, the dot product is the sum of the products of the corresponding entries of the two sequences of numbers
      # https://www.geeksforgeeks.org/doolittle-algorithm-lu-decomposition/?ref=lbp
      # Now take care of U element starting second col toward end of the matrix in the row.
      U[z, i] = (A[z, i] - np.dot(L[z, :z], U[:z, i]))
      if(debug): print("----------------------------")  
      if(debug): print("in U loop i:", i , " z:", z)
      if(debug): print(A)
      if(debug): print(U)
      if(debug): print(L) 
      
    if(debug): print("----------------------------------------------")
    for k in range(z+1, n):
      # Now take care of L elements starting second row towards bottom of the matrix 
      L[k, z] = (A[k, z] - np.dot(L[k, :z], U[:z, z])) / U[z, z]
      if(debug): print("----------------------------")    
      if(debug): print("in L loop k:", k , " z:", z)
      if(debug): print(A)
      if(debug): print(U)
      if(debug): print(L) 

  return (L, U)


def crout_1(A):
    """
    Returns the lower-triangular matrix L and the unit upper-triangular
    matrix U such that L*U = the given matrix A.
    The computation uses Crout's Algorithm to perform LU decomposition on A.
    INPUT:
        - A: list; the square matrix to decompose
    OUTPUT:
        - list; the matrix L followed by the matrix U
    """
    # This is Crout's Algorithm.
    n = len(A)
    L = [[0] * n for i in range(n)]
    U = [[0] * n for i in range(n)]
    for j in range(n):
        U[j][j] = 1             # set the j,j-th entry of U to 1
        for i in range(j, n):  # starting at L[j][j], solve j-th column of L
            alpha = float(A[i][j])
            for k in range(j):
                alpha -= L[i][k]*U[k][j]
            L[i][j] = alpha
        for i in range(j+1, n):# starting at U[j][j+1], solve j-th row of U
            tempU = float(A[j][i])
            for k in range(j):
                tempU -= L[j][k]*U[k][i]
            if int(L[j][j]) == 0:
                L[j][j] = e-40
            U[j][i] = tempU/L[j][j]
    return [L, U]


def crout(A):
    """
    Returns the lower-triangular matrix L and the unit upper-triangular
    matrix U such that L*U = the given matrix A.
    Uses Crout's Algorithm to perform LU decomposition on A.
    INPUT:
        - A: The square matrix to decompose
    OUTPUT:
        - The matrix L followed by the matrix U
    """     
    n = len(A)
    L = np.zeros([n, n], dtype = float) 
    U = np.zeros([n, n], dtype = float)
    
   
    for z in range(n):
        
        U[z,z] = 1  # Diagonal elements are all ones for Lower trangular matrix         
        
        # for-loop starting at L[j][j] in order to solve the j-th column of L
        
        for j in range(z,n):

            # Declaring a temporary L to store values and insert them later in the L matrix

            temporary_L = float(A[j,z]) # 1. Take A[j,z]
            
            for k in range(z):
                
                temporary_L -= L[j,k]*U[k,z] # 2. Keep subtracting L[j,k]*U[k,z] 
                
            L[j,z] = temporary_L # 3. Assign to L[j,z]
            
        # for-loop starting at U[j][j+1] in order to solve the j-th row of U
        
        for j in range(z+1, n):
            
            # Declaring a temporary U to store values and insert them later in the U matrix
            
            temporary_U = float(A[z,j]) # 1. Take A[z,j]
            
            for k in range(z):
                
                temporary_U -= L[z,k]*U[k,j] # 2. Keep subtracting L[z,k]*U[k,j] 
                
            U[z,j] = temporary_U / L[z,z] # 3. Assign to  U[z,j]
    
    
    return (L, U)


def solveusing(algorithm_used,A, b):

    L, U = algorithm_used(A)
    
    print("L = " + str(L) + "\n")
    print("U = " + str(U) + "\n")
   
    y = forward_substitution(L,b)
    x = backward_substitution(U,y)
        
    return x


# Printing out the results

# Check if the matrix is square
if A.shape[0] != A.shape[1]:
    print('Input argument is not a square matrix.')
    exit()

print("---------------------------------------------------------------------------------------")

print("Objective 1: Write a program to test if the matrix A is strictly diagonally dominant.")

# Check if the matrix is diagonally dominant 
if(is_strictly_diagonally_dominant_np(A)):
    print("Matrix A is diagonally dominant")
else:
    print("Matrix A is not diagonally dominant")

print("---------------------------------------------------------------------------------------")

print("Objective 2: LU Decomposition: Write a program to obtain a solution of (*) using both methods of Doolittle")
print("and Crout. Your program in each case should list L and U; followed by the solution for x")

print("")
print("Make the Matrix diagonally dominant")
print("--------------------------------------------")
#Pivot to get diogonially dominant matrix 
A, piv = pivot_matrix(A) 
b = b[piv]
print(A)
print("--------------------------------------------")

# For Doolittle 
#print("Using Doolittle's algorithm:" + "\n")
#print("------------------------------------------------")
#print( "x = " + str(solveusing(doolittle,A,b)) + "\n" )

# For Crout 
print("\n" + "Using Crout's algorithm:" + "\n")
print("------------------------------------------------")
print( "x = " + str(solveusing(crout,A,b)) + "\n")