function [b,d,e,c]=ThomasGivens(c,d,e,b);
% THOMASGIVENS solves a tridiagonal system of linear equations
%   [b,d,e,c]=ThomasGivens(c,d,e,b) solves a tridiagonal linear system
%   using Givens rotations. The coefficient matrix is
%   A=diag(c,-1)+diag(d)+diag(e,1), and the right hand side b is
%   overwritten with the solution. The R factor is also returned,
%   R=diag(d)+diag(e,1)+diag(c,2).

n=length(d);
e(n)=0;
for i=1: n-1                                         % elimination
  if c(i)~=0 
    t=d(i)/c(i); si=1/sqrt(1+t*t); co=t*si;
    d(i)=d(i)*co+c(i)*si; h=e(i);
    e(i)=h*co+d(i+1)*si; d(i+1)=-h*si+d(i+1)*co;
    c(i)=e(i+1)*si; e(i+1)=e(i+1)*co;
    h=b(i); b(i)=h*co+b(i+1)*si;
    b(i+1)=-h*si+b(i+1)*co;
  end;
end;
b(n)=b(n)/d(n);                                      % backsubstitution
b(n-1)=(b(n-1)-e(n-1)*b(n))/d(n-1); 
for i=n-2:-1:1,
  b(i)=(b(i)-e(i)*b(i+1)-c(i)*b(i+2))/d(i);
end;
