%% QR-Verfahren
% Wir betrachten folgende Matrizen (siehe Bespiele aus dem Buch):
% A_1=\pmatrix{1 & 1 & 4  & 2 \cr 0 &7 &3 & 20\cr 2 & 7 & 3 & 12 \cr 2 & 28 & 15 & 20},  A_2=\pmatrix{17 & 24 & 1 & 8 &15 \cr
%  23 & 5 & 7 &14 & 16 \cr
%  4 & 6 & 13 &20 & 22 \cr
%  10 & 12 & 19 & 21 & 3 \cr
%  11 & 18 &25 & 2 & 9}, A_3=\pmatrix{2 & 3 & 4 & 5& 6 \cr
%  4 & 4& 5& 6 & 7 \cr
%  0 & 3 & 6 & 7 & 8 \cr
%  0 & 0& 2 & 8 & 9 \cr
%  0 & 0 & 0 & 1 & 10},
% A_4={\rm tridiag}(-1,2,-1) \in \mathbb{R}^{n\times n},   A_5: (n\times n)-Hilbertmatrix
% Diese fünf Matrizen haben alle ausschließlich reelle Eigenwerte.
% Bei den Beispielen 4 und 5 muss auch die Dimension festgelegt werden.
% Es wird eine der obenstehenden Matrizen gewählt:

choice=1;
n=10; % nur für Beispiele 4 und 5 relevant
A=Makematrix(choice,n)
Atest=A;

% Die Matrizen A_3, A_4 haben bereits obere Hessenberggestalt. Die Matrizen A_1, ~A_2 und A_5 müssen erst über eine orthogonale Ähnlichkeitstransformation auf Hessenbergform transformiert werden (siehe Buch). Dazu wird die Matlab-Routine hess verwendet:

if ((choice==1)||(choice==2)||(choice==5))
    disp('Hessenbergmatrix');
    A=hess(A)
    disp('------------------------------------');
    disp('------------------------------------');
end

%
%  Das QR-Verfahren wird auf die obere Hessenberg-Matrix A angewendet.
% Die Methode hat folgende Struktur:
% Bestimmung eines Verschiebungsparameters: Zur Bestimmung des Verschiebungsparameters \mu wird eine im Buch beschriebene Methode verwendet.
% QR-Iteration mit Verschiebung: Zur Bestimmung der QR-Zerlegung wird der Einfachheit halber die Matlab-Routine qr verwendet. Diese Vorgehensweise ist nicht effizient, weil die Hessenberg-Gestalt der Matrix nicht berücksichtigt wird (besser wäre eine Elimination der Nebendiagonaleinträge mit Givens-Rotationen).
% Deflation: Gestartet wird mit der n \times n-Matrix A_n:=A. Nach einigen QR-Iterationen ist der Nebendiagonaleintrag a_{n-1,n} hinreichend klein (|a_{n-1,n}| \leq 10^{-14}) und die linke obere (n-1)\times(n-1)-Teilmatrix A_{n-1} wird weiter bearbeitet, usw., bis eine 2 \times 2obere Dreiecksmatrix übrig bleibt.
% Ausgabe:
% Falls n \leq 10 gilt, wird  die m \times m-Matrix A_m (m=n,n-1,\ldots,2), mit Eintrag  (A_m)_{m,m-1} etwa Null, gezeigt (auf diese Matrix wird Deflation angewendet).
% Eigenwertapproximation: Annäherung eines Eigenwertes ( = Eintrag (A_m)_{m,m}).
% AnzahlQRIter: die benötigte Anzahl der QR-Iterationen, um die  Eigenwertapproximation (A_m)_{m,m} zu bestimmen.
% Schießlich werden alle (nach Größe sortierten) Eigenwertapproximationen \tilde \lambda_i, 1 \leq i \leq n, gezeigt.
% Zur Kontrolle wird abschließend die MATLAB-Routine eig zur Bestimmung der Eigenwerte der Matrix A verwendet. Die mit eig berechneten Eigenwerte \hat \lambda_iwerden mit den Eigenwerten \tilde \lambda_i verglichen,  und die maximale   Abweichung MaxAbweichung:=\max_{1 \leq i \leq n} |\tilde \lambda_i - \hat \lambda_i| wird berechnet.
%% Aktivität:
% Wählen Sie unterschiedliche Bespielmatrizen und variieren Sie die Dimension n. Wie verhalten sich die QR-Iterationszahlen AnzahlQRIter in jedem Deflationsschritt? Hängen diese Iterationszahlen (wesentlich) von der Dimension n ab?

n=size(A,1);
Eigenwert=zeros(n,1);
disp('QR-Verfahren:')
for  m=n:-1:2                   %Deflationsschleife
    fprintf('Matrix nach wiederholter QR-Iteration. Dimension der Matrix ist m=%d ',m);
    A=A(1:m,1:m);
    I=eye(m);
    AnzahlQRIter=0;
    while ((abs(A(m,m-1)) > 10^(-14)) && (AnzahlQRIter <50))
        AnzahlQRIter=AnzahlQRIter+1;
        mu=Shift(A(m-1:m,m-1:m)); % Shiftparameter wird bestimmt
        [Q,R]=qr(A-mu*I);         % QR-Zerlegung wird berechnet 
        A=R*Q+mu*I;               % QR-Iteration
    end
    %% Ausgabe von Ergebnissen
    if (n <= 10) 
        A 
    end
    AnzahlQRIter
    Eigenwert(m)=A(m,m);
    fprintf('Eigenwertapproximation:\t %8.4e\n', A(m,m));
    disp('------------------------------------');
end
Eigenwert(1)=A(1,1);
fprintf('Eigenwertapproximation:\t %8.4e\n', A(1,1));
disp('------------------------------------');
disp('------------------------------------');
fprintf('\nAlle Eigenwertapproximationen: \n');
format shorte;
sort(Eigenwert)
disp('------------------------------------');
fprintf('\nMaxAbweichung: \n');
TestGenauigkeit(Atest,Eigenwert)

%

function [A]=Makematrix(k,n);
  if (k==1)
      A=[1 1 4 2;
         0 7   3 20;
         2 7  3 12;
         2 28 15 20];
  end
  if (k==2)
      A=[17 24 1 8  15;
         23 5  7 14 16;
         4  6 13 20 22;
        10 12 19 21 3 ;
        11 18 25 2  9];
  end
  if (k==3)
      A=[2 3 4 5 6;
         4 4 5 6 7;
         0 3 6 7 8;
         0 0 2 8 9;
         0 0 0 1 10];
  end
   if (k==4)
     A=full(gallery('tridiag',n,-1,2,-1));
   end
   if (k==5)
     A=hilb(n);
   end
end

function [mu]=Shift(B)
     d=det(B);
     tr=trace(B);
     eig1=0.5*(tr+sqrt(tr^2-4*d)); eig2=0.5*(tr-sqrt(tr^2-4*d));
     realeig1=real(eig1); realeig2=real(eig2);
     if (abs(realeig1 -B(2,2)) <= abs(realeig2 -B(2,2)))
         mu=realeig1;
     else 
         mu=realeig2;
     end
end

function [error]=TestGenauigkeit(A,Eigenwert)
    n=size(A,1);
    E=eig(A);
    error=norm(sort(E)-sort(Eigenwert),inf);
end
    