Compare the result of Gaussian elimination with the output of numpy.linalg.solve

Question

In the following code I have implemented Gaussian elimination with partial pivoting for a general square linear system Ax=b. I have tested my code and it produced the right output. I have used it to solve Ax=b where A is a random 100x100 matrix and b is a random 100x1 vector.**

However now I am looking for some help with comparing my solution against the solution obtained by using numpy.linalg.solve. How to add this comparison to my code?

import numpy as np
def GEPP(A, b, doPricing = True):
    '''
    Gaussian elimination with partial pivoting.
    input: A is an n x n numpy matrix
           b is an n x 1 numpy array
    output: x is the solution of Ax=b 
        with the entries permuted in 
        accordance with the pivoting 
        done by the algorithm
    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)

    # k represents the current pivot row. Since GE traverses the matrix in the 

    # upper right triangle, we also use k for indicating the k-th diagonal 

    # column index.

    # Elimination

    for k in range(n-1):

        if doPricing:

            # Pivot

            maxindex = abs(A[k:,k]).argmax() + k

            if A[maxindex, k] == 0:


                raise ValueError("Matrix is singular.")

            # Swap

            if maxindex != k:

                A[[k,maxindex]] = A[[maxindex, k]]

                b[[k,maxindex]] = b[[maxindex, k]]

        else:

            if A[k, k] == 0:

                raise ValueError("Pivot element is zero. Try setting doPricing to True.")

       #Eliminate

       for row in range(k+1, n):

           multiplier = A[row,k]/A[k,k]

           A[row, k:] = A[row, k:] - multiplier*A[k, k:]

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

    # Back Substitution

    x = np.zeros(n)

    for k in range(n-1, -1, -1):

        x[k] = (b[k] - np.dot(A[k,k+1:],x[k+1:]))/A[k,k]

    return x



if __name__ == "__main__":
    A = np.round(np.random.rand(100, 100)*10)
    b =  np.round(np.random.rand(100)*10)
    print (GEPP(np.copy(A), np.copy(b), doPricing = False))

Answer 1

To compare the two solutions, use np.allclose which verifies that two arrays agree, element by element, within some reasonably margins. One should not expect exact equality in floating point arithmetics.

my_solution = GEPP(np.copy(A), np.copy(b), doPricing=False)
numpy_solution = np.linalg.solve(A, b)
if np.allclose(my_solution, numpy_solution):
    print("Yay")
else:
    print("NOOOOO")

This prints "Yay".