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