function [x_min,f_min,fcteval,compteur,iteration,temps,comptons,drapeau] = trmgould(x0,R,prob)

% Call: [x_min,f_min,fcteval,compteur,iteration,temps,comptons,drapeau] = trmgould(x0,R,prob)
% This is the trust region method found in Gould et al. [1999]
% This algorithm is used to minimize an unconstrained function. The function's
% objective should be available using the file 'objective.m' and the function's
% objective, gradient and Hessian, using the file 'objgradhess.m'. This
% function uses the subroutine 'newtrust.m' to solve the trust region
% subproblem
% 
% INPUT VARIABLES:
% 
% x0 ... initial estimate of a minimizer of the function
% R  ... initial radius of the trust region
% prob ... number of the problem to be solve (used to access the files
%           objective.m and objgradhess.m)
%
% OUTPUT VARIABLES:
%
% x_min   ... an approximate of a minimizer of the function
% f_min   ... objective at x_min. This approximates the minimum of the function
% fcteval ... number of fonctions evaluations
% compteur .. number of matrix-vector multiplications
% iteration . number of iterations done by the trust region method 
% temps   ... time taken to find the approximate solution
% comptons .. relative distance of the optimal Lagrange multipliers of the 
%             trust region subproblems to the minimum eigenvalue of the
%             Hessians
% drapeau ... vector of integers (0,1,2) used to distinguish the different
%             cases for the trust region subproblemes

% GLOBAL VARIABLES
global A D counter % global variables used by the newtrust4b.m subroutine 


global A counter x

b = ' '; % define blanks
n = size(x0,1);


temps = cputime;

% Constant used to update the trust region radius (i.e. evaluate the
% performance of the trust region step
n1 = 0.01;
n2 = 0.95;
gam1 = 0.5;
gam2 = 2;
s = R;

x = x0; % x is the current approximation to a minimizer of the function

[fopt,g,A] = objgradhess(x,prob);; % fopt is the objective evaluated at x (i.e.
                                   % the current approximation to the minimum
                                   % of the function, g is the gradient and A
                                   % is the Hessian.

% initiaze variables
fcteval = 1;
iteration = 1;
compteur = 0;

  disp(['iteration',b,b,b,b,'time elapsed so far',b,b,b,b,'f(x)',b,b,b,b,b,b,b,b,'norm of the gradient']);
  disp([b,b,b,b,num2str(iteration),b,b,b,b,b,b,b,b,b,b,b,b,b,num2str(cputime-temps),b,b,b,b,b,b,b,b,b,b,b,b,b,num2str(fopt),b,b,b,b,b,b,b,b,b,b,b,num2str(norm(g)),b,b,b,b]);

if (norm(g)/(1+abs(fopt)) < 0.00001) % relative error on the norm of the
                                     % gradient 

   % Initial approximate is good enough
   x_min = x;
   f_min = fopt;
   return; % END OF ALGORITHM
end

while (norm(g)/(abs(fopt)+1) > 0.00001)

epsilon = 10^(-6);
dgaptol = min(epsilon,10^(-4)*norm(g)); % relative duality gap. This is used
                                        % so that initially we don't ask a lot
                                        % of precision for the trust region
                                        % subproblems, but as we get close to
                                        % the solution (gradient is small), we
                                        % do.
  % SOLVE THE TRUST REGION SUBPROBLEM
  [d,qopt,lopt,lambdamin,iter,drapeau(iteration)] = newtrust(-g,s,dgaptol);

  % d is the step direction
  % qopt is the predicted value of the objective at x+d (predicted by the 
  % quadratic model)
  % lopt is the Lagranga multiplier of the optimal solution to the trust
  % region subproblem.
  % lambdamin is the minimum eigenvalue of the current matrix A (the Hessian at
  % x).
  % drapeau(iteration) is an integer that tells us which case occurred in the
  % trust region subproblem (hard case, case 1 and 2 and easy case).

  comptons(iteration) = (lambdamin-lopt)/(abs(lambdamin) + 1); % relative 
                        % distance of the optimal Lagrange multiplier to the
                        % minimum eigenvalue of A (Hessian at x). This is just
                        % used for the purpose of analysing the iterates.;

  compteur = compteur + counter; % counter is the number of matrix-vector
                                 % multiplication done in the last trust region
                                 % subproblem. So update compteur, the total
                                 % number of matrix-vector multiplications


  f_predicted = fopt + 1/2*qopt;
  f_real = objective(x+d,prob);



  r = (fopt - f_real)/(fopt - f_predicted); % r is a measure of our improvement
                                            % r=1 is good, r=0 is bad
  if (r > n1)
     % improvement of the objective function is good enough so take a step
     x = x+d;
     [fopt,g,A] = objgradhess(x,prob);
     fcteval = fcteval + 1;
  end

  if (r > n2)
    % improvement of the objective function is great, increase the trust region
    % radius
    s = gam2*s;
  end

  if (r < n1)
    % improvement of the objective function is not satisfactory. Reduce the 
    % trust region radius
    s = gam1*s;
  end

  disp('new values');
  r
  s

  iteration = iteration + 1; % one more iteration has been done


  % display info on the iteration
  disp(['iteration',b,b,b,b,'time elapsed so far',b,b,b,b,'f(x)',b,b,b,b,b,b,b,b,'norm of the gradient']);
  disp([b,b,b,b,num2str(iteration),b,b,b,b,b,b,b,b,b,b,b,b,b,num2str(cputime-temps),b,b,b,b,b,b,b,b,b,b,b,b,b,num2str(fopt),b,b,b,b,b,b,b,b,b,b,b,num2str(norm(g)),b,b,b,b]);

end

% return the approximation to the optimal solutions
f_min = fopt;
x_min = x;

temps = cputime - temps; % temps = time to solve the problem