Gauss elimination with pivoting in python

Question

This is my code with gauss algorithm and it seems work in "def MEG", but in "def MEGPivoting" when i try to exchanges my rows with the gauss algorithm with pivoting and i print the new matrix, it printed the seems matrix like the original. Why? Is there a problem with the for cicle or it's another problem? why the new matrix is the copy of the first?

import numpy as np
import time

#Eliminazione di Gausss senza pivoting
def MEG(A,b):
    #calcolo dimensione vettore
    n = len(A)
    #ciclo fino alla penultima colonna (dove non eliminiamo nulla)
    for j in range(n-1):
        #ciclo per le righe (partendo dalla seconda)
            for i in range(j+1, n):
                #calcolo del coeff m
                m = A[i][j]/A[j][j]
                A[i][j] = 0
                for k in range(j+1,n):
                    #formula Gauss
                    A[i][k] = A[i][k] - m*A[j][k]
                    #anche per i termini b
                b[i] = b[i] - m*b[j]
    return A
 
#Costruisco l' algoritmo
def MEGPivoting(A,b):
    n = len(A)
    for j in range(1,n-1):
        #individuo l' elemento pivot
        amax = abs(A[j][j])
        imax = j
        for i in range (j+1, n):
            if abs(A[i][j]>amax):
                amax = abs(A[i][j])
                imax = i
        #eventuale scambio di riga
        if imax > j:
            for k in range (j, n):
            #scambio
                A[j][k],A[imax][k] = A[imax][k],A[j][k]
                b[j],b[imax]=b[imax],b[j]

        print('matrice con riga scambiate')
        print(A)
        
        #metodo di eliminazione di Gauss
        MEG(A,b)
        return  A 
    
#Creo la matrice
n = 5
M = 10
A = np.random.rand(n,n)*2*M - M 
print('Matrice = ')
print(A)
xsol = np.ones(n) #vettore incognite tutti 1
inizio = time.time()
b = np.dot(A,xsol) #vero vettore incognite ricavato da A*xsol
x = MEGPivoting(A,b)
fine = time.time()
#tempo totale
tempo = fine-inizio
print('\n\nSoluzione con pivoting = ')
print(x)
print('\n\nTempo impiegato: %f'%tempo)

Answer 1

Python does the operations in a kind of pass by address, so when you pass the array A to the MEGPivoting() method you have to do a copy operation inside. At the beginning of the method, when calculating n = len(A) also try to insert the command A = np.copy(A) . This will serve precisely to make a copy of the matrix, so that when you carry out the operations the original matrix A will not be modified, the operations will be performed on its copy. You also need to add it to the MEG() method. Obviously you must also make a copy of the vector of the known terms b with the command b = np.copy(b).

The code should be something like this:

#Eliminazione di Gausss senza pivoting
def MEG(A,b):
   A = np.copy(A) # copia della matrice A
   b = np.copy(b) # copia del vettore dei termini noti b
   # calcolo dimensione vettore
   n = len(A)
   #ciclo fino alla penultima colonna (dove non eliminiamo nulla)
   for j in range(n-1):
       #ciclo per le righe (partendo dalla seconda)
           for i in range(j+1, n):
               #calcolo del coeff m
               m = A[i][j]/A[j][j]
               A[i][j] = 0
               for k in range(j+1,n):
                  #formula Gauss
                  A[i][k] = A[i][k] - m*A[j][k]
                  #anche per i termini b
               b[i] = b[i] - m*b[j]
    return A


def MEGPivoting(A,b):
   A = np.copy(A) # copia della matrice A
   b = np.copy(b) # copia del vettore dei termini noti b
   n = len(A)
   for j in range(1,n-1):
       # pivot
       amax = abs(A[j][j])
       imax = j
       for i in range (j+1, n):
           if abs(A[i][j]>amax):
               amax = abs(A[i][j])
               imax = i
       # scambio di riga
       if imax > j:
           for k in range (j, n):
           #scambio
               A[j][k],A[imax][k] = A[imax][k],A[j][k]
               b[j],b[imax]=b[imax],b[j]

       print('matrice con riga scambiate')
       print(A)
    
       # di eliminazione di Gauss
       MEG(A,b)
   return  A