 function A=piv(Ain,row,col)
%function A=piv(Ain,row,col)
%   pivot element (row,col)
[m,n]=size(Ain);
if Ain(row,col)<=0,
   disp('Warning: pivot element not positive');
end
arow=Ain(row,:);
pivot=Ain(row,col);
cols=setdiff(1:n,col);
rows=setdiff(1:m,row);
A=zeros(size(Ain));
for i=rows
for j=cols
    A(i,j)=Ain(i,j)-Ain(i,col)*Ain(row,j)/pivot;
end
end
A(row,:)=arow/pivot;
