function [l,u,p]=ge3piv(a)
%routine to triangulate matrix a into upper and lower triangular matrices l and u.
% permutation matrix p used for row pivoting.
n=length(a);
l=zeros(n);
p=eye(n);
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)
        t1=p(k,:) ; %store coeffts of equation k
        p(k,:)=p(k+ind-1,:); %copy over to k, coeffts from equation (k+ind-1)
        p(k+ind-1,:)=t1; %restore coeffts of equation k to equation (k+ind-1)
    end
   for i=k+1:n
      v(i)=a(i,k)/a(k,k);
      a(i,k)=v(i);
      for j=k+1:n
         a(i,j)=a(i,j)-v(i)*a(k,j);
      end
   end
end
u=triu(a);
l=tril(a,-1)+diag(ones(n,1));
