From oleg@nrcbsa.bio.nrc.ca Thu May 23 08:52:24 1991 Return-Path: Received: from nrcnet0.nrc.ca by CS.UTK.EDU with SMTP (5.61++/2.5.1s-UTK) id AA26046; Thu, 23 May 91 08:52:12 -0400 Message-Id: <9105231251.AA28692@ nrcbsa.bio.nrc.ca> Date: Thu, 23 May 91 08:51:36 EDT From: Dr. Oleg Keselyov To: dongarra@cs.utk.edu Subject: HL_vector.shar Status: RO #--------------------------------CUT HERE------------------------------------- #! /bin/sh # # This is a shell archive. Save this into a file, edit it # and delete all lines above this comment. Then give this # file to sh by executing the command 'sh file'. The files # will be extracted into the current directory owned by # you with default permissions. # The archive contains a package of high level vector operations. Involves the # Aitken-Lagrange interpolation over the table of uniform or # arbitrary mesh, Hook-Jeevse local multidimensional minimizer. # The contents of the archive # # Source code for the programs being submitted # ali.c Aitken-Lagrange interpolation over the table of # uniform or arbitrary mesh # hjmin.c Hook-Jeevse local minimum finder # # Validation routines # valin.c Verify alin() function (Arbitrary mesh interpolation) # valis.c Verify alis() function (Uniform mesh interpolation) # vhjmin.c Verify Hook-Jeevse minimizer # # Verification protocols # valin.dat On Silicon Graphics IRIS 4D/280S # valis.dat # vhjmin.dat # valin.dat.vax On VAX 6000 station # valis.dat.vax # vhjmin.dat.vax echo 'x - ali.c' sed 's/^X//' << '________This_Is_The_END________' > ali.c X/* X ************************************************************************ X * C math library X * Aitken-Lagrange interpolation X * X * function ALIS - interpolate function value for a given argument value X * using function values tabulated over the uniform grid X * X * Input X * double alis(q,x0,s,y) X * double q Argument value specified X * double x0 Argument for the 1. grid X * double s Grid mesh, >0 X * VECTOR y Vector of function values X * tabulated at points X * x0 + s*(i-V_lwb(v)) X * Vector has to contain at X * least 2 elements X * X * function ALIN - interpolate function value for a given argument value X * using function values tabulated over the non-uniform grid X * X * Input X * double alin(q,x,y) X * const double q Argument value specified X * const VECTOR x Vector of argument values X * const VECTOR y Vector of function values X * tabulated at points x[i] X * Vector has to contain at X * least 2 elements X * Output X * Both functions return the interpolated function value at point q X * Interpolation is terminated either X * - if the difference between two successive interpolated X * values is absolutely less than EPSILON X * - if the absolute value of this difference stops X * diminishing X * - after (N-1) steps, N being the no. of elements in vector y X * X * Algorithm X * Aitken scheme of Lagrange interpolation X * Algorithm is described in a book X * "Mathematical software for computers", Institute of X * Mathematics of the Belorussian Academy of Sciences, X * Minsk, 1974, issue #4 X * p. 146 (description of ALI, DALI subroutines) X * p. 180 (description of ATSE, DATSE subroutines) X * X ************************************************************************ X */ X X X#include "assert.h" X#include "vector.h" X#include "math.h" X#include X X X/* X * Package private data X */ X Xstatic VECTOR Arg; /* [1:n] Arranged table of arguments */ Xstatic VECTOR Val; /* [1:n] Arranged table of function values*/ Xstatic double Q; /* Argument the function is to be */ X /* interpolated at */ X X Xstatic void allocate(n) Xconst int n; /* No. of points arrays must be allocated at*/ X{ X if( Arg != (VECTOR)0 ) /* Vectors were allocated at previous calls*/ X if( V_no_elems(Arg) == n ) X return; /* Vectors were allocated for a proper no. */ X else X { X V_free(Arg); X V_free(Val); X } X Arg = V_new(1,n); X Val = V_new(1,n); X} X X X/* X * Aitken - Lagrange process X * X * Arg and Val tables are assumed to be arranged in the proper way X * (see root module) X * X */ X Xstatic double aitken() X{ X register double valp = *V_elem(Val,1); /* Best approximation found*/ X register double diffp = HUGE_VAL; /* abs(valp-prev. to valp)*/ X register int i,j; X X/* V_print(Arg,"Arg - interpolation nodes"); V_print(Val,"Val"); */ X /* Compute j-th row of the Aitken scheme and */ X /* place it in the Val array */ X for(j=2; j<=V_upb(Val); j++) X { X register double argj = *V_elem(Arg,j); X register V$elem_type *valj = V_elem(Val,j); X register double diff; X for(i=1; i<=j-1; i++) X *valj = ( *V_elem(Val,i)*(Q-argj) - *valj * (Q-*V_elem(Arg,i)) )/ X (*V_elem(Arg,i) - argj); X X/* printf("\nValj = %g, valp = %g, argj=%g",*valj,valp,argj); */ X diff = fabs( *valj - valp ); X X if( j>2 && diff == 0 ) /* Exact result has been achieved*/ X break; X X if( j>4 && diff > diffp ) /* Difference stops diminishing */ X break; /* after the 4. step */ X X valp = *valj; diffp = diff; X } X X return valp; X} X X X X/* X *======================================================================= X * Function ALIS root module X * X * Arranges input data in the form suitable for aitken and call the latter X * X * Algorithm X * Fill in vectors arg and val[i] = y(arg[i]) so that distance abs(x-arg[i]) X * would increase as i increases. X * X */ X Xdouble alis(q,x0,s,y) Xdouble q; /* Argument value specified */ Xdouble x0; /* Argument for the 1. grid */ Xdouble s; /* Grid mesh */ XVECTOR y; /* Vector of function values */ X{ X int n = V_no_elems(y); /* No. of grids */ X X assure( n > 1, "too few function values, DIM(y) must be > 1"); X assure( s > 0, "grid mesh has to be a positive number"); X allocate(n); X Q = q; X X /* Preparation step */ X /* Selection is done by computing the subscript */ X /* js of that argument, which is next to x. */ X /* Afterwards neighbouring argument values are */ X /* tested and selected in the above sense. */ X X { X int js = (int)( (q-x0)/s + 1.5 ); /* Index for a point */ X /* x0+s*i closest to q */ X register int dir; /* Direction to the next grid */ X /* close to q */ X register int jcurr, jleft, jright; X register int i; X X if( js<1 ) /* Test for the case of extrapo-*/ X js = 1; /* lation to the left end */ X else if( js > n ) X js = n; /* or to the right end */ X X dir = ( q > x0 + (js-1)*s ? 1 : -1 ); X X jcurr = js; jleft = js; jright = js; X for(i=1; i<=n; ++i) /* Pick up elements x0+s*i */ X { /* in the neighbourhood of x */ X *V_elem(Arg,i) = x0 + (jcurr-1)*s; X *V_elem(Val,i) = *V_elem(y,jcurr-1+V_lwb(y)); X if( jright >= n ) X dir = -1; X if( jleft <= 1 ) X dir = 1; X if( dir > 0 ) X { X jcurr = ++jright; X dir = -1; X } X else X { X jcurr = --jleft; X dir = 1; X } X } X } X X return aitken(); X} X X X X X/* X * Function that defines the order while sorting 'indices' array X * with qsort X * Function gets pointers to two elements in indices array, i.e. X * to the two indices for the Arg[i] and Arg[j] elements. X * Function returns -1, 0, or +1 depending on the fact if X * abs(Arg[i]-Q) is less, equal, or greater than abs(Arg[j]-Q) X */ X Xstatic int compare_indices(ip,jp) Xconst int * ip; Xconst int * jp; X{ X register double delta = fabs(*V_elem(Arg,*ip) - Q) - X fabs(*V_elem(Arg,*jp) - Q); X if( delta < 0 ) X return -1; X else if(delta == 0) X return 0; X else X return 1; X} X X X/* X *======================================================================= X * Function ALIN root module X * X * Arranges input data in the form suitable for aitken and call the latter X * X * Algorithm X * Fill in vectors arg and val[i] = y(arg[i]) so that distance abs(q-arg[i]) X * would increase as i increases. X * X */ X Xdouble alin(q,x,y) Xconst double q; /* Argument value specified */ Xconst VECTOR x; /* Vector of argument values */ Xconst VECTOR y; /* Vector of function values */ X{ X int n = V_no_elems(x); /* No. of grids */ X X V_are_compatible(x,y); X assure( n > 1, "too few function values, DIM(x) must be > 1"); X allocate(n); X Q = q; X /* Preparation step */ X /* Selection is done by sorting x,y arrays */ X /* in the way mentioned above. Fisrt an array */ X /* of indices is created and sorted, then arg, */ X /* val arrays are filled in with indices sorted */ X X { X int * indices = calloc(V_no_elems(Arg)+V_lwb(Arg),sizeof(int)); X register int i; X X for(i=V_lwb(Arg); i<=V_upb(Arg); i++)/* 1..n correspond to the */ X { /* unsorted x,y */ X indices[i] = i; X *V_elem(Arg,i) = *V_elem(x,V_lwb(x)+i-V_lwb(Arg)); X } X X /* Sort indices in the way that */ X /* |q-x[ind[i]]| < |q-x[ind[j]]|*/ X /* for all j>i */ X qsort(&indices[V_lwb(Arg)],n,sizeof(int),compare_indices); X X for(i=V_lwb(Arg); i<=V_upb(Arg); i++) X { X register int ind = indices[i]; X *V_elem(Arg,i) = *V_elem(x,V_lwb(x)+ind-V_lwb(Arg)); X *V_elem(Val,i) = *V_elem(y,V_lwb(y)+ind-V_lwb(Arg)); X } X X free(indices); X } X X return aitken(); X} ________This_Is_The_END________ if test `wc -l < ali.c` -ne 281; then echo 'shar: ali.c was damaged during transit (should have had 281 lines)' fi echo 'x - hjmin.c' sed 's/^X//' << '________This_Is_The_END________' > hjmin.c X/* X ************************************************************************ X * C math library X * function HJMIN - find local minimum of a given function X * by the Hook-Jeevse method X * X * Input X * double hjmin(b,h0,funct) X * VECTOR b On input contains the initial X * guess to the minimum location X * On return has the vector X * specifying for the minimum location X * const VECTOR h0 Initial values for the steps along X * each directions X * double f(const VECTOR x) Procedure to compute a function X * value at the specified point X * X * Output X * Hjmin returns the function value at the point of minimum X * X * Algorithm X * Hook-Jeevse method of direct search for a function minimum X * The method is of the 0. order (i.e. requiring no gradient computation) X * See X * B.Bondi. Methods of optimization. An Introduction - M., X * "Radio i sviaz", 1988 - 127 p. (in Russian) X * X * X * Notes X * static VECTORs h and bas are used as work arrays. They are not X * destroyed (freed) on exit so that next call to Hjmin could use X * them if they are still appropriate for that call. X * X ************************************************************************ X */ X X X#include "assert.h" X#include "math.h" X#include "vector.h" X X/* X * Static data X */ X Xstatic VECTOR h; /* Steps along axes */ Xstatic VECTOR bas; /* Base point */ Xstatic double f_bas; /* Function value at it */ Xstatic double (*f)(const VECTOR x); /* Function being minimized */ X X X/* X * Examine the function f in the vicinity of the bas point X * by making tentative steps fro/back along the every coordinate. X * When the function is found to decrease, the point bas is correspondingly X * updated. Examination() returns the minimal function value found in X * the investigated region. X * X */ Xstatic double examination() X{ X register V_pointer(ph,h); X register V_pointer(pbas,bas); X register int i; X X /* Perform step along a coordinate */ X for(i=V_lwb(bas); i<=V_upb(bas); i++,ph++,pbas++) X { X register double basi_old = *pbas; /* Old coordinate value */ X register double f_new; X X if( *pbas = basi_old + *ph, f_new = (*f)(bas), f_new < f_bas ) X f_bas = f_new; /* Step caused f to decrease, OK*/ X else if( *pbas = basi_old - *ph, f_new = (*f)(bas), f_new < f_bas ) X f_bas = f_new; X else X *pbas = basi_old; /* No fall was found along this coord*/ X } X return f_bas; X} X X X X/* X * Root module X */ X Xdouble hjmin(b,h0,ff) XVECTOR b; Xconst VECTOR h0; Xdouble (*ff)(const VECTOR x); X{ X double f_min; /* Min function value found */ X /* beforehand */ X X const double tau = 10*EPSILON; /* Termination criterium */ X const double threshold = 1e-8; /* Threshold for the function */ X /* decay to be treated as */ X /* significant */ X const double step_reduce_factor = 10; X register int i; X X f = ff; X V_are_compatible(b,h0); X X /* Allocate h and bas if needed */ X if( h == (VECTOR)0 || V_lwb(h) != V_lwb(h0) || V_upb(h) != V_upb(h0) ) X { /* h, bas nust be (re)allocated */ X if( h != (VECTOR)0 ) X { X V_free(h); X V_free(bas); X } X h = V_new(V_lwb(b),V_upb(b)); X bas = V_new(V_lwb(b),V_upb(b)); X } X V_assign(h,h0); X V_assign(bas,b); f_bas = f_min = (*f)(b); X X X for(;;) /* Main iteration loop */ X { /* b is a next approximation to min */ X if( examination() < f_min - threshold ) X { /* Function falled markedly */ X do /* Perform search on pattern */ X { /* bas - pattern, b - new min appr*/ X register V_pointer(pbas,bas); X register V_pointer(pmin,b); X for(i=V_lwb(bas); i<=V_upb(bas); i++) X { X register double t = (*pbas - *pmin) + *pbas; X *pmin++ = *pbas; X *pbas++ = t; X } X f_min = f_bas; f_bas = (*f)(bas); X } /* Continue search until f doesn't*/ X while( examination() < f_min - threshold ); /* decrease */ X V_assign(bas,b); f_bas = f_min; /* Save the best approx found */ X } X else /* Function didn't fall markedly*/ X { /* upon examination near bas */ X register int success = 1; /* Try to reduce the step */ X register V_pointer(ph,h); X register V_pointer(pbas,bas); X for(i=V_lwb(h); i<=V_upb(h); i++) X { X *ph /= step_reduce_factor; X success &= ( *ph++ /(1 + fabs(*pbas++)) < tau ); X } X if( success ) X return f_min; X } X X } X} ________This_Is_The_END________ if test `wc -l < hjmin.c` -ne 155; then echo 'shar: hjmin.c was damaged during transit (should have had 155 lines)' fi echo 'x - valin.c' sed 's/^X//' << '________This_Is_The_END________' > valin.c X/* X * Verify ALIN functions X */ X X#include "vector.h" X#include "math.h" X Xstatic double f(x) /* Function involved in the test */ Xdouble x; X{ X return sin(x) * exp(-x/10); X} X Xstatic VECTOR Nodes; X Xstatic void test(msg,q,ia,ib) /* Test interpolation to q with */ Xconst char * msg; Xdouble q; /* x[ia,ib], y[ia:ib] */ Xint ia,ib; X{ X VECTOR x = V_new(ia,ib); X VECTOR y = V_new(ia,ib); X register int i; X X for(i=V_lwb(y); i<=V_upb(y); i++) X *V_elem(y,i) = f( (*V_elem(x,i) = *V_elem(Nodes,i)) ); X X printf("\n\n%s\nNodes within the range [%d:%d] are used in interpolation\n", X msg,ia,ib); X X printf("\nInterpolated value to x= %g is\t%g",q,alin(q,x,y)); X printf("\nExact function value at this point\t%g\n",f(q)); X X V_free(x); X V_free(y); X} X Xmain() X{ X printf("\n\n\t\tVerify ALIN procedure for a function sin(x) * exp(-x/10)\n"); X X { X register int i; X Nodes = V_new(2,51); X for(i=V_lwb(Nodes); i<=11; i++) X *V_elem(Nodes,i) = (i-1.)/10; X for(i=12; i<=V_upb(Nodes); i++) X *V_elem(Nodes,i) = (i-6.)/5; X X printf("\t\tGrids 0.1 .. 1 have the mesh 0.1, and 1..9 have the mesh 0.2"); X } X X test("Extrapolation to the left",0.0,2,11); X test("Extrapolation to the right",8.0,2,41); X test("Interpolation just to the 16. node",2.0,2,31); X test("Interpolation to middle of table",1.55,2,50); X X} ________This_Is_The_END________ if test `wc -l < valin.c` -ne 58; then echo 'shar: valin.c was damaged during transit (should have had 58 lines)' fi echo 'x - valis.c' sed 's/^X//' << '________This_Is_The_END________' > valis.c X/* X * Verify ALIS functions X */ X X#include "vector.h" X#include "math.h" X Xstatic double s = 0.1; X Xstatic double f(x) /* Function involved in the test */ Xdouble x; X{ X return exp(-x); X} X Xstatic void test(msg,x,ia,ib) /* Test interpolation to x with */ Xconst char * msg; Xdouble x; /* y[ia:ib] */ Xint ia,ib; X{ X VECTOR y = V_new(ia,ib); X register int i; X X for(i=V_lwb(y); i<=V_upb(y); i++) X *V_elem(y,i) = f(s*i); X X printf("\n\n%s\nThe function being tabulated over [%g:%g]\n", X msg,ia*s,ib*s); X X printf("\nInterpolated value to x= %g is\t%g",x,alis(x,ia*s,s,y)); X printf("\nExact function value at this point\t%g\n",f(x)); X} X Xmain() X{ X printf("\n\n\t\tVerify ALIS procedure for a function exp(-x)\n"); X printf("\t\t\tTabulated unifomly with the grid mesh of %g\n\n",s); X X test("Extrapolation to the left",0.0,1,10); X test("Extrapolation to the right",30.5*s,0,28); X test("Interpolation just to the 10. node",10*s,5,15); X test("Interpolation to middle of table",5.6*s,1,10); X X} ________This_Is_The_END________ if test `wc -l < valis.c` -ne 44; then echo 'shar: valis.c was damaged during transit (should have had 44 lines)' fi echo 'x - valin.dat' sed 's/^X//' << '________This_Is_The_END________' > valin.dat X X X Verify ALIN procedure for a function sin(x) * exp(-x/10) X Grids 0.1 .. 1 have the mesh 0.1, and 1..9 have the mesh 0.2 X XExtrapolation to the left XNodes within the range [2:11] are used in interpolation X XInterpolated value to x= 0 is 8.51708e-07 XExact function value at this point 0 X X XExtrapolation to the right XNodes within the range [2:41] are used in interpolation X XInterpolated value to x= 8 is 0.443167 XExact function value at this point 0.444547 X X XInterpolation just to the 16. node XNodes within the range [2:31] are used in interpolation X XInterpolated value to x= 2 is 0.74447 XExact function value at this point 0.74447 X X XInterpolation to middle of table XNodes within the range [2:50] are used in interpolation X XInterpolated value to x= 1.55 is 0.85623 XExact function value at this point 0.85623 ________This_Is_The_END________ if test `wc -l < valin.dat` -ne 31; then echo 'shar: valin.dat was damaged during transit (should have had 31 lines)' fi echo 'x - valis.dat' sed 's/^X//' << '________This_Is_The_END________' > valis.dat X X X Verify ALIS procedure for a function exp(-x) X Tabulated unifomly with the grid mesh of 0.1 X X X XExtrapolation to the left XThe function being tabulated over [0.1:1] X XInterpolated value to x= 0 is 0.999998 XExact function value at this point 1 X X XExtrapolation to the right XThe function being tabulated over [0:2.8] X XInterpolated value to x= 3.05 is 0.0473609 XExact function value at this point 0.0473589 X X XInterpolation just to the 10. node XThe function being tabulated over [0.5:1.5] X XInterpolated value to x= 1 is 0.367879 XExact function value at this point 0.367879 X X XInterpolation to middle of table XThe function being tabulated over [0.1:1] X XInterpolated value to x= 0.56 is 0.571209 XExact function value at this point 0.571209 ________This_Is_The_END________ if test `wc -l < valis.dat` -ne 33; then echo 'shar: valis.dat was damaged during transit (should have had 33 lines)' fi echo 'x - valin.dat.vax' sed 's/^X//' << '________This_Is_The_END________' > valin.dat.vax X X X Verify ALIN procedure for a function sin(x) * exp(-x/10) X Grids 0.1 .. 1 have the mesh 0.1, and 1..9 have the mesh 0.2 X XExtrapolation to the left XNodes within the range [2:11] are used in interpolation X X Xvector [1:10] - Arg - interpolation nodes X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X X X Xvector [1:10] - Val X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X9.884006e-02 0.194735 0.286786 0.374149 0.456044 0.53176 0.600665 0.662203 0.715907 0.761394 X X XValj = 2.944693e-03, valp = 9.884006e-02, argj=0.2 XValj = -8.998138e-04, valp = 2.944693e-03, argj=0.3 XValj = -5.620640e-05, valp = -8.998138e-04, argj=0.4 XValj = 7.443349e-06, valp = -5.620640e-05, argj=0.5 XValj = 1.169593e-06, valp = 7.443349e-06, argj=0.6 XValj = 8.517077e-07, valp = 1.169593e-06, argj=0.7 XValj = 1.696237e-06, valp = 8.517077e-07, argj=0.8 XInterpolated value to x= 0 is 8.517077e-07 XExact function value at this point 0 X X XExtrapolation to the right XNodes within the range [2:41] are used in interpolation X X Xvector [1:40] - Arg - interpolation nodes X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X 7 6.8 6.6 6.4 6.2 6 5.8 5.6 5.4 5.2 X X 11 12 13 14 15 16 17 18 19 20 X------------------------------------------------------------------------------- X 5 4.8 4.6 4.4 4.2 4 3.8 3.6 3.4 3.2 X X 21 22 23 24 25 26 27 28 29 30 X------------------------------------------------------------------------------- X 3 2.8 2.6 2.4 2.2 2 1.8 1.6 1.4 1.2 X X 31 32 33 34 35 36 37 38 39 40 X------------------------------------------------------------------------------- X 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 X X X Xvector [1:40] - Val X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X0.32625 0.250326 0.161021 6.145556e-02 -4.469759e-02 -0.153346 -0.26013 -0.360585 -0.450327 -0.525232 X X 11 12 13 14 15 16 17 18 19 20 X------------------------------------------------------------------------------- X-0.581617 -0.61641 -0.627301 -0.612866 -0.572666 -0.5073 -0.418426 -0.308736 -0.181887 -4.238836e-02 X X 21 22 23 24 25 26 27 28 29 30 X------------------------------------------------------------------------------- X0.104544 0.253179 0.397478 0.531338 0.648834 0.74447 0.813426 0.85178 0.856709 0.826645 X X 31 32 33 34 35 36 37 38 39 40 X------------------------------------------------------------------------------- X0.761394 0.715907 0.662203 0.600665 0.53176 0.456044 0.374149 0.286786 0.194735 9.884006e-02 X X XValj = 0.705868, valp = 0.32625, argj=6.8 XValj = 0.505138, valp = 0.705868, argj=6.6 XValj = 0.395847, valp = 0.505138, argj=6.4 XValj = 0.434259, valp = 0.395847, argj=6.2 XValj = 0.450398, valp = 0.434259, argj=6 XValj = 0.445535, valp = 0.450398, argj=5.8 XValj = 0.443449, valp = 0.445535, argj=5.6 XValj = 0.443167, valp = 0.443449, argj=5.4 XValj = 0.441998, valp = 0.443167, argj=5.2 XInterpolated value to x= 8 is 0.443167 XExact function value at this point 0.444547 X X XInterpolation just to the 16. node XNodes within the range [2:31] are used in interpolation X X Xvector [1:30] - Arg - interpolation nodes X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X 2 2.2 1.8 1.6 2.4 2.6 1.4 2.8 1.2 1 X X 11 12 13 14 15 16 17 18 19 20 X------------------------------------------------------------------------------- X 3 0.9 0.8 3.2 0.7 0.6 3.4 0.5 3.6 0.4 X X 21 22 23 24 25 26 27 28 29 30 X------------------------------------------------------------------------------- X 0.3 3.8 0.2 0.1 4 4.2 4.4 4.6 4.8 5 X X X Xvector [1:30] - Val X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X0.74447 0.648834 0.813426 0.85178 0.531338 0.397478 0.856709 0.253179 0.826645 0.761394 X X 11 12 13 14 15 16 17 18 19 20 X------------------------------------------------------------------------------- X0.104544 0.715907 0.662203 -4.238836e-02 0.600665 0.53176 -0.181887 0.456044 -0.308736 0.374149 X X 21 22 23 24 25 26 27 28 29 30 X------------------------------------------------------------------------------- X0.286786 -0.418426 0.194735 9.884006e-02 -0.5073 -0.572666 -0.612866 -0.627301 -0.61641 -0.581617 X X XValj = 0.74447, valp = 0.74447, argj=2.2 XValj = 0.74447, valp = 0.74447, argj=1.8 XInterpolated value to x= 2 is 0.74447 XExact function value at this point 0.74447 X X XInterpolation to middle of table XNodes within the range [2:50] are used in interpolation X X Xvector [1:49] - Arg - interpolation nodes X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X 1.6 1.4 1.8 1.2 2 1 0.9 2.2 0.8 0.7 X X 11 12 13 14 15 16 17 18 19 20 X------------------------------------------------------------------------------- X 2.4 0.6 2.6 0.5 0.4 2.8 0.3 0.2 0.1 3 X X 21 22 23 24 25 26 27 28 29 30 X------------------------------------------------------------------------------- X 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 X X 31 32 33 34 35 36 37 38 39 40 X------------------------------------------------------------------------------- X 5.2 5.4 5.6 5.8 6 6.2 6.4 6.6 6.8 7 X X 41 42 43 44 45 46 47 48 49 X------------------------------------------------------------------------------- X 7.2 7.4 7.6 7.8 8 8.2 8.4 8.6 8.8 X X X Xvector [1:49] - Val X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X0.85178 0.856709 0.813426 0.826645 0.74447 0.761394 0.715907 0.648834 0.662203 0.600665 X X 11 12 13 14 15 16 17 18 19 20 X------------------------------------------------------------------------------- X0.531338 0.53176 0.397478 0.456044 0.374149 0.253179 0.286786 0.194735 9.884006e-02 0.104544 X X 21 22 23 24 25 26 27 28 29 30 X------------------------------------------------------------------------------- X-4.238836e-02 -0.181887 -0.308736 -0.418426 -0.5073 -0.572666 -0.612866 -0.627301 -0.61641 -0.581617 X X 31 32 33 34 35 36 37 38 39 40 X------------------------------------------------------------------------------- X-0.525232 -0.450327 -0.360585 -0.26013 -0.153346 -4.469759e-02 6.145556e-02 0.161021 0.250326 0.32625 X X 41 42 43 44 45 46 47 48 49 X------------------------------------------------------------------------------- X0.38632 0.428786 0.452664 0.457738 0.444547 0.414328 0.368939 0.310769 0.242614 X X XValj = 0.853013, valp = 0.85178, argj=1.4 XValj = 0.856146, valp = 0.853013, argj=1.8 XValj = 0.856207, valp = 0.856146, argj=1.2 XValj = 0.856229, valp = 0.856207, argj=2 XValj = 0.85623, valp = 0.856229, argj=1 XValj = 0.85623, valp = 0.85623, argj=0.9 XInterpolated value to x= 1.55 is 0.85623 XExact function value at this point 0.85623 ________This_Is_The_END________ if test `wc -l < valin.dat.vax` -ne 191; then echo 'shar: valin.dat.vax was damaged during transit (should have had 191 lines)' fi echo 'x - valis.dat.vax' sed 's/^X//' << '________This_Is_The_END________' > valis.dat.vax X X X Verify ALIS procedure for a function exp(-x) X Tabulated unifomly with the grid mesh of 0.1 X X X XExtrapolation to the left XThe function being tabulated over [0.1:1] X X Xvector [1:10] - Arg - interpolation nodes X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X X X Xvector [1:10] - Val X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X0.904837 0.818731 0.740818 0.67032 0.606531 0.548812 0.496585 0.449329 0.40657 0.367879 X X XValj = 0.990944, valp = 0.904837, argj=0.2 XValj = 0.999138, valp = 0.990944, argj=0.3 XValj = 0.999918, valp = 0.999138, argj=0.4 XValj = 0.999992, valp = 0.999918, argj=0.5 XValj = 0.999998, valp = 0.999992, argj=0.6 XValj = 0.999998, valp = 0.999998, argj=0.7 XInterpolated value to x= 0 is 0.999998 XExact function value at this point 1 X X XExtrapolation to the right XThe function being tabulated over [0:2.8] X X Xvector [1:29] - Arg - interpolation nodes X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X 2.8 2.7 2.6 2.5 2.4 2.3 2.2 2.1 2 1.9 X X 11 12 13 14 15 16 17 18 19 20 X------------------------------------------------------------------------------- X 1.8 1.7 1.6 1.5 1.4 1.3 1.2 1.1 1 0.9 X X 21 22 23 24 25 26 27 28 29 X------------------------------------------------------------------------------- X 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 X X X Xvector [1:29] - Val X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X6.081006e-02 6.720551e-02 7.427358e-02 8.208500e-02 9.071796e-02 0.100259 0.110803 0.122456 0.135335 0.149569 X X 11 12 13 14 15 16 17 18 19 20 X------------------------------------------------------------------------------- X0.165299 0.182684 0.201897 0.22313 0.246597 0.272532 0.301194 0.332871 0.367879 0.40657 X X 21 22 23 24 25 26 27 28 29 X------------------------------------------------------------------------------- X0.449329 0.496585 0.548812 0.606531 0.67032 0.740818 0.818731 0.904837 1 X X XValj = 4.482143e-02, valp = 6.081006e-02, argj=2.7 XValj = 4.776407e-02, valp = 4.482143e-02, argj=2.6 XValj = 4.729968e-02, valp = 4.776407e-02, argj=2.5 XValj = 4.736654e-02, valp = 4.729968e-02, argj=2.4 XValj = 4.735743e-02, valp = 4.736654e-02, argj=2.3 XValj = 4.736087e-02, valp = 4.735743e-02, argj=2.2 XValj = 4.737027e-02, valp = 4.736087e-02, argj=2.1 XInterpolated value to x= 3.05 is 4.736087e-02 XExact function value at this point 4.735892e-02 X X XInterpolation just to the 10. node XThe function being tabulated over [0.5:1.5] X X Xvector [1:11] - Arg - interpolation nodes X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X 1 0.9 1.1 0.8 1.2 0.7 1.3 0.6 1.4 0.5 X X 11 X------------------------------------------------------------------------------- X 1.5 X X X Xvector [1:11] - Val X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X0.367879 0.40657 0.332871 0.449329 0.301194 0.496585 0.272532 0.548812 0.246597 0.606531 X X 11 X------------------------------------------------------------------------------- X0.22313 X X XValj = 0.367879, valp = 0.367879, argj=0.9 XValj = 0.367879, valp = 0.367879, argj=1.1 XInterpolated value to x= 1 is 0.367879 XExact function value at this point 0.367879 X X XInterpolation to middle of table XThe function being tabulated over [0.1:1] X X Xvector [1:10] - Arg - interpolation nodes X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X 0.6 0.5 0.7 0.4 0.8 0.3 0.9 0.2 1 0.1 X X X Xvector [1:10] - Val X X 1 2 3 4 5 6 7 8 9 10 X------------------------------------------------------------------------------- X0.548812 0.606531 0.496585 0.67032 0.449329 0.740818 0.40657 0.818731 0.367879 0.904837 X X XValj = 0.571899, valp = 0.548812, argj=0.5 XValj = 0.57124, valp = 0.571899, argj=0.7 XValj = 0.571208, valp = 0.57124, argj=0.4 XValj = 0.571209, valp = 0.571208, argj=0.8 XValj = 0.571209, valp = 0.571209, argj=0.3 XInterpolated value to x= 0.56 is 0.571209 XExact function value at this point 0.571209 ________This_Is_The_END________ if test `wc -l < valis.dat.vax` -ne 141; then echo 'shar: valis.dat.vax was damaged during transit (should have had 141 lines)' fi echo 'x - vhjmin.c' sed 's/^X//' << '________This_Is_The_END________' > vhjmin.c X/* Verify HJMIN function */ X X#include "math.h" X#include "vector.h" X Xdouble hjmin(VECTOR b, const VECTOR h0, double (*f)(const VECTOR x)); X Xstatic int iter_count; X X /* Simplified vector print */ Xstatic void pr_vector(v) XVECTOR v; X{ X register int i; X for(i=V_lwb(v); i<=V_upb(v); i++) X printf("%9.3g ",*V_elem(v,i)); X} X X X/* Rosenbroke function */ X Xdouble f1(v) Xconst VECTOR v; X{ X register double x1 = *V_elem(v,1); X register double x2 = *V_elem(v,2); X iter_count++; X return 100*powi(x2 - x1*x1, 2) + powi(1 - x1, 2); X} X Xvoid test1() X{ X const int n = 2; X VECTOR b0 = V_build(1,n,-1.2,1.,"END"); /* Initial guess */ X VECTOR b = V_new(1,n); /* Min location found */ X VECTOR bm = V_build(1,n,1.0,1.0,"END"); /* Exact min location */ X VECTOR h0 = V_build(1,n,10.,10.,"END"); /* Initial step */ X X V_assign(b,b0); X iter_count = 0; X printf("\n\n\tRosenbroke function\n"); X printf("\n\n\tf = 100*(x2-x1^2)^2 + (1-x1)^2\n\n"); X printf("\nInitial guess b0 = "); pr_vector(b0); X printf("\nFunction value at it f0 = %g",f1(b0)); X printf("\nInitial steps h0 = "); pr_vector(h0); X printf("\n"); X printf("\nMinimum f value found f = %g", hjmin(b,h0,f1)); X printf("\n at b = "); pr_vector(b); X printf("\nExact min location bm = "); pr_vector(bm); X printf("\nNo. of iterations ni = %d\n",iter_count); X X} X X X X/* Bocks function */ X Xdouble f2(v) Xconst VECTOR v; X{ X register double x1 = *V_elem(v,1); X register double x2 = *V_elem(v,2); X iter_count++; X return powi( exp(-1./10) - exp(-x1/10) + exp(-10./10) - exp(-x2/10), 2); X} X Xvoid test2() X{ X const int n = 2; X VECTOR b0 = V_build(1,n,0.,0.,"END"); /* Initial guess */ X VECTOR b = V_new(1,n); /* Min location found */ X VECTOR bm = V_build(1,n,10.0,1.0,"END"); /* Exact min location */ X VECTOR h0 = V_build(1,n,10.,10.,"END"); /* Initial step */ X X V_assign(b,b0); X iter_count = 0; X printf("\n\n\tBocks function\n"); X printf("\n\n\tf = [ exp(-1/10) - exp(-x1/10) + exp(-10/10) -exp(-x2/10) ]^2\n\n"); X printf("\nInitial guess b0 = "); pr_vector(b0); X printf("\nFunction value at it f0 = %g",f2(b0)); X printf("\nInitial steps h0 = "); pr_vector(h0); X printf("\n"); X printf("\nMinimum f value found f = %g", hjmin(b,h0,f2)); X printf("\n at b = "); pr_vector(b); X printf("\nExact min location bm = "); pr_vector(bm); X printf("\nNo. of iterations ni = %d\n",iter_count); X X} X X X X/* Mile & Cuntrell function */ X Xdouble f3(v) Xconst VECTOR v; X{ X register double x1 = *V_elem(v,1); X register double x2 = *V_elem(v,2); X register double x3 = *V_elem(v,3); X register double x4 = *V_elem(v,4); X iter_count++; X return powi( exp(x1)-x2, 4) + 100*powi(x2-x3,6) + powi(atan(x3-x4),4) + X powi(x1,8); X} X Xvoid test3() X{ X const int n = 4; X VECTOR b0 = V_build(1,n,1.,2.,2.,2.,"END"); /* Initial guess */ X VECTOR b = V_new(1,n); /* Min location found */ X VECTOR bm = V_build(1,n,0.,1.,1.,1.,"END"); /* Exact min location */ X VECTOR h0 = V_build(1,n,10.,10.,10.,10.,"END"); /* Initial step */ X X V_assign(b,b0); X iter_count = 0; X printf("\n\n\tMile & Cuntrell function\n"); X printf("\n\n\tf = [ exp(x1)-x2 ]^4 +100(x2-x3)^6 + atan(x3-x4)^4 + x1^8\n\n"); X printf("\nInitial guess b0 = "); pr_vector(b0); X printf("\nFunction value at it f0 = %g",f3(b0)); X printf("\nInitial steps h0 = "); pr_vector(h0); X printf("\n"); X printf("\nMinimum f value found f = %g", hjmin(b,h0,f3)); X printf("\n at b = "); pr_vector(b); X printf("\nExact min location bm = "); pr_vector(bm); X printf("\nNo. of iterations ni = %d\n",iter_count); X X} X X/* Powell function */ X Xdouble f4(v) Xconst VECTOR v; X{ X register double x1 = *V_elem(v,1); X register double x2 = *V_elem(v,2); X register double x3 = *V_elem(v,3); X register double x4 = *V_elem(v,4); X iter_count++; X return powi( x1+10*x2, 2) + 5*powi(x3-x4,2) + powi(x2-2*x3,4) + X 10*powi(x1-x4,4); X} X Xvoid test4() X{ X const int n = 4; X VECTOR b0 = V_build(1,n,3.,-1.,0.,1.,"END"); /* Initial guess */ X VECTOR b = V_new(1,n); /* Min location found */ X VECTOR bm = V_build(1,n,0.,0.,0.,0.,"END"); /* Exact min location */ X VECTOR h0 = V_build(1,n,10.,10.,10.,10.,"END"); /* Initial step */ X X V_assign(b,b0); X iter_count = 0; X printf("\n\n\tPowell function\n"); X printf("\n\n\tf = (x1+10*x2)^2 + 5(x3-x4)^2 + (x2-2x3)^4 + 10(x1-x4)^4\n\n"); X printf("\nInitial guess b0 = "); pr_vector(b0); X printf("\nFunction value at it f0 = %g",f4(b0)); X printf("\nInitial steps h0 = "); pr_vector(h0); X printf("\n"); X printf("\nMinimum f value found f = %g", hjmin(b,h0,f4)); X printf("\n at b = "); pr_vector(b); X printf("\nExact min location bm = "); pr_vector(bm); X printf("\nNo. of iterations ni = %d\n",iter_count); X X} X X Xmain() X{ X printf("\n\n\n\tVerify HJMIN program to minimize a function value\n"); X test1(); X test2(); X test3(); X test4(); X} ________This_Is_The_END________ if test `wc -l < vhjmin.c` -ne 174; then echo 'shar: vhjmin.c was damaged during transit (should have had 174 lines)' fi echo 'x - vhjmin.dat' sed 's/^X//' << '________This_Is_The_END________' > vhjmin.dat X X X X Verify HJMIN program to minimize a function value X X X Rosenbroke function X X X f = 100*(x2-x1^2)^2 + (1-x1)^2 X X XInitial guess b0 = -1.2 1 XFunction value at it f0 = 24.2 XInitial steps h0 = 10 10 X XMinimum f value found f = 6.97611e-11 X at b = 1 1 XExact min location bm = 1 1 XNo. of iterations ni = 270 X X X Bocks function X X X f = [ exp(-1/10) - exp(-x1/10) + exp(-10/10) -exp(-x2/10) ]^2 X X XInitial guess b0 = 0 0 XFunction value at it f0 = 0.528941 XInitial steps h0 = 10 10 X XMinimum f value found f = 0 X at b = 10 1 XExact min location bm = 10 1 XNo. of iterations ni = 89 X X X Mile & Cuntrell function X X X f = [ exp(x1)-x2 ]^4 +100(x2-x3)^6 + atan(x3-x4)^4 + x1^8 X X XInitial guess b0 = 1 2 2 2 XFunction value at it f0 = 1.26618 XInitial steps h0 = 10 10 10 10 X XMinimum f value found f = 4.1179e-27 X at b = 1.49e-08 1 1 1 XExact min location bm = 0 1 1 1 XNo. of iterations ni = 227 X X X Powell function X X X f = (x1+10*x2)^2 + 5(x3-x4)^2 + (x2-2x3)^4 + 10(x1-x4)^4 X X XInitial guess b0 = 3 -1 0 1 XFunction value at it f0 = 215 XInitial steps h0 = 10 10 10 10 X XMinimum f value found f = 0 X at b = 0 0 0 0 XExact min location bm = 0 0 0 0 XNo. of iterations ni = 154 ________This_Is_The_END________ if test `wc -l < vhjmin.dat` -ne 68; then echo 'shar: vhjmin.dat was damaged during transit (should have had 68 lines)' fi echo 'x - vhjmin.dat.vax' sed 's/^X//' << '________This_Is_The_END________' > vhjmin.dat.vax X X X X Verify HJMIN program to minimize a function value X X X Rosenbroke function X X X f = 100*(x2-x1^2)^2 + (1-x1)^2 X X XInitial guess b0 = -1.2 1 XFunction value at it f0 = 24.2 XInitial steps h0 = 10 10 X XMinimum f value found f = 2.402913e-10 X at b = 1 1 XExact min location bm = 1 1 XNo. of iterations ni = 272 X X X Bocks function X X X f = [ exp(-1/10) - exp(-x1/10) + exp(-10/10) -exp(-x2/10) ]^2 X X XInitial guess b0 = 0 0 XFunction value at it f0 = 0.528941 XInitial steps h0 = 10 10 X XMinimum f value found f = 0 X at b = 10 1 XExact min location bm = 10 1 XNo. of iterations ni = 93 X X X Mile & Cuntrell function X X X f = [ exp(x1)-x2 ]^4 +100(x2-x3)^6 + atan(x3-x4)^4 + x1^8 X X XInitial guess b0 = 1 2 2 2 XFunction value at it f0 = 1.26618 XInitial steps h0 = 10 10 10 10 X XMinimum f value found f = 7.888609e-27 X at b = 5.960e-08 1 1 1 XExact min location bm = 0 1 1 1 XNo. of iterations ni = 235 X X X Powell function X X X f = (x1+10*x2)^2 + 5(x3-x4)^2 + (x2-2x3)^4 + 10(x1-x4)^4 X X XInitial guess b0 = 3 -1 0 1 XFunction value at it f0 = 215 XInitial steps h0 = 10 10 10 10 X XMinimum f value found f = 0 X at b = 0 0 0 0 XExact min location bm = 0 0 0 0 XNo. of iterations ni = 162 ________This_Is_The_END________ if test `wc -l < vhjmin.dat.vax` -ne 68; then echo 'shar: vhjmin.dat.vax was damaged during transit (should have had 68 lines)' fi .