home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / magazine / drdobbs / 1986 / 03 / metric.mar < prev    next >
Text File  |  1986-03-31  |  52KB  |  2,180 lines
  1. /* GLOBAL.H                      08/18/85
  2.  *
  3.  * Global include file for the Variable Metric Minimizer program.
  4.  *
  5.  * (c) Copyright 1985, Billybob Software.  All rights reserved.
  6.  * This program may be reproduced for personal, non-profit use only.
  7.  */
  8.  
  9. #define        C80
  10.  
  11. /*
  12.  * remove the above line if you are not using C/80
  13.  */
  14.  
  15. #ifdef    C80
  16.  
  17. /* 
  18.  * required for the Software Toolworks C/80 Version 3.0 compiler
  19.  */
  20.  
  21. #define        double        float
  22. #include    "fprintf.h"
  23. #include    "scanf.h"
  24.  
  25. /*
  26.  * the following is the <math.h> for C/80
  27.  */
  28.  
  29. float    sin() , cos() , atan() , sqrt() , exp() ,
  30.     pow() , pow10() , ln() , log() , fabs() , atof()  ;
  31.  
  32. #else
  33.  
  34. /*
  35.  * other compilers
  36.  */
  37.  
  38. #include    <stdio.h>
  39. #include    <math.h>
  40. #include    <ctype.h>
  41.  
  42. /*
  43.  * ctype.h may be needed for an isspace used in read.c
  44.  * some systems need it, some don't! 
  45.  * check your math.h for fabs and atof declarations!
  46.  */
  47.  
  48. #endif
  49.  
  50. double    dot() ;
  51. struct    bound
  52. {
  53.     int fl ;
  54.     double up ;
  55.     double lo ;
  56.     double mi ;
  57. };
  58.  
  59. struct    xydata
  60. {
  61.     double x ;
  62.     double y ;
  63. };
  64.  
  65. #define        BOUND        struct bound
  66. #define        DATA        struct xydata
  67. #define     VECMAX        10
  68. #define        MATMAX        VECMAX*VECMAX
  69. #define        NEPS        2
  70. #define        STDES        0
  71. #define        DFP        1
  72. #define        BFGS        2
  73. #define        ALLDONE        -1
  74. #define        NEGEDM        -3000
  75. #define        BADMIN        -2000
  76. #define        BADLS        -1000
  77. #define        MAXITERS    0
  78. #define        OK        0
  79. #define        EDM1        1000
  80. #define        EDM2        2000
  81. #define        TOOSML        3000
  82.  
  83. Listing TWO
  84. /* VMM.C                     11/04/85
  85.  *
  86.  * The main driver routine, function library routine
  87.  *             and all library functions for vmm.
  88.  * (c) Copyright 1985, Billybob Software.  All rights reserved.
  89.  * This program may be reproduced for personal, non-profit use only.
  90.  */
  91.  
  92. #include    "global.h"
  93.  
  94. #define        DMAX    50    /* Size of experimental data array    */
  95. #define        NFUN    3    /* Num of functions in function library    */
  96.  
  97. /************************************************************************/
  98.  
  99. main()
  100. {
  101.     static    int    nparam ;    /* number of parameters        */
  102.     static    double    const[VECMAX] ;    /* constants used in function    */
  103.     static    double    param[VECMAX] ;    /* parameters to be found    */
  104.     static    double    dstep[VECMAX] ; /* step sizes for derivatives    */
  105.     static    double    epsilon[VECMAX];/* stopping criteria        */
  106.     static    BOUND    limit[VECMAX] ;    /* limits if constrained    */
  107.     static    char    *pname[VECMAX] ;/* parameter names         */
  108.     static    double    g0 ;        /* expected value of minimum    */
  109.     static    int    debug ;        /* debug level, 0 if off    */
  110.     static    int    itermin ;    /* minimum number of iterations    */
  111.     static    int    iterlim ;    /* maximum number of iterations    */
  112.     static    int    nreset ;    /* reset y after this many iters*/
  113.     static    DATA    exp[DMAX] ;    /* xy experimental data        */
  114.     static    int    npoints ;    /* number of data points    */
  115.     static    double    *(*fun)() ;    /* ptr to function to minimize    */
  116.     static    int    iters ;        /* number of iterations it took    */
  117.     static    int    retcode ;    /* return code from minimization*/
  118.     static    double    gmin ;        /* the value at the minimum    */
  119.     static     int    method ;    /* method used to update y    */
  120.     static  int    control ;    
  121.     int    i , j ;
  122.  
  123.     for ( i=1 ; ; ++i )
  124.     {
  125.         raz ( &nparam , &itermin , &iterlim , &nreset , &debug ,
  126.             limit , dstep , param , &g0 , epsilon , &method ) ;
  127.  
  128.          if( (control = reader ( &fun , const , &nparam , param , 
  129.                 limit ,    dstep , pname , epsilon , &itermin ,
  130.                 &iterlim , &nreset , &g0 , &debug , exp ,
  131.                 &npoints , DMAX , &method)) == ALLDONE ) 
  132.             break ;
  133.         
  134.         if ( control )
  135.             continue ;
  136.  
  137.         if ( dfault ( nparam , &itermin , &iterlim , &nreset ,
  138.                 epsilon , limit , dstep , debug ) )
  139.             continue ;
  140.  
  141.         retcode = minim ( nparam , fun , param , dstep , limit ,
  142.                 g0 , itermin , iterlim , epsilon , &gmin , 
  143.                 &iters , debug , nreset , exp , npoints ,
  144.                 const , method ) ;
  145.  
  146.         printf("\t***************************\n") ;
  147.         printf("\t* results from dataset %2d *\n", i ) ;
  148.         printf("\t***************************\n") ;
  149.         printf("stopped after %1d iterations, ", iters    ) ;
  150.  
  151.         switch ( retcode )
  152.         {
  153.         case NEGEDM : printf("negative e.d.m.\n")    ; break ;
  154.         case BADMIN : printf("new min > last min\n") ; break ;
  155.         case BADLS  : printf("line search failed\n") ; break ;
  156.         case MAXITERS: printf("maximum iterations\n"); break ;
  157.         case EDM1   : printf("e.d.m. within tol.\n") ; break ;
  158.         case EDM2   : printf("e.d.m. within tol.\n") ; break ;
  159.         case TOOSML : printf("change too small\n")   ; break ;
  160.         default     : printf("unknown reason!!!\n")  ; break ;
  161.         }
  162.  
  163.         printf("value of the minimum was %11.4e\n", gmin ) ;
  164.         printf("parameter values at minimum:\n") ;
  165.  
  166.         for (j=0 ; j<nparam ; ++j) 
  167.             printf("%15s  %11.4e\n" , pname[j] , param[j]);
  168.  
  169.         printf("*********************************************\n") ;
  170.     }
  171. }
  172.  
  173. /************************************************************************
  174.  *    Function library. Note that NFUN is #defined at the top of this
  175.  *    module.
  176.  */
  177.  
  178. funlib ( funname , fun , np , nc , pname )
  179. char    funname[] ;        /* name of function to be minimized    */
  180. int    *fun ;            /* address of function (non-portable)    */
  181. int    *np ;            /* number of parameters in the function */
  182. int    *nc ;            /* number of constants in the function    */
  183. char    *pname[] ;        /* parameter names vector        */
  184. {
  185.     /*
  186.      * use funname to return fun and np, nc from a list
  187.      * for each new function in the library, you must:
  188.      *    (1) declare the function "double" and link it in the load module
  189.      *    (2) declare the function below under "functions available"
  190.      *    (3) add an entry in the function table,including parameter names
  191.      *    (4) add an entry in the switch table in the code
  192.      *    (5) update the value of NFUN
  193.      *
  194.      * Available functions are:
  195.      */
  196.  
  197.     double    cohen() ;
  198.     double    sinus() ;
  199.     double    rosen() ;
  200.  
  201.     /* 
  202.      * function table
  203.      */
  204.  
  205.     static struct
  206.     {
  207.         char *name ;         /* function name in funname call         */
  208.         int  nnpp ;         /* number of parameters in the function */
  209.         int  nncc ;         /* number of constants in the function  */
  210.         char *ppnn[VECMAX]; /* list of names for the parameters         */
  211.  
  212.     } fctlib[NFUN]    = {
  213.  
  214. /*  0  */    { "cohen" ,  2, 0 , "alpha" , "beta"         } ,
  215. /*  1  */    { "sine"  ,  4, 0 , "P0" , "P1" , "P2" , "P3"    } ,
  216. /*  2  */    { "rosen" , 10, 5 , "x1" , "x2" , "x3" , "x4" , "x5" ,
  217.                          "x6" , "x7" , "x8" , "x9" , "x10" }     
  218.     } ;
  219.  
  220.     int    j,  k;
  221.  
  222.     for ( j=0 ; j<NFUN ; ++j ) 
  223.         if ( ! strcmp( funname , fctlib[j].name ) )
  224.             break ;
  225.  
  226.     if (j==NFUN)
  227.         return( 1 ) ;
  228.  
  229.     *np = fctlib[j].nnpp ;
  230.     *nc = fctlib[j].nncc ;
  231.  
  232.     for (k=0 ; k < *np ; ++k)
  233.         pname[k] = fctlib[j].ppnn[k]  ;
  234.  
  235.     switch( j )
  236.     {
  237.     case 0 : *fun = cohen ; break ;
  238.     case 1 : *fun = sinus ; break ;
  239.     case 2 : *fun = rosen ; break ;
  240.     }
  241.  
  242.     return( 0 ) ;
  243. }    
  244.  
  245. /************************************************************************/
  246.  
  247. double    cohen ( const , xx , data , npoints , user )
  248. double    const[] ;    /* vector of constants used by this fcn        */
  249. double    xx[] ;        /* the current values of the vector x        */
  250. DATA    data[] ;    /* dummy                    */
  251. int     npoints ;    /* dummy                    */
  252. int    user ;        /* dummy, reserved for each user!        */
  253. {
  254.     /*
  255.      * computes the function to be minimized 
  256.      * 
  257.      * this is the test function given in Cohen, page 280
  258.      * note that there are no constants for this fcn --
  259.      * the vector const is included for compatability with other fcns
  260.      *
  261.      * note program calls function with user = 0
  262.      */
  263.  
  264.     double  t1 , t2 , t3 , t4 ;
  265.  
  266.     t1 = (xx[1] - xx[0]*xx[0]) ;
  267.     t2 = t1 * t1 ;
  268.     t3 = (xx[0] * (1.0 - xx[1]) + 1.0) ;
  269.     t4 = t3 * t3 ;
  270.     return ( t2 + t4 ) ;
  271. }
  272.  
  273. /************************************************************************/
  274.  
  275. double sinus ( const , xx , data , npoints , user )
  276. double     const[] ;    /* vector of constants for this fcn        */
  277. double    xx[] ;        /* parameter vector                */
  278. DATA    data[] ;    /* data from file                */
  279. int    npoints ;    /* number of data points            */
  280. int    user ;        /* user parameter                */
  281. {
  282.     /*
  283.      * extract parameters for polynomial fit to sine function
  284.      */
  285.  
  286.     double  t1 , t2 , t3 , t4 ;
  287.     int    i;
  288.     
  289.     for ( i=0 , t4 = 0.0 ; i<npoints ; i++ )
  290.     {
  291.         t1 = data[i].x * data[i].x ;
  292.         t2 = data[i].x *
  293.             (xx[0] + t1*(xx[1] + t1*(xx[2] + t1*xx[3]) )) ;
  294.         t3 = (data[i].y - t2)/data[i].y ;
  295.         t4 += t3 * t3 ;
  296.     }
  297.  
  298.     return( t4 /= (double)(npoints-4) ) ;
  299. }
  300.  
  301. /************************************************************************/
  302.  
  303. double    rosen ( const , xx , data , npoints , user )
  304. double    const[] ;    /* vector of constants used by this fcn        */
  305. double    xx[] ;        /* the current values of the vector x        */
  306. DATA    data[] ;    /* dummy                    */
  307. int     npoints ;    /* dummy                    */
  308. int    user ;        /* dummy, reserved for each user!        */
  309. {
  310.     /*
  311.      * Rosenbrock's test function
  312.      */
  313.  
  314.     double  t1 , t2 , t3 ;
  315.     int    l , m ;
  316.  
  317.     for ( l=0 , m=0 , t3 = 0.0 ; l<10 ; l += 2 , ++m )
  318.     {
  319.         t1 = xx[l] ;
  320.         t2 = xx[l+1] ;
  321.         t3 += const[m]*(t2-t1*t1)*(t2-t1*t1) + (1.0-t1)*(1.0-t1) ;
  322.     }
  323.     return( t3 ) ;
  324. }
  325.  
  326. /* LISTING THREE
  327.  * MINIM.C                        11/04/85
  328.  *
  329.  * The main minimization routine
  330.  * (c) Copyright 1985, Billybob Software.  All rights reserved.
  331.  * This program may be reproduced for personal, non-profit use only.
  332.  */
  333.  
  334. #include    "global.h"
  335.  
  336. /***************************************************************************/
  337.  
  338. minim ( n , fun , xextern , xstep , xlim , g0 , itermin , iterlim , epsilon , 
  339.     gmin , iters , debug , nreset , data , npoints , const , method )
  340.  
  341. int    n ;        /* number of parameters in the fit        */
  342. double    (*fun)() ;    /* pointer to the function to minimize        */
  343. double    xextern[] ;    /* vector containing the starting values of
  344.                the parameters; returns values at minimum    */
  345. double    xstep[] ;    /* step sizes on parameters for deriv calcs    */
  346. BOUND    xlim[] ;    /* limits on the vector x            */
  347. double    g0 ;        /* expected value of the minimum of the fcn    */
  348. int    itermin ;    /* minimum number of iterations            */
  349. int    iterlim ;    /* limit on the number of iterations        */
  350. double     epsilon[] ;    /* iteration control parameters            */
  351. double    *gmin ;        /* value of minimum for returned vector x    */
  352. int     *iters ;    /* number of iterations to converge        */
  353. int    debug ;        /* flag = 0 for no print, >0 for debug print    */
  354. int    nreset ;    /* reset limit for the y matrix            */
  355. DATA    data[] ;    /* the data                    */
  356. int    npoints ;    /* number of data points            */
  357. double    const[] ;    /* vector of constants required by fcn        */
  358. int    method ;    /* update method for the matrix y        */
  359. {
  360.     /*
  361.      * variable metric minimization algorithm
  362.      * Reference (using DFP method) is:     Numerical Analysis
  363.      *                    A.M. Cohen et. al.
  364.      *                    Halstead Press
  365.      *                    John Wiley and Sons
  366.      *                    1973
  367.      *                    pages 276-280
  368.      * See also original papers by Fletcher and Powell
  369.      */
  370.  
  371.     int    i , j , rc ;
  372.     static    int    ycount ;    /* counter for resetting the Y matrix     */
  373.     static    double    y[MATMAX] ; /* the Y matrix                  */
  374.     static    double    a[MATMAX] ; /* the A and B matrices are calculated to */
  375.     static    double    b[MATMAX] ; /* update the Y matrix at each iteration  */
  376.     static    double    xintern[VECMAX];/* parameters in internal coordinates */
  377.     static    double    zint[VECMAX] ;    /* gradients in internal coordinates  */
  378.     static    double    znewi[VECMAX] ;    /* same for the next iteration           */
  379.     static    double    sigma[VECMAX] ;    /* changes vector              */
  380.     static    double    s[VECMAX] ;    /* changes direction vector          */
  381.     static    double    xi[VECMAX] ;    /* change in gradient vector between  */
  382.                      * iterations                  */
  383.     static    double    xold[VECMAX] ;    /* holds previous values of x            */
  384.     static    double    alpha ;        /* change scaling parameter          */
  385.     static    double    g ;        /* value of fcn for a given vector x  */
  386.     static    double    gnew ;        /* likewise for the next iteration    */
  387.     double    dot() ;
  388.  
  389.     /*
  390.      * get started
  391.      */
  392.  
  393.     g = (*fun)( const , xextern , data , npoints , 0 ) ;
  394.  
  395.     gozinta( n , xintern , xextern , xlim ) ;
  396.     derivs ( n , fun , xintern , xstep , xlim , zint , data , 
  397.                         npoints , const , debug);
  398.     /*
  399.      * main loop
  400.      */
  401.  
  402.     for ( i = 1 , ycount = nreset ; i <= iterlim ; ++i , ++ycount )
  403.     {
  404.         if (ycount == nreset)
  405.         {
  406.             reset( n , y );
  407.             ycount = 0 ;
  408.         }
  409.         if (debug)
  410.         {
  411.             printf("\t\t++++++++++++++++++++++\n") ;
  412.             printf("\t\t+ Begin iteration %2d +\n", i ) ;
  413.             printf("\t\t++++++++++++++++++++++\n") ;
  414.             printf("external parameter values are:\n") ;
  415.             vout( n , xextern ) ;
  416.             printf("internal parameter values are:\n") ;
  417.             vout( n , xintern ) ;    
  418.             printf("internal gradient vector is:\n") ;
  419.             vout( n , zint ) ;
  420.             printf("approximate inverse hessian matrix is:\n") ;
  421.             mout( n , y ) ;
  422.             printf("the current minimum is --->%11.4e<---\n", g ) ;
  423.         }
  424.  
  425.         matvec ( n , y , zint , s ) ;
  426.         for (j=0 ; j<n ; ++j)
  427.         {
  428.             xold[j] = xintern[j] ;
  429.             s[j] *= -1.0 ;
  430.         }
  431.  
  432.         if (debug>1)
  433.         {
  434.             printf("vector s is:\n");    vout( n , s ) ;
  435.         }
  436.  
  437.         /*
  438.          * do line search and quit if it fails
  439.          * line search was successful if getalpha returns OK ( zero ) 
  440.          */
  441.  
  442.         if ( rc = getalpha ( n , fun , xintern , xstep , xlim , zint , 
  443.              s , g ,g0 , &alpha , debug , data , npoints , const ) )
  444.                         break ;
  445.         /*
  446.          * update parameters
  447.          */
  448.  
  449.         for (j=0 ; j<n ; ++j)
  450.         {
  451.             sigma[j] = alpha * s[j] ;
  452.             xintern[j] += sigma[j];
  453.         }
  454.         if (debug>1)
  455.         {
  456.             printf("vector sigma is:\n") ;    vout( n , sigma ) ;
  457.         }
  458.  
  459.         /*
  460.          * get values for next iteration
  461.          */
  462.  
  463.         gozouta ( n , xintern , xextern , xlim );
  464.         gnew = (*fun)( const , xextern , data , npoints , 0 ) ;
  465.  
  466.         derivs (n , fun , xintern , xstep , xlim , znewi ,
  467.                     data , npoints , const , debug ) ;
  468.  
  469.         for (j=0 ; j<n ; ++j)
  470.             xi[j] = znewi[j] - zint[j] ;
  471.  
  472.         if (debug>1)
  473.         {
  474.             printf("orthogonality satisfied to %11.4e\n",
  475.                     dot ( n , znewi , sigma ) ) ;
  476.             printf("vector xi is:\n") ;
  477.             vout( n , xi ) ;
  478.         }
  479.  
  480.         /*
  481.          * decide whether to continue or not
  482.          * continue if decide returns OK ( zero )
  483.          */
  484.  
  485.         if ( rc = decide( g , gnew , n , znewi , sigma , y ,
  486.                     epsilon , itermin , i , debug ) )
  487.             break ;
  488.         /*
  489.          * update all other quantities for the next pass
  490.          * three methods available:
  491.          * steepest descent (slowest and included only for completeness)
  492.          * Davidon Fletcher Powell (default, empirically the favorite)
  493.          * Broyden Fletcher Goldfarb Shanno (thought to be better by
  494.          * some)
  495.          * The method is selected by keyword on input.
  496.          */
  497.  
  498.         if ( method == STDES )
  499.             ; /* nothing to do, y is not updated */
  500.  
  501.         else if ( method == DFP )
  502.         {
  503.             dfpa ( n , sigma , xi , a , debug ) ;
  504.             dfpb ( n ,   y   , xi , b , debug ) ;
  505.             gety ( n ,   a   ,  b , y , debug ) ;
  506.         }
  507.         else if ( method == BFGS )
  508.         {
  509.             bfgsa ( n , y , sigma , xi , a , debug ) ;
  510.             bfgsb ( n , y , sigma , xi , b , debug ) ;
  511.             gety  ( n , a ,   b   ,  y ,     debug ) ;
  512.         }
  513.  
  514.         for (j=0 ; j<n ; ++j)
  515.             zint[j] = znewi[j] ;
  516.  
  517.         gozouta ( n , xintern , xextern , xlim ) ;
  518.         g = gnew ;
  519.  
  520.         /*
  521.          * all done with this pass
  522.          */
  523.     }
  524.  
  525.     /*
  526.      * outside of the main loop here
  527.      */
  528.  
  529.     *iters = ( i > iterlim ) ? iterlim : i ;
  530.  
  531.     /*
  532.      * first handle the case of abnormal exits
  533.      * these have returned rc negative above, from either getalpha or decide
  534.      */
  535.  
  536.     if ( rc < 0 )
  537.     {
  538.         for (j=0; j<n; ++j)
  539.             xintern[j] = xold[j] ;
  540.         gozouta ( n , xintern , xextern , xlim ) ;
  541.         *gmin = g ;
  542.     }
  543.     else
  544.     {
  545.         /*
  546.          * regular exit and fall through These represent positive or
  547.          * zero values of rc; zero means max iterations.
  548.          */
  549.  
  550.         gozouta ( n , xintern , xextern , xlim ) ;
  551.         *gmin = gnew ;
  552.     }
  553.  
  554.     if (debug)
  555.     {
  556.         printf("\t\t++++++++++++++++++++++++\n") ;
  557.         printf("\t\t+ exiting minimization +\n") ;
  558.         printf("\t\t++++++++++++++++++++++++\n") ;
  559.         printf("\tnumber of iterations was: %11d\n",*iters) ;
  560.         printf("\t          return code is: %11d\n",rc) ;
  561.         printf("\t              minimum is: %11.4e\n", *gmin) ;
  562.         printf("best parameter values:\n") ;
  563.         vout( n , xextern ) ;
  564.         printf("last internal gradient vector:\n") ;
  565.         vout( n , znewi ) ;    
  566.     }
  567.     return( rc ) ;
  568. }
  569.  
  570. /* LISTING FOUR
  571.  * AUXL.C                      11/04/85
  572.  *
  573.  * The line search and decision routines
  574.  * (c) Copyright 1985, Billybob Software.  All rights reserved.
  575.  * This program may be reproduced for personal, non-profit use only.
  576.  */
  577.  
  578. #include    "global.h"
  579.  
  580. /************************************************************************/
  581.  
  582. getalpha ( n , fun , xx , xstep , xlim , z , s , g , g0 ,
  583.          alpha , debug , data , npoints , const )
  584. int    n ;        /* number of parameters                */
  585. double    (*fun)() ;    /* pointer to the function to minimize        */
  586. double    xx[] ;        /* parameter vector                */
  587. double  xstep[] ;    /* stepsize for parameter derivative calcs    */
  588. BOUND    xlim[] ;    /* limits on parameters                */
  589. double    z[] ;        /* gradient vector                */
  590. double  s[] ;        /* changes direction vector            */
  591. double    g ;        /* value of fcn for the given vector xx        */
  592. double  g0 ;        /* expected value of the minimum of the fcn    */
  593. double    *alpha ;    /* scaling parameter we are trying to find    */
  594. int    debug ;        /* flag = 0 for no print, >0 for debug print    */
  595. DATA    data[] ;    /* the data                    */
  596. int    npoints ;    /* number of data points            */
  597. double  const[] ;    /* constants vector needed by fcn        */
  598. {
  599.     /*
  600.      * find alpha according to Cohen, pp 279-280
  601.      */
  602.  
  603.     int    i ;
  604.     static    double    u ;        /* see Cohen            */
  605.     static    double    eta ;        /* see Cohen            */
  606.     static    double  t[VECMAX] ;    /* see Cohen            */
  607.     static    double  text[VECMAX] ;  /* same, in external coordinates    */
  608.     static    double    zt[VECMAX] ;    /* the z for fcn evaluated at x = t */
  609.     static    double    gt ;        /* the g that results from fcn at   */
  610.                     /*            x = t        */
  611.     static    double  nu ;        /* see Cohen            */
  612.     static    double    d ;        /* see Cohen            */
  613.     static    double  w ;        /* see Cohen            */
  614.     static    double    temp1 ;
  615.     static    double    temp2 ;
  616.  
  617.     if (debug>1) {
  618.         printf("\t\t+++++++++++++++++++++++++\n") ;
  619.         printf("\t\t+ line search procedure +\n") ;
  620.         printf("\t\t+++++++++++++++++++++++++\n") ;
  621.     }
  622.  
  623.     u = dot ( n , z , s ) ;
  624.     eta = -2.0 * ( g - g0 ) / u ;
  625.  
  626.     if (eta > 1.0)
  627.         eta = 1.0 ;
  628.  
  629.     for (i=0 ; i<n ; ++i)
  630.         t[i] = xx[i] + eta * s[i] ;
  631.  
  632.     gozouta ( n , t , text , xlim ) ;
  633.  
  634.     gt = (*fun)( const , text , data , npoints , 0 ) ;
  635.     derivs ( n, fun, t, xstep, xlim, zt, data, npoints, const, debug );
  636.     nu = dot ( n , zt , s ) ;
  637.  
  638.     if (debug>1)
  639.     {
  640.         printf("u and nu are: %11.4e %11.4e\n", u , nu ) ;
  641.         printf("g and gt are: %11.4e %11.4e\n", g , gt ) ;
  642.         printf("eta is: %11.4e\n", eta ) ;
  643.         if (debug>2)
  644.         {
  645.             printf("vector t is:\n") ;
  646.             vout( n , t ) ;
  647.             printf("vector zt is:\n") ;
  648.             vout( n , zt ) ;
  649.         }        
  650.     }
  651.  
  652.     d = 3*(g - gt)/eta + u + nu ;
  653.  
  654.     if ( (temp1 = d*d) < (temp2 = u*nu) )
  655.     {
  656.         if (debug>1)
  657.         {
  658.             printf("cubic interpolation has failed, ") ;
  659.             printf("about to take sqrt of negative number!\n") ;
  660.             printf("d squared is %11.4e ", temp1 ) ;
  661.             printf("and u*nu is %11.4e\n", temp2 ) ;
  662.             printf("difference is %11.4e\n", temp1-temp2 ) ;
  663.         }
  664.         if ( u < 0.0  &&  nu < 0.0  &&  gt < g )
  665.         {
  666.             *alpha = eta ;
  667.             if (debug>1)
  668.                 printf("return with alpha = eta\n") ;
  669.             return( OK ) ;
  670.         }
  671.         else
  672.         {
  673.             *alpha = 0.0 ;
  674.             if (debug>1)
  675.                 printf("line search has failed\n") ;
  676.             return( BADLS ) ;
  677.         }
  678.     }
  679.     else
  680.     {
  681.         w = sqrt(temp1 - temp2) ;
  682.         *alpha = eta * ( 1.0 - (nu+w-d)/(nu-u+2.0*w) ) ;
  683.         if (debug>1)
  684.             printf("alpha returned is %11.4e\n", *alpha ) ;
  685.         return( OK ) ;
  686.     }
  687. }
  688.  
  689. /***************************************************************************/
  690.  
  691. decide ( g , gnew , n , znewi , sigma , y , 
  692.         epsilon , itermin , iterations , debug )
  693.  
  694. double    g ;            /* old value of minimum            */
  695. double    gnew ;            /* new value of minimum            */
  696. int    n ;            /* number of parameters            */
  697. double    znewi[] ;        /* vector of new gradients        */
  698. double    sigma[] ;        /* change vector for parameters        */
  699. double    y[] ;            /* the transformation matrix        */
  700. double    epsilon[] ;        /* cutoff criteria vector        */
  701. int    itermin ;        /* minimum number of iterations        */
  702. int    iterations ;        /* current iteration number        */
  703. int     debug ;            /* debug flag                */
  704. {
  705.     /*
  706.      * makes the decision to continue iterating or not
  707.      * returns a zero if we want to keep going
  708.      * returns a positive number if we have a normal stopping reason
  709.      * returns a negative number if we have a catastrophic stopper
  710.      */
  711.  
  712.     static    double    edm ;
  713.     static    double    temp[VECMAX] ;
  714.  
  715.     if (debug>1)
  716.     {
  717.         printf("\t\t+++++++++++++++++++++++++\n") ;
  718.         printf("\t\t+ decision making logic +\n") ;
  719.         printf("\t\t+++++++++++++++++++++++++\n") ;
  720.         printf("old minimum was: %11.4e, new minimum is %11.4e\n",
  721.                  g , gnew ) ;
  722.         printf("ratio new/old is: %11.4e\n", gnew/g ) ;
  723.     }    
  724.  
  725.     /*
  726.      * exit if the new "minimum" is greater than or equal to the old minimum
  727.      */
  728.  
  729.     if (gnew >= g)
  730.     {
  731.         if (debug>1)
  732.             printf("STOP!!!   new minimum is not lower!\n") ;
  733.         return( BADMIN ) ;
  734.     }
  735.  
  736.     /*
  737.      * if the edm (estimated distance to the minimum) is negative we have
  738.      * a catastrophic problem and should stop
  739.      */
  740.  
  741.     matvec( n , y , znewi , temp ) ;
  742.       edm = dot( n , znewi , temp ) ;
  743.  
  744.     if( edm < 0.0 )
  745.     {
  746.         if (debug>1)
  747.             printf("STOP!!!   edm is negative = %11.4e\n", edm);
  748.         return( NEGEDM ) ;
  749.     }
  750.  
  751.     /*
  752.      * now look at normal exits if the minimum number of iterations has
  753.      * been accomplished
  754.      */
  755.  
  756.     if (iterations < itermin) {
  757.         if (debug>1) printf("-->keep going, too few iterations\n") ;
  758.         return( OK ) ;
  759.     }
  760.  
  761.     /*
  762.      * edm test: two ways to calculate, stop if either satisfies
  763.      */
  764.  
  765.     if (edm < epsilon[0])
  766.     {
  767.         if (debug>1)
  768.             printf("STOP, close enough, edm = %11.4e\n", edm);
  769.         return( EDM1 ) ;
  770.     }
  771.  
  772.     edm = dot( n , sigma , sigma ) ;
  773.  
  774.     if (edm < epsilon[0])
  775.     {
  776.         if (debug>1)
  777.             printf("STOP, close enough, edm = %11.4e\n", edm);
  778.         return( EDM2 ) ;
  779.     }
  780.  
  781.     /*
  782.      * % change in g less than "something" means that approach is too slow
  783.      */
  784.  
  785.     if ( g-gnew < epsilon[1]*g )
  786.     {
  787.         if (debug>1) 
  788.             printf("STOP, too slow, fractional change = %11.4e\n",
  789.                     (g-gnew)/g ) ;
  790.         return( TOOSML ) ;
  791.     }
  792.  
  793.     /*
  794.      * fall through case
  795.      * note that case of maximum number of iterations is handled in the
  796.      * calling program
  797.      */
  798.  
  799.     if (debug>1)
  800.         printf("--->keep going, no stoppers found...\n") ;
  801.  
  802.     return( OK ) ;
  803. }
  804.  
  805. /* LISTING FIVE
  806.  * UP.C                     11/04/85
  807.  *
  808.  * Routines for updating the matrix y
  809.  * (c) Copyright 1985, Billybob Software.  All rights reserved.
  810.  * This program may be reproduced for personal, non-profit use only.
  811.  */
  812.  
  813. #include    "global.h"
  814.  
  815. /*-----------------------------------------------------------------*/
  816.  
  817. dfpa ( n , sigma , xi , a , debug )
  818. int    n ;        /* number of parameters                */
  819. double    *sigma ;    /* changes vector                */
  820. double    *xi ;        /* change in gradient vector            */
  821. double    *a ;        /* the result                    */ 
  822. int    debug ;        /* flag = 0 for no print, >0 for debug print    */
  823. {
  824.     /*
  825.      * Compute the matrix A which is used to correct Y
  826.      * Davidon Fletcher Powell method
  827.      */
  828.  
  829.     int    i ;
  830.     double    norm ;
  831.     double    *t ;
  832.  
  833.     t = a ;
  834.     cross( n , sigma , sigma , a ) ;
  835.     norm = 1.0 / dot( n , sigma , xi ) ;
  836.     i = n*n ;    
  837.  
  838.     while( i-- )
  839.         *a++ *= norm ;
  840.  
  841.     if (debug>1)
  842.     {
  843.         printf("AAAAA correcting matrix    AAAAA is:\n") ;
  844.         mout( n , t ) ;
  845.     }
  846. }
  847.  
  848. /*-----------------------------------------------------------------*/
  849.  
  850. dfpb ( n , y , xi , b , debug )
  851. int    n ;        /* number of parameters                */
  852. double    *y ;        /* current Y matrix                */
  853. double    *xi ;        /* change in gradient vector            */
  854. double    *b ;        /* the result                    */
  855. int    debug ;        /* flag = 0 if no print, >0 for debug print    */
  856. {
  857.     /*
  858.      * compute the matrix B which is used to update Y
  859.      * Davidon Fletcher Powell method
  860.      */
  861.  
  862.     int    i ;
  863.     static    double    temp[VECMAX] ;
  864.     double    *t ;
  865.     double    norm ;
  866.  
  867.     t = b ;
  868.     matvec( n , y , xi , temp ) ;
  869.     norm = - 1.0 / dot( n , temp , xi ) ;
  870.     cross( n , temp , temp , b ) ;
  871.     i = n*n ;
  872.  
  873.     while( i-- )
  874.         *b++ *= norm ; 
  875.  
  876.     if (debug>1)
  877.     {
  878.         printf("BBBBB correcting matrix    BBBBB is:\n") ;
  879.         mout( n , t ) ;
  880.     }
  881. }
  882.  
  883. /*-----------------------------------------------------------------*/
  884.  
  885. bfgsa ( n , y , sigma , xi , a , debug )
  886. int    n ;        /* number of parameters                */
  887. double    *y ;        /* current Y matrix                */
  888. double    *sigma ;    /* changes vector                */
  889. double    *xi ;        /* change in gradient vector            */
  890. double    *a ;        /* the result                    */ 
  891. int    debug ;        /* flag = 0 for no print, >0 for debug print    */
  892. {
  893.     /*
  894.      * compute the matrix A which is used to correct Y
  895.      * Broyden Fletcher Goldfarb Shanno method
  896.      */
  897.  
  898.     int    i , j ;
  899.     static    double    temp[VECMAX] ;
  900.     static     double    tmpa[MATMAX] ;
  901.     double    norm ;
  902.  
  903.     norm = 1.0 / dot ( n , sigma , xi ) ;
  904.     matvec( n , y , xi , temp ) ;
  905.  
  906.     for (j=0 ; j<n ; ++j)
  907.         temp[j] = sigma[j] - temp[j] ;
  908.  
  909.     cross( n , temp , sigma , tmpa ) ;
  910.     cross( n , sigma , temp ,  a   ) ;
  911.     i = n*n ;    
  912.  
  913.     for (j=0 ; j<i ; ++j)
  914.         a[j] = norm * ( a[j] + tmpa[j] ) ;
  915.  
  916.     if (debug>1)
  917.     {
  918.         printf("AAAAA correcting matrix    AAAAA is:\n") ;
  919.         mout( n , a ) ;
  920.     }
  921. }
  922.  
  923. /*-----------------------------------------------------------------*/
  924.  
  925. bfgsb ( n , y , sigma , xi , b , debug )
  926. int    n ;        /* number of parameters                */
  927. double    *y ;        /* current Y matrix                */
  928. double    *sigma ;    /* changes vector                */
  929. double    *xi ;        /* change in gradient vector            */
  930. double    *b ;        /* the result                    */
  931. int    debug ;        /* flag = 0 if no print, >0 for debug print    */
  932. {
  933.     /*
  934.      * compute the matrix B which is used to update Y
  935.      * Broyden Fletcher Goldfarb Shanno method
  936.      */
  937.  
  938.     int    i ;
  939.     static    double    temp[VECMAX] ;
  940.     double    *t ;
  941.     double    norm ;
  942.  
  943.     t = b ;
  944.     norm = 1.0 / dot( n , sigma , xi ) ;
  945.     norm = - norm * norm ;    
  946.     matvec( n , y , xi , temp ) ;
  947.  
  948.     for (i=0 ; i<n ; ++i)
  949.         temp[i] = sigma[i] - temp[i] ;
  950.  
  951.     norm *= dot( n , temp , xi ) ;
  952.     cross( n , sigma , sigma , b ) ;         
  953.     i = n*n ;
  954.  
  955.     while( i-- )
  956.         *b++ *= norm ; 
  957.  
  958.     if (debug>1)
  959.     {
  960.         printf("BBBBB correcting matrix    BBBBB is:\n") ;
  961.         mout( n , t ) ;
  962.     }
  963. }
  964.  
  965. /*-----------------------------------------------------------------*/
  966.  
  967. gety ( n , a , b , y , debug )
  968. int    n ;        
  969. double    *a ;
  970. double    *b ;
  971. double    *y ;
  972. int    debug ;
  973. {
  974.     /*
  975.      * use matrices A and B to correct Y
  976.      */
  977.     int    i ;
  978.     double    *t ;
  979.  
  980.     t = y ;
  981.     i = n*n ;
  982.  
  983.     while( i-- )
  984.         *y++ += *a++ + *b++ ;
  985.  
  986.     if (debug>1)
  987.     {
  988.         printf("YYYYY new matrix YYYYY is:\n") ;
  989.         mout( n , t ) ;
  990.     }
  991. }
  992.  
  993.  * LISTING SIX
  994.  * UTIL.C                         11/04/85
  995.  *
  996.  * Contains miscellaneous routines needed for variable metric minimization 
  997.  * program, including matrix and vector utilities
  998.  * (c) Copyright 1985, Billybob Software.  All rights reserved.
  999.  * This program may be reproduced for personal, non-profit use only.
  1000.  */
  1001.  
  1002. #include    "global.h"
  1003. #define        INALINE        6
  1004.  
  1005. /*-----------------------------------------------------------------*/
  1006.  
  1007. gozinta( n , xintern , xextern , xlim )
  1008. int    n ;
  1009. double    xintern[] ;
  1010. double    xextern[] ;
  1011. BOUND    xlim[] ;
  1012. {
  1013.     /*
  1014.      * transform all external coordinates to internal coordinates
  1015.      */
  1016.  
  1017.     int    i ;
  1018.  
  1019.     for ( i=0 ; i<n ; ++i )
  1020.         tointern( i , xintern , xextern , xlim ) ;
  1021. }
  1022.  
  1023. /*-----------------------------------------------------------------*/
  1024.  
  1025. gozouta( n , xintern , xextern , xlim )
  1026. int    n ;
  1027. double    xintern[] ;
  1028. double    xextern[] ;
  1029. BOUND    xlim[] ;
  1030. {
  1031.     /*
  1032.      * transform all internal coordinates to external coordinates
  1033.      */
  1034.  
  1035.     int    i ;
  1036.  
  1037.     for ( i=0 ; i<n ; ++i )
  1038.         toextern( i , xintern , xextern , xlim ) ;
  1039. }
  1040.  
  1041. /*-----------------------------------------------------------------*/
  1042.  
  1043. tointern ( j , xintern , xextern , xlim )
  1044. int    j ;            /* parameter to transform        */
  1045. double    xintern[] ;        /* internal coordinates            */
  1046. double    xextern[] ;        /* external coordinates            */
  1047. BOUND    xlim[] ;        /* limits on the parameters        */
  1048. {
  1049.     /*
  1050.      * transform to unbounded "internal" coordinates used by
  1051.      * minimization program.
  1052.      */
  1053.  
  1054.     double    y ;
  1055.  
  1056.     if( xlim[j].fl )
  1057.     {
  1058.         if ( xlim[j].mi == 0.0 )
  1059.             xintern[j] = xextern[j] ;
  1060.         else
  1061.         {
  1062.             y = ( (xextern[j]-xlim[j].lo) / xlim[j].mi ) - 1.0 ;
  1063.             xintern[j] = atan( y / sqrt( 1.0 - y*y ) ) ;
  1064.         }
  1065.     }
  1066.     else
  1067.         xintern[j] = xextern[j] ;
  1068. }
  1069.  
  1070. /*-----------------------------------------------------------------*/
  1071.  
  1072. toextern ( j , xintern , xextern , xlim )
  1073. int    j ;            /* parameter to transform        */
  1074. double    xintern[] ;        /* internal coordinates            */
  1075. double    xextern[] ;         /* external coordinates            */
  1076. BOUND    xlim[] ;        /* limits on the parameters        */
  1077. {
  1078.     /*
  1079.      * transforms to bounded "external" coordinates known to the real world
  1080.      */
  1081.  
  1082.     if( xlim[j].fl )
  1083.         xextern[j] = xlim[j].lo + xlim[j].mi *
  1084.                         ( sin(xintern[j]) + 1.0 ) ;
  1085.     else
  1086.         xextern[j] = xintern[j] ;
  1087. }
  1088.  
  1089. /*-----------------------------------------------------------------*/
  1090.  
  1091. derivs ( n , fun , xintern , xstep , xlim , z , data , m , const , debug )
  1092. int    n ;            /* number of parameters            */
  1093. double    (*fun)() ;        /* pointer to the function to minimize    */  
  1094. double    xintern[] ;        /* coordinates vector            */
  1095. double    xstep[] ;         /* step size to take on each component    */
  1096. BOUND    xlim[] ;        /* limit vector                */
  1097. double    z[] ;            /* derivative vector            */
  1098. DATA    data[] ;        /* the data                 */
  1099. int    m ;            /* number of data points        */
  1100. double    const[] ;        /* constants required by fcn        */
  1101. int    debug ;            /* debug flag                */
  1102. {
  1103.     /*
  1104.      * compute the derivatives of the vector x using a simple
  1105.      * finite difference.
  1106.      */
  1107.  
  1108.     int     i ;
  1109.     static    double    xextern[VECMAX];/* throwaway vector in ext coords */
  1110.     static    double    f1 , f2 ;    /* values of fcn at +/- stepsize  */
  1111.     static    double    xtemp ;    
  1112.     static    double    eps ;
  1113.  
  1114.     gozouta( n , xintern , xextern , xlim ) ;
  1115.     for ( i=0 ; i<n ; ++i )
  1116.     {
  1117.         xtemp        = xintern[i] ;
  1118.         eps        = fabs( xtemp ) * xstep[i] ;
  1119.         xintern[i] = xtemp + eps ;
  1120.  
  1121.         toextern( i , xintern , xextern , xlim ) ;
  1122.         f1 = (*fun)( const , xextern , data , m , 0 ) ;
  1123.         xintern[i] = xtemp - eps ;
  1124.  
  1125.         toextern( i , xintern , xextern , xlim ) ;
  1126.         f2 = (*fun)( const , xextern , data , m , 0 ) ;
  1127.         z[i] = (f1 - f2) / ( 2.0 * eps ) ;
  1128.  
  1129.         if (debug>2)
  1130.                 printf("i,f1,f2,step,deriv %d %12e %12e %12e %12e\n",
  1131.                     i, f1, f2, eps , z[i] ) ;
  1132.         xintern[i] = xtemp ;
  1133.         toextern( i , xintern , xextern , xlim ) ;
  1134.     }
  1135.  
  1136.     if (debug==2)
  1137.     {
  1138.         printf("derivatives are:\n");
  1139.         vout( n , z ) ;
  1140.     }
  1141. }
  1142.  
  1143. /*-----------------------------------------------------------------*/
  1144.  
  1145. raz ( nparam , itermin , iterlim , nreset , debug , limit , dstep ,
  1146.         param , g0 , epsilon , method )
  1147. int    *nparam ;        /* number of parameters        */
  1148. int    *itermin ;        /* minimum number of iterations    */
  1149. int    *iterlim ;        /* maximum number of iterations    */
  1150. int    *nreset ;        /* reset y after nreset iters    */
  1151. int    *debug ;        /* debug flag            */
  1152. BOUND    limit[] ;        /* limit vector            */
  1153. double    dstep[] ;        /* step size for derivatives    */
  1154. double    param[] ;        /* initial values of parameters    */
  1155. double    *g0 ;            /* expected value of minimum    */
  1156. double    epsilon[] ;        /* stopping criteria vector    */
  1157. int    *method ;        /* update method for y        */
  1158. {
  1159.     /*
  1160.      * reinitialize important parameters before reading, so that we
  1161.      * can later test to see if they were set in the dataset
  1162.      */
  1163.  
  1164.     int    j ;
  1165.  
  1166.     *debug  = *nparam = *itermin = *iterlim = *nreset = 0 ;
  1167.     *g0     = 0.0 ;
  1168.     *method = DFP ;
  1169.  
  1170.     for ( j=0 ; j<VECMAX ; ++j )
  1171.     {
  1172.         limit[j].fl = 0   ;
  1173.         dstep[j]    = 0.0 ;
  1174.         epsilon[j]  = 0.0 ;
  1175.         param[j]    = 0.0 ;
  1176.     }
  1177. }
  1178.  
  1179. /*-----------------------------------------------------------------*/
  1180.  
  1181. dfault ( n , itermin , iterlim , nreset , epsilon , xlim , dstep , debug )
  1182. int    n ;            /* number of parameters            */
  1183. int    *itermin ;        /* minimum number of iterations        */
  1184. int    *iterlim ;        /* maximum number of iterations        */
  1185. int    *nreset ;        /* reset y matrix after this many iters */
  1186. double    epsilon[] ;        /* stopping criteria vector        */
  1187. BOUND    xlim[] ;        /* limit vector for parameters         */
  1188. double    dstep[] ;        /* step sizes for derivative calc    */
  1189. int    debug ;            /* flag for debug            */
  1190. {
  1191.     /*
  1192.      * this routine sets defaults if parameters not set with input data 
  1193.      */
  1194.  
  1195.     int    j ;
  1196.  
  1197.     if (debug)
  1198.     {
  1199.         printf("\t\t+++++++++++++++++++\n") ;
  1200.         printf("\t\t+ initializations +\n") ;
  1201.         printf("\t\t+++++++++++++++++++\n") ;
  1202.     }
  1203.  
  1204.     /*
  1205.      * check that the number of parameters has been set; fatal error if not!
  1206.      */
  1207.  
  1208.     if ( !n )
  1209.         return( -1 ) ;
  1210.  
  1211.     /*
  1212.      * set up other defaults as required -- these depend on knowing n
  1213.      */
  1214.  
  1215.     if ( !*itermin )
  1216.         *itermin = n ;
  1217.  
  1218.     if ( !*iterlim )
  1219.         *iterlim = 2*n ;
  1220.  
  1221.     if ( !*nreset  )
  1222.         *nreset = (3*n)/2 ;
  1223.  
  1224.     for (j=0 ; j<n ; ++j)
  1225.         if (dstep[j] == 0.0)
  1226.             dstep[j] = 0.01 ;
  1227.     if (debug)
  1228.     {
  1229.         printf("min iters = %3d, max iters = %3d, ", *itermin ,
  1230.                         *iterlim ) ;
  1231.         printf("reset y every %3d iterations\n", *nreset ) ;
  1232.         printf("step sizes for derivatives:\n") ;  vout( n , dstep ) ;
  1233.     }
  1234.  
  1235.     if ( epsilon[0] == 0.0 )
  1236.     {
  1237.         epsilon[0] = 1.0e-06 ;
  1238.         epsilon[1] = 0.001 ;
  1239.         if (debug) printf("epsilons set to defaults:\n") ; 
  1240.     }
  1241.     else if (debug)
  1242.         printf("epsilons from input data:\n") ;
  1243.  
  1244.     if (debug)
  1245.         vout( NEPS , epsilon ) ;
  1246.  
  1247.     /*
  1248.      * note that if you add other stopping criteria (more elements in
  1249.      * the epsilon vector) you will have to modify this code
  1250.      *
  1251.      * the following defaults were set in raz and remain if not changed
  1252.      * by the data read in reader:
  1253.      *    starting values of parameters:    0.0
  1254.      *    expected value of minimum:         0.0
  1255.      *    constrained/unconstrained:    unconstrained (all)
  1256.      *    parameter names:        blank
  1257.      *    updating method for y :        Davidon-Fletcher-Powell (DFP)
  1258.      *
  1259.      * variables used for intern<-->extern conversions:
  1260.      */
  1261.  
  1262.     for (j=0 ; j<n ; j++) 
  1263.         if (xlim[j].fl)
  1264.             xlim[j].mi = (xlim[j].up - xlim[j].lo) / 2.0 ;
  1265.  
  1266.     return( 0 ) ;
  1267. }
  1268.  
  1269. /*-----------------------------------------------------------------*/
  1270.  
  1271. vout ( n , a )
  1272. int    n ;
  1273. double    *a ;
  1274. {
  1275.     /* 
  1276.      * output the floating point vector a with n components 
  1277.      * INALINE values to a line, indent succeeding lines appropriately
  1278.      */
  1279.  
  1280.     praline( n , a , 0 ) ;
  1281.     putchar('\n') ;
  1282. }
  1283.  
  1284. /*-----------------------------------------------------------------*/
  1285.  
  1286. praline( n , a , indent )
  1287. int    n ;
  1288. double    *a ;
  1289. int    indent ;
  1290. {
  1291.     /*
  1292.      * print out as many lines as required, indenting as we go
  1293.      */
  1294.  
  1295.     int    i ;
  1296.  
  1297.     if ( indent )
  1298.     {
  1299.         putchar('\n') ;
  1300.         for ( i=indent ; i ; --i )
  1301.             putchar(' ');
  1302.     }
  1303.  
  1304.     for ( i=INALINE ;  i && n  ;  --i , --n ) 
  1305.         printf("%11.4e " , *a++ ) ;
  1306.  
  1307.     if ( !n )
  1308.         return ;
  1309.  
  1310.     praline( n , a , ++indent ) ;
  1311. }
  1312.  
  1313. /*-----------------------------------------------------------------*/
  1314.  
  1315. mout ( n , a )
  1316. int    n ;
  1317. double    *a ;
  1318. {
  1319.     /* 
  1320.      * output the floating point matrix a with n by n components 
  1321.      */
  1322.  
  1323.     int    i ;
  1324.     double    *p ;
  1325.  
  1326.     for (i=n , p=a ; i ; --i , p += n )
  1327.         vout( n , p ) ;
  1328. }
  1329.  
  1330. /*-----------------------------------------------------------------*/
  1331.  
  1332. double dot( n , a , b )
  1333. int    n ;
  1334. double    *a , *b ;
  1335. {
  1336.     /*
  1337.      * dot product of the vectors a and b, n components each
  1338.      */
  1339.  
  1340.     double    c ;
  1341.  
  1342.     c = 0.0 ;
  1343.     while( n-- )
  1344.         c += *a++ * *b++ ;
  1345.     return( c ) ;
  1346. }
  1347.  
  1348. /*-----------------------------------------------------------------*/
  1349.  
  1350. reset ( n , y )
  1351. int    n ;
  1352. double    *y ;
  1353. {
  1354.     /*
  1355.      * set the n by n matrix y equal to the identity matrix
  1356.      */
  1357.  
  1358.     int    i, j;
  1359.  
  1360.     for( i = n-1, j = n ; i ;  --i, j = n )
  1361.     {
  1362.         *y++ = 1.0 ;
  1363.         while ( j-- )
  1364.             *y++ = 0.0 ;
  1365.     }
  1366.     *y = 1.0 ;
  1367. }
  1368.  
  1369. /*-----------------------------------------------------------------*/
  1370.  
  1371. cross( n , v1 , v2 , m )
  1372. int    n ;
  1373. double    *v1 ;
  1374. double    *v2 ; 
  1375. double    *m ;
  1376. {
  1377.     /* 
  1378.      * product of two vectors of dimension n yielding an n by n matrix
  1379.      */
  1380.  
  1381.     int     i , j ;
  1382.     double    *p ;
  1383.  
  1384.     for( i=n , p=v2 ; i ; --i , ++v1 , p=v2 )
  1385.         for( j=n ; j ; --j )
  1386.             *m++ = *v1 * *p++ ;
  1387. }
  1388.  
  1389. /*-----------------------------------------------------------------*/
  1390.  
  1391. matvec( n , m , v1 , v2 )
  1392. int    n ;
  1393. double    *m ;
  1394. double    *v1 ;
  1395. double    *v2 ;
  1396. {
  1397.     /*
  1398.      * product of an n by n matrix and a column vector with n components
  1399.      * yielding a second n component vector:  m * v1 = v2
  1400.      */
  1401.     int    i , j ;
  1402.     double    *p ;
  1403.  
  1404.     for ( i=n , p=v1 ; i ; --i , ++v2 , p=v1 ) 
  1405.         for ( j=n , *v2=0.0 ; j ; --j )
  1406.             *v2 += *m++ * *p++ ;
  1407. }
  1408.  
  1409. /* LISTING SEVEN
  1410.  * READ.C                         12/05/85
  1411.  *
  1412.  * reads standard input to get appropriate parameters
  1413.  * echoes the data as it is read
  1414.  * (c) Copyright 1985, Billybob Software.  All rights reserved.
  1415.  * This program may be reproduced for personal, non-profit use only.
  1416.  */
  1417.  
  1418. #include    "global.h"
  1419. #define        LEN        130
  1420. #define        MAXKEY        19    
  1421. #define        DEBUGON        if (*debug)
  1422. #define        LF        putchar('\n')
  1423.  
  1424. /************************************************************************/
  1425.  
  1426. reader ( fun , const , nparam , param , limit , dstep , pname ,
  1427.              epsilon , itermin , iterlim , nreset , g0 , debug , 
  1428.              data , npoints , npmax , method )
  1429.  
  1430. int    *fun ;        /* address of function (non-portable)        */
  1431. double    const[] ;    /* values of constants used in fcn        */
  1432. int    *nparam ;    /* number of parameters to be found by fitting    */
  1433. double    param[] ;    /* starting values of parameters        */
  1434. BOUND    limit[] ;    /* on-off flag, upper and lower limits          */
  1435. double    dstep[] ;    /* step sizes used to calculate derivatives    */
  1436. char    *pname[] ;    /* parameter names                */
  1437. double    epsilon[] ;    /* vector of stopping criteria for iteration    */
  1438. int    *itermin ;    /* minimum number of iterations required    */
  1439. int    *iterlim ;    /* maximum number of iterations allowed        */
  1440. int    *nreset ;    /* # of iters to reset y matrix         */
  1441. double    *g0;        /* expected value of minimum            */
  1442. int    *debug ;    /* debug mode                    */
  1443. DATA    data[] ;    /* data from input file                */
  1444. int    *npoints ;    /* number of data points             */
  1445. int    npmax ;        /* maximum number of data points allowed    */
  1446. int    *method ;    /* dfp or bfgs updating method            */
  1447. {
  1448.     int    nconst ;    /* number of constants used in fcn    */
  1449.     char    instring[LEN] ;
  1450.     char    *restofline ;
  1451.     char    *firstword ; 
  1452.     static  char    funname[] = "               " ;
  1453.     static    char    *key[MAXKEY] =
  1454.     {
  1455.         "bfgs" ,        /*       0     */
  1456.         "constants" ,        /*       1    */
  1457.         "debug" ,        /*       2    */
  1458.         "derivsteps" ,        /*       3    */
  1459.         "dfp" ,            /*       4    */
  1460.         "end" ,            /*       5    */
  1461.         "epsilons" ,        /*       6    */
  1462.         "exmin" ,        /*       7      */
  1463.         "funname" ,        /*       8    */    
  1464.         "iterlim" ,        /*       9    */    
  1465.         "itermin" ,        /*     10    */
  1466.         "limitflags" ,        /*     11    */
  1467.         "lowerlimits" ,        /*     12    */
  1468.         "newdata" ,        /*    13    */
  1469.         "next" ,        /*     14    */
  1470.         "pstart" ,        /*     15    */
  1471.         "reset" ,        /*     16    */
  1472.         "sd" ,            /*      17    */
  1473.         "upperlimits"        /*     18    */
  1474.     } ;
  1475.  
  1476.     /*
  1477.      * note above order; keywords are tested below in same order...
  1478.      * remember this when adding keywords! 
  1479.      */
  1480.  
  1481.     int    i , j , ncmd , ncom , bad ;
  1482.  
  1483. #ifdef    C80
  1484.     int chan ;
  1485. #else
  1486.     FILE    *chan ;
  1487.     FILE    *fopen() ;
  1488.     char    *fgets() ;
  1489. #endif
  1490.  
  1491.     /*
  1492.      * loop until a "next" or "end" command is picked up
  1493.      */
  1494.  
  1495.     ncmd = ncom = bad = 0 ;
  1496.     for (;;)
  1497.     {
  1498.  
  1499. #ifdef    C80
  1500.         if ( ! getline( instring , LEN ) )
  1501.             return( ALLDONE ) ;
  1502.  
  1503.         printf("%s\n", instring) ;
  1504.  
  1505.         /*
  1506.          * C-80 library routine getline returns the line stored into 
  1507.          * "instring" with the null at the end, but with the '\n'
  1508.          * stripped off.
  1509.          */
  1510.  
  1511. #else
  1512.         if (fgets( instring , LEN , stdin ) == NULL)
  1513.             return(ALLDONE);
  1514.  
  1515.         printf("%s", instring) ;
  1516. #endif
  1517.  
  1518.         restofline = instring ;
  1519.         getnext( &restofline , &firstword ) ;
  1520.  
  1521.         for ( j=0 ; j<MAXKEY ; ++j ) 
  1522.             if ( ! strcmp( firstword , key[j] ) )
  1523.                 break ;
  1524.  
  1525.         switch (j)
  1526.         {
  1527. /* bfgs        */
  1528.         case 0:
  1529.             ++ncmd ;
  1530.             DEBUGON printf("bfgs update requested\n") ;
  1531.             *method = BFGS ;
  1532.             break ;
  1533. /* constants    */
  1534.         case 1:
  1535.             ++ncmd ;
  1536.             DEBUGON printf("%d constants in %s:\n", nconst ,
  1537.                         funname );
  1538.             for (i=0 ; i < nconst ; ++i)
  1539.             {
  1540.                 getnext( &restofline , &firstword ) ;
  1541.                 if( check( firstword ) )
  1542.                 {
  1543.                     ++bad; 
  1544.                     break;
  1545.                 } 
  1546.                 const[i] = (double) atof( firstword ) ;
  1547.                 DEBUGON printf("%e  ", const[i]);
  1548.             }
  1549.             DEBUGON LF;
  1550.             break ;
  1551. /* debug    */
  1552.         case 2:
  1553.             ++ncmd ;    
  1554.             getnext( &restofline , &firstword ) ;
  1555.             if ( check( firstword ) )
  1556.             {
  1557.                 ++bad;
  1558.                 break;
  1559.             } 
  1560.  
  1561.             *debug = atoi( firstword ) ;
  1562.             DEBUGON printf("debug level is %d\n", *debug);
  1563.             break ;
  1564. /* derivsteps    */        
  1565.         case 3:
  1566.             ++ncmd ;
  1567.             DEBUGON printf("derivsteps are:\n");
  1568.             for (i=0 ; i < *nparam ; ++i)
  1569.             {
  1570.                 getnext( &restofline , &firstword ) ;    
  1571.                 if ( check( firstword ) )
  1572.                 {
  1573.                     ++bad;
  1574.                     break;
  1575.                 } 
  1576.                 dstep[i] = (double) atof( firstword ) ;
  1577.                 DEBUGON printf("%e  ", dstep[i]);
  1578.             }
  1579.             DEBUGON LF;
  1580.             break ;
  1581. /* dfp        */
  1582.         case 4:
  1583.             ++ncmd ;
  1584.             DEBUGON printf("dfp update requested\n") ;
  1585.             *method = DFP ;
  1586.             break ;
  1587. /* end        */        
  1588.         case 5:
  1589.             DEBUGON printf("end of all data\n") ;
  1590.             return( ALLDONE ) ;
  1591. /* epsilons    */
  1592.         case 6:
  1593.             ++ncmd ;
  1594.             DEBUGON printf("%d epsilons for cutoffs:\n", 
  1595.                                 NEPS );
  1596.             for (i=0 ; i < NEPS ; ++i)
  1597.             {
  1598.                 getnext( &restofline , &firstword ) ;
  1599.                 if( check( firstword ) )
  1600.                 {
  1601.                     ++bad;
  1602.                     break;
  1603.                 } 
  1604.                 epsilon[i] = (double) atof( firstword );
  1605.                 DEBUGON printf("%e  ", epsilon[i]);
  1606.             }
  1607.             DEBUGON LF;
  1608.             break ;
  1609. /* exmin    */
  1610.         case 7:
  1611.             ++ncmd ;
  1612.             getnext( &restofline , &firstword ) ;
  1613.             if ( check( firstword ) )
  1614.             {
  1615.                 ++bad;
  1616.                 break;
  1617.             } 
  1618.             *g0 = (double) atof( firstword ) ;
  1619.             DEBUGON printf("expected minimum is %11.4e\n", *g0);
  1620.             break ;
  1621. /* funname    */
  1622.         case 8:
  1623.             ++ncmd ;
  1624.             getnext( &restofline , &firstword ) ;
  1625.             if( check( firstword ) )
  1626.             {
  1627.                 ++bad;
  1628.                 break;
  1629.             } 
  1630.             strcpy( funname , firstword ) ;
  1631.             DEBUGON printf("function to minimize is %s\n",
  1632.                             funname ) ;
  1633.  
  1634.             if ( funlib( funname, fun, nparam, &nconst, pname))
  1635.             {
  1636.                 DEBUGON printf("function not in library!\n") ;
  1637.                 ++bad ;
  1638.             }
  1639.             break ;
  1640. /* iterlim    */
  1641.         case 9:
  1642.             ++ncmd ;
  1643.             getnext( &restofline , &firstword ) ; 
  1644.             if( check( firstword ))
  1645.             {
  1646.                 ++bad;
  1647.                 break;
  1648.             } 
  1649.             *iterlim = atoi( firstword ) ;
  1650.             DEBUGON printf("iterlim is %d\n", *iterlim);
  1651.             break ;
  1652. /* itermin    */
  1653.         case 10:
  1654.             ++ncmd ;
  1655.             getnext( &restofline , &firstword ) ;
  1656.             if( check( firstword ) )
  1657.             {
  1658.                 ++bad;
  1659.                 break;
  1660.             } 
  1661.             *itermin = atoi( firstword ) ;
  1662.             DEBUGON printf("itermin is %d\n", *itermin);
  1663.             break ;
  1664. /* limitflags    */
  1665.         case 11:
  1666.             ++ncmd ;
  1667.             DEBUGON printf("limitflags are:\n") ;
  1668.             for (i=0 ; i < *nparam ; ++i)
  1669.             {
  1670.                 getnext( &restofline , &firstword ) ;
  1671.                 if ( check( firstword ) )
  1672.                 {
  1673.                     ++bad;
  1674.                     break;
  1675.                 } 
  1676.                 limit[i].fl = atoi( firstword ) ;
  1677.                 DEBUGON printf("%d  ", limit[i].fl);
  1678.             }
  1679.             DEBUGON LF;
  1680.             break ;
  1681. /* lowerlimits    */
  1682.         case 12:
  1683.             ++ncmd ;
  1684.             DEBUGON printf("lowerlimits are:\n") ;
  1685.             for (i=0 ; i < *nparam ; ++i)
  1686.             {
  1687.                 getnext( &restofline , &firstword ) ;
  1688.                 if ( check( firstword ) )
  1689.                 {
  1690.                     ++bad;
  1691.                     break;
  1692.                 } 
  1693.                 limit[i].lo = (double) atof( firstword ) ;
  1694.                 DEBUGON printf("%e  ", limit[i].lo);
  1695.             }
  1696.             DEBUGON LF;
  1697.             break ;
  1698. /* newdata    */
  1699.         case 13:
  1700.             ++ncmd ;
  1701.             getnext( &restofline , &firstword ) ;
  1702.             DEBUGON printf("datafile requested was %s\n" ,
  1703.                             firstword ) ;
  1704.             if ( check( firstword ) )
  1705.             {
  1706.                 ++bad;
  1707.                 break;
  1708.             } 
  1709.  
  1710.             chan = fopen( firstword , "r");
  1711.  
  1712.             if ( !chan ) 
  1713.             {
  1714.                 printf("datafile can't be opened!\n") ;
  1715.                 return( ALLDONE ) ;
  1716.             }
  1717.  
  1718.             fscanf( chan , "%d" , npoints ) ;
  1719.             DEBUGON printf("%d datapoints in file\n", *npoints);
  1720.  
  1721.             if (*npoints > npmax)
  1722.             {
  1723.                 ++bad ;
  1724.                 printf("more data than allowed!\n") ;
  1725.                 printf("%d datapoints in file\n", *npoints);
  1726.                 printf("%d points allowed for\n", npmax ) ;
  1727.                 break ;
  1728.             }
  1729.  
  1730.             for (i=0 ; i<*npoints ; ++i)
  1731.             {
  1732. #ifdef    C80
  1733.                 fscanf(chan, "%f %f", &data[i].x, &data[i].y);
  1734. #else
  1735.                 fscanf(chan,"%lf %lf",&data[i].x, &data[i].y);
  1736. #endif
  1737.                 DEBUGON printf("%3d %f %f\n" , i+1 ,
  1738.                     data[i].x , data[i].y ) ;
  1739.             }
  1740.             fclose( chan ) ;
  1741.             break ;
  1742. /* next        */        
  1743.         case 14:
  1744.             ++ncmd ;
  1745.             DEBUGON printf("end of this data set\n");
  1746.             printf("%2d command lines read\n", ncmd) ;
  1747.             printf("%2d comment lines read\n", ncom) ;
  1748.  
  1749.             if( ! *nparam )
  1750.                 ++bad ; 
  1751.  
  1752.             return( bad ) ;
  1753. /* pstart    */        
  1754.         case 15:
  1755.             ++ncmd ;
  1756.             DEBUGON printf("starting values are:\n" ) ;
  1757.             for (i=0 ; i < *nparam ; ++i) 
  1758.             {
  1759.                 getnext( &restofline , &firstword ) ;
  1760.                 if ( check( firstword ) )
  1761.                 {
  1762.                     ++bad;
  1763.                     break;
  1764.                 } 
  1765.                 param[i] = (double) atof( firstword ) ;
  1766.                 DEBUGON printf("%e  ", param[i]);
  1767.             }
  1768.             DEBUGON LF;
  1769.             break ;
  1770. /* reset    */                
  1771.         case 16:
  1772.             ++ncmd ;
  1773.             getnext( &restofline , &firstword ) ;
  1774.             if( check( firstword ) )
  1775.             {
  1776.                 ++bad;
  1777.                 break;
  1778.             } 
  1779.             *nreset = atoi( firstword ) ;
  1780.             DEBUGON printf("reset is %d\n", *nreset);
  1781.             break ;    
  1782. /* sd        */
  1783.         case 17:
  1784.             ++ncmd ;
  1785.             DEBUGON printf("steepest descent update requested\n") ;
  1786.             *method = STDES ;
  1787.             break ;
  1788. /* upperlimits    */        
  1789.         case 18:
  1790.             ++ncmd ;
  1791.             DEBUGON printf("upperlimits are:\n");
  1792.             for (i=0 ; i < *nparam ; ++i) 
  1793.             {
  1794.                 getnext( &restofline , &firstword ) ;
  1795.                 if ( check( firstword ) )
  1796.                 {
  1797.                     ++bad;
  1798.                     break;
  1799.                 } 
  1800.                 limit[i].up = (double) atof( firstword ) ;
  1801.                 DEBUGON printf("%e  ", limit[i].up);
  1802.             }
  1803.             DEBUGON LF;
  1804.             break ;
  1805.  
  1806.             /* 
  1807.              * anything else is a comment; throw away and go
  1808.              * to next line.
  1809.              */
  1810.  
  1811.         default:
  1812.             ++ncom ;
  1813.             DEBUGON printf("***%s*** taken as comment\n", 
  1814.                             firstword ) ;
  1815.         }
  1816.     }
  1817. }
  1818.  
  1819. /*************************************************************************/
  1820.  
  1821. getnext( string , next )
  1822. char    **string ;            
  1823. char    **next ;            
  1824. {
  1825.     /*
  1826.      * splits the input string "string" into two pieces:
  1827.      *    "next" contains the first word with no leading or trailing blanks
  1828.      *    "string" then contains the rest of the line
  1829.      */
  1830.  
  1831.     int    length ;
  1832.     char    *p ;
  1833.  
  1834.     length = strlen( p = *string ) ;
  1835.     while( length-- > 0  &&  isspace( *p++ ) )
  1836.         ;
  1837.  
  1838.     *next = --p ;
  1839.  
  1840.     while( length-- >= 0 && !isspace( *++p ) )
  1841.         ;
  1842.  
  1843.     *p = '\0' ;
  1844.     *string = ++p ;
  1845. }
  1846. /*************************************************************************/
  1847.  
  1848. check( string )
  1849. char *string ;
  1850. {
  1851.     /*
  1852.      * checks if a string starts with a blank
  1853.      * if the string was obtained with getnext, this means the string 
  1854.      * is blank unconditionally report this error
  1855.      */    
  1856.  
  1857.     if ( ! *string )
  1858.     {
  1859.         printf("Unexpected blank encountered!\n") ;
  1860.         return( 1 ) ;
  1861.     }
  1862.     else
  1863.         return( 0 ) ;
  1864. }
  1865.  
  1866.  
  1867. LISTING EIGHT
  1868. Test input for the VMM program
  1869. File name is "a1.dat"
  1870. Sample input using the Cohen function
  1871. Demonstrates using the limit flags and detail debug printout
  1872. funname        cohen
  1873. pstart        1.0    1.0
  1874. limitflags    1    1
  1875. lowerlimits    0.5    0.75
  1876. upperlimits    2.0    3.0
  1877. iterlim        10
  1878. debug        2
  1879. next
  1880. This case tests a function which uses "experimental" data
  1881. Note the non-default epsilons for this case
  1882. funname        sine
  1883. pstart        1.57    -0.65    0.08    -0.005
  1884. epsilons    1.0e-14    1.0e-6
  1885. newdata        b:sine.dat
  1886. iterlim        12
  1887. reset        8
  1888. debug        1
  1889. next
  1890. The following cases use the Rosenbrock function
  1891. These correspond to the benchmarks...
  1892. Case 1
  1893. funname    rosen
  1894. pstart    -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0
  1895. constants  100.  100.  100.  100.  100.
  1896. sd
  1897. next
  1898. Case 2
  1899. funname    rosen
  1900. pstart    -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0
  1901. constants  10.  10.  10.  10.  10.
  1902. sd
  1903. next
  1904. Case 3
  1905. funname    rosen
  1906. pstart    -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0
  1907. constants  1.  1.  1.  1.  1.
  1908. sd
  1909. next
  1910. Case 4
  1911. funname    rosen
  1912. pstart  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8
  1913. constants  100.  100.  100.  100.  100.
  1914. sd
  1915. next
  1916. Case 5
  1917. funname    rosen
  1918. pstart  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8
  1919. constants  10.  10.  10.  10.  10.
  1920. sd
  1921. next
  1922. Case 6
  1923. funname    rosen
  1924. pstart  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8
  1925. constants  1.  1.  1.  1.  1.
  1926. sd
  1927. next
  1928. Case 7
  1929. funname    rosen
  1930. pstart    -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0
  1931. constants  100.  100.  100.  100.  100.
  1932. dfp
  1933. reset  5
  1934. next
  1935. Case 8
  1936. funname    rosen
  1937. pstart    -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0
  1938. constants  10.  10.  10.  10.  10.
  1939. dfp
  1940. reset  5
  1941. next
  1942. Case 9
  1943. funname    rosen
  1944. pstart    -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0
  1945. constants  1.  1.  1.  1.  1.
  1946. dfp
  1947. reset  5
  1948. next
  1949. Case 10
  1950. funname    rosen
  1951. pstart  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8
  1952. constants  100.  100.  100.  100.  100.
  1953. dfp
  1954. reset  5
  1955. next
  1956. Case 11
  1957. funname    rosen
  1958. pstart  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8
  1959. constants  10.  10.  10.  10.  10.
  1960. dfp
  1961. reset  5
  1962. next
  1963. Case 12
  1964. funname    rosen
  1965. pstart  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8
  1966. constants  1.  1.  1.  1.  1.
  1967. dfp
  1968. reset  5
  1969. next
  1970. Case 13
  1971. funname    rosen
  1972. pstart    -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0
  1973. constants  100.  100.  100.  100.  100.
  1974. bfgs
  1975. reset  5
  1976. next
  1977. Case 14
  1978. funname    rosen
  1979. pstart    -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0
  1980. constants  10.  10.  10.  10.  10.
  1981. bfgs
  1982. reset  5
  1983. next
  1984. Case 15
  1985. funname    rosen
  1986. pstart    -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0  -1.2  1.0
  1987. constants  1.  1.  1.  1.  1.
  1988. bfgs
  1989. reset  5
  1990. next
  1991. Case 16
  1992. funname    rosen
  1993. pstart  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8
  1994. constants  100.  100.  100.  100.  100.
  1995. bfgs
  1996. reset  5
  1997. next
  1998. Case 17
  1999. funname    rosen
  2000. pstart  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8
  2001. constants  10.  10.  10.  10.  10.
  2002. bfgs
  2003. reset  5
  2004. next
  2005. Case 18
  2006. funname    rosen
  2007. pstart  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8  1.2  0.8
  2008. constants  1.  1.  1.  1.  1.
  2009. bfgs
  2010. reset  5
  2011. next
  2012. end
  2013.  
  2014.  
  2015. /* LISTING NINE
  2016.  * The way it is done in vmm
  2017.  * This works and will be portable
  2018.  * so long as addresses are the 
  2019.  * same size as ints.
  2020.  * Fine on 16 bit systems with "small
  2021.  *    model" e.g. 64K programs
  2022.  */
  2023.  
  2024. #define        double        float    
  2025. #include    "fprintf.h"
  2026. #include    "scanf.h"
  2027. /*---------------------------------*/
  2028. main()
  2029. {
  2030. double    *(*fun)() ;
  2031. double    a , b ;
  2032.  
  2033.     while(1) {
  2034.         readit( &fun , &a , &b ) ;
  2035.         doit( fun , a , b ) ;
  2036.     }
  2037. }
  2038. /*----------------------------------*/
  2039. readit( fun , a , b )
  2040. int    *fun ;
  2041. double     *a ;
  2042. double     *b ;
  2043. {
  2044. char    string[10] ;
  2045.  
  2046.     scanf("%s %f %f", string , a , b ) ;
  2047.     
  2048.     funlib( fun , string ) ;
  2049. }
  2050. /*------------------------------------*/
  2051. funlib( fun , string )
  2052. int    *fun ;
  2053. char    string[] ;
  2054. {
  2055. double    sum() ;
  2056. double  mul() ;
  2057.     
  2058.     if ( !strcmp( string, "add" ) ) *fun = sum ;
  2059.     else if ( !strcmp( string, "mul" ) ) *fun = mul ;
  2060.     else exit(0);
  2061. }
  2062. /*---------------------------------------*/
  2063. doit( fun , a , b )
  2064. double    (*fun)() ;
  2065. double    a ;
  2066. double    b ;
  2067. {
  2068. double    g ;
  2069.     g = (*fun)( a , b ) ;
  2070.     printf("%f is the result\n" , g ) ;
  2071. }
  2072. /*----------------------------------------*/
  2073. double    sum( a , b ) 
  2074. double     a ;
  2075. double    b ;
  2076. {
  2077.     printf("%f + %f = %f\n" , a , b , a+b ) ;
  2078.     return( a + b ) ;
  2079. }
  2080. /*----------------------------------------*/
  2081. double    mul( a , b )
  2082. double     a ;
  2083. double    b ;
  2084. {
  2085.     printf("%f * %f = %f\n" , a , b , a*b ) ;
  2086.     return( a * b ) ;
  2087. }
  2088.  
  2089.  
  2090. /* LISTING TEN
  2091.  * The correct and most portable way to pass an object which is a
  2092.  * pointer to a function
  2093.  */
  2094.     
  2095. #include    <stdio.h>
  2096. /*---------------------------------*/
  2097. main()
  2098. {
  2099.     double    (*fun)() ;        /* fun is a pointer to a function
  2100.                      *  which returns a double
  2101.                      */
  2102.     double    (*readit())() ;        /* readit is a function that
  2103.                      * returns a pointer to a function
  2104.                      * that returns a double
  2105.                      */
  2106.     double    a , b ;
  2107.  
  2108.     while(1) {
  2109.         fun = readit( &a , &b ) ;
  2110.         doit( fun , a , b ) ;
  2111.     }
  2112. }
  2113.  
  2114. /*----------------------------------*/
  2115.  
  2116. double    (*readit( a , b ))()    /* note position of arguments    */
  2117. double *a ;
  2118. double *b ;
  2119. {
  2120.     double    (*fun)() ;
  2121.     double    (*funlib())() ;        /*  funlib is a function which 
  2122.                      *  returns a pointer to a function
  2123.                      *  which returns a double
  2124.                      */
  2125.     char    string[10] ;
  2126.  
  2127.     scanf("%s %lf %lf", string , a , b ) ;
  2128.     
  2129.     fun = funlib( string ) ;
  2130.     return( fun ) ;
  2131. }
  2132.  
  2133. /*------------------------------------*/
  2134.  
  2135. double    (*funlib( string ))()    /* note position of arguments  */
  2136. char    string[] ;
  2137. {
  2138.     double    (*fun)() ;
  2139.     double    sum() ;
  2140.     double  mul() ;
  2141.         
  2142.     if ( !strcmp( string, "add" ) )
  2143.         fun = sum ;
  2144.     else if ( !strcmp( string, "mul" ) )
  2145.         fun = mul ;
  2146.     else
  2147.         exit(0);
  2148.  
  2149.     return( fun ) ;
  2150. }
  2151.  
  2152. /*---------------------------------------*/
  2153. doit( fun , a , b )
  2154. double    (*fun)() ;
  2155. double    a ;
  2156. double    b ;
  2157. {
  2158.     double    g ;
  2159.     g = (*fun)( a , b ) ;
  2160.     printf("%f is the result\n" , g ) ;
  2161. }
  2162. /*----------------------------------------*/
  2163. double    sum( a , b ) 
  2164. double     a ;
  2165. double    b ;
  2166. {
  2167.     printf("%f + %f = %f\n" , a , b , a+b ) ;
  2168.     return( a + b ) ;
  2169. }
  2170. /*----------------------------------------*/
  2171. double    mul( a , b )
  2172. double     a ;
  2173. double    b ;
  2174. {
  2175.     printf("%f * %f = %f\n" , a , b , a*b ) ;
  2176.     return( a * b ) ;
  2177. }
  2178.  
  2179. *END OF LISTINGS
  2180.