function [y, t]=onestep0(Ax,x,p,rhos)
%
%  It perform line search to solve
%
%            inf_t rho(x+t*p), 
%
%  where
%
%       rho(x)=x' A x/(x' x), Rayleigh quotient of a Hermitian matrix A
%
% 
% RCL 6/27/2009
%---------------------------------------------------------------
%
%  Input 
%
%        Ax     n-by-1; it is A*x
%        x     n-vector, x'*x=1
%        p     n-vector, p'*p=1
%        rhos  [rho(x)  rho(p)]
%
%  Output
%
%        y     x+t*p, in case t=infty, y=p instead
%        t     optimal argument
%              = 0    inf_t rho(x+t*p)=rho(x):      no improvement
%
%AA=[x'; p']*A*[x p]; 
tmp=Ax'*p; 
AA=[rhos(1)  tmp;
    conj(tmp) rhos(2)];
%BB=[x'; p']*[x p];
tmp=x'*p;
BB=[ 1    tmp;
    conj(tmp) 1];

% a=AA(2,2)*(BB(1,2)+BB(2,1))-(AA(1,2)+AA(2,1))*BB(2,2);
tmp1 = AA(2,2)*(BB(1,2)+BB(2,1)); tmp2 = (AA(1,2)+AA(2,1))*BB(2,2);
a=tmp1-tmp2;
if abs(a) <= 4*eps*(abs(tmp1)+abs(tmp2))
   a = 0;
end

% b=2*( AA(2,2)*BB(1,1)-AA(1,1)*BB(2,2) );
tmp3 = AA(2,2); tmp4 = AA(1,1); 
b = 2*(tmp3-tmp4);
if abs(b) <= 8*eps*(abs(tmp3)+abs(tmp4))
   b = 0;
end

% c=(AA(1,2)+AA(2,1))*BB(1,1)-AA(1,1)*(BB(1,2)+BB(2,1));
tmp5 = (AA(1,2)+AA(2,1)); tmp6 = AA(1,1)*(BB(1,2)+BB(2,1));
c = tmp5 - tmp6;
if abs(c) <= 4*eps*(abs(tmp5)+abs(tmp6))
   c = 0;
end

if a ~= 0,
   info = 0; tmp = b*b-4*a*c;
   if b <= 0,
      t = (-b+sqrt(abs(tmp)))/(2*a);
   else
      t = -2*c/(b+sqrt(abs(tmp)));
   end
   y=x+t*p;
else  % a == 0
   if  b == 0 % a = b == 0, too
       y = x; t = 0;      % no improvement
   elseif b < 0   % a = 0, b < 0
       t=1/0; y=p;
   else % a = 0, b>0
      t=-c/b; y= x+t*p;
   end 
end
 
