function  [X, info, perf] = conj_grad(fun,par, x0, opts)
%CONJ_GRAD  Conjugate gradient method for nonlinear optimization:
%  Find  xm = argmin{f(x)} , where  x  is an n-vector and the scalar
%  function  f  with gradient  g  (with elements  g(i) = Df/Dx_i )
%  must be given by a MATLAB function with with declaration
%            function  [f, g] = fun(x, par)
%  par  holds parameters of the function.  It may be dummy.  
%  
%  Call
%      [X, info {, perf}] = conj_grad(fun,par, x0, opts)
%
%  Input parameters
%  fun  :  String with the name of the function.
%  par  :  Parameters of the function.  May be empty.
%  x0   :  Starting guess for  x .
%  opts :  Vector with 9 elements:
%          opts(1) :  Choice of CG method:
%              opts(1) = 1 :  Fletcher - Reeves, 
%                otherwise :  Polak - Ribiere 
%          opts(2) :  Choice of line search method:
%              opts(2) = 1 :  Exact l. s., otherwise : Soft l.s.
%          opts(3) :  Upper bound on initial step
%          opts(4:6)  used in stopping criteria:
%              ||g||_inf <= opts(4)                     or 
%              ||dx||_2 <= opts(5)*(opts(5) + ||x||_2)  or
%              no. of function evaluations exceeds  opts(6) . 
%          opts(7:9)  used as linesearch parameters lpar(3:5) .
%              See LINESEARCH .
%          Any illegal element in  opts  is replaced by its
%          default value :
%              opts(1:6) = [2  2  1  1e-4*||g(x0)||_inf  1e-6  100]
%          and
%              opts(2)=1:  opts(7:9) = [1e-4  1e-6  10]
%              opts(2)=2:  opts(7:9) = [1e-1  1e-2  10]
%
%  Output parameters
%  X    :  If  perf  is present, then array, holding the iterates
%          columnwise.  Otherwise, computed solution vector.
%  info :  Performance information, vector with 6 elements:
%          info(1:3) = final values of 
%              [f(x)  ||g||_inf  ||dx||_2] 
%          info(4:5) = no. of iteration steps and evaluations of (f,g)
%          info(6) = 1 :  Stopped by small gradient
%                    2 :  Stopped by small x-step
%                    3 :  Stopped by  opts(6) .  
%  perf :  (optional). If present, then array, holding 
%            perf(1:2,:) = values of  f(x) and ||g||_inf
%            perf(3:5,:) = Line search info:  values of  
%                          alpha, phi'(alpha), no. fct. evals.
%
%  Use of other MATLAB functions :  LINESEARCH .
%
%  Method
%  See Chapter 4 in  P.E. Frandsen, K. Jonasson, H.B. Nielsen, 
%  O. Tingleff: "Unconstrained Optimization", IMM, DTU.  1999.
   verbose = 0;
%  Hans Bruun Nielsen,  IMM, DTU.  99.08.05-09 / 08.23 / 12.07
    alfa=-1;
   %  Check call 
   
   [x  n  f  g  opts] = check(fun,par,x0,opts);
   %  Finish initialization
   k = 1;   kmax = opts(6);
   ng = norm(g, inf);   n2g = norm(g);   gg = n2g^2;
   gam = 0;   h = zeros(size(x));   nh = 0;   neval = 1;
   Trace = nargout > 2;
   if  Trace
         X = x(:)*ones(1,kmax+1);
         perf = [f; ng; zeros(3,1)]*ones(1,kmax+1);
       end
   found = ng <= opts(4);
   nit=0;
   
   while  ~found
     %  Previous values
     hpr = h;   nhpr = nh;  ggpr = gg; nit=nit+1;
     if  opts(1) == 2,   gpr = g; end
     %  New direction and scaling
     k = k+1;   h = gam*h - g;   nh = norm(h);
     %  Check descent
     gh = dot(g,h);     
     if  dot(g,h) >= -1e-3 * n2g * nh    % Not descent
%       disp('Flipping Descent');
       h = -g;   nh = n2g;
     end
     % Scale, aiming at  alpha = 1
     if  k > 2,  sc = .9*alfa*nhpr/nh;
     else    % First step.  Try not to move too far
       sc = opts(3)/32/nh;
     end
     %  Line search
     lspar = [32  opts([2 7:8])  min(kmax-neval, opts(9))];

     [al2  f2  g2  dval2] = linesearch(fun,par,x,f,g, sc*h, lspar);
     if al2==0
       h = -g;
       [al  f  g  dval] = linesearch(fun,par,x,f,g, sc*h, lspar); 
     else
       al = al2; f = f2; g = g2; dval = dval2;
     end
	   
     alfa = sc*al;   ng = norm(g, inf);   n2g = norm(g);  gg = n2g^2;
     %  Update  gamma  and  x
     if  opts(1) == 1,  gam = gg/ggpr;
     else,              gam = dot(g-gpr,g)/ggpr; end
     x = x + alfa*h;  
%     x = x + jitter_mag/nit*2*(0.5-rand(size(x))).*x;
     neval = neval + dval;
     if  Trace
           X(:,k) = x(:);   perf(1:2,k) = [f; ng];
           perf(3:5,k-1) = [alfa; dot(h,g); dval]; end
     f_hist(k) = full(f);
     %  Check stopping criteria
     if      ng <= opts(4)
         found = 1;
     elseif  alfa*nh <= opts(5)*(opts(5) + norm(x)) 
         found = 2;
     elseif  neval >= kmax
         found = 3;
     end
     if verbose
       fprintf('It=%d Val=%g NEval=%d\n',nit,full(f),neval);
     end
   end  % iteration

%   fprintf('It=%d Val=%g NEval=%d\n',nit,full(f),neval)
   %  Set return values
   if  Trace
     X = X(:,1:k);   perf = perf(:,1:k);
   else,  X = x;  end
   info = [f  ng  alfa*nh  k-1  neval  found];

% ==========  auxiliary function  =================================

function  [x,n, f,g, opts] = check(fun,par,x0,opts0)
%  Check function call
   sx = size(x0);   n = max(sx);
   if  (min(sx) > 1)
       error('x0  should be a vector'), end
   x = x0;   [f g] = feval(fun,x,par);
   sf = size(f);   sg = size(g);
   if  any(sf-1) | ~isreal(f)
       error('f  should be a real valued scalar'), end
   if  any(sg - sx)
       error('g  should be a vector of the same type as  x'), end

   so = size(opts0);
   if  (min(so) ~= 1) | (max(so) < 9) | any(~isreal(opts0(1:9)))
       error('opts  should be a real valued vector of length 9'), end
   opts = opts0(1:9);   opts = opts(:).';
   i = find(opts(1:2) > 2);
   if  length(i)    % Set default values
     d = [2  2];   opts(i) = d(i);
   end
   i = find(opts(1:6) <= 0);
   if  length(i)    % Set default values
%     d = [2  2  1  1e-4*norm(g, inf)  1e-6  10000];
     d = [2  2  1  1e-4*norm(g, inf)  1e-6  5000];
     opts(i) = d(i);
   end
   i = find(opts(7:9) <= 0);
   if  length(i)    % Set default values
     d = [1e-6  1e-6  10; 1e-1  1e-2  5];
     %      d = [1e-7  1e-7  10; 1e-5  1e-5  100];
     opts(6+i) = d(opts(2),i);
   end

return