function [x,fp,info]=uncmin(x0,n,func,fscale)
%UNCMIN  Unconstrained minimization.
%        [x,f,info]=uncmin(x0,n,func)  solves the following problem:
%
%                          min  f(x)                           (*)
%                           x
%        where x0 is the initial guess for the optimal x which is a
%        vector of length n and where  func is the name (in string
%        syntax) of the .m file containing code for the function 
%        f in (*).
%

%        References: 
%        [1] Kahaner, Moler, and Nash, NUMERICAL METHODS AND SOFTWARE,
%        [2] Gill, Murray, and Wright, PRACTICAL OPTIMIZATION

%        Ellen McGrattan,  3-6-89
%        Revised  6-18-89

%-----------------------------------------------------------------------
%                      USER-DEFINED PARAMETERS
%-----------------------------------------------------------------------

epsm=2.22e-16;                                 % machine epsilon
maxits=500;                                    % maximum iterations
gradtl=epsm^(1/3);                             % tolerance for gradient
steptl=sqrt(epsm)                              % tolerance for step length
typsiz=max(abs(x0),ones(n,1))                  % typical x size
eval(['fscale=abs(',func,'(x0));'])            % typical function size
fscale=1.0;                    
iexp=1;                                        % expense of evaluation


p=zeros(n,1);
itncnt=0;
iretcd=-1;
sx=ones(n,1)./typsiz

stpsiz=sqrt(sum((x0.^2).*(sx.^2)))
stepmx=max(1e+3*stpsiz,1e+3)
ndigit=fix(-log10(epsm))

iagflg=0;
iahflg=0;

rnf=max(10^(-ndigit),epsm)
if iagflg;
  eval(['[f,g]=',func,'(x0);'])
else;
  eval(['f=',func,'(x0);'])
  for i=1:n;
    stepsz=zeros(n,1);
    stepsz(i)=sqrt(rnf)*max(abs(x0(i)),1/sx(i));
    eval(['g(i,1)=(',func,'(x0+stepsz)-f)/stepsz(i);'])
  end;
end;
[info,icscmx]=optstp(n,x0,f,g,x0,itncnt,0,gradtl,steptl,sx,fscale,maxits, ..
                     iretcd,0);
if info;
  x=x0;
  fp=f;
  return;
end;      
if iexp;
  A=diag(sx);
else;
%
% EVALUATE ANALYTIC OR FINITE DIFFERENCE HESSIAN AND CHECK ANALYTIC
% HESSIAN IF REQUESTED (ONLY IF USER-SUPPLIED ANALYTIC HESSIAN
% ROUTINE D2FCN FILLS ONLY LOWER TRIANGULAR PART AND DIAGONAL OF A).
%
  disp('WARNING: where is the analytical or numerical hessian for iexp=0?')
end;
%
% ITERATION
% ---------

while info==0;
  itncnt=itncnt+1;
%
% SOLVE FOR NEWTON STEP:  AP=-G
%
  p=-(tril(A)*tril(A)')\g;
  [x,fp,p,mxtake,iretcd]=lnsrch(n,x0,f,g,p,func,stepmx,steptl,sx);
  while (iretcd==1 & iagflg==0);
    p=-(tril(A)*tril(A)')\g;
    [x,fp,p,mxtake,iretcd]=lnsrch(n,x0,f,g,p,func,stepmx,steptl,sx);
    iagflg=-1;
    for i=1:n;
      stepsz=zeros(n,1);
      stepsz(i)=rnf^(1/3)*max(abs(x0(i)),1./sx(i));
      eval(['g(i,1)=(',func,'(x0+stepsz)-',func,'(x0-stepsz))/(2*stepsz(i));'])
    end;
  end;

  p=x-x0;
  if iagflg==-1;
    for i=1:n;
      stepsz=zeros(n,1);
      stepsz(i)=rnf^(1/3)*max(abs(x(i)),1./sx(i));
      eval(['gp(i,1)=(',func,'(x+stepsz)-',func,'(x-stepsz))/(2*stepsz(i));'])
    end;
  end;
  if iagflg==0;
    for i=1:n;
      stepsz=zeros(n,1);
      stepsz(i)=sqrt(rnf)*max(abs(x(i)),1/sx(i));
      eval(['gp(i,1)=(',func,'(x+stepsz)-fp)/stepsz(i);'])
    end;
  end;
  if iagflg==1;
    eval(['[tem,gp]=',func,'(x);']);
  end;
  [info,icscmx]=optstp(n,x,fp,gp,x0,itncnt,icscmx,gradtl,steptl,sx, ..
                       fscale,maxits,iretcd,mxtake);
  if info==0;
    if ~iexp;
      disp('WARNING: what about iexp=0 case?')
    else;
      [A,noupdt]=secfac(n,x0,g,A,x,gp,epsm,itncnt,rnf,iagflg,noupdt);
      if any( imag(eig(A)) );
        disp('A does not have real eigs')
        [imag(eig(A)),imag(eig(triu(A)+triu(A,1)'))]
        A
      end;
    end;
    f=fp;
    x0=x;
    g=gp;
  end;
end;
%
% TERMINATION
% -----------
% RESET XPLS,FPLS,GPLS,  IF PREVIOUS ITERATE SOLUTION
%
if info==3;
  fp=f;
  x=x0;
  gp=g;
end;
