function [V,H] = qriterRS(A);
%
%
%
%   Performs QR iteration in real arithmetic
%
%           Shifts are selected from "unwanted Ritz values
%           from the leading m by m block.
%
%   Input:  A -- a real n by n symmetric matrix
%
%   Output: V -- a real n by n orthogonal matrix
%
%           H -- a real n by n  tridiagonal matrix
%   
%                AV = VH
%
%           The leading k by k block contains the algebraically 
%           largest eigenvalues of A.
%
%           Here m = 20, k = 6
%
%           Note:  In real life, you would take advantage of symmetry.
%                  This code does not do that.
% 
%   D.C. Sorensen
%   2 Mar 2000
%
%
    format short e
    [n1,n1] = size(A);
    v = randn(n1,1);
%
%   Reduce A to tridiagonal form with full Arnoldi process
%

    [V,H] = ArnoldiC(A,n1,v);
    t1 = eig(A);
    [t2,ii]  = sort(-real(t1)); 
    t1 = t1(ii);

    m = 20; k = 6;
    t1 = t1(1:m);
%
%   Find eigenvalues of leading m by m block 
%   and sort acording to largest real part
%
    mu = eig(H(1:m,1:m));
    [t,ii] = sort(-real(mu));
    mu = mu(ii);
    iter = 1; ynorm = 1;
    while( iter < 25 & ynorm > 1000*eps)
       iter = iter + 1;
       m = 20; k = 6;
       j = m;
       while(j > k),
%
%            Apply unwanted eigenvalues as shifts with QR steps
%

             [V,H] = qrstep(V,H,mu(j),1,n1);
             j = j-1;
       end

%
%      Display some convergence information
%
       Hvals = H(1:j,1:j);
       [Y,D] = eig(Hvals);
       y = abs(Y(j,:)*H(j+1,j));
       y = y'
       ynorm = norm(y(1:j-1));
       Niter = iter
%
%      Compute eigenvalues of leading m by m block 
%
       Hvals = H(1:m,1:m);
       [Y,D] = eig(Hvals);
%
%      Sort the eigenvalues according to largest real part
%

       mu = diag(D);
       [t,ii] = sort(-real(mu));
       mu = mu(ii);
%
%      Display some convergence information
%

       yvec = abs(Y(m,:)*H(m+1,m));
       yvec = yvec(ii)';
       y_diff_mu_lam = [yvec abs(abs(mu) - abs(t1)) abs(mu) abs(t1)]
%
%      ---------------------Graphics to Display Convergence-----------------
%
          n = m;
          x = ones(n*(n+1)/2 , 1); y = x;
          phd = zeros(n-1,1);
          phy = zeros(n-1,1);
   
          m = 16;
          colors = bluemap(m);
          x = [-1  -1  n+2 n+2];
          y = [-1  n+2 -1  n+2];
   
          ph1 = plot(x,y,'o');
          axis tight;
          set(ph1,'MarkeredgeColor',[1 1 1])
          hold on;
          indx = 1;
          Hmax = max(abs(diag(Hvals,-1)));
          
          x = [1:n];
          y = [n:-1:1];
   
          for j = 1:n-1, 
              jj = min(n-1,j);
              i = min(n,jj+1);
              gval = max(eps, abs(Hvals(i,jj)/Hmax));
              gval = -log(gval);
              gval = max(1,ceil(gval));
              gval = min(16,gval);
              gval = 16 - gval + 1;
   
              yval = max(eps, abs(yvec(jj)));
              yval = -log(yval);
              yval = max(1,ceil(yval));
              yval = min(16,yval);
              yval = 16 - yval + 1;
   
              ii = n-i+1; jj = j; 
              phy(j) = plot(jj,0,'o');
              phd(j) = plot([jj,jj+1],[ii, ii+1],'o');
              axis tight;
              set(phd(j),'Color',colors(m,:))
              set(phd(j),'MarkerSize',9)
              set(phd(j),'MarkeredgeColor',[1 1 1])
   
              set(phy(j),'Color',colors(m,:))
              set(phy(j),'MarkerSize',9)
              set(phy(j),'MarkeredgeColor',[1 1 1])
   
              if(abs(Hvals(i,j)) > 0),
                 set(phd(j),'MarkerFaceColor', colors(gval,:));
              end
              if(abs(yvec(j)) > 100*eps),
                 set(phy(j),'MarkerFaceColor', colors(yval,:));
              end
          end
          x = [x n]; y = [y 0];
          ph = plot(x,y,'o');
          axis tight;
          set(ph,'Color',colors(m,:))
          set(ph,'MarkerSize',9)
          set(ph,'MarkeredgeColor',[1 1 1])
          set(ph,'MarkerFaceColor',colors(m,:))
          x1 = [ .7 k+.3];
          y1 = [n-k-.3   n+.3];
          x2 = [ .7  .7];
          x3 = [ k+.3  k+.3];
          y2 = [n-k-.3   n-k-.3];
          y3 = [n+.3   n+.3];
          plot(x1,y2,'m--',x1,y3,'m--',x2,y1,'m--',x3,y1,'m--');
          axis tight;
          pause(.2)
          hold off
       end