function   Aout = h95v(A);
%%  function   A = h95v(A);
%%  To obtain a matrix in Hessenberg form similar to the matrix A
%%  by repeated application of the Householder transformation
%%  INPUT
%%    A - a nxn matrix  (does not need to be symmetric)
%%  OUTPUT
%%    Aout -  a matrix with similar to A and with zeros below its
%%            subdiagonal. A and Aout have the same eigenvalues.
%%  Verbose mode gives the Householder transformation and the Aout for
%%    each stage

[n,m]=size(A);Aout=A;
for k=1:n-2; %%Step 1
   %% Step 2
   q=norm(Aout(k+1:n,k),2);  %% first set up the choice of alpha
   %% Step 3
   if Aout(k+1,k)==0
      alpha=q;
   else
      alpha=-q*sign(Aout(k+1,k));
   end;
   %% Step 4
   rsq=alpha.^2-alpha*Aout(k+1,k); %% find 2r^2
   %% Step 5
   v(k,1)=0;  
   v(k+1,1)=Aout(k+1,k)-alpha;     
   v(k+2:n,1)=Aout(k+2:n,k);  %% note w=v/2r
   fprintf('the Householder transform at step %g is\n',k)
   householder=eye(n)-v*v'/rsq
   pause
   u(1:n,1)=(1/rsq)*Aout(1:n,k+1:n)*v(k+1:n,1); %%Step 6
   y(1,1:n)=(1/rsq)*v(k+1:n,1)'*Aout(k+1:n,1:n); %%Step 6
   prod= v(k+1:n,1)'*u(k+1:n,1); %%Step 7
   z(1:n,1)=u(1:n,1)-prod*v(1:n,1)/(2*rsq); %%Step 8
   for l=k+1:n; %Step 9
      Aout(1:k,l)=Aout(1:k,l)-u(1:k,1)*v(l,1);
      Aout(l,1:k)=Aout(l,1:k)-y(1,1:k)*v(l,1);
      Aout(k+1:n,l)=Aout(k+1:n,l)-2*v(l,1)*z(k+1:n,1)-y(1,l)*v(k+1:n,1);%Step 11
      Aout(k+1:n,l)=Aout(k+1:n,l)+v(l,1)*u(k+1:n,1);
   end;
Aout
pause
   clear u v y z;
end;

