%
% function [y,G,Psi,a]=hypgiv(n,p,q,x,Phi)
%
% =========================================================================
%
% Computation of the hyperbolic Givens transformation matrix
% 
% =========================================================================
% input data: n ... dimension of the vector x,
%      p,q (p<q)... components of x
%                   (x(p) will be changed, x(q) will be annihilated)
%             x ... the reduced vector
%           Phi ... the given diagonal matrix
%                   (the diagonal entries of Phi are +1 or -1)
% output data:
%             y ... vector x after the transformation
%             G ... matrix of hyperbolic Givens rotation
%                   (i.e. y=G*x)
%           Psi ... diagonal matrix Phi after the transformation
%                   (i.e. Psi=G'*Phi*G,
%                    G is hyperexcgange or hypernormal
%                         with respect to Phi)
%             a ... parameter:
%                   a=0 - the rotation is skipped
%                         (the p-th component of x is nonzero,
%                          the q-th one is already zero)
%                   a=1 - the classical Givens rotation
%                         (i.e. Phi(p,p)=Phi(q,q))
%                   a=2 - hyperbolic Givens rotation, |x(p)|>|x(q)|
%                   a=3 - hyperbolic Givens rotation, |x(p)|<|x(q)|
%                   a=4 - hyperbolic energy is zero,
%                         |x(p)|=|x(q)| and Phi(q,q)=-Phi(p,p)
%                   a=5 - the components are exchanged
%                         (x(p) is zero, x(q) nonzero)
%
%
function [y,G,Psi,a]=hypgiv(n,p,q,x,Phi)
y=x; G=eye(n); Psi=Phi;
%
if abs(x(p)*x(q))>eps
    if sign(Phi(p,p))==sign(Phi(q,q))
	 a=1;
	 ro=sqrt(abs(x(p))*abs(x(p))+abs(x(q))*abs(x(q)));
	 ro1=1/ro;
	 G(p,p)=abs(x(p))*ro1;
	 G(q,q)=G(p,p);
	 G(p,q)=(x(p)/abs(x(p)))*x(q)'*ro1;
	 G(q,p)=-G(p,q)';
	 y(p)=x(p)*ro/abs(x(p));
	 y(q)=0;
     else
         if abs(x(p))-abs(x(q))>eps
	   a=2;
	   ro=abs(x(p))*sqrt(1-(abs(x(q))*abs(x(q)))/(abs(x(p))*abs(x(p))));
	   G(p,p)=abs(x(p))/ro;
	   G(q,q)=G(p,p);
	   G(p,q)=-(x(p)/abs(x(p)))*(x(q)'/ro);
	   G(q,p)=G(p,q)';
	   y(p)=x(p)*ro/abs(x(p));
	   y(q)=0;
	  else
           if abs(x(q))-abs(x(p))>eps
	      a=3;
	      ro=abs(x(q))*sqrt(1-(abs(x(p))*abs(x(p)))/(abs(x(q))*abs(x(q))));
	      G(p,q)=abs(x(q))/ro;
	      G(q,p)=G(p,q);
	      G(p,p)=-(x(q)/abs(x(q)))*(x(p)'/ro);
	      G(q,q)=G(p,p)';
	      y(p)=x(q)*ro/abs(x(q));
	      y(q)=0;
	      Psi=G'*Phi*G;
	    else u=-x(p)*x(q)'/(abs(x(q))*abs(x(q)));
		  a=4;
		  G(p,p)=2;
		  G(p,q)=u;
		  G(q,p)=u'*sqrt(2);
                  G(q,q)=sqrt(2);
		  Psi=G'*Phi*G;
		  y(p)=x(p);
		  y(q)=0;
           end
	 end
      end
    elseif abs(x(q))>eps
	  a=5;
	  G(p,p)=0;
	  G(p,q)=1;
	  G(q,p)=1;
	  G(q,q)=0;
	  y(p)=x(q);
	  y(q)=0;
	 Psi=G'*Phi*G;
    else a=0;
end