function [Q,H,nqr]=scriptqralg(A,demo)
% [Q,H,nqr]=scriptqralg(A)
% [Q,H,nqr]=scriptqralg(A,demo)
%
% Berechne reelle Schurform zu A.
%
% --------------------------------------------------------------------------------------------
% Input
% A           reelle nxn Matrix
% demo        demo=0  wirkt als wenn man [Q,H]=scriptqralg(A) aufruft
%             demo~=0 erfordert das h"aufige Tastendr"ucken, damit die Demo weiterl"auft
%
% Output
% Q, R       Schurform A=Q*R*Q'
% nqr        Anzahl QR-Iterationsschritte

if nargin==1
   demo=0;
end
[Q,H]=hess(A);

figure(1)
spy(H,'*k')
title('Hessenbergform')
if demo
   pause
else
   pause(0.5)
end

n=size(A,1);
k=0;
l=n;
flag=-1;
nqr=0;
while flag
   % Deflation
   m=0;
   E=sparse(n,n);
   for i=l-1:-1:k+1
       if abs(H(i+1,i))<=eps*(abs(H(i,i))+abs(H(i+1,i+1)))
	  H(i+1,i)=0;
	  E(i+1,i)=1;
	  if m<=1
	     m=0;
	     l=i;
	  else % m>1
	     k=i;
	  end
       else
	  m=m+1;
       end
   end
   if m<=1 
      m=0;
      for i=k:-1:1
	  if H(i+1,i)==0 
	     if m<=1
		m=0;
		l=i;
	     else % m>1
		k=i;
	     end
	  else
	     m=m+1;
	  end
      end
   end
   if m<=1
      flag=0;
   end
   figure(1)
   spy(H,'*k')
   hold on
   spy(E,'or')
   hold off
   title('Deflation')
   pause(0.5)

   figure(1)
   spy(H,'*k')
   title('Deflation')
   if demo
      pause
   else
      pause(0.5)
   end
   
   if flag
      % QR-Iterationsschritt, impliziter Doppelshift
      % H(k+1:l,k+1:l) mindestens 3x3
      t=H(l-1,l-1)+H(l,l);
      d=det(H(l-1:l,l-1:l));
      
      b=H(k+1:k+3,k+1:k+2)*H(k+1:k+2,k+1);
      b(1:2)=b(1:2)-t*H(k+1:k+2,k+1); 
      b(1)=b(1)+d;
      [Q1,R1]=qr(b);
      H(k+1:k+3,k+1:n)=Q1'*H(k+1:k+3,k+1:n);
      H(1:min(k+4,l),k+1:k+3)=H(1:min(k+4,l),k+1:k+3)*Q1;
      Q(:,k+1:k+3)=Q(:,k+1:k+3)*Q1;

      figure(1)
      spy(H,'*k')
      hold on
      E=sparse(n,n);
      E(k+1:k+3,k+4:n)=1;
      spy(E,'*b')
      E=sparse(n,n);
      E(1:k,k+1:k+3)=1;
      if l-k>3
	 E(k+4,k+3)=1;
      end
      spy(E,'*r')
      E=sparse(n,n);
      E(k+1:k+3,k+1:k+3)=1;
      spy(E,'*m')
      E=sparse(n,n);
      E(k+3,k+1)=1;
      if l-k>3
	 E(k+4,k+1:k+2)=1;
      end
      spy(E,'*g')
      spy(E,'og')
      hold off
      title('Ein Schritt QR-Zerlegung')
      if demo
	 pause
      else
	 pause(0.5)
      end
      
      for i=k+1:l-3
	  [Q2,R2]=qr(H(i+1:i+3,i));
	  
	  H(i+1:i+3,i:n)=Q2'*H(i+1:i+3,i:n);
	  H(1:min(i+4,l),i+1:i+3)=H(1:min(i+4,l),i+1:i+3)*Q2;
	  H(i+2:i+3,i)=0;
	  Q(:,i+1:i+3)=Q(:,i+1:i+3)*Q2;

	  figure(1)
	  spy(H,'*k')
	  hold on
	  E=sparse(n,n);
	  E(i+1:i+3,i+3:n)=1;
	  E(i+1,i)=1;
	  spy(E,'*b')
	  E=sparse(n,n);
	  E(1:i,i+1:i+3)=1;
	  if l-i>3
	     E(i+4,i+3)=1;
	  end
	  spy(E,'*r')
	  E=sparse(n,n);
	  E(i+1:i+3,i+1:i+3)=1;
	  spy(E,'*m')
	  E=sparse(n,n);
	  E(i+3,i+1)=1;
	  if l-i>3
	     E(i+4,i+1:i+2)=1;
	  end
	  spy(E,'*g')
	  spy(E,'og')
	  hold off
	  title(['Buckelschieben, Iterationsschritt ' num2str(nqr+1)])
	  if demo
	     pause(0.5)
	  else
	     pause(0.01)
	  end
      end

      [Q3,R3]=qr(H(l-1:l,l-2));
	  
      H(l-1:l,l-2:n)=Q3'*H(l-1:l,l-2:n);
      H(1:l,l-1:l)=H(1:l,l-1:l)*Q3;
      H(l,l-2)=0;
      Q(:,l-1:l)=Q(:,l-1:l)*Q3;

      figure(1)
      spy(H,'*k')
      hold on
      E=sparse(n,n);
      E(l-1:l,l+1:n)=1;
      E(l-1,l-2)=1;
      spy(E,'*b')
      E=sparse(n,n);
      E(1:l-2,l-1:l)=1;
      spy(E,'*r')
      E=sparse(n,n);
      E(l-1:l,l-1:l)=1;
      spy(E,'*m')
      E=sparse(n,n);
      hold off
      title(['Buckelschieben, Iterationsschritt ' num2str(nqr+1)])
      if demo
	 pause
      else
	 pause(0.5)
      end
      nqr=nqr+1;
   end
end