function UD=ad_gauss(AD,la)
%performs a Gaussian elimination on compressed diagonal form of a square matrix
%AD contains diagonals columnwise, la is a pointer to the diagonals
%UD is the upper triangular matrix on compressed diagonal form.
[n,m]=size(AD);
UD=AD;
ma=max(la);
js=find(la==ma)+1;
is=length(la)-js+1;
for k=1:n-1
   for i=k+1:min(n,k+is)
      for j=i:min(n,k+ma)
         UD(i,j-i+1)=UD(i,j-i+1)-(UD(k,js+i-k-1)/UD(k,1))*UD(k,j-k+1);
      end
      for j=k+1:i-1
         UD(j,js+i-j-1)=UD(j,js+i-j-1)-(UD(k,js+i-k-1)/UD(k,1))*UD(k,j-k+1);
      end
   end
end
UD=UD(1:n,1:1+ma);
