ADMB Documentation  11.1.1913
 All Classes Files Functions Variables Typedefs Friends Defines
newfmin.cpp
Go to the documentation of this file.
00001 /*
00002  * $Id: newfmin.cpp 1906 2014-04-17 20:52:49Z johnoel $
00003  *
00004  * Author: Unknown
00005  * Copyright (c) 2009-2012 ADMB Foundation
00006  *
00007  * This file was originally written in FORTRAN II by an unknown author.
00008  * In the 1980s, it was ported to C and C++ and extensively modified by
00009  * David Fournier.
00010  *
00011  */
00017 #include <fvar.hpp>
00018 #if defined(_WIN32)
00019   #include <windows.h>
00020   #include <admodel.h>
00021 #endif
00022 #if defined(__BORLANDC__)
00023   #include <signal.h>
00024 #endif
00025 #ifdef __ZTC__
00026   #include <conio.h>
00027 #endif
00028 #include <fvar.hpp>
00029 extern int ctlc_flag;
00030 
00031 #if defined (__WAT32__) || defined (_MSC_VER)
00032   #include <conio.h>
00033 #endif
00034 #if defined(__CYGWIN__)
00035   int kbhit(void)
00036   {
00037     return 0;
00038   }
00039 #endif
00040 #ifdef __ZTC__
00041   #include <iostream.hpp>
00042   #include <disp.h>
00043   #define endl "\n"
00044 #endif
00045   #include <signal.h>
00046   extern "C" void onintr(int k)
00047   {
00048     signal(SIGINT, exit_handler);
00049     ctlc_flag = 1;
00050     if (ad_printf)
00051       (*ad_printf)(
00052 "\npress q to quit or c to invoke derivative checker or s to stop optimizing: "
00053       );
00054   }
00055 #ifdef __NDPX__
00056   #include <iostream.hxx>
00057 #endif
00058 
00059 #if defined (_MSC_VER)
00060   void __cdecl clrscr();
00061 #else
00062   extern "C" void clrscr();
00063   #define getch getchar
00064 #endif
00065 
00066 #include <math.h>
00067 #include <stdlib.h>
00068 #include <stdio.h>
00069 #include <ctype.h>
00070 
00071 int* pointer_to_phase = 0;
00072 
00073 double dafsqrt( double x );
00074   void tracing_message(int traceflag,const char *s);
00075   void tracing_message(int traceflag,const char *s,int *pn);
00076   void tracing_message(int traceflag,const char *s,double *pd);
00077   void tracing_message(int traceflag,const char *s,double d);
00078 
00083   void tracing_message(int traceflag,const char *s)
00084   {
00085     if (traceflag)
00086     {
00087       ofstream ofs;
00088       ofs.open("adtrace",ios::app);
00089       ofs << s << endl;
00090       ofs.close();
00091     }
00092   }
00093 
00098   void tracing_message(int traceflag,const char *s,int *pn)
00099   {
00100     if (traceflag)
00101     {
00102       ofstream ofs;
00103       ofs.open("adtrace",ios::app);
00104       ofs << s << " " << *pn << endl;
00105       ofs.close();
00106     }
00107   }
00108 
00113   void tracing_message(int traceflag,const char *s,double *pd)
00114   {
00115     if (traceflag)
00116     {
00117       ofstream ofs;
00118       ofs.open("adtrace",ios::app);
00119       ofs << s << " " << *pd << endl;
00120       ofs.close();
00121     }
00122   }
00123 
00128   void tracing_message(int traceflag,const char *s,double d)
00129   {
00130     if (traceflag)
00131     {
00132       ofstream ofs;
00133       ofs.open("adtrace",ios::app);
00134       ofs << s << " " << d << endl;
00135       ofs.close();
00136     }
00137   }
00138 int log_values_switch=0;
00139 ofstream logstream("fmin.log");
00140 
00145 void print_values(const double& f, const dvector & x,const dvector& g)
00146 {
00147   logstream << setprecision(13) << f << endl;
00148   logstream << setprecision(13) << x << endl;
00149   logstream << setprecision(13) << g << endl;
00150 }
00151 adtimer* pfmintime = 0;
00152 extern int traceflag;
00153 //#pragma warn -sig
00154 
00155 #ifdef _MSC_VER
00156 BOOL CtrlHandler( DWORD fdwCtrlType )
00157 {
00158   if (fdwCtrlType == CTRL_C_EVENT)
00159   {
00160     ctlc_flag = 1;
00161     if (ad_printf)
00162       (*ad_printf)("\npress q to quit or c to invoke derivative checker: ");
00163     return true;
00164   }
00165   return false;
00166 }
00167 #endif
00168 
00215 void fmm::fmin(const double& _f, const dvector &_x, const dvector& _g)
00216 {
00217   if (log_values_switch)
00218   {
00219     print_values(_f,_x,_g);
00220   }
00221   if (pfmintime==0) pfmintime=new adtimer;
00222   tracing_message(traceflag,"A3");
00223 
00224   /* Remember gradient and function values
00225      resulted from previous function evaluation */
00226   dvector& g=(dvector&) _g;
00227   double& f=(double&) _f;
00228 
00229   /* Create local vector x as a pointer to the argument vector _x */
00230   independent_variables& x= (independent_variables&) _x;
00231     #ifdef DIAG
00232       cout << "On entry to fmin: " << *this << endl;
00233     #endif
00234   tracing_message(traceflag,"A4");
00235 
00236 #ifdef _MSC_VER
00237   SetConsoleCtrlHandler((PHANDLER_ROUTINE)CtrlHandler, true);
00238 #endif
00239 
00240   /* Check the value of control variable ireturn:
00241         -1 (exit status)
00242          0 (initialization of function minimizer)
00243          1 (call1 - x update and convergence check)
00244          2 (call2 - line search and Hessian update)
00245          >=3 (derivative check)
00246   */
00247 #if !defined (_MSC_VER)
00248   #if defined( __SUN__) && !(defined __GNU__)
00249     #if defined( __HP__)
00250         if (ireturn <= 0 )
00251         {
00252           signal(SIGINT, &onintr);
00253         }
00254     #else
00255         if (ireturn <= 0 )
00256         {
00257           signal(SIGINT, (SIG_PF)&onintr);
00258         }
00259     #endif
00260   #endif
00261 #endif
00262 #if defined( __GNU__) || defined (__BORLANDC__)
00263       if (ireturn <= 0 )
00264       {
00265         signal(SIGINT, &onintr);
00266       }
00267 #endif
00268 #ifdef __ZTC__
00269       if (ireturn <= 0 )
00270       {
00271         if (disp_inited == 0)
00272         {
00273           disp_open();
00274           disp_usebios();
00275         }
00276       }
00277 #endif
00278   tracing_message(traceflag,"A5");
00279       if (ireturn ==1 && dcheck_flag ==0)
00280       {
00281         ireturn = 3;
00282       }
00283       if (ireturn >= 3)
00284       {
00285          /* Entering derivative check */
00286          derch(f, x, g, n, ireturn);
00287          return;
00288       }
00289       if (ireturn == 1) goto call1;
00290       if (ireturn == 2) goto call2;
00291 
00292      /* we are here because ireturn=0
00293         Start initializing function minimizer variables */
00294       fbest=1.e+100;
00295   tracing_message(traceflag,"A6");
00296 
00297       /* allocate Hessian
00298          h - object of class dfsdmat, the memory is allocated
00299              only for elements of lower triagonal matrix */
00300       if (!h) h.allocate(n);
00301 
00302       /* initialize w, the dvector of size 4*n:
00303          w(1,n) contains gradient vector in the point x_new = x_old + alpha*p_k
00304          w(n+1,2n) contains direction of linear search, i.e. transpose(p_k)
00305          w(2n+1,4n) are used for Hessian updating terms */
00306       w.initialize();
00307 
00308       /* alpha - the Newton step size */
00309       alpha=1.0; /* full Newton step at the beginning */
00310       ihflag=0;
00311 
00312       /* validity check for dimensions and indexing of
00313          independent vector and gradient vector
00314          Note, this function will work correctly only if
00315          indices start at 1  */
00316      if (n==0)
00317      {
00318        cerr << "Error -- the number of active parameters"
00319          " fmin must be > 0\n";
00320        ad_exit(1);
00321      }
00322   tracing_message(traceflag,"A7");
00323      if (x.indexmin() !=1)
00324      {
00325        cerr << "Error -- minimum valid index"
00326          " for independent_variables in fmin must be 1\n"
00327         << " it is " << x.indexmin() << "\n";
00328         ad_exit(1);
00329      }
00330      if (x.size() <n)
00331      {
00332        cerr << "Error -- the size of the independent_variables"
00333         " which is " << x.size() << " must be >= " << n << "\n"
00334         << " the number of independent variables in fmin\n";
00335         ad_exit(1);
00336      }
00337   tracing_message(traceflag,"A8");
00338      if (g.indexmin() !=1)
00339      {
00340        cerr << "Error -- minimum valid index"
00341          " for the gradient vector in fmin must be 1\n"
00342         << " it is " << g.indexmin() << "\n";
00343         ad_exit(1);
00344      }
00345      if (g.size() <n)
00346      {
00347        cerr << "Error -- the size of the gradient vector"
00348         " which is " << g.size() << " must be >=\n"
00349         << " the number of independent variables in fmin\n";
00350         ad_exit(1);
00351      }
00352      /* end of validity check */
00353   tracing_message(traceflag,"A9");
00354 
00355      /* xx is reserved for the updated values of independent variables,
00356         at the moment put the current values */
00357      for (i=1; i<=n; i++)
00358            xx.elem(i)=x.elem(i);
00359   tracing_message(traceflag,"A10");
00360 
00361       /* itn - iteration counter */
00362       itn=0;
00363       /* icc - obsolete calls counter? */
00364       icc=0;
00365 
00366        /* initialize funval vector,
00367           it will contain last 10 function values (10-th is the most recent)
00368           needed to stop minimization in case if f(1)-f(10)<min_improve  */
00369        for (i=1; i< 11; i++)
00370           funval.elem(i)=0.;
00371   tracing_message(traceflag,"A11");
00372       ihang = 0;  /* flag, =1 when function minimizer is not making progress */
00373       llog=1;
00374       np=n+1;
00375       n1=n-1;
00376       nn=n*np/2;
00377       is=n;
00378       iu=n;
00379       iv=n+n;
00380       ib=iv+n;
00381       iexit=0;    /* exit code */
00382   tracing_message(traceflag,"A12");
00383 
00384       /* Initialize hessian to the unit matrix */
00385       h.elem(1,1) = 1;
00386       for (i=2; i<=n; i++)
00387       {
00388         for ( j=1; j<i; j++)
00389         {
00390            h.elem(i,j)=0;
00391         }
00392         h.elem(i,i)=1;
00393       }
00394   tracing_message(traceflag,"A13");
00395 
00396       /* dmin - minimal element of hessian used to
00397          verify its positive definiteness */
00398       dmin=h.elem(1,1);
00399       for ( i=2; i<=n; i++)
00400       {
00401          if(h.elem(i,i)<dmin)
00402             dmin=h.elem(i,i);
00403       }
00404       if (dmin <= 0.)    /* hessian is not positive definite? */
00405          goto label7020; /* exit */
00406       if(dfn == 0.)
00407          z = 0.0;
00408   tracing_message(traceflag,"A14");
00409       for (i=1; i<=n; i++)
00410       {
00411         xsave.elem(i)=x.elem(i);
00412         x.elem(i)=xx.elem(i);
00413       }
00414       ireturn=1; /* upon next entry will go to call1 */
00415   tracing_message(traceflag,"A15");
00416       if (h.disk_save())
00417       {
00418         cout << "starting hessian save" << endl;
00419         h.save();
00420         cout << "finished hessian save" << endl;
00421       }
00422   tracing_message(traceflag,"A16");
00423       return;
00424       /* End of initialization */
00425   call1: /* first line search step and x update */
00426   tracing_message(traceflag,"A17");
00427       if (h.disk_save())
00428       {
00429         cout << "starting hessian restore" << endl;
00430         h.restore();
00431         cout << "finished hessian restore" << endl;
00432       }
00433   tracing_message(traceflag,"A18");
00434       for (i=1; i<=n; i++)
00435       {
00436         x.elem(i)=xsave.elem(i);
00437       }
00438       ireturn=3;
00439   tracing_message(traceflag,"A19");
00440       {
00441       }
00442       for ( i=1; i<=n; i++)
00443          gbest.elem(i)=g.elem(i);
00444   tracing_message(traceflag,"A20");
00445       funval.elem(10) = f;
00446       df=dfn;
00447       if(dfn == 0.0)
00448          df = f - z;
00449       if(dfn < 0.0)
00450          df=fabs(df * f);
00451       if(df <= 0.0)
00452          df=1.;
00453 label20: /* check for convergence */
00454       ic=0; /* counter for number of calls */
00455       iconv = 1; /* convergence flag: 1 - convergence, 2 - not yet */
00456       for ( i=1; i<=9; i++)
00457          funval.elem(i)= funval.elem(i+1);
00458       funval.elem(10) = f;
00459       /* check if function value is improving */
00460       if ( itn>15 && fabs( funval.elem(1)-funval.elem(10))< min_improve )
00461          ihang = 1;
00462       gmax = 0;
00463       /* satisfy convergence criterion? */
00464       for ( i=1; i<=n; i++)
00465       {
00466         if(fabs(g.elem(i)) > crit) iconv = 2;
00467         if(fabs(g.elem(i)) > fabs(gmax) ) gmax = g.elem(i);
00468       }
00469       /* exit if either convergence or no improvement has been achieved
00470          during last 10 iterations */
00471       if( iconv == 1 || ihang == 1)
00472       {
00473          ireturn=-1;
00474          goto label92;
00475       }
00476       /* otherwise proceed to the Newton step (label21) */
00477       if(iprint == 0)
00478          goto label21; /* without printing */
00479       if(itn == 0)
00480          goto label7000; /* printing Initial statistics first */
00481 #if defined(USE_DDOUBLE)
00482 #undef double
00483       if(fmod(double(itn),double(iprint)) != 0)
00484          goto label21;
00485 #define double dd_real
00486 #else
00487       if(fmod(double(itn),double(iprint)) != 0)
00488          goto label21;
00489 #endif
00490       if (llog) goto label7010;
00491 #if defined (_MSC_VER) && !defined (__WAT32__)
00492         if (!scroll_flag) clrscr();
00493 #endif
00494 label7003: /* Printing table header */
00495       if (iprint>0)
00496       {
00497         if (ad_printf)
00498         {
00499           (*ad_printf)("%d variables; iteration %ld; function evaluation %ld",
00500            n, itn, ifn);
00501           if (pointer_to_phase)
00502           {
00503             (*ad_printf)("; phase %d", *pointer_to_phase);
00504           }
00505           (*ad_printf)(
00506            "\nFunction value %15.7le; maximum gradient component mag %12.4le\n",
00507 #if defined(USE_DDOUBLE)
00508   #undef double
00509               double(f), double(gmax));
00510   #define double dd_real
00511 #else
00512               f, gmax);
00513 #endif
00514         }
00515       }
00516 /*label7002:*/
00517       /* Printing Statistics table */
00518       if(iprint>0)
00519       {
00520         fmmdisp(x, g, n, this->scroll_flag,noprintx);
00521       }
00522 label21 : /* Calculating Newton step */
00523       /* found good search direction, increment iteration number */
00524       itn=itn+1;
00525       for (i=1; i<=n; i++)
00526          x.elem(i)=xx.elem(i);
00527       w.elem(1)=-g.elem(1);
00528       pfmintime->get_elapsed_time_and_reset();
00529 
00530       /* solving system of linear equations H_(k+1) * (x_(k+1)-x(k)) = -g_k
00531          to get next search direction
00532          p_k = (x_(k+1)-x(k)) = - inv(H_(k+1)) * g_k */
00533       for (i=2; i<=n; i++)
00534       {
00535         i1=i-1;
00536         z=-g.elem(i);
00537         double * pd=&(h.elem(i,1));
00538         double * pw=&(w.elem(1));
00539         for (j=1; j<=i1; j++)
00540         {
00541           z-=*pd++ * *pw++;
00542         }
00543         w.elem(i)=z;
00544       }
00545       w.elem(is+n)=w.elem(n)/h.elem(n,n);
00546       {
00547         dvector tmp(1,n);
00548         tmp.initialize();
00549         for (i=1; i<=n1; i++)
00550         {
00551           j=i;
00552           double * pd=&(h.elem(n-j+1,n-1));
00553           double qd=w.elem(is+np-j);
00554           double * pt=&(tmp(1));
00555           for (int ii=1; ii<=n1; ii++)
00556           {
00557             *pt++ +=*pd-- * qd;
00558           }
00559           w.elem(is+n-i)=w.elem(n-i)/h.elem(n-i,n-i)-tmp(i);
00560         }
00561       }/* w(n+1,2n) now contains search direction
00562           with current Hessian approximation */
00563       gs=0.0;
00564       for (i=1; i<=n; i++)
00565          gs+=w.elem(is+i)*g.elem(i);/* gs = -inv(H_k)*g_k*df(x_k+alpha_k*p_k) */
00566       iexit=2;
00567       if(gs >= 0.0)
00568          goto label92; /* exit with error */
00569       gso=gs;
00570       /* compute new step length 0 < alpha <= 1 */
00571       alpha=-2.0*df/gs;
00572       if(alpha > 1.0)
00573         alpha=1.0;
00574       df=f;
00575       tot=0.0;
00576 label30: /* Taking a step, updating x */
00577       iexit=3;
00578       if (ialph) { goto label92; } /* exit if step size too small */
00579       /* exit if number of function evaluations exceeded maxfn */
00580       if( ifn >= maxfn)
00581       {
00582         maxfn_flag=1;
00583         goto label92;
00584       }
00585       else
00586       {
00587         maxfn_flag=0;
00588         iexit=1;
00589       }
00590       if(quit_flag) goto label92;
00591       for (i=1; i<=n; i++)
00592       {
00593         /* w(n+1,2n) has the next direction to go */
00594         z=alpha*w.elem(is+i);
00595         /* new independent vector values */
00596         xx.elem(i)+=z;
00597       }
00598       for (i=1; i<=n; i++)
00599       { /* save previous values and update x return value */
00600         xsave.elem(i)=x.elem(i);
00601         gsave.elem(i)=g.elem(i);
00602         x.elem(i)=xx.elem(i);
00603         fsave = f;
00604       }
00605       fsave = f;
00606       ireturn=2;
00607       if (h.disk_save())
00608       {
00609         cout << "starting hessian save" << endl;
00610         h.save();
00611         cout << "finished hessian save" << endl;
00612       }
00613       return; // end of call1
00614   call2: /* Line search and Hessian update */
00615       if (h.disk_save())
00616       {
00617         cout << "starting hessian restore" << endl;
00618         h.restore();
00619         cout << "finished hessian restore" << endl;
00620       }
00621       /* restore x_k, g(x_k) and g(x_k+alpha*p_k) */
00622       for (i=1; i<=n; i++)
00623       {
00624         x.elem(i)=xsave.elem(i); //x_k
00625         w.elem(i)=g.elem(i);     //g(x_k+alpha*p_k)
00626         g.elem(i)=gsave.elem(i); //g(x_k)
00627       }
00628       fy = f;
00629       f = fsave; /* now fy is a new function value, f is the old one */
00630       ireturn=-1;
00631       /* remember current best function value, gradient and parameter values */
00632       if (fy <= fbest)
00633       {
00634         fbest=fy;
00635         for (i=1; i<=n; i++)
00636         {
00637           x.elem(i)=xx.elem(i);
00638           gbest.elem(i)=w.elem(i);
00639         }
00640       }
00641       /* what to do if CTRL-C keys were pressed */
00642       if (use_control_c==1)
00643       {
00644 #if defined(__BORLANDC__)
00645          if ( kbhit() || ctlc_flag|| ifn == dcheck_flag )
00646 #elif defined(_MSC_VER)
00647          //if ( kbhit() || ifn == dcheck_flag )
00648          if ( _kbhit() || ctlc_flag || ifn == dcheck_flag )
00649 #else
00650          if(ctlc_flag || ifn == dcheck_flag )
00651 #endif
00652          {
00653             int c=0;
00654             if (ifn != dcheck_flag)
00655             {
00656             #if defined(__DJGPP__)
00657               c = toupper(getxkey());
00658             #else
00659               c = toupper(getch());
00660             #endif
00661             }
00662             else
00663               c='C';
00664             if ( c == 'C')
00665             {
00666               for (i=1; i<=n; i++)
00667               {
00668                 x.elem(i)=xx.elem(i);
00669               }
00670               ireturn = 3;
00671               derch(f, x , w, n, ireturn);
00672               return;
00673             }
00674             else if(c=='S')
00675             {
00676               //set convergence criteria to something high to stop now
00677               crit=100000.0;
00678               return;
00679             }
00680             else
00681             {
00682               if ( c == 'Q'|| c == 'N')
00683               {
00684                 quit_flag=c;
00685                 goto label92;
00686               }
00687               else
00688               {
00689                 quit_flag=0;
00690               }
00691             }
00692          }
00693        }
00694        if (quit_flag)
00695        {
00696          if (quit_flag==1) quit_flag='Q';
00697          if (quit_flag==2) quit_flag='N';
00698          goto label92;
00699        }
00700       /* Otherwise, continue */
00701        /* Note, icc seems to be unused */
00702        icc+=1;
00703        if( icc >= 5)
00704           icc=0;
00705       /* ic - counter of calls without x update */
00706       ic++;
00707       if( ic >imax)
00708       {
00709          if (iprint>0)
00710          {
00711            if (ad_printf)
00712              (*ad_printf)("  ic > imax  in fminim is answer attained ?\n");
00713            fmmdisp(x, g, n, this->scroll_flag,noprintx);
00714          }
00715          ihflag=1;
00716          ihang=1;
00717          goto label92;
00718       }
00719       /* new function value was passed to fmin, will increment ifn */
00720       ifn++;
00721 
00722       gys=0.0;
00723 
00724       /* gys = transpose(p_k) * df(x_k+alpha_k*p_k) */
00725       for (i=1; i<= n; i++)
00726          gys+=w.elem(i)*w.elem(is+i);
00727 
00728       /* bad step; unless modified by the user, fringe default = 0 */
00729       if(fy>f+fringe)
00730       {
00731          goto label40; /* backtrack */
00732       }
00733       /* Otherwise verify if the search step length satisfies
00734          strong Wolfe condition */
00735       if(fabs(gys/gso)<=.9)
00736          goto label50; /* proceed to Hessian update */
00737 /* or slightly modified constant in Wolfe condition for the number of
00738  calls > 4 */
00739       if(fabs(gys/gso)<=.95 && ic > 4)
00740          goto label50; /* proceed to Hessian update */
00741       if(gys>0.0) /* not a descent direction */
00742          goto label40; /* backtrack */
00743       tot+=alpha;
00744       z=10.0;
00745       if(gs<gys)
00746          z=gys/(gs-gys);
00747       if(z>10.0)
00748          z=10.0;
00749       /* increase step length */
00750       alpha=alpha*z;
00751       if (alpha == 0.)
00752       {
00753          ialph=1;
00754          #ifdef __ZTC__
00755          if (ireturn <= 0)
00756          {
00757            disp_close();
00758          }
00759          #endif
00760          return;
00761       }
00762       f=fy;
00763       gs=gys;
00764       goto label30; /* update x with new alpha */
00765 label40: /* new step is not acceptable, stepping back and
00766             start backtracking along the Newton direction
00767             trying a smaller value of alpha */
00768       for (i=1;i<=n;i++)
00769          xx.elem(i)-=alpha*w.elem(is+i);
00770       if (alpha == 0.)
00771       {
00772         ialph=1;
00773         return;
00774       }
00775       /* compute new alpha */
00776       z=3.0*(f-fy)/alpha+gys+gs;
00777       zz=dafsqrt(z*z-gs*gys);
00778       z=1.0-(gys+zz-z)/(2.0*zz+gys-gs);
00779       if (fabs(fy-1.e+95) < 1.e-66)
00780       {
00781         alpha*=.001;
00782       }
00783       else
00784       {
00785         if (z>10.0)
00786         {
00787           cout << "large z" << z << endl;
00788           z=10.0;
00789         }
00790         alpha=alpha*z;
00791       }
00792       if (alpha == 0.)
00793       {
00794          ialph=1;
00795         if (ialph)
00796         {
00797           if (ad_printf) (*ad_printf)("\nFunction minimizer: Step size"
00798             "  too small -- ialph=1");
00799         }
00800          return;
00801       }
00802       goto label30; /* update x with new alpha */
00803 label50: /* compute Hessian updating terms */
00804       alpha+=tot;
00805       f=fy;
00806       df-=f;
00807       dgs=gys-gso;
00808       xxlink=1;
00809       if(dgs+alpha*gso>0.0)
00810          goto label52;
00811       for (i=1;i<=n;i++)
00812          w.elem(iu+i)=w.elem(i)-g.elem(i);
00813       /* now w(n+1,2n) = df(x_k+alpha_k*p_k)-df(x_k) */
00814       sig=1.0/(alpha*dgs);
00815       goto label70;
00816 label52: /* compute Hessian updating terms */
00817       zz=alpha/(dgs-alpha*gso);
00818       z=dgs*zz-1.0;
00819       for (i=1;i<=n;i++)
00820          w.elem(iu+i)=z*g.elem(i)+w.elem(i);
00821       sig=1.0/(zz*dgs*dgs);
00822       goto label70;
00823 label60: /* compute Hessian updating terms */
00824       xxlink=2;
00825       for (i=1;i<=n;i++)
00826          w.elem(iu+i)=g.elem(i);
00827       if(dgs+alpha*gso>0.0)
00828          goto label62;
00829       sig=1.0/gso;
00830       goto  label70;
00831 label62: /* compute Hessian updating terms */
00832       sig=-zz;
00833       goto label70;
00834 label65: /* save in g the gradient df(x_k+alpha*p_k) */
00835       for (i=1;i<=n;i++)
00836          g.elem(i)=w.elem(i);
00837       goto  label20; //convergence check
00838 label70:  // Hessian update
00839       w.elem(iv+1)=w.elem(iu+1);
00840       pfmintime->get_elapsed_time_and_reset();
00841       for (i=2;i<=n;i++)
00842       {
00843          i1=i-1;
00844          z=w.elem(iu+i);
00845          double * pd=&(h.elem(i,1));
00846          double * pw=&(w.elem(iv+1));
00847          for (j=1;j<=i1;j++)
00848          {
00849            z-=*pd++ * *pw++;
00850          }
00851          w.elem(iv+i)=z;
00852       }
00853       pfmintime->get_elapsed_time_and_reset();
00854       for (i=1;i<=n;i++)
00855       {  /* BFGS updating formula */
00856          z=h.elem(i,i)+sig*w.elem(iv+i)*w.elem(iv+i);
00857          if(z <= 0.0)
00858             z=dmin;
00859          if(z<dmin)
00860             dmin=z;
00861          h.elem(i,i)=z;
00862          w.elem(ib+i)=w.elem(iv+i)*sig/z;
00863          sig-=w.elem(ib+i)*w.elem(ib+i)*z;
00864        }
00865       for (j=2;j<=n;j++)
00866       {
00867          double * pd=&(h.elem(j,1));
00868          double * qd=&(w.elem(iu+j));
00869          double * rd=&(w.elem(iv+1));
00870          for (i=1;i<j;i++)
00871          {
00872             *qd-=*pd * *rd++;
00873             *pd++ +=w.elem(ib+i)* *qd;
00874          }
00875       }
00876       if (xxlink == 1) goto label60;
00877       if (xxlink == 2) goto label65;
00878 /*label90:*/
00879       for (i=1;i<=n;i++)
00880          g.elem(i)=w.elem(i);
00881 label92: /* Exit with error */
00882       if (iprint>0)
00883       {
00884         if (ialph)
00885         {
00886           if (ad_printf)
00887            (*ad_printf)("\nFunction minimizer: Step size too small -- ialph=1");
00888         }
00889         if (ihang == 1)
00890         {
00891           if (ad_printf)
00892             (*ad_printf)(
00893            "Function minimizer not making progress ... is minimum attained?\n");
00894 #if defined(USE_DDOUBLE)
00895 #undef double
00896            if (ad_printf)
00897            (*ad_printf)("Minimprove criterion = %12.4le\n",double(min_improve));
00898 #define double dd_real
00899 #else
00900            if (ad_printf)
00901              (*ad_printf)("Minimprove criterion = %12.4le\n",min_improve);
00902 #endif
00903         }
00904       }
00905       if(iexit == 2)
00906       {
00907         if (iprint>0)
00908         {
00909           if (ad_printf)
00910             (*ad_printf)("*** grad transpose times delta x greater >= 0\n"
00911            " --- convergence critera may be too strict\n");
00912           ireturn=-1;
00913         }
00914       }
00915 #if defined (_MSC_VER) && !defined (__WAT32__)
00916         if (scroll_flag == 0) clrscr();
00917 #endif
00918       if (maxfn_flag == 1)
00919       {
00920         if (iprint>0)
00921         {
00922           if (ad_printf)
00923             (*ad_printf)("Maximum number of function evaluations exceeded");
00924         }
00925       }
00926       if (iprint>0)
00927       {
00928         if (quit_flag == 'Q')
00929           if (ad_printf) (*ad_printf)("User initiated interrupt");
00930       }
00931       if(iprint == 0) goto label777;
00932       if (ad_printf) (*ad_printf)("\n - final statistics:\n");
00933       if (ad_printf)
00934         (*ad_printf)("%d variables; iteration %ld; function evaluation %ld\n",
00935                        n, itn, ifn);
00936 #if defined(USE_DDOUBLE)
00937 #undef double
00938       if (ad_printf)
00939         (*ad_printf)(
00940              "Function value %12.4le; maximum gradient component mag %12.4le\n",
00941              double(f), double(gmax));
00942       if (ad_printf)
00943         (*ad_printf)(
00944           "Exit code = %ld;  converg criter %12.4le\n",iexit,double(crit));
00945 #define double dd_real
00946 #else
00947       if (ad_printf)
00948         (*ad_printf)(
00949           "Function value %12.4le; maximum gradient component mag %12.4le\n",
00950           f, gmax);
00951       if (ad_printf)
00952         (*ad_printf)(
00953           "Exit code = %ld;  converg criter %12.4le\n",iexit,crit);
00954 #endif
00955       fmmdisp(x, g, n, this->scroll_flag,noprintx);
00956 label777: /* Printing final Hessian approximation */
00957          if (ireturn <= 0)
00958          #ifdef DIAG
00959            if (ad_printf) (*ad_printf)("Final values of h in fmin:\n");
00960            cout << h << "\n";
00961          #endif
00962          #ifdef __ZTC__
00963          {
00964            disp_close();
00965          }
00966          #endif
00967          return;
00968 label7000:/* Printing Initial statistics */
00969       if (iprint>0)
00970       {
00971 #if defined (_MSC_VER) && !defined (__WAT32__)
00972         if (!scroll_flag) clrscr();
00973 #endif
00974         if (ad_printf) (*ad_printf)("\nInitial statistics: ");
00975       }
00976       goto label7003;
00977 label7010:/* Printing Intermediate statistics */
00978    if (iprint>0)
00979    {
00980 #if defined (_MSC_VER) && !defined (__WAT32__)
00981      if (!scroll_flag) clrscr();
00982 #endif
00983      if (ad_printf) (*ad_printf)("\nIntermediate statistics: ");
00984    }
00985    llog=0;
00986    goto label7003;
00987 label7020:/* Exis because Hessian is not positive definite */
00988   if (iprint > 0)
00989   {
00990     if (ad_printf) (*ad_printf)("*** hessian not positive definite\n");
00991   }
00992 #ifdef __ZTC__
00993   if (ireturn <= 0)
00994   {
00995     disp_close();
00996   }
00997 #endif
00998 }
01004 double dafsqrt(double x)
01005 {
01006   return x > 0 ? sqrt(x) : 0.0;
01007 }