% Input variables:
%   xk    - current iterate
%   dk    - current descent direction
%   t0    - initial step size parameter
%   gamma - factor to enlarge t in loop (A.1)
%   alpha - factor for the Armijo rule
%   rho   - factor for the gradient condition
%
% Version May 25, 2009 (Stefan Volkwein)
function tk = wolfe_powell(fun,grad_fun,xk,dk,t0,gamma,alpha, rho)

% Check input parameters t0, gamma, alpha, and rho:
  if t0<=0
      error('Initial step size parameter must be positive!');
  end
  if gamma<=1
      error('gamma must be larger than 1!');
  end
  if alpha<=0
      error('alpha must be positive!');
  end
  if rho<=alpha
      error('rho must be larger than alpha!');
  end
  if rho>=1
      error('rho must be smaller than 1!');
  end
  
% Compute phi(t) for t=0 and for t=t0:
  phi0=fun(xk);
  phit=fun(xk+t0*dk);
  
% Compute the gradient of phi(t) for t=0 and t=t0:
  grad_phi0=grad_fun(xk)'*dk;
  grad_phit=grad_fun(xk+t0*dk)'*dk;
  
% Check descent property for dk:
  if (grad_phi0>=-eps)
      error('dk is not a descent direction!');
  end

% Initialize the step size parameter:
  t=t0;
  
% Loop (A.1):
  flag_A1=1;
  while flag_A1
      if (phit>=phi0+alpha*t*grad_phi0)
          % Armijo rule is violated.
          % Fix interval [a,b] containing Wolfe-Powell step size parameter:
            a=0;
            b=t;
          % Leave loop (A.1) and goto (B.0):
            flag_A1=0;
            flag_B0=1;
      else
          % Armijo rule holds.
          if grad_phit>=rho*grad_phi0
              % Gradient condition holds, i.e., Wolfe-Powell conditions are
              % satisfied. Outputs are t=t_k, phi(t), and phi'(t):
                tk=t;
                %fk=phit;
                %grad_fk=grad_phit;
              % Leave loop (A.1) and do not goto (B.0):
                flag_A1=0;
                flag_B0=0;
          else
              % Gradient condition is violated. Enlarge the interval [0,t]
              % containing a Wolfe-Powell stepsize:
                t=gamma*t;
          end
      end
  end
  
  if flag_B0
      % In [0,t]=[a,b] there is at least one t satisfying the Wolfe-Powell
      % conditions:
        t1=a;
        t2=b;
        
      % Loop (B.1):
        flag_B1=1;
        while flag_B1
            % Choose trial step size parameter:
            quadratic=2*(phit-phi0-grad_phi0*(t2-t1));
            if quadratic>eps
                t=t1+(-grad_phi0*(t2-t1)^2)/quadratic; 
            else
                % Choose mean value:
                t=(t1+t2)/2;
            end
            % Evaluate the function phi at t:
              phit=fun(xk+t*dk);
              if phit>=phi0+alpha*t*grad_phi0
                  % Armijo rule is violated.
                  % Fix interval [a,b] containing Wolfe-Powell step size
                  % parameter: t1 is not changed.
                    t2=t;
              else
                  % Armijo rule holds.
                  % Evaluate the gradient of phi at t:
                    grad_phit=grad_fun(xk+t*dk)'*dk;
                    if grad_phit>=rho*grad_phi0
                        % Gradient condition holds, i.e., Wolfe-Powell
                        % conditions are satisfied. Outputs are t=t_k,
                        % phi(t), and phi'(t):
                          tk=t;
                          %fk=phit;
                          %grad_fk=grad_phit;
                        % Leave loop (B.1):
                          flag_B1=0;
                    else
                        % Gradient condition is violated.
                        % Fix interval [t1,t2] containing Wolfe-Powell step
                        % size parameter: t2 is not changed.
                          t1=t;
                    end
              end
        end
  end