function  [ Q1, V1 ] = RefinEig( Q, V, B )
%  [Q, V] = RefinEig( Q, V, B )  applies iterative refinement to try
%  to improve the accuracy of the eigenvectors  Q  and eigevalues  V ,
%  of a square matrix  B ,  previously delivered by  [Q  V] = eig(B) .
%  Prudence requires a check that the residual  [B, Q]*[Q; -V]  be not
%  much bigger after  RefinEig  than before,  since  RefinEig  could
%  worsen rather than improve accuracy if the eigenvalues are not well
%  separated or if the eigenvectors are too nearly linearly dependent.
%  Sometimes repetition of  [Q, V] = RefinEig( Q, V, B )  pays off.
%  This version of  RefinEig  works with  MATLAB  versions  3.5 & 4.2 .
%    RefinEig  works best on  ix86-based PCs  and  680x0-based Macs,
%    whose floating-point arithmetics accumulate matrix products with
%    11  extra bits of precision.         (C) W. Kahan,  18 Jan. 1996

I = eye(size(Q)) ;  u = diag(I) ;  n = length(u) ;
v = diag(V) ;  V = diag(v) ; % ...  ensures that  V  is diagonal.
if  all(all(B == B')) ,  % ...  deal first with  Hermitian  B .
   dQ = [Q', I]*[Q; -I]*0.5 ;  d = min([ norm(dQ,1), norm(dQ,inf)])^2 ;
   if  (0.25 - d) ~= 0.25 ,  % ... which is not uncommon, ...
			   if  (0.125 - d*d) == 0.125,  %  ...  the likliest case  ...
							   dQ = dQ - dQ*(1.5*dQ - 2.5*dQ*dQ) ;
							 else   % ...  Q  requires full re-orthogonalization:
          [dQ, dZ, Q] = svd(Q) ;  Q = dQ*Q' ;
          dQ = [Q', I]*[Q; -I]*0.5 ; end, end % ... small dQ = (Q'*Q-I)/2 .
   dQ = Q*dQ ;  % ...  P = Q - dQ  is accurately unitary.
   dH = [B, Q, dQ, B]*[Q; -V; V; -dQ] ; % ... = B*(Q-dQ) - (Q-dQ)*V .
   dH = [-dQ', Q']*[dH;dH] ;  dH = (dH+dH')*0.5 ; % ... = P'*(B*P-P*V) .
   dh = diag(dH) ;  y = v*0.5 ;  yt = y' ;  dx = dh*0.5 ;  dxt = dx' ;
   dZ = (dH-diag(dh)) ./ ((y(:,u') - yt(u,:)) + (dx(:,u') - dxt(u,:)) + I) ;
   dZ = dZ ./( sqrt(real(dZ.*conj(dZ)) + 1) + 1 ) ;
   dZ = dZ ./( sqrt(real(dZ.*conj(dZ)) + 1) + 1 ) ;  K = ~finite(dZ) ;
   if  any(any(K))
         a = 1/(1 + sqrt(2)) ;  A = triu(K)*a ;  A = A - A' ;
         dZ(K) = sign(dH(K)).*A(K) ;  end
   dX = dZ*dH ;  dX = ((dX' + dX) - dX*dZ) + dH ;  dx = real(diag(dX)) ;
   A = real(conj(dZ).*dZ) ;  a = norm(A, 1) ; % ...  A = abs(dZ).^2 .
   if (8192 + a*a) == 8192 ,   % ...  the normal case
            dv = dx - A*dx ; % ...  dv = (I+A)\dx .
      else
            [L, U] = lu(I+A) ;
            dv = U\(L\dx) ;  dv = dv - U\(L\([I, -I, A]*[dv; dx; dv])) ;
            k = ~finite(dv) ;
            if any(k),  dv(k) = dx(k) ;  end,  end
			dX = dX + dZ*diag(dv)*dZ ;  dX = (dX+dX')*0.5 ; 
   y = 2*v ;  yt = y' ;  dvt = dv' ;
   dZ = ( y(:,u') - yt(u,:) ) + ( dv(:,u') - dvt(u,:) ) + I ;
   dZ = ( dX - diag(diag(dX))) ./ dZ ;   a = norm(dZ, 1) ;
   if (8192 + a*a) == 8192 ,   % ...  the normal case
          dX = (dZ - dZ*dZ)*2 ;  %  ...  = 2*((I+dZ)\dZ) .
      else
          K = ~( abs(dZ) < 1024 ) ; % ... Will this work on  PC-Matlab 3.5 ?
          if  any(any(K))
               A = triu(K)*1024 ;  A = A - A' ;
               dZ(K) = sign(dX(K)).*A(K) ;  end
          dX = ( (I + dZ)\dZ )*2 ;  end
   dQ = [Q, -dQ, dQ]*[dX; dX; I] ;  Q1 = Q - dQ ;  V1 = diag(v + dv) ;
   return,  end  %  Subsequent lines deal with non-Hermitian  B .
dZ = [B, Q]*[Q; -V] ;  [L, U] = lu(Q) ;  dC = U\(L\dZ) ;  % ...  residual.
dC = dC - U\(L\([Q, -I]*[dC; dZ])) ;  %  ...  iteratively refines  dC .
dv = diag(dC) ;  % ...  1st order approx'n to eigenvalue corrections.
vt = v.' ;  dvt = dv.' ;
S = 0.5*( (vt(u,:) - v(:,u')) + (dvt(u,:) - dv(:,u')) ) ;  % ... skew.
dZ = dC - diag(dv) ;  Y = sqrt( dZ.*dZ.' + S.*S ) ;  % ... symmetric.
K = real(Y.*conj(S)) ;  K = (K < 0) + triu(K==0, 1) ; % ... 0's & 1's .
Y(K) = -Y(K) ;  % ... Turns  Y  skew.
S = S + Y ;
%  Because  MATLAB  bungles division of complex numbers by  Inf ,  division
%  by zero cannot be prevented simply by setting  S = S + Inf*(~S) .
dZ = dZ ./ (S+I) ; % ... exact for diag-sum of  2-by-2 's,  else  1st order.
K = ~finite(dZ) ;  if any(any(K)),  dZ(K) = I(K) ;  end %  NaN & Inf --> 0 .
dC = dC + dC*dZ ;  dvt = diag(dC).' ;  ...  2nd order approximation.
dV = diag(dvt) ;  S = ( vt(u,:) - v(:,u') ) + ( dvt(u,:) - dV ) + I ;
dZ =  (dC - dV) ./ S ;   % ...  2nd order approximation.
K = ~finite(dZ) ;  if any(any(K)) ,  dZ(K) = I(K) ;  end
dZ = Q*dZ ;
for  j = 1:n , % ... Normalize column  j  of  Q1  just as  eig  should.
  q = Q(:,j) ;  z = dZ(:,j) ; ...  q+z  is improved eigenvector.
  q1 = [q', -1]*[q; 1] ; % ... = q'*q - 1 .
  d = (z'*z + 2*real(q'*z)) + q1 ; % ... = (q+z)'*(q+z)^2 - 1.
  d = d/(1 + sqrt(1+d)) ;  dq = ( z - q*d )/(1 + d) ;
  if  dq'*dq*64 < 1 + q1 , % ... additive normalization is best:
        Q1(:,j) = q + dq ;
     else ,  % ... normalization by division may be more accurate:
        q = q+z ;  Q1(:,j) = q/sqrt( q'*q ) ;  end,  end  % ... dq, j
V1 = diag(v + dvt.') ;  %  End  RefinEig.m