[x, ns]=sweep(s,b)

function [x, ns]=sweep(s,b)
%SWEEP double balayage de Cholevski : solution de s*x=b pour s tridiagonale
%      [x, ns]=sweep(s,b)
% s    est la matrice du systme i.e.
%         soit une matrice carre creuse ou non
%         soit une matrice 3xn les colonnes 1, 2 et 3 sont alors
%         respectivement la sous diagonale s(1,1) == 0
%                        la diagonale principale
%                        la sur diagonale  s(3,n) == 0
%      remarque une matrice 3x3 est considre tre sous la deuxime forme
%      des que s(1,1) == 0.
%
% b    est une matrice de seconds membres, 1 par colonne.
%
% x    solutions du systme de mme taille que b
%
% ns   optionnel contient la matrice sous forme carre creuse ou non
%
%    Cette mthode est prfrable  l'oprateur \ si n > 500

if length(s)==0
   error('systme vide');
end;

[n,m] = size(s);
sdf = ( (n==m) & (n==3) & (s(1,1)==0) )  | (n~=m)
if sdf
   i = [1 2:n   ,1:n  1:n-1 3]';
   j = [3 1:n-1, 1:n, 2:n   1]';
   s = sparse( i, j, s(:));
   %  Ceci est equivalent 
   %  s = diag(s(2:n,1),-1) + diag(s(:,2))+  diag(s(2:n,3),+1)
   %  mais produit une matrice creuse.
end;
if nargout == 2
   ns = s;
end;
if abs(s(1,1)) == 0
    disp('chec balayage : pivot nul en premire ligne.');
end;
if length(s) == 1
   x = b./s;
   return;
end;

[n,m] = size(b);
x = zeros(size(b));
G = zeros(size(b));
B = s(1,1);
x(1,:) = b(1,:)./B;
for j = 2:n
   G(j) = s(j-1,j)./B;
   B = s(j,j)-s(j-1,j).*G(j);
   if B == 0
      disp(['chec balayage : pivot nul en ligne ', num2str(j)]);
   end;
   x(j,:) = (b(j,:)-s(j,j-1).*x(j-1))./B;
end;
for j=n-1:-1:1
   x(j,:)=x(j,:)-G(j+1).*x(j+1,:);
end;