A=[1 2 1
   2 2 3
   1 3 6];
%A = [3 -2 1 4; -2 -6 2 -1; 1 2 -2 5; 4 -1 5 -7];

W = diag([1 1 1 ]);
for i = 1:11
    Q = diag([1 1 1]);
    [p,q] = maxmat(tril(A,-1)+triu(A,1))
    eta = (A(q,q)-A(p,p))/(2*A(p,q))
    t = sign(eta)/(abs(eta)+sqrt(eta^2+1))
    c = 1/sqrt(t^2+1)
    s = t*c
    Q(p,p) = c;
    Q(q,q) = c;
    Q(p,q) = s;
    Q(q,p) = -s;
    disp(Q)
    A = Q'*A*Q
    W = W*Q
end

%x3 = W(:,3)./W(3,3);