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