# -*- coding: UTF-8 -*-
#!/usr/bin/python3
import numpy as np

def GENP(A, b):
    '''
    Gaussian elimination with no pivoting.
    % input: A is an n x n nonsingular matrix
    %        b is an n x 1 vector
    % output: x is the solution of Ax=b.
    % post-condition: A and b have been modified. 
    '''
    n =  len(A)
    if b.size != n:
        raise ValueError("Invalid argument: incompatible sizes between A & b.", b.size, n)

    for pivot_row in range(n-1):
        for row in range(pivot_row+1, n):
            multiplier = A[row][pivot_row]/A[pivot_row][pivot_row]

            A[row][pivot_row] = multiplier
            for col in range(pivot_row + 1, n):
                A[row][col] = A[row][col] - multiplier*A[pivot_row][col]

            b[row] = b[row] - multiplier*b[pivot_row]

    #print(A)
    #print(b)
    x = np.zeros(n)
    k = n-1
    x[k] = b[k]/A[k,k]
    while k >= 0:
        x[k] = (b[k] - np.dot(A[k,k+1:],x[k+1:]))/A[k,k]
        k = k-1
    return x

def GEPP(A, b):
    '''
    Gaussian elimination with partial pivoting.
    % input: A is an n x n nonsingular matrix
    %        b is an n x 1 vector
    % output: x is the solution of Ax=b.
    % post-condition: A and b have been modified. 
    '''
    n =  len(A)
    if b.size != n:
        raise ValueError("Invalid argument: incompatible sizes between A & b.", b.size, n)

    for k in range(n-1):
        maxindex = abs(A[k:,k]).argmax() + k
        if A[maxindex, k] == 0:
            raise ValueError("Matrix is singular.")

        if maxindex != k:
            A[[k,maxindex]] = A[[maxindex, k]]
            b[[k,maxindex]] = b[[maxindex, k]]
        for row in range(k+1, n):
            multiplier = A[row][k]/A[k][k]

            A[row][k] = multiplier
            for col in range(k + 1, n):
                A[row][col] = A[row][col] - multiplier*A[k][col]

            b[row] = b[row] - multiplier*b[k]
    #print(A)
    #print(b)
    x = np.zeros(n)
    k = n-1
    x[k] = b[k]/A[k,k]
    while k >= 0:
        x[k] = (b[k] - np.dot(A[k,k+1:],x[k+1:]))/A[k,k]
        k = k-1
    return x

if __name__ == "__main__":
    '''
    A = np.array([[1.,-1.,1.,-1.],[1.,0.,0.,0.],[1.,1.,1.,1.],[1.,2.,4.,8.]])
    b = np.array([[14.],[4.],[2.],[2.]])    
    '''
    #generate matrix A
    A = np.zeros((84,84))
    A[0][0]=6
    A[0][1]=1
    A[83][82]=8
    A[83][83]=6
    col = 0
    for row in range(1, 83):
        A[row][col]  =8
        A[row][col+1]=6
        A[row][col+2]=1
        col=col+1
    #print(A)

    #generate matrix A
    b = np.zeros((84, 1))
    b[0]=7
    b[83]=14
    for row in range(1, 83):
        b[row]=15
    #print(b)

    print("Gaussian elimination with no pivoting.")
    print(GENP(np.copy(A), np.copy(b)))

    print("Gaussian elimination with partial pivoting.")
    print(GEPP(A,b))