#!/usr/bin/env python
#
# file: $ENGR_2011/lecture_02.py
#
# A simple Python program that demonstrates Gaussian elimination and rank.
#

# import system modules
#
import numpy as np
from numpy.linalg import matrix_rank
import os
import sys

# function: gauss_elim
#
# arguments:
#  a: the input matrix (a numpy 2D array)
#
# returns:
#  x: the solution as a numpy vector
#
def gauss_elim(a):

    # create space for the output array
    #
    n = len(a)
    x = np.zeros(n)

    # Applying Gauss Elimination
    for i in range(n):
        if a[i][i] == 0.0:
            sys.exit('Divide by zero detected!')
        
        for j in range(i+1, n):
            ratio = a[j][i]/a[i][i]
        
            for k in range(n+1):
                a[j][k] = a[j][k] - ratio * a[i][k]

    # Back Substitution
    x[n-1] = a[n-1][n]/a[n-1][n-1]

    for i in range(n-2,-1,-1):
        x[i] = a[i][n]
    
        for j in range(i+1,n):
            x[i] = x[i] - a[i][j]*x[j]
            x[i] = x[i]/a[i][i]

    # exit gracefully
    #
    return x

# function: main
#
def main(argv):

    # create an augmented matrix
    #
    a = np.array([[1, -2, 1, 0], [2, 1, -3, 0], [4, -7, 1, -1]])
    print("The input matrix is:")
    print(a)

    x = gauss_elim(a)

    # Displaying solution
    #
    print("The solution is:")
    print(x)

    # Also compute and display the rank of a matrix
    #
    b = np.array([[1, -2, 1, 0], [2, 1, -3, 0], [4, -7, 1, -1]])
    print("The input matrix is:")
    print(b)
    print("The rank is:", matrix_rank(b))

    # exit gracefully
    #
    return True

# begin gracefully
#
if __name__ == '__main__':
    main(sys.argv[0:])

#
# end of file