function [lambda,V]=tridiewev(a,b)
% [lambda,V]=tridiewev(a,b);
% computes all eigenvalues lambda and the corresponding 
% eigenvectors 
% V(:,i) of a symmetric unreduced tridiagonal matrix
% with main diagonal a and subdiagonal b
% using bisection and a rudimentary form a invers iteration
% namely solving 
% P(T-lambda(i)*I) =L*R
% R*V(:,i)= R(k,k)*(0,...,0,1)' and V(k,i)=1 if R(k,k)=0
% where R(k,k) is the first under the smallest pivot, seen from the
% bottom
% 
% 
  n=length(a);
  lambda=zeros(n,1);
  V=zeros(n);
  b(n)=0;
  for i=1:n-1
      if b(i)==0
          error('tridiewev: subdiagonal element=0');
      end
  end
  siz=size(a);
  if siz(1)==1
      a=a';
  end
  siz=size(b);
  if siz(1)==1
      b=b';
  end
  q=zeros(n,1);
  A=zeros(n,4);
  % gerschgorin circles
  low=min([a(1)-abs(b(1)),min(a(2:n-1)-abs(b(1:n-2))-abs(b(2:n-1))), ...
          a(n)-abs(b(n-1))]);
  up=max([a(1)+abs(b(1)),max(a(2:n-1)+abs(b(1:n-2))+abs(b(2:n-1))), ...
          a(n)+abs(b(n-1))]);
  % iteration
  for i=1:n
      %determine lambda(i) 
      if i>1
          low=lambda(i-1);
      end
      lowl=low;
      upl=up;
      while abs(upl-lowl) > eps*(max(abs(upl),abs(lowl))+eps)
          mu=(lowl+upl)/2;
          count=0;
          q(1)=a(1)-mu;
          if q(1)==0 
              q(1)=eps;
          end
          if q(1) <0 
              count=1;
          end
          for k=2:n
              q(k)=a(k)-mu-b(k-1)^2/q(k-1);
              if k<n && q(k)==0
                  q(k)=eps;
              end
              if q(k) < 0
                  count=count+1;
              end
          end
          if count < i
              lowl=mu;
          else
              upl=mu;
          end
      end
      lambda(i)=mu;
      % set matrix A, banded storage scheme
      A(1,1)=0; % undefined value
      A(2:n,1)=b(1:n-1);
      A(1:n,2)=a-mu;
      A(1:n-1,3)=b(1:n-1);
      A(n,3)=0;
      A(1:n,4)=zeros(n,1); 
      %additional  superdiagonal in 
      % case of row interchange
      % elimination with pivots for a tridiagonal, 
      % written in full form
      % for i=1:n-1
      %  if abs(A(i+1,i)) > abs(A(i,i)) 
      %    for k=0:2
      %       term=A(i+1,i+k);
      %       A(i+1,i+k)=A(i,i+k);
      %       A(i,i+k)=term;
      %    end
      %   end
      %   A(i+1,i)=A(i+1,i)/A(i,i);
      %   for k=1:2
      %     A(i+1,i+k)=A(i+1,i+k)-A(i+1,i)*A(i,i+k);
      %   end
      % end
      % indextransformation A(i,j)->  A(i,2-(i-j)) , j=i-1,i,i+1,i+2
      for j=1:n-1
          if abs(A(j+1,1)) > abs(A(j,2))
             for k=0:2
                 term=A(j+1,k+1);
                 A(j+1,k+1)=A(j,k+2);
                 A(j,k+2)=term;
             end
         end
         A(j+1,1)=A(j+1,1)/A(j,2);
         for k=1:2
             A(j+1,k+1)=A(j+1,k+1)-A(j+1,1)*A(j,k+2);
         end
      end
      % upper triangular matrix contained in the three columns
      % A(i,j),j=2,4
      k=n;
      pivmin=abs(A(n,2));
      for j=n-1:-1:1
          if abs(A(j,2)) < pivmin 
              k=j;
              pivmin=abs(A(j,2));
          end
      end
      V(k,i)=1;
      if k>1 
        V(k-1,i)=-A(k-1,3)/A(k-1,2);
      end
      for j=k-2:-1:1
         V(j,i)=-(A(j,3)*V(j+1,i)+A(j,4)*V(j+2,i))/A(j,2);
      end
      V(:,i)=V(:,i)/norm(V(:,i));
  end