function [lambda,fxlambda,gxlambda,count]=wolfes(ffunc,gfunc,n,x,s,fx,gx,c1,c2,alpha,beta)
%INPUT: Only the first five parameter is neccesary.
%fx and gx are the function value and gradient at point x. It's for
%reduction of unnecssary recalculation of these values in applications
%alpha has to be positive and less than 1 if present.
%beta has to be greater than 1 if present.
%
%OUTPUT:
%(1) This line search always stops for differentiable function ffunc. If
%ffunc is continous, lambda is a step length which satisfies Wolfe
%condition; otherwise, this line search may not stop, and even it stops,
%lambda may not satisfy Wolfe condition since it may not exist.
%(2) Whether ffunc is continous or not, when this line search stops,
%fxlambda and gxlambda are always the function value and the gradient of at
%step length lambda, and count is the total number of function evaluations
%including gradient evaluations.
%(3) If s is an ascending direction, lambda returned will be negtive. If s
%is orthogonal to gradient at x, lambda will be 0.
%

lambda=0; count=0;
if nargin<11,beta=1.5; end
if nargin<10,alpha=0.5; end
if nargin<9,c2=0.9; end
if nargin<8,c1=0.1; end
if nargin<7,gx=feval(gfunc,x); count=count+n; end
if nargin<6,fx=feval(ffunc,x); count=count+1; end
fxlambda=fx; gxlambda=gx;

%modulate s according to whether s is a descent direction
glambda0=s'*gx;
ascending=0;
if glambda0==0
    return;
elseif glambda0>0
    ascending=1;
    s=-s;
    glambda0=-glambda0;
end

%determine a step length lambda such that f() goes below line c1
lambda=1;
fxlambda=feval(ffunc,x+lambda*s); count=count+1;
while fxlambda>fx+c1*lambda*glambda0
    lambda=alpha*lambda;
fxlambda=feval(ffunc,x+lambda*s); count=count+1;
end
gxlambda=feval(gfunc,x+lambda*s); count=count+n;
if lambda==0 | fxlambda==fx, if ascending, lambda=-lambda; end; return; end
glambda=s'*gxlambda;
if glambda>=c2*glambda0
    if ascending, lambda=-lambda; end %unmodulate s
    return;
end

%determine an interval [a,b] where there must be a point satisfies Wolfe
% condition
if fxlambda>=fx+c2*lambda*glambda0
    a=0; b=lambda;
else
    a=lambda;
    while fxlambda<fx+c2*lambda*glambda0
        lambda=beta*lambda;
fxlambda=feval(ffunc,x+lambda*s); count=count+1;
end
    if isinf(fxlambda) | isinf(lambda)
        gxlambda=feval(gfunc,x+lambda*s); count=count+n;
        if ascending, lambda=-lambda; end %unmodulate s
        return;
    end
    p=a; q=lambda;
    while 1
        lambda=(p+q)/2;
fxlambda=feval(ffunc,x+lambda*s); count=count+1;
        if lambda==p | lambda==q
            %if this happens, function f() must be non-continous and there
            %is no point which satisfies Wolfe condition. we choose p as
            %lambda and return, which is below line c2.
            if lambda==q
                lambda=p;
fxlambda=feval(ffunc,x+lambda*s); count=count+1;
            end
            gxlambda=feval(gfunc,x+lambda*s); count=count+n;
            if ascending, lambda=-lambda; end %unmodulate s
            return;
        end
        if fxlambda>fx+c1*lambda*glambda0
            q=lambda;
        elseif fxlambda<fx+c2*lambda*glambda0
            p=lambda;
        else
            break
        end
    end
    b=lambda;
end
%find a Wolfe point which satisfies Wolfe condition
while 1
    lambda=(b+a)/2;
fxlambda=feval(ffunc,x+lambda*s);
gxlambda=feval(gfunc,x+lambda*s);
    count=count+1+n;
glambda=s'*gxlambda;
if glambda>=c2*glambda0
        if ascending, lambda=-lambda; end %unmodulate s
        return;
    end
    if lambda==a | lambda==b
        %if this happens, function f() must be non-continous and there is
        %no point which satisfies Wolfe condition. we choose a as lambda
        %and return, which is below line c2.
        if lambda==b
            lambda=a;
fxlambda=feval(ffunc,x+lambda*s);
gxlambda=feval(gfunc,x+lambda*s);
            count=count+1+n;
        end
        if ascending, lambda=-lambda; end %unmodulate s
        return;
    end
    if fxlambda>=fx+c2*lambda*glambda0
        b=lambda;
    else
        a=lambda;
    end
end

return