function [V,D,it]=Jacobi_val_12_2(A)
% Cálculo por Jacobi de valores y vectores propios de una MATRIZ SIMÉTRICA A
%
 tol=sqrt(eps)*norm(A,'fro');
 D=A; n=length(A); V=eye(n);
 [m1,p]=max(triu(abs(D),1));  % En (p,q) elemento mayor valor no en diagonal
 [~, q]=max(m1);              %    Posición fila máximo valor no cero en L(A)
 p=p(q); it=0;                %    Posición columna máximo anterior
 while off(D)>tol             % Proceos iterativo; necesita rutina off mia
    J=Jac_Rot(D,p,q);         % Se hacen cero Dpq y Dqp (p debe ser < q)
    D([p q],:)=J'*D([p q],:);
    D(:,[p q])=D(:,[p q])*J;
    V(:,[p q])=V(:,[p q])*J;
    [m1,p]=max(triu(abs(D),1));
    [~, q]=max(m1);
    p=p(q);
    it=it+1;
 end
 [D,I]=sort(diag(D)); V=V(:,I);
end

function J=Jac_Rot(A,b,d)
% Cálculo de rotación de Jacobi 2x2 en la matriz A.
%
 if A(b,d)~=0
    tau=(A(d,d)-A(b,b))/2/A(b,d);
    if tau>=0
       t=1/(tau+sqrt(1+tau^2)); 
    else
       t=-1/(-tau+sqrt(1+tau^2));
    end
    c=1/sqrt(1+t^2);
    s=c*t;
 else
    c=1; s=0; 
 end
 J=[c s; -s c];           % Igual que Givens
end

function a=off(A) 
% Calcula off de la matriz cuadrada A: raiz cuadrada de la suma 
%    de los cuadrados de los coefieicntes de A no en la
%    diagonal principal; también sqrt(sum(sum(triu(a.^2,1)))).
%
 n=length(A); a=0;
 for k=1:n-1
    a=a+sum(A(k,k+1:n).^2); 
 end
 a=sqrt(a);
end
