function [lambda,v,resnorm,flag,history,historyend]=rqi(a,b,x0,itermax)
% [lambda,v,resnorm,flag]=rqi(a,b,x0) 
% performs the rqi - iteration beginning with
% initial vector x0 for a unreduced symmetric tridiagonal 
% matrix with diagonal a(1:n) and b(1:n-1)
% should b contain more than n-1 elements those get destroyed
% [a(1),b(1),0,...,0;b(1),a(2),b(2),0,....a(n)];
% lambda and  v are the final eigenvalue obtained and 
% the corresponding eigenvector, resnorm is the 2-norm
% of A*v-lambda*v. 
% v has length one.
% iteration is terminated if either A-lambda*I is
% singular to working precision (with flag=1) or 
% if more than itermax steps are required for this 
% with flag=0
% history(1:historyend,1:2) stores the rayleighquotient
% and the last pivot value of the total iteration
% 
% 
  n=length(a);
  if n==1 
    lambda=a(1);
    v(1)=1;
    resnorm=0;
    flag=1;
    return;
  end
  x=x0;
  siz=size(x);
  if siz(1)==1 
    x=x';
  end
  xnorm=norm(x);
  if xnorm == 0 
    error('rqi: initial vector is zero ');
  end
  x=x/xnorm;
  y=zeros(n,1);
  if length(b) < n-1 
    error('rqi: b as not the proper length');
  end
  b(n)=0;
  for i=1:n-1
      if b(i)==0
          error('rqi: 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
  history=zeros(itermax,2);
  A=zeros(n,4);
  pivotbnd=eps*(norm(a,inf)+2*norm(b,inf));
  % iteration
  flag=0;
  iter=0;
  while iter < itermax
      iter=iter+1;
      %determine rho(i) 
      xs=[x(2:n);0];
      rho=sum(a.*(x.^2)+2*b.*x.*xs); % xnorm=1
      history(iter,1)=rho;
      % 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-rho;
      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))
          % do a row interchange
             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
             term=x(j+1);
             x(j+1)=x(j);
             x(j)=term;
         end
         A(j+1,1)=A(j+1,1)/A(j,2);  % the multiplier
         for k=1:2
             A(j+1,k+1)=A(j+1,k+1)-A(j+1,1)*A(j,k+2);
         end
         x(j+1)=x(j+1)-A(j+1,1)*x(j);
      end
      % upper triangular matrix contained in the three columns
      % A(i,j),j=2,4 , i=1,..,n
      history(iter,2)=A(n,2);
      if abs(A(n,2)) <= pivotbnd
        y(n)=1;
        y(n-1)=-A(n-1,3)/A(n-1,2);
        for j=n-2:-1:1
          y(j)=-(A(j,3)*y(j+1)+A(j,4)*y(j+2))/A(j,2);
        end
        y=y/norm(y);
        flag=1;
        x=y;
        historyend=iter;
        break
      else
        y(n)=x(n)/A(n,2);
        y(n-1)=(x(n-1)-A(n-1,3)*y(n))/A(n-1,2);
        for j=n-2:-1:1
          y(j)=-(A(j,3)*y(j+1)+A(j,4)*y(j+2))/A(j,2);
        end
        y=y/norm(y);
      end
      x=y;        
  end
  if flag==0
      historyend=iter-1;
  end
  v=x;
  lambda=rho;
  y(1)=a(1)*v(1)+b(1)*v(2)-rho*v(1);
  y(n)=a(n)*v(n)+b(n-1)*v(n-1)-rho*v(n);
  for i=2:n-1
    y(i)=b(i-1)*v(i-1)+a(i)*v(i)+b(i)*v(i+1)-rho*v(i);
  end
  resnorm=norm(y);