function x = solLR(L,R,b)

 
% x = solLR(L,R,b)
% 
% Risolve il sistema
%
%  (L*R)*x = b
%
% dove L e' una matrice triangolare inferiore con diagonale unitaria
%      R e' una matrice triangolare superiore

                  
[n,m]=size(L);

if(n~=m), 
   error('La matrice deve essere quadrata'),
end

% Risolve il sistema triangolare inferiore

x = b(:);

for i=2:n

    x(i)=x(i)-L(i,1:i-1)*x(1:i-1);

end


% Risolve il sistema triangolare superiore

if (R(n,n)==0),
   error('La matrice R e'' singolare: impossibile proseguire'),
end

x(n) = x(n)/R(n,n); 

for i=n-1:-1:1

   if (R(i,i)==0),
      error('La matrice R e'' singolare: impossibile proseguire'),
   end

   x(i) = ( x(i) - R(i,i+1:n)*x(i+1:n) ) / R(i,i);

end    

return
