PROC (4) = dfpmin2(p,ftol,maxsec,maxit,&fct,&dfct,&linm);
/* PURPOSE: minimizes a function called FCT using the
   Broyden-Fletcher-Goldfarb-Shanno (BFGS) method.The procedure is tailored
   after DFPMIN in Numerical Recipes by Bo Honore and very slightly modified
   by in Aug96 by Alan G. Isaac. (Modifications are noted; original is in archive.)

This code is provided, without performance guarantees,
for public non-commercial use.

   FORMAT: {p,fret,iter,tim}= dfpmin2(p,ftol,maxsec,maxit,&fct,&dfct,&linm);
example: {p,fret,iter,tim}=dfpmin2(.5|2,.001,1000,200,&myfct,&gradp,&linmin2);

   INPUT

   p      (nx1)    starting vector
   ftol   (1x1)    fractional convergence tolerance in the function value
   maxsec (1x1)    maximum number of seconds allowed for running program
   maxit  (1x1)    maximum number of iterations allowed
   fct    proc     function we want to minimize
   dfct   proc     routine calculating the gradient of the function
                   (agi: now must return a row vector. E.g., use proc gradp.)
   linm   proc     routine performing line minimization
                   (agi: proc linmin2 is included below and may be used)

   OUTPUT

   p      (nx1)    minimizing vector
   fret   (1x1)    value of function at minimium,p
   iter   (1x1)    iterations taken
   tim    (1x1)    seconds elapsed

*/
    local tim,date1,fp,xi,g,hessin,hdg,dg,iter,eps,fret,n,tol,fac,fad, fae,
      fct:proc,dfct:proc,linm:proc;

    tim=0;
    date1=date;
    eps=0.0000000001;
    tol=0.00001;
    n=rows(p);
    fp=fct(p);
    g=dfct(&fct,p)';  @agi: change to use gradp@
    xi=-g;
    hessin=eye(n);
    iter=1;
    do while iter<=maxit;
    {p,fret}=linm(p,xi,tol,&fct);
    tim=ethsec(date1,date)/100;
    /* agi: monitor() was not provided despite the line
    call monitor(p,fret,iter,tim);
           so I've added a monitoring subroutine
    */
    gosub monitor;
    if ((2*abs(fp-fret)).le ftol*(abs(fp)+abs(fret)+eps));
    RETP(p,fret,iter,tim);
    endif;
    if (tim .ge maxsec); "Maximum number of seconds exceeded";
    RETP(p,fret,iter,tim);
    endif;
    fp=fct(p);
    dg=g;
    fret=fct(p);
    g=dfct(&fct,p)';             @agi: again, change to use gradp@
    dg=g-dg;
    hdg=hessin*dg;
    fac=dg'*xi;
    fae=dg'*hdg;
    if (fac==0); RETP(p,fret,iter,tim); else; fac=1/fac; endif;
    if (fae==0); RETP(p,fret,iter,tim); else; fad=1/fae; endif;
    dg=fac*xi-fad*hdg;
    hessin=hessin+fac*xi*xi'-fad*hdg*hdg'+fae*dg*dg';
     xi=-hessin*g;
    @agi: the next line added to allow escaping the program@
     if (iter-1)/10==trunc((iter-1)/10);gosub brk;endif;
    iter=iter+1;
    endo; "Maximum number of iterations exceeded";
    RETP(p,fret,iter,tim);
/* begin subroutines (agi: these are additions)*/
monitor:
     scroll 1|1|15|80|5|7;
     locate 1,1;
     print "iteration: "$+ftocv(iter,1,0);
     print "seconds elapsed: "$+ftocv(tim,1,0);
     print "current value of function: "$+ftocv(fret,1,5);
     print "current parameter vector (transposed): ";
     p';
     print;
     pause(0.3);
return;

brk:
    locate 24,1;
    printdos "\27[7m";
    printdos "Press [Esc] to exit program.";
    printdos "\27[0m";
    if key==27;end;endif;
return;

ENDP;


PROC (2) = linmin2(pii,xi,tol,&fct);

/*
   This routine performs an ad hoc n-dimensional Golden Section Search for
   the minimum of a function called FCT. Includes very small changes by agi.

   INPUT

   pii   (nx1)    vector of starting values   @agi: change name here and below@
   xi    (nx1)    direction vector
   tol   (1x1)    tolerance
   fct   proc     function we want to minimize

   OUTPUT

   p     (nx1)    minimizing vector
   f     (nx1)    value of fct at minimum,p
*/

    local f1,f2,fm,f,p1,p2,pm,p,xd,x1,x2,d1,d2,fx1,fx2,itry,fct:proc;
    d1=0.618; d2=0.382;
    p1=pii;   @agi: named pii instead of p@
    p2=pii+xi;
    f1=fct(p1);
    itry=0;
  be:
    itry=itry+1;
    if (itry.eq 10); "linmin itry=10"; stop; endif;
    f2=fct(p2);
    if(f1.eq f2); p2=p2+xi; goto be; endif;
    if (f1.le f2);
    pm=p1; fm=f1;
    f1=f2; p1=p2;
    else;
    pm=p2; fm=f2;
    endif;
    xd=pm-p1;
    p2=pm+xd;
     f2=fct(p2);
        local ctr;ctr=0; @agi:add ctr; used for escaping program@
    do until (f2.gt fm);
      xd=2*xd;
      p2=pm+xd;
      f2=fct(p2);
        @agi: next line added to allow escape from program@
      if ctr/10==trunc(ctr/10);gosub brk;endif;ctr=ctr+1;
    endo;
/* p1'; pm'; p2'; */
    stit:
    x1=d1*p1+d2*p2;
    x2=d2*p1+d1*p2;
    fx1=fct(x1);
    fx2=fct(x2);
    if (fx1.le fx2);
    p2=x2;
    else;
    p1=x1;
    endif;
    if (sqrt(sumc((p1-p2).^2)).ge tol); goto stit; endif;
    if (fx1.lt fx2);
    p=x1;
    else;
    p=x2;
    endif;
    f=fct(p);
RETP(p,f);
end;
brk:
    locate 24,1;
    printdos "\27[7m";
    printdos "Press [Esc] to exit program.";
    printdos "\27[0m";
    if key==27;end;endif;
    locate 1,1;
return;
ENDP;
