title: generating random matrices in python
body: <p>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. However now I am trying to do the following but I am not quite sure how to code it, looking for some help with this! </p>
  
  <p><strong>I want to test my implementation by solving Ax=b where A is a random 100x100 matrix and b is a random 100x1 vector.</strong></p>
  
  <p>In my code I have put in the matrices </p>
  
  <pre><code>A = np.array([[3.,2.,-4.],[2.,3.,3.],[5.,-3.,1.]])
  
  b =  np.array([[3.],[15.],[14.]])
  </code></pre>
  
  <p>and gotten the following correct output:</p>
  
  <pre><code>[3. 1. 2.]
  [3. 1. 2.]
  </code></pre>
  
  <p><strong>but now how do I change it to generate the random matrices?</strong></p>
  
  <p>here is my code below:</p>
  
  <pre><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 &amp; 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.array([[3.,2.,-4.],[2.,3.,3.],[5.,-3.,1.]])
  
      b =  np.array([[3.],[15.],[14.]])
  
      print (GEPP(np.copy(A), np.copy(b), doPricing = False))
  
      print (GEPP(A,b))
  </code></pre>