title: Implementing Gaussian elimination with partial pivoting
body: <p>I am writing a program to implement Gaussian elimination with partial pivoting in MATLAB. I created an integer array to store the interchange of rows, instead of directly exchanging the rows.</p>
  
  <p>However, I could not obtain the correct result and I could not figure out the problem. How should I modify my code to get the right answer?</p>
  
  <pre><code>A=[1, -2, 1; 2, 1, -3; 4, -7, 1];
  b=[0, 5, 1];
  b=b';
  
  function x = GE(A,b)
  
  [m,n]= size(A); 
  if m ~= n
    disp('Not a square');
  end
  
  for p=1:n
     array(p)=p;
  endfor
  
  A = [A,b];      
  
  %elimination
  for i = 1:n-1
      pivot = i;
      %select pivot
      for j = i+1:n
          if  abs(A(array(i),i)) &lt; abs(A(array(i),i)) %row interchange
              temp = array(i);
              array(i) = array(j);
              array(j) = temp;
          end
      end
  
      while (pivot &lt;= n &amp;&amp; A(pivot,i)== 0)
          pivot = pivot+1;
      end
      if pivot &gt; n
          disp('No unique solution');
          break
      else
          if pivot &gt; i
              tem = array(i);
              array(i) = pivot
              pivot= tem;
          end
      end
  
      for j = i+1:n
          m = -A(array(j),i)/A(array(i),i);
          for k = i+1:n+1 
              A(array(j),k) = A(array(j),k) + m*A(array(i),k);
          end
      end
  end
  
  if A(n,n) == 0
      disp('No unique solution');
      return
  end
  
  %backward substitution
  x(n) = A(array(n),n+1)/A(array(n),n);
  for i = n - 1:-1:1
      sum = 0;
      for j = i+1:n
          sum = sum + A(array(i),j)*x(j);
      end
      x(i) = (A(array(i),n+1) - sum)/A(array(i),i);
  end
  
  endfunction
  </code></pre>