title: How to implement LU decomposition with partial pivoting in Python?
body: <p>I want to implement my own LU decomposition P,L,U = my_lu(A), so that given a matrix A, computes the LU decomposition with partial pivoting. But I only know how to do it without pivoting.
  Can anyone help to do the partial pivoting?</p>
  
  <pre><code>def lu(A):
  
      import numpy as np
  
      # Return an error if matrix is not square
      if not A.shape[0]==A.shape[1]:
          raise ValueError("Input matrix must be square")
  
      n = A.shape[0] 
  
      L = np.zeros((n,n),dtype='float64') 
      U = np.zeros((n,n),dtype='float64') 
      U[:] = A 
      np.fill_diagonal(L,1) # fill the diagonal of L with 1
  
      for i in range(n-1):
          for j in range(i+1,n):
              L[j,i] = U[j,i]/U[i,i]
              U[j,i:] = U[j,i:]-L[j,i]*U[i,i:]
              U[j,i] = 0
      return (L,U)
  </code></pre>