function [V,H] = qriterR(A);
%
%   Performs QR iteration in real arithmetic
%
%   Input:  A -- a real n by n matrix
%
%   Output: V -- a real n by n orthogonal matrix
%
%           H -- a real n by n upper triangular matrix
%   
%                AV = VH
%                        
%                leading k by k block of H has eigenvalues of
%                largest real part
%
% 
%   D.C. Sorensen
%   2 Mar 2000
%
    format short e
    [n1,n1] = size(A);
%
%   Reduce A to Hessenberg form  AV = VH
%
    v = randn(n1,1);
    [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);
%
%            Decrease j by 2 if m_j imaginary and by 1 if real
%
             if (abs(imag(mu(j))) > 0),
                j = j-2;
             else
                j = j-1;
             end
       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);
       mu = diag(D);
%
%      Sort the eigenvalues according to largest real part
%
       [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)));
          for j = 1:n, 
              for i = 1:j,
                 ii = n-i+1; jj = n-j+1;
                 x(indx) = ii; y(indx) = j;
                 indx = indx + 1;
              end
              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,ii,'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
          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