function P1()
clc
randn('seed',0)
j = 60;
C = randn(j,j);
u = randn(j,1);
disp('You:')
[Z,H,alpha] = SpecialHess(C,u);
errZ = norm(Z'*Z - eye(j,j));
errH = norm(H - triu((Z'*C*Z),-1));
erru = norm(Z'*u - [zeros(j-1,1);alpha]);
fprintf('alpha                 = %20.15f\n',alpha)
fprintf('norm(Z''*Z - I)        = %20.15f\n',errZ)
fprintf('norm(Z''*A*Z - H)      = %20.15f\n',errH)
fprintf('norm(Z''*u - alpha*ej) = %20.15f\n',erru)
disp(' ')
disp('CVL:')
[Z,H,alpha] = SpecialHessCVL(C,u);
errZ = norm(Z'*Z - eye(j,j));
errH = norm(H - triu((Z'*C*Z),-1));
erru = norm(Z'*u - [zeros(j-1,1);alpha]);
fprintf('alpha                 = %20.15f\n',alpha)
fprintf('norm(Z''*Z - I)        = %20.15f\n',errZ)
fprintf('norm(Z''*A*Z - H)      = %20.15f\n',errH)
fprintf('norm(Z''*u - alpha*ej) = %20.15f\n',erru)



function [Z,H,alpha] = SpecialHessCVL(C,u)
% C is a j-by-j matrix and u is a nonzero column j-vector
% Z is a j-by-j orthogonal matrix with the property that Z'*CZ  = H
%   is upper Hessenberg, Z'*u = alpha*I(:,j) where I = eye(j,j).

% Suppose P'*u = alpha*I(:,1). 
% If M is a j-by-j matrix and we use hess to compute the decomposition 
% then Q'*(P'*M*P)Q = H1, then (P*Q)'*u = alpha*I(:,1). If E = I(:,j:-1:1)
%  then              
%                  (PQE)'M(PQE) = E*H1*E is lower hessenberg and 
%  and
%                  (PQE)'*u = alpha*I(:,j)
% If M = C' then
%                  (PQE)*C*(PQE) = E*H1'*E is upper hessenberg
% and
%                  (PQE)'*u = alpha*I(:,j)
% Thus, Z = PQE.

  j = length(u);
% Compute P = I - 2vv' so P*u = alpha*I(:,1)
  v = u; 
  if u(1)>0
      v(1) = v(1) + norm(u); v = v/norm(v); alpha = -norm(u);
  else
      v(1) = v(1) - norm(u); v = v/norm(v); alpha = norm(u);
  end
% Compute C1 = P'*C'*P..
  C1 = C'; C1 = C1 - (2*v)*(v'*C1); C1 = C1 - (C1*v)*(2*v)';
% Compute Hessenberg decomposition of C1...
  [Q,H1] = hess(C1);
  H = H1'; H = H(j:-1:1,j:-1:1);
% Form Z..
  Z = Q - (2*v)*(v'*Q);  %P*Q
  Z = Z(:,j:-1:1);       %P*Q*E
  
 
% A more elegant solution is to apply hess to [0 u'; u C']
% and recognize that Q(2:j+1,2:j+1) and H(2:j,2:j), suitably permuted,
% have what you want.

