addpath('../routines/');

% Größe der Testmatrizen
v=10:10:50;
% Anzahl Testmatrizen pro Dimension
anz=5;


j=1;
fl=zeros(2,length(v));
it=zeros(1,length(v));
fehler=zeros(2,length(v));
abstand=zeros(1,length(v));

fprintf(' n |      flops      |Fehler, Zerlegung| max. Abstand, Eigenw.\n');
fprintf('   |  QRALG |  schur |  QRALG |  schur |                      \n');
fprintf('--------------------------------------------------------------\n');
for n=v
    for i=1:anz
        A=randn(n); 

	
	% QR-Algorithmus, implizite Variante, Francis-Shift
	flops(0);
	[m,V,H]=QRALG(A,'IFRANCIS');
	fl(1,j)=fl(1,j)+flops; 
	it(j)=it(j)+m;
	fehler(1,j)=fehler(1,j)+norm(V*H*V'-A,1)/norm(A,1);

	% Eigenwerte aus Blockdiagonalmatrix bestimmen
	ewH=zeros(n,1);
	i=1;
	while i<=n
	      if i<n
		 if H(i+1,i)==0
		    ewH(i)=H(i,i);
		    i=i+1;
		 else
		    ewH(i:i+1)=eig(H(i:i+1,i:i+1));
		    i=i+2;
		 end
	      else
		 ewH(n)=H(n,n);
		 i=i+1;
	      end
        end % while i
    
	
	% Schur-Zerlegung aus MATLAB
	flops(0);
	[U,T]=schur(A);
	fl(2,j)=fl(2,j)+flops; 
	fehler(2,j)=fehler(2,j)+norm(U*T*U'-A,1)/norm(A,1);
    
	% Eigenwerte aus Blockdiagonalmatrix bestimmen
	ewT=zeros(n,1);
	i=1;
	while i<=n
	      if i<n
		 if T(i+1,i)==0
		    ewT(i)=T(i,i);
		    i=i+1;
		 else
		    ewT(i:i+1)=eig(T(i:i+1,i:i+1));
		    i=i+2;
		 end
	      else
		 ewT(n)=T(n,n);
		 i=i+1;
	      end
        end % while i
    

	% beste Paarung der Eigenwerte ewH, ewT
	for i=1:n
	    [a,k]=min(abs(ewT(i)-ewH(i:n)));
	    k=k+i-1;
	    ewH([i k])=ewH([k i]);
	end % i
	abstand(j)=abstand(j)+norm(ewT-ewH,'inf')/max(abs(ewT));
    
    end % for i
    fl(:,j)=fl(:,j)/anz;
    fehler(:,j)=fehler(:,j)/anz;
    abstand(j)=abstand(j)/anz;
    it(j)=it(j)/anz;
    fprintf('%3d|%8.1e|%8.1e|%8.1e|%8.1e| %8.1e\n',n,...
	    fl(1,j),fl(2,j),fehler(1,j),fehler(2,j),abstand(j));
    j=j+1;
      
end % for n

figure(1);
loglog(v, fl(1,:), '+', v, fl(2,:),'o');
title('Rechenzeitvergleich QRALG/schur');
xlabel('n, `+` QRALG, `o` schur');
ylabel('flops');

figure(2);
loglog(v, fehler(1,:), '+', v, fehler(2,:),'o');
title('Vergleich der Residuum QRALG/schur');
xlabel('n, `+` QRALG, `o` schur');
ylabel('Fehler');

figure(3);
loglog(v, abstand(1,:), 'xr');
title('Abstand der berechneten Eigenwerte');
xlabel('n');
ylabel('maximaler Abstand');
