function P = Houses(A, c, id)
if nargin<1
   help houses; return
end
[m n]= size(A);
%  Generate the Householder matrix P=P_j
%
%  Usage:  P=houses(A,c, id)    for a n x n matrix A 
%                where 1<= c <= m is an integer (for column index) 
%                    id=1 for A=QR purpose,       p<=n-1 is fine.
%                                        (A can be m x n for id=1)
%                and id=2 for Householder method P'AP (NB c<=n-2)
% E.g.
%  A1 = rand(6,4);  A_orig = A1; %%%%%%%%% Case of id=1
%    P = houses(A1,1,1), A1=P*A1, Q=P;
%    P = houses(A1,2,1), A1=P*A1, Q=Q*P;
%    P = houses(A1,3,1), A1=P*A1, Q=Q*P;
%    P = houses(A1,4,1), R =P*A1, Q=Q*P, A=A_orig
%  A = rand(5); A=A+A'; A_orig=A; %%%%%%%%% Case of id=2
%    P = houses(A,1,2), A=P*A*P, Q=P;
%    P = houses(A,2,2), A=P*A*P, Q=Q*P;
%    P = houses(A,3,2), T=P*A*P, Q=Q*P, A = A_orig,  spy( abs(A)>1.0e-12 )
if nargin<3
   help houses; error('Need to input c and id ...')
end
if id ==2 
  if n ~= m, disp('Error in Houses.m: Matrix not square !'),
          Size_Input_Matrix = size(A),return,end
end
  if c>=n, disp(['This col ' num2str(c) ' is not required']); end
if c>n | c<1
   help houses; error('Check col no "c" ...');
end
   if size(A)==size(A')
   s = norm(A-A',1)>eps;     %% 1 for umsymme pr 0 for symmetry
   else
   s=1;
   end
if s==0, S='T'; else, S='H'; end

if id==2
   P = eye(n);   % For large n, speye(n) may be better... 
   j = c+1; 
else
   P = eye(m);   % For large m, speye(m) may be better... 
   j = c;
end;

      x=A(j:m,c); 
     s = norm(x); v=[]; 
     k = 1:j;
 v(k,1)=zeros(j,1);
   sig = sign(A(j,c));
v(j,1) = A(j,c)+sig*s ;
     k = j+1:m;
v(k,1) = A(k,c);

P = P - 2*(v*v')/(v'*v);
if j>=n-1 | (j>=n-1 & id==2)
  fprintf('houses.m --------- Last transform %d  ', c)
else
  fprintf('houses.m --------- ---- Col %d ', c)
end
  if id==2
   fprintf('of P^T*A*P = %s or A=P*%s*P^T\n',S,S) 
   else
   fprintf('of P^T*A = R or A=P*R\n')
  end
