% reduce to Hessenberg

function A = hessenberg(A)

[m, m1] = size(A);


for k=1:m-2
    
    x = A(k+1:m,k);
    v = x;
    v(1) = v(1) + sign(x(1))*norm(x);
    v = v/norm(v);
    
    A(k+1:m,k:m) = A(k+1:m,k:m) - 2*v*(v'*A(k+1:m,k:m));
    
    A(1:m,k+1:m) = A(1:m,k+1:m) - 2*(A(1:m,k+1:m)*v)*v';
    
end

