function [u,b]=ge2piv(a,b)
%routine to triangulate matrix a into upper triangular matrix u
%and modify loadvector b accordingly.
%with row pivoting.
n=length(a);
for k=1:n-1
  %apply row pivoting if necessary to equation k and below
    [pivot,ind]=max(abs(a(k:n,k))); % determine row with max value on col k
    if ind>1   %ind=index of row with max coefft relative to current pivot
        t1=a(k,:) ; %store coeffts of equation k
        a(k,:)=a(k+ind-1,:); %copy over to k, coeffts from equation (k+ind-1)
        a(k+ind-1,:)=t1; %restore coeffts of equation k to equation (k+ind-1)
        %ditto for the load vector
        t2=b(k);
        b(k)=b(k+ind-1);
        b(k+ind-1)=t2;
    end
   for i=k+1:n
      v(i)=a(i,k)/a(k,k);
      b(i)=b(i)-v(i)*b(k);
      for j=k+1:n
         a(i,j)=a(i,j)-v(i)*a(k,j);
      end
   end
end
u=triu(a);
