function [T,S] = pvlform(H)
% [T,S]=PVLFORM(H) computes the PVL form of the matrix H
% H: Hamiltonian matrix
% T: PVL form of H
% S: orthogonal symplectic matrix such that H=STS'
n = size(H,1)/2;
T = H;
S = eye(2*n);
for k = 1:n-1
    U = symplectic_red_naive(T(:,k),k+1);
    T = U*T*U';
    S = S*U';