function P = Givens(A, p,q, id)
% _____________________________________________________________________
% Obtain the Givens matrix P=P(p,q) for A=PAP (2-sided id=2) Hessenberg
%                                    or A=P*A (1-sided id=1) QR factors
% Usage :   P=givens(A,p,q, id)   for id=2 (Both side) or id=1 QR
% E.g.
%  A1 = rand(3,2);  A_orig = A1; %%%%%%%%% Case of id=1
%    P = givens(A1,1,2,1), A1=P*A1, Q=P;
%    P = givens(A1,1,3,1), A1=P*A1, Q=P;
%    P = givens(A1,2,3,1), R =P*A1, Q=Q*P, A=A_orig
%  A = rand(4); A=A+A'; A_orig=A; %%%%%%%%% Case of id=2
%    P = givens(A,2,3,2), A=P*A*P, Q=P;
%    P = givens(A,2,4,2), A=P*A*P, Q=P;
%    P = givens(A,3,4,2), T=P*A*P, Q=Q*P; A = A_orig,  spy( abs(A)>1.0e-12 )

if nargin<1, help givens, return, end
m = size(A,1); P = eye(m);
if id==1, j = p; else
          j = p-1;
   if j<1,p,j, error(['p = ' num2str(p) ' is too small']); end, end
if p>=q, q, error(['p = ' num2str(p) ' should be smaller than q']); end
   t = norm(imag(A),1)>eps;  %% 1 for complex or 0 for real only
   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

ratio = A(q,j)/A(p,j); %% Lower/Upper

  if id==2 & p==1, error('		Wrong p q pairs'); end
  if q-p<=0, error('		Wrong p q pairs'); end

if t==0  %% Real case
theta = atan( ratio );
 sina = sin(theta); cost = cos(theta);  sinb=sina;

else %% Complex case
 bot = sqrt( 1+norm(ratio)^2 );
 sina = conj(ratio)/bot;       sinb = ratio /bot;  
  cost = 1/bot;
end

            P(p,p)=cost; P(p,q)=sina; P(q,p)=sinb; P(q,q)=-cost;
if p>=m-1
  fprintf('givens.m --------- Last Step (%d,%d) for column %d ', p,q,j)
  if id==2
   fprintf('of Q^T*A*Q = %s or A=Q*%s*Q^T\n',S,S)
   else
   fprintf('of Q^T*A = R or A=Q*R\n')
  end
else
  fprintf('givens.m --------- Step (%d,%d) for column %d ', p,q,j)
  if id==2
   fprintf('of Q^T*A*Q = %s or A=Q*%s*Q^T\n',S,S)
   else
   fprintf('of Q^T*A = R or A=Q*R\n')
  end
end
