Implementing a partial pivot If condition for Gaussian Elimination in Python-CodePudding

I've figured out how to code a function such that you can row reduce and solve linear algebra problems, the only issue I'm running into is setting up the if condition that does the partial pivoting when the constant that is used to row reduce is equal to 0. I've attempted it below and the logic makes sense but am struggling to understand why my solution doesn't work.

import numpy as np

def gaussElim(A,B):
    M = np.concatenate((A,B), axis=1)# Combines the two matrices (assuming they have the same ammount of rows and matrix B has been )
    nr, nc  = M.shape
    for r in range (nr):
        const = M[r][r]
        if const == 0: # **This is the condition that is tripping me up**
            for i in range (nr-1):
                M[r][i]=M[r 1][i]
                M[r 1][i] = M[r][i]
                const = M[r][r]
        for c in range (r,nc):
            M[r][c] = M[r][c]/const
        for rr in range(nr):
            if rr!= r :
                const = M[rr][r]
                for c in range(r,nc):
                    M[rr][c] = M[rr][c] -  const * M[r][c]
    return M[:, nc-1]



Mrx = np.array([ [1.0,3,2,4,3,1], [-4,0,3,2,3,4], [3,-1,3,2,2,5], [3,3,12,2,-
6,-4], [-1,-2,-3,7,6,4], [7,5,0,0,4,2] ]) 
Rhs = np.array([[ 4, 5, 6, 10, 6, -8 ]]) # This is a row vectorv 
RhsT = Rhs.T # Rhs.T is transpose of Rhs and a column vector 
S = gaussElim(Mrx,RhsT)  
print(S)  
 

A1 = np.array([[2.0,1,-1],[2,1,-2],[1,-1,1]]) 
b1 = np.array([[1.0],[-2],[2]]) 
S1 = gaussElim(A1,b1) 
print (S1)
x = np.linalg.solve(A1,b1)
print(x)

Yes I have looked at other people's solutions to Gaussian elimination but I want to understand why more specifically my solution to partial pivoting isn't working. Thank you!

Imputing A1, b1 into my function gives me

[1, 0.8, 1.8]

The correct answer is

[1, 2, 3]

The print statements were so I could see if my function was working as intended.

CodePudding user response：

One issue is this code:

            for i in range (nr-1):
                M[r][i]=M[r 1][i]   # line 1
                M[r 1][i] = M[r][i] # line 2
                const = M[r][r]     # line 3

The line I've commented as line 2 simply undoes the work of line 1. If you're attempting to swap values, try replacing lines 1 and 2 with this:

                M[r][i], M[r 1][i] = M[r 1][i], M[r][i]

... or, if you prefer to really be explicit about the swap:

                temp = M[r][i]
                M[r][i]=M[r 1][i]   # line 1
                M[r 1][i] = temp    # line 2

A second (probably benign) issue in your code is that line 3 above (const = M[r][r]) does not need to be inside the loop as it currently is, and I believe you can dedent one level with no change in the ultimate result.

A third (benign) issue is that M[r][c] = M[r][c]/const can be simplified to M[r][c] /= const, assuming the original arrays A and B (and thus M) are floats.

Similarly, M[rr][c] = M[rr][c] - const * M[r][c] can be simplified to M[rr][c] -= const * M[r][c].

Putting it all together (most importantly, correcting the first issue above):

import numpy as np

def gaussElim(A,B):
    M = np.concatenate((A,B), axis=1)# Combines the two matrices (assuming they have the same ammount of rows and matrix B has been )
    nr, nc  = M.shape
    for r in range (nr):
        const = M[r][r]
        if const == 0: # **This is the condition that is tripping me up**
            for i in range (nr-1):
                M[r][i], M[r 1][i] = M[r 1][i], M[r][i]
            const = M[r][r]
        for c in range (r,nc):
            M[r][c] /= const
        for rr in range(nr):
            if rr!= r :
                const = M[rr][r]
                for c in range(r,nc):
                    M[rr][c] -= const * M[r][c]
    return M[:, nc-1]



Mrx = np.array([ [1.0,3,2,4,3,1], [-4,0,3,2,3,4], [3,-1,3,2,2,5], [3,3,12,2,-
6,-4], [-1,-2,-3,7,6,4], [7,5,0,0,4,2] ]) 
Rhs = np.array([[ 4, 5, 6, 10, 6, -8 ]]) # This is a row vectorv 
RhsT = Rhs.T # Rhs.T is transpose of Rhs and a column vector 
S = gaussElim(Mrx,RhsT)  
print(S)  
 

A1 = np.array([[2.0,1,-1],[2,1,-2],[1,-1,1]]) 
b1 = np.array([[1.0],[-2],[2]]) 
S1 = gaussElim(A1,b1) 
print (S1)
x = np.linalg.solve(A1,b1)
print(x)

Output:

[-0.97124946  2.88607869 -1.80198876  3.47492434 -6.16688284  4.51794207]
[0.33333333 1.33333333 1.        ]
[[1.]
 [2.]
 [3.]]