function x=EliminationBandMatrix(p,q,B,b);
% ELIMINATIONBANDMATRIX solves a linear system with a banded matrix
%   x=EliminationBandMatrix(p,q,B,b); solves the banded linear system
%   Ax=b with partial pivoting. The columns of B contain the nonzero
%   diagonals of the matrix A. The first q columns of B contain the
%   lower diagonals of A (augmented by leading zeros) the remaining
%   columns of B contain the diagonal of A and the p upper diagonals,
%   augmented by trailing zeros. The vector b contains the right-hand
%   side.

n=length(B);
B=[B,zeros(n,q)];                          % augment B with q columns
normb=norm(B,1);
for i=1:n
  maximum=0;                               % search pivot
  for k=i:min(i+q,n)
    if abs(B(k,i-k+q+1))>maximum,
      kmax=k; maximum=abs(B(k,i-k+q+1));
    end
  end
  if maximum<1e-14*normb;                  % only small pivots
    error('matrix is singular')
  end
  if i~=kmax                               % interchange rows
    h=B(kmax,i-kmax+q+1:min(n,i+2*q+p-kmax+1));
    B(kmax,i-kmax+q+1:min(n,i+2*q+p-kmax+1))=B(i,q+1:min(n,2*q+p+1));
    B(i,q+1:min(n,2*q+p+1))=h;
    h=b(kmax); b(kmax)=b(i); b(i)=h;
  end
  for k=i+1:min(n,i+q)                     % elimination step
    B(k,i-k+q+1)=B(k,i-k+q+1)/B(i,q+1);
  end 
  for k=i+1:min(n,i+q)
    b(k)=b(k)-B(k,i-k+q+1)*b(i);
    for j=i+1:min(n,i+p+q)
      B(k,j-k+q+1)=B(k,j-k+q+1)-B(k,i-k+q+1)*B(i,j-i+q+1);
    end
  end
end
for i=n:-1:1                               % back substitution 
  s=b(i);
  for j=i+1:min(n,i+q+p)
    s=s-B(i,j-i+q+1)*x(j);
  end 
  x(i)=s/B(i,q+1);
end 
x=x(:);
