 function x = gaussp (A, b)
%-----------------------------
% function x = gaussp(A, b) 
% solves A x=b by Gaussian 
% elimination with partial pivoting
%-----------------------------
 n = size(A,1) ;
 A = [A,b] 
 pause 
for k=1:n-1 
  [t,ip] = max(abs(A(k:n,k))); 
  ip = ip+k-1;
  disp('swapping:')
  disp([ip, k]) 
%% swap if needed -  
  if (ip ~= k)
    temp = A(k,k:n+1) ;
    A(k,k:n+1) = A(ip,k:n+1); 
    A(ip,k:n+1) = temp;
    A
   end
%%
  for i=k+1:n
      piv = A(i,k) / A(k,k) ;
%%    A(i,k+1:n+1)=A(i,k+1:n+1)-piv*A(k,k+1:n+1);
      A(i,k:n+1)=A(i,k:n+1)-piv*A(k,k:n+1);
  end
 A
 pause 
 end 
 x = backsolv(A,A(:,n+1)); 
