clear; clc;

%% Erzeugen einer SPD Matrix
N  = 2;
N1 = N+1;
n  = N1*N1;
d0 = -ones(n-N1,1);
d1 = 4*ones(n,1);
for (i=1:N1)
    d1(i*N1) = 3;
end
d2 = -ones(n-1,1);
for (j=1:N1-1)
    d2(j*N1) = 0;
end

D  = diag(d1) + diag(d2,1) + diag(d2,-1) + diag(d0,N1) + diag(d0,-N1);
A1 = N^2*D + eye(N1^2);                                       % A1 symmetrisch

for i=1:n                                                   % Alle Unterdeterminanten positiv 
    det( A1(1:i, 1:i));                                       % -> Matrix positiv definit  
    if (det (A1(1:i,1:i)) > 0)
        'true';
    else
        'false';
    end
end

%% Matrix auf tridiagonale Form bringen
    
H = A1;
Q = eye(n);

for (k=1:n-2)
    u        = zeros(n-k,1);
    u(1:n-k) = H(k+1:n,k);
    t        = norm(u,inf);
    u        = u./t;
    t        = sign(u(1))*norm(u);
    u(1)     = u(1) + t;
    t        = t*u(1);
    P        = eye(n);
    for (i=1:n-k)
        for (j=1:n-k)
            P(k+i,k+j) = P(k+i,k+j) - u(i)*u(j)/t;
        end
    end
    Q = P*Q;
    H = P*H*P';
 end
 Q = Q';
 QHQ = Q*H*Q';            % Vergleiche Matrizen A1 und QHQ                 

norm(A1-QHQ)

 A = H;

  %% QR Algorithmus zur Berechnung der Eigenwerte

 a = zeros(n,1);
 for (i=1:n)
    a(i) = A(i,i);      % Hauptdiagonale von A
 end

 b = zeros(n-1,1);
 for (i=1:n-1)
    b(i) = A(i,i+1);    % Nebendiagonale(n) von A
 end

 kmax = 10000000;
 x = zeros(n,1);
 y = x;
 g = x;
 c = x;
 s = x;
 d = x;
 e = x;
 for k=1:kmax
    x = a(1);
    y = b(1);
    for j=1:(n-1)
        g(j)   = sqrt(x(j)^2 + b(j)^2);
        c(j)   = x(j) / g(j);
        s(j)   = b(j) / g(j);
        d(j)   = c(j)*y(j) + s(j)*a(j+1);
        x(j+1) = -s(j)*y(j) + c(j)*a(j+1);
        if j ~= (n-1)
            e(j)   = s(j)*b(j+1);
            y(j+1) = c(j)*b(j+1);
        end
    end
    g(n) = x(n);
    a(1) = s(1)*d(1) + c(1)*g(1);
    b(1) = s(1)*g(2);
    for (j=2:n-1)
        a(j) = s(j)*d(j) + c(j)*c(j-1)*g(j);
        b(j) = s(j)*g(j+1);
    end
    a(n) = c(n-1)*g(n);
    if (norm(b) < 0.00001)
        break
    end
 end

 eig_A = eig(A);       % Vergleiche selbst berechnete Werte mit Eigenwertfunktion von Matlab

norm(sort(a(:))-sort(eig_A(:)))
