function [L,U] = luband(A,p,q)
% Algorithm of Crout: LU-factorization of band matrix A
% p = upper bandwidth, q = lower bandwidth 
n = size(A,1); if (size(A,2) ~= n), error('n ~= m'); end
L = eye(n); U = zeros(n,n); 
for k=1:n
  % Compute row of U (1+q nonzero entries only) 
  for j=k:min(n,k+q)
      U(k,j) = A(k,j) - L(k,1:k-1)*U(1:k-1,j); end
  % Compute column of L (p nonzero off-diagonal entries only) 
  for i=k+1:min(n,k+p) 
      L(i,k) = (A(i,k) - L(i,1:k-1)*U(1:k-1,k)) /U(k,k); end
end
        