Gaussian elimination in numerical

Question

In the following code I have implemented Gaussian elimination without partial pivoting for a general square linear system Ax = b.

However I am looking for some help with implementing the following two requirements,

1) I want to make sure that my function terminates if a zero pivot is encountered.

2) And I want to test it by solving Ax = b where A is a random 100x100 matrix and b is a random 100x1 vector.

Looking for some help with adding these two things into the code I have produced! thanks!

 def linearsolver(A,b):
  n = len(A)
  M = A

  i = 0
  for x in M:
   x.append(b[i])
   i += 1

  for k in range(n):
   for i in range(k,n):
     if abs(M[i][k]) > abs(M[k][k]):
        M[k], M[i] = M[i],M[k]
     else:
        pass

   for j in range(k+1,n):
       q = float(M[j][k]) / M[k][k]
       for m in range(k, n+1):
          M[j][m] -=  q * M[k][m]

  x = [0 for i in range(n)]

  x[n-1] =float(M[n-1][n])/M[n-1][n-1]
  for i in range (n-1,-1,-1):
    z = 0
    for j in range(i+1,n):
        z = z  + float(M[i][j])*x[j]
    x[i] = float(M[i][n] - z)/M[i][i]
  print (x)

Answer 1

I'm answering this here anyways. People will say you should answer things on Stackoverflow if it is primarily coding, however, mathematical questions on StackOverflow are instantly overlooked.

Gaussian Elimination without Pivoting

import numpy as np
import math


def forward_elimination(A, b, n):
    """
    Calculates the forward part of Gaussian elimination.
    """
    for row in range(0, n-1):
        for i in range(row+1, n):
            factor = A[i,row] / A[row,row]
            for j in range(row, n):
                A[i,j] = A[i,j] - factor * A[row,j]

            b[i] = b[i] - factor * b[row]

        print('A = \n%s and b = %s' % (A,b))
    return A, b

def back_substitution(a, b, n):
    """"
    Does back substitution, returns the Gauss result.
    """
    x = np.zeros((n,1))
    x[n-1] = b[n-1] / a[n-1, n-1]
    for row in range(n-2, -1, -1):
        sums = b[row]
        for j in range(row+1, n):
            sums = sums - a[row,j] * x[j]
        x[row] = sums / a[row,row]
    return x

def gauss(A, b):
    """
    This function performs Gauss elimination without pivoting.
    """
    n = A.shape[0]

    # Check for zero diagonal elements
    if any(np.diag(A)==0):
        raise ZeroDivisionError(('Division by zero will occur; '
                                  'pivoting currently not supported'))

    A, b = forward_elimination(A, b, n)
    return back_substitution(A, b, n)

Generates random matrix with given condition number

def gen_rand_matrix(m,n,kappa):
#
#
    k= min(m,n)
    (Q1,R1) = np.linalg.qr(np.random.rand(m,m))
    (Q2,R2) = np.linalg.qr(np.random.rand(n,n))

    U = Q1[0:m,0:k]
    V = Q2[0:k,0:n]
    j=k-1

    l = kappa**(1/j)
    x1 = np.arange(0,j+1)
    x2 = np.flip(x1,axis=0)
    sing = np.power(l,x2)
    S = np.diag(sing)
    Vt = np.transpose(V)

    M = np.dot(np.dot(U,S),Vt)
    return M

This part here is testing

kappa = 1/(math.exp(1e-3)-1)
m=100

A  = gen_rand_matrix(m,m,kappa)
x = np.random.rand(m,1)
b = np.dot(A,x)
xhat1 = gauss(A,b)

xerror1 = np.linalg.norm(xhat1-x)/np.linalg.norm(x)

Matlab Implementation Of Partial Pivoting From a Numerical Methods Book

There is an implementation here