% tridia.m
% fd, 10 novembre 2006. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  resolution de systeme lineaire pour une matrice tridiagonale 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function xx=tridia(n,aa,bb,cc,sm) 
% parametres d entree  n    ordre du systeme lineaire 
%                      aa   diagonale de la matrice  
%                      bb   diagonale du dessus  
%                      cc   diagonale du dessous  
%                      sm   second membre  
% parametres de sortie xx   solution du systeme lineaire  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
al = aa ; be = bb ; ga = cc ;    %  utilitaires matriciels 
yy = sm ;                        %  un tableau tampon   
%
%                        decomposition L*U
al(1,1:1) = aa(1,1:1) ; 
for j=2:n ; 
be(1,j:j) = bb(1,j:j) ; 
ga(1,j:j) = cc(1,j:j) / al(1,j-1:j-1) ; 
al(1,j:j) = aa(1,j:j) -  be(1,j:j) * ga(1,j:j) ; 
end ; 
%                        descente
yy(1,1:1) = sm(1,1:1) ; 
for j=2:n ; 
yy(1,j:j) = sm(1,j:j) - ga(1,j:j) * yy(1,j-1:j-1) ;
end ; 
%                        remontee
xx(1,n:n) = yy(1,n:n) / al(1,n:n) ; 
for k=1:n-1 ; j = n-k ; 
xx(1,j:j) = ( yy(1,j:j) - be(1,j+1:j+1) * xx(1,j+1:j+1)) / al(1,j:j) ; 
end ; 
resi = sm ; 
for j=1:n ; resi(1,j:j) = resi(1,j:j) -  aa(1,j:j) * xx(1,j:j) ; end ; 
for j=1:n-1 ; 
resi(1,j:j) = resi(1,j:j) -  bb(1,j+1:j+1) * xx(1,j+1:j+1) ; end ;
for j=2:n ; 
resi(1,j:j) = resi(1,j:j) - cc(1,j:j) * xx(1,j-1:j-1) ; end ;
zero = 0 ; 
for j=1:n ; zero = max (zero,  resi(1,j:j)**2 ) ; end ; 
zero = zero /  n ; zero = sqrt (zero) ; 
if (zero > 10**(-12)) zero , end ; 
 

 