clear all;
close all;
% Extrem primitive Version des Jacobi-Verfahrens,
% nur zur Veranschaulichung des Prinzips.
n=5; 
A=2*rand(n,n)-1; % Zufallsmatrix 
A=A+A'          % symmetrisieren
ei=eig(A)       % wir merken uns das MATLAB-Ergebnis
Fnormsq=sum(sum(A.*A)') % Frobeniusnorm zum Quadrat, invariant
NA=Fnormsq-sum(diag(A).*diag(A)) % Nebendiagonal-Quadratsumme
while NA>1.0e-12 % Jacobi-Schleife 
    for j=1:n    % ueber alle unteren Elemente der Matrix A 
        for i=1:j-1;
            % ab hier GENAU wie im Buch, S. 281
            tau=(A(i,i)-A(j,j))/sqrt(4*A(i,j)^2+(A(i,i)-A(j,j))^2);
            c=sqrt((1+tau)/2);
            s=sign(A(i,j))*sqrt((1-tau)/2);
            % ab hier SEHR primitiv, mit kompletten Matrizen
            B=eye(n);
            B(i,i)=c;
            B(j,j)=c;
            B(i,j)=s;
            B(j,i)=-s;
            NA=NA-2*A(i,j)^2; % Update der Quadratsumme     
            A=B*A*B';         % Transformation
        end 
    end
end
A
NA