import numpy as np
import copy
# 
def gauss(A,b):
    A = np.copy(A)
    b = np.copy(b)
    n = np.size(A,0)

    for k in range(n-1):
        for i in range(k+1,n):
            piv = A[i,k]/A[k,k]
            for j in range(k+1,n):
                A[i,j] -= piv*A[k,j]
            b[i] -= piv*b[k]
    print (A)
    print (b)
    # -- note: in-place back-solve b ovrteritten
    backsolve_in(A,b)
    return b

# pivoting 
def gaussp(A,b):
    A = np.copy(A)
    b = np.copy(b)
    n = np.size(A,0)
    temp = np.zeros(n)
    
    for k in range(n-1):
        [t, ip] = maxidx(A, k, k)
        if (ip != k):
            #print(f'swapping rows {k:3d} & {ip:3d}')
            print("swapping rows %2d & %2d" % (k,ip))
            temp[k:] = A[k,k:]
            A[k,k:]  = A[ip,k:] 
            A[ip,k:] = temp[k:]
            t = b[k] 
            b[k] = b[ip]
            b[ip]  = t
        for i in range(k+1,n):
            piv = A[i,k]/A[k,k]
            for j in range(k+1,n):
                A[i,j] -= piv*A[k,j]
            b[i] -= piv*b[k]
    print (A)
    print (b)
    # -- note: in-place back-solve b ovrteritten
    backsolve_in(A,b)
    return b

def gaussLU(A):
    n = np.size(A,0)
    L = np.identity(n) 
    U = np.copy(A) 

    for k in range(n-1):
        for i in range(k+1,n):
            piv = U[i,k]/U[k,k]
            L[i,k] = piv
            for j in range(k+1,n):
                U[i,j] -= piv*U[k,j]
            U[i,k] = 0.0
    # -- Done -- return L and U 
    return (L,U)




def backsolve_in(A,x):
    # inplace back-solve. on return  Rhs x <- sol.
    # uses only upper triang. part of A.
    n = np.size(A,0)
    for i in range(n-1,-1,-1):
        t = x[i]
        for j in range(i+1,n):
            t -= A[i,j]*x[j]
        x[i] = t/A[i,i]
        
def maxidx(A, k, l):
# max entry in column l - starting in row k     
    m,n = np.shape(A)
    imax = -1;
    tmax = -np.Inf
    for i in range(k,m):
        t = np.abs(A[i,l])
        if (t>tmax):
            tmax = t
            imax = i
    return(tmax,imax)