function [H,V]=IFRANCIS(H,V,f,l)
% [H,V]=IFRANCIS(A,V,f,l)
% H=IFRANCIS(A,f,l)
%
% QR-Schritt mit implizitem Francis-Shift für reelle Hessenbergmatrizen.
% Dabei wird eine (l-f+1)x(l-f+1)-Matrix Q aus dem QR-Schritt für H(f:l,f:l)
% bestimmt. Mit Hilfe dieser Matrix wird H mit seinem neuen Wert H_neu=W'*H*W
% überschrieben. Dabei gilt W(l:f,l:f)=Q und sonst ist W überall die Identität.
% Bei Bedarf wird V mit V_neu=V*W aufdatiert.
% 
% -----------------------------------------------------------------------------
%
% Input
% -----
% l,f           1<=f<=l<=n. Block H(f:l,f:l) von H der durch QR-Schritt
%               bearbeitet wird 
% 
% -----------------------------------------------------------------------------
%
% Input/Output
% ------------
% H             nxn Obere Hessenbergmatrix
% V             optional, Transformationsmatrix


% Bei drei Parametern werden nur die ersten drei (H,V,f) übergeben. Wir
% interpretieren sie aber als (H,f,l)
if nargin==3
   l=f;
   f=V;
   clear V;
end % if

n=size(H,1);
v=f:l;


% -----------------------------------------------------------------------------
% impliziter Francis-Shift
% Berechne Spur und Determinante des unteren 2x2 Blocks
% Spur der Matrix H(l-1:l,l-1:l)
s=H(l-1,l-1)+H(l,l);
% Determinante der Matrix H(l-1:l,l-1:l)
t=H(l-1,l-1)*H(l,l)-H(l-1,l)*H(l,l-1);
% -----------------------------------------------------------------------------

% Nur bei Blöcken der Dimension >=3 macht der Francis-Shift Sinn
if f+2<=l
   % Berechne die erste Spalte von (H-lambda I)(H-lambda' I)==(H^2-sH+tI), 
   % diese hat nur drei Eintraege x,y,z
   x=H(f,f)*H(f,f)+H(f,f+1)*H(f+1,f)-s*H(f,f)+t;
   y=H(f+1,f)*H(f,f)+H(f+1,f+1)*H(f+1,f)-s*H(f+1,f);
   z=H(f+2,f+1)*H(f+1,f);
   
   % Bulge-chasing 
   % Alle bis auf auf den den letzten Schritt bearbeiten einen 2x2 Bulge
   v=f:n;
   w=1:min(f+3,l);
   for k=f-1:l-3
       [h,beta,mu]=HOUSE([x;y;z]);
       q=k+1:k+3;
       H(q,v)=ROWHOUSE(H(q,v),h,beta);
       H(w,q)=COLHOUSE(H(w,q),h,beta);

       if nargin==4
	  V(:,q)=COLHOUSE(V(:,q),h);	  
       end % if

       x=H(k+2,k+1);
       y=H(k+3,k+1);
       if k<l-3
	  z=H(k+4,k+1);
       end % if

       if k>=f
	  v(1)=[];
       end % if
       if k+5<=l
	  w=[w k+5];
       end
   end % for k
   % Letzter Schritt, nur noch ein einzelner Eintrag
   k=l-2;
   [c,s]=GIVENS(x,y);
   q=k+1:k+2;
   H(q,v)=ROWROT(H(q,v),c,s);
   H(w,q)=COLROT(H(w,q),c,s);

   if nargin==4
      V(:,q)=COLROT(V(:,q),c,s);
   end % if
   
elseif f+1==l 
   % Möglicherweise reelle Eigenwerte, versuche billigeren Wilkinson-Shift
   r=s*s-4*t;
   if r>=0
      if H(f,f)<=H(l,l)
         mu=(s+sqrt(r))/2;
      else
         mu=(s-sqrt(r))/2;
      end % if-else
      % Berechne die erste Spalte von (H-mu I)
      % diese hat nur zwei Eintraege x,y
      x=H(f,f)-mu;
      y=H(f+1,f);

      
      % 1 Schritt Bulge-chasing 
      [c,s]=GIVENS(x,y);
      v=f:n;
      w=1:l;
      q=f:l;
      H(q,v)=ROWROT(H(q,v),c,s);
      H(w,q)=COLROT(H(w,q),c,s);
      if nargin==4
	 V(:,q)=COLROT(V(:,q),c,s);
      end % if
   end % if r>=0
end % if-elseif f+2<=l







