function [T,U,V]=urv(H)
% [T,U,V]=URV(H) computes the URV decomposition of the Hamiltonian matrix H
% H: a Hamiltonian matrix
% T: URV factor
% U,V: orthogonal symplectic matrices such that H=UTV'
n = size(H,1)/2;
U = eye(2*n);
V = eye(2*n);
T = H;
for k=1:n-1
    S = symplectic_red_naive(T(:,k),k);
    U = U*S';
    T = S*T;
    v = T(n+k,:)';
    S = symplectic_red_naive([v(n+1:2*n);-v(1:n)],k+1); % S(Jv,k+1)
    V = V*S';
    T = T*S';
end
S = symplectic_red_naive(T(:,n),n);
U = U*S';
T = S*T;