function [V,R,ritz] = Iram_Cnv(A,eigA,which,k,m,f);
%
%
%   Input:  A     -- an n by n matrix
%
%           eigA  -- eigenvalues of A
%
%           which -- a string indicating which eigenvalues are
%                      wanted.  Options are:
%
%                         'LR', 'SR', 'LM','LA','SA'
%
%           k    -- a positive integer 
% 
%           m    -- a positive integer  k < m << n
%
%                     Recommended: m > 2k
%   
%           f    -- an nonzero n vector
%                  
%                     If f does not appear, f <- randn(n,1) will occur
% 
%
%   Output: V    -- an n by k orthogonal matrix
%
%           R    -- a k by k quasi upper triangular matrix
%
%           ritz -- a vector containing Ritz estimates for the 
%                   k Ritz values (eigenvalues of R);
%
% 
%           with   norm(AV - VR)  ~  norm(ritz) 
%
%           Hence, a partial real Schur form of order k
%
%
%   D.C. Sorensen
%   2 March 2000
%  
    nmax = 30;
    ph1 = zeros(nmax,1);
    ph2 = zeros(nmax,1);

    plot(real(eigA),imag(eigA),'k+');
    hold on;

    [n,n] = size(A);
    if (nargin < 5), f = randn(n,1); end
    
    V = zeros(n,m); H = zeros(m,m);

    tol = sqrt(eps); ritz = 1; nitr = 0;

    ko = 1; 

    while (norm(ritz) > tol & nitr < nmax),
       nitr = nitr+1;

       [V,H,f] = Arnold(A,V,H,f,ko,m);

       [mu, ritz]= select_shifts(H,which);

       Q = eye(m);
       j = m;
       while(j > k)
%
%            Apply unwanted eigenvalues as shifts with QR steps
%
             [Q,H] = qrstep(Q,H,mu(j),1,m);
%
%            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
       ko = j;
       
       yvec = norm(f)*ritz;
       ritz = norm(f)*ritz(1:ko);
       f = V*Q(:,ko+1)*H(ko+1,ko) + f*Q(m,ko);
       V(:,1:ko) = V*Q(:,1:ko);

       Ritz_est = ritz'
       Nitrs = nitr

             

%
%      ---------------------Graphics to Display Convergence-----------------
%
          
   
          ph1(nitr) = plot(real(mu(1:ko)), imag(mu(1:ko)),'o');
          set(ph1(nitr),'MarkeredgeColor',[1 1 1])

          ph2(nitr) = plot(real(mu(ko+1:m)), imag(mu(ko+1:m)),'o');
          set(ph2(nitr),'MarkeredgeColor',[1 1 1])
   end

   mlt = 4;
   colors = bluemap(mlt*nitr);

   for j = 1:nitr
      set(ph1(j),'Color',colors((mlt-1)*nitr+j,:))
%     set(ph1(j),'Color',colors(mlt*nitr-2,:))
      set(ph1(j),'MarkerSize',7)
      set(ph1(j),'MarkeredgeColor',[1 1 1])
      set(ph1(j),'MarkerFaceColor', colors((mlt-1)*nitr+j,:));
%     set(ph1(j),'MarkerFaceColor', colors(mlt*nitr-2,:));

      set(ph2(j),'Color',colors(j,:))
      set(ph2(j),'MarkerSize',7)
      set(ph2(j),'MarkeredgeColor',[1 1 1])
      set(ph2(j),'MarkerFaceColor', colors(j,:));
      pause(.5)
   end
%
%  Change basis to partial Schur form
%

   [Q,R] = schur(H(1:ko,1:ko));
   V = V(:,1:ko)*Q;
   mu = eig(R);
   ph = plot(real(mu), imag(mu),'ro');
   set(ph,'MarkerSize',7)

   hold off