%
%   Hyperbolic Givens transformation -
%  a special treatment in the case of zero hyperbolic energy
%
%  explanation:
%  suppose that the hyperbolic energy of the reduced vector is in some
%  step of the algorithm equal to zero. In this case the procedure
%  hypgiv gives the matrix Psi with Psi(p,p)=\pm 2. In the next
%  step, i.e. annihilating the next component of x, we call instead 
%  of hypgiv the procedure zerogiv, which gives again
%  the diagonal matrix Psi with all diagonal elements \pm 1
%  (except the case a=4), see the Table 3
%
%  [y,G,Psi,a]=zerogiv(n,p,q,x,Phi)
%
%  input data:
%          n ... dimension of x
%  p,q (p<q) ... components of x for reduction
%                (we will annihilate x(q))
%          x ... the reduced vector
%        Phi ... the given diagonal matrix with Phi(p,p)=\pm 2,
%                others diagonal entries of Phi are \pm 1
%  output data:
%          y ... vector x after the transformation
%          G ... matrix for reduction (see Table 3), G*x=y
%        Psi ... matrix Phi after the transformation,
%                Psi has on diagonal only \pm 1 (except the case a=4),
%                Psi=G'*Phi*G
%          a ... parameter, the same meaning as in procedure hypgiv
%                (if a=4 the hyperbolic
%                 energy of the reduced vector is again zero)
%
%
function [y,G,Psi,a]= zerogiv(n,p,q,x,Phi)
y=x; Psi=Phi; G=eye(n);a=0;

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)=x(p)'*ro1/sqrt(2);
    G(p,q)=x(q)'*ro1/sqrt(2);
    G(q,p)=ro1*x(q);
    G(q,q)=-x(p)*ro1;
    y(p)=ro/sqrt(2);
    y(q)=0;
else e=abs(x(p))-abs(x(q));
     if abs(e)< eps
	a=4;
	u=-x(p)/x(q);
        G(p,p)=sqrt(2);
        G(p,q)=u/sqrt(2);
	G(q,p)=u'*sqrt(2);
        y(p)=x(p)/sqrt(2);
	y(q)=0;
      else if e>0
	    a=2;
	    ro=sqrt(abs(x(p))*abs(x(p))-abs(x(q))*abs(x(q)));
	    z=1;
           else
	    a=3;
	   ro=sqrt(abs(x(q))*abs(x(q))-abs(x(p))*abs(x(p)));
	   z=-1;
           end
           ro1=1/ro;
           G(p,p)=ro1*x(p)'/sqrt(2);
           G(p,q)=-ro1*x(q)'/sqrt(2);
           G(q,p)=ro1*x(q);
	   G(q,q)=-ro1*x(p);
           y(p)=z*ro/sqrt(2);
	   y(q)=0;
       end
end
Psi=G'*Phi*G;