cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
optimize_subplx.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 #include "cddefines.h"
4 #include "optimize.h"
5 
6 STATIC void calcc(long,realnum*,long,long,int,realnum[]);
7 STATIC double cdasum(long,realnum[],long);
8 STATIC void cdaxpy(long,double,realnum[],long,realnum[],long);
9 STATIC void cdcopy(long,realnum[],long,realnum[],long);
10 STATIC void csscal(long,double,realnum[],long);
11 STATIC double dist(long,realnum[],realnum[]);
12 STATIC void evalf(long,long[],realnum[],long,realnum[],realnum*,long*);
13 STATIC void fstats(double,long,int);
14 STATIC void newpt(long,double,realnum[],realnum[],int,realnum[],int*);
15 STATIC void order(long,realnum[],long*,long*,long*);
16 STATIC void partx(long,long[],realnum[],long*,long[]);
17 STATIC void setstp(long,long,realnum[],realnum[]);
18 STATIC void simplx(long,realnum[],long,long[],long,int,realnum[],
19  realnum*,long*,realnum*,realnum[],long*);
20 STATIC void sortd(long,realnum[],long[]);
21 STATIC void start(long,realnum[],realnum[],long,long[],realnum*,int*);
22 STATIC void subopt(long);
23 
24 /*lint -e801 goto is deprecated */
25 /*lint -e64 type mismatch */
26 /*********************************************************************
27  *
28  * NB this is much the original coding and does not use strong typing
29  *gjf mod: to avoid conflict with exemplar parallel version of lapack,
30  *dasum changed to cdasum
31  *daxpy changed to cdaxpy
32  *dcopy changed to cdcopy
33  *sscal changed to csscal
34  *
35  *********************************************************************** */
36 struct t_usubc {
38  beta,
39  gamma,
40  delta,
41  psi,
42  omega;
43  long int nsmin,
44  nsmax,
45  irepl,
46  ifxsw;
48  fstop;
49  long int nfstop,
50  nfxe;
52  ftest;
53  int minf,
54  initx,
55  newx;
56  } usubc;
57 struct t_isubc {
59  sfstop,
60  sfbest;
61  int IntNew;
62  } isubc;
63 
64 void optimize_subplex(long int n,
65  double tol,
66  long int maxnfe,
67  long int mode,
68  realnum scale[],
69  realnum x[],
70  realnum *fx,
71  long int *nfe,
72  realnum work[],
73  long int iwork[],
74  long int *iflag)
75 {
76  static int cmode;
77  static long int i,
78  i_,
79  ifsptr,
80  ins,
81  insfnl,
82  insptr,
83  ipptr,
84  isptr,
85  istep,
86  istptr,
87  nnn,
88  ns,
89  nsubs;
90  static realnum dum,
91  scl[1],
92  sfx,
93  xpscl,
94  xstop,
95  xstop2;
96 
97  static const realnum bnsfac[2][3] = { {-1.0f,-2.0f,0.0f}, {1.0f,0.0f,2.0f} };
98 
99  DEBUG_ENTRY( "optimize_subplex()" );
100 
101  /* Coded by Tom Rowan
102  * Department of Computer Sciences
103  * University of Texas at Austin
104  *
105  * Jason Ferguson:
106  * Changed on 8/4/94, in order to incorparate into cloudy
107  * changes made are: double precision to real,
108  * a 2-D function f(n,x) into 1-D f(x),
109  * and the termination check.
110  *
111  * optimize_subplex uses the subplex method to solve unconstrained
112  * optimization problems. The method is well suited for
113  * optimizing objective functions that are noisy or are
114  * discontinuous at the solution.
115  *
116  * optimize_subplex sets default optimization options by calling the
117  * subroutine subopt. The user can override these defaults
118  * by calling subopt prior to calling optimize_subplex, changing the
119  * appropriate common variables, and setting the value of
120  * mode as indicated below.
121  *
122  * By default, optimize_subplex performs minimization.
123  *
124  * input
125  *
126  * ffff - user supplied function f(n,x) to be optimized,
127  * declared external in calling routine
128  * always uses optimize_func - change 97 dec 8
129  *
130  * n - problem dimension
131  *
132  * tol - relative error tolerance for x (tol .ge. 0.)
133  *
134  * maxnfe - maximum number of function evaluations
135  *
136  * mode - integer mode switch with binary expansion
137  * (bit 1) (bit 0) :
138  * bit 0 = 0 : first call to optimize_subplex
139  * = 1 : continuation of previous call
140  * bit 1 = 0 : use default options
141  * = 1 : user set options
142  *
143  * scale - scale and initial stepsizes for corresponding
144  * components of x
145  * (If scale(1) .lt. 0.,
146  * ABS(scale(1)) is used for all components of x,
147  * and scale(2),...,scale(n) are not referenced.)
148  *
149  * x - starting guess for optimum
150  *
151  * work - real work array of dimension .ge.
152  * 2*n + nsmax*(nsmax+4) + 1
153  * (nsmax is set in subroutine subopt.
154  * default: nsmax = MIN(5,n))
155  *
156  * iwork - integer work array of dimension .ge.
157  * n + INT(n/nsmin)
158  * (nsmin is set in subroutine subopt.
159  * default: nsmin = MIN(2,n))
160  *
161  * output
162  *
163  * x - computed optimum
164  *
165  * fx - value of f at x
166  *
167  * nfe - number of function evaluations
168  *
169  * iflag - error flag
170  * = -2 : invalid input
171  * = -1 : maxnfe exceeded
172  * = 0 : tol satisfied
173  * = 1 : limit of machine precision
174  * = 2 : fstop reached (fstop usage is determined
175  * by values of options minf, nfstop, and
176  * irepl. default: f(x) not tested against
177  * fstop)
178  * iflag should not be reset between calls to
179  * optimize_subplex.
180  *
181  * common
182  * */
183 
184  if( (mode%2) == 0 )
185  {
186 
187  /* first call, check input
188  * */
189  if( n < 1 )
190  goto L_120;
191  if( tol < 0.0f )
192  goto L_120;
193  if( maxnfe < 1 )
194  goto L_120;
195  if( scale[0] >= 0.0f )
196  {
197  for( i=1; i <= n; i++ )
198  {
199  i_ = i - 1;
200  xpscl = x[i_] + scale[i_];
201  if( fp_equal( xpscl, x[i_] ) )
202  goto L_120;
203  }
204  }
205  else
206  {
207  scl[0] = (realnum)fabs(scale[0]);
208  for( i=1; i <= n; i++ )
209  {
210  i_ = i - 1;
211  xpscl = x[i_] + scl[0];
212  if( fp_equal( xpscl, x[i_] ) )
213  goto L_120;
214  }
215  }
216  if( ((mode/2)%2) == 0 )
217  {
218  subopt(n);
219  }
220  else if( usubc.alpha <= 0.0f )
221  {
222  goto L_120;
223  }
224  else if( usubc.beta <= 0.0f || usubc.beta >= 1.0f )
225  {
226  goto L_120;
227  }
228  else if( usubc.gamma <= 1.0f )
229  {
230  goto L_120;
231  }
232  else if( usubc.delta <= 0.0f || usubc.delta >= 1.0f )
233  {
234  goto L_120;
235  }
236  else if( usubc.psi <= 0.0f || usubc.psi >= 1.0f )
237  {
238  goto L_120;
239  }
240  else if( usubc.omega <= 0.0f || usubc.omega >= 1.0f )
241  {
242  goto L_120;
243  }
244  else if( (usubc.nsmin < 1 || usubc.nsmax < usubc.nsmin) ||
245  n < usubc.nsmax )
246  {
247  goto L_120;
248  }
249  else if( n < ((n - 1)/usubc.nsmax + 1)*usubc.nsmin )
250  {
251  goto L_120;
252  }
253  else if( usubc.irepl < 0 || usubc.irepl > 2 )
254  {
255  goto L_120;
256  }
257  else if( usubc.ifxsw < 1 || usubc.ifxsw > 3 )
258  {
259  goto L_120;
260  }
261  else if( usubc.bonus < 0.0f )
262  {
263  goto L_120;
264  }
265  else if( usubc.nfstop < 0 )
266  {
267  goto L_120;
268  }
269 
270  /* initialization
271  * */
272  istptr = n + 1;
273  isptr = istptr + n;
274  ifsptr = isptr + usubc.nsmax*(usubc.nsmax + 3);
275  insptr = n + 1;
276  if( scale[0] > 0.0f )
277  {
278  cdcopy(n,scale,1,work,1);
279  cdcopy(n,scale,1,&work[istptr-1],1);
280  }
281  else
282  {
283  cdcopy(n,scl,0,work,1);
284  cdcopy(n,scl,0,&work[istptr-1],1);
285  }
286  for( i=1; i <= n; i++ )
287  {
288  i_ = i - 1;
289  iwork[i_] = i;
290  }
291  *nfe = 0;
292  usubc.nfxe = 1;
293  if( usubc.irepl == 0 )
294  {
295  isubc.fbonus = 0.0f;
296  }
297  else if( usubc.minf )
298  {
299  isubc.fbonus = bnsfac[0][usubc.ifxsw-1]*usubc.bonus;
300  }
301  else
302  {
303  isubc.fbonus = bnsfac[1][usubc.ifxsw-1]*usubc.bonus;
304  }
305  if( usubc.nfstop == 0 )
306  {
307  isubc.sfstop = 0.0f;
308  }
309  else if( usubc.minf )
310  {
312  }
313  else
314  {
315  isubc.sfstop = -usubc.fstop;
316  }
317  usubc.ftest = 0.0f;
318  cmode = false;
319  isubc.IntNew = true;
320  usubc.initx = true;
321  nnn = 0;
322  evalf(nnn,iwork,(realnum*)&dum,n,x,&sfx,nfe);
323  usubc.initx = false;
324 
325  /* continuation of previous call
326  * */
327  }
328  else if( *iflag == 2 )
329  {
330  if( usubc.minf )
331  {
333  }
334  else
335  {
336  isubc.sfstop = -usubc.fstop;
337  }
338  cmode = true;
339  goto L_70;
340  }
341  else if( *iflag == -1 )
342  {
343  cmode = true;
344  goto L_70;
345  }
346  else if( *iflag == 0 )
347  {
348  cmode = false;
349  goto L_90;
350  }
351  else
352  {
353  return;
354  }
355 
356  /* subplex loop
357  * */
358 L_40:
359 
360  for( i=0; i < n; i++ )
361  {
362  work[i] = (realnum)fabs(work[i]);
363  }
364 
365  sortd(n,work,iwork);
366  partx(n,iwork,work,&nsubs,&iwork[insptr-1]);
367  cdcopy(n,x,1,work,1);
368  ins = insptr;
369  insfnl = insptr + nsubs - 1;
370  ipptr = 1;
371 
372  /* simplex loop
373  * */
374 
375 L_60:
376  ns = iwork[ins-1];
377 
378 L_70:
379  simplx(n,&work[istptr-1],ns,&iwork[ipptr-1],maxnfe,cmode,x,&sfx,
380  nfe,&work[isptr-1],&work[ifsptr-1],iflag);
381 
382  cmode = false;
383  if( *iflag != 0 )
384  goto L_121;
385 
386  if( ins < insfnl )
387  {
388  ins += 1;
389  ipptr += ns;
390  goto L_60;
391  }
392 
393  /* end simplex loop
394  * */
395  for( i=0; i < n; i++ )
396  {
397  work[i] = x[i] - work[i];
398  }
399 
400  /* check termination
401  * */
402 L_90:
403  /* new stop criterion: calculate distance between
404  * previous and new best model, if distance is
405  * smaller than tol return. */
406  istep = istptr-1;
407  xstop = 0.f;
408  xstop2 = 0.f;
409  for( i=0; i < n; i++ )
410  {
411  realnum ttemp;
412  xstop += POW2(work[i]);
413  ttemp = (realnum)fabs(work[istep]);
414  /* chng to avoid side effects
415  * xstop2 = MAX2(xstop2,(realnum)fabs(work[istep]));*/
416  xstop2 = MAX2( xstop2 , ttemp );
417  istep++;
418  }
419 
420  if( sqrt(xstop) > tol || xstop2*usubc.psi > (realnum)tol )
421  {
422  setstp(nsubs,n,work,&work[istptr-1]);
423  goto L_40;
424  }
425 
426  /* end subplex loop
427  * */
428 
429  *iflag = 0;
430 L_121:
431  if( usubc.minf )
432  {
433  *fx = sfx;
434  }
435  else
436  {
437  *fx = -sfx;
438  }
439  return;
440 
441  /* invalid input
442  * */
443 L_120:
444 
445  *iflag = -2;
446  return;
447 }
448 /********************************************************************* */
449 STATIC void subopt(long int n)
450 {
451 
452  DEBUG_ENTRY( "subopt()" );
453 
454 
455  /* Coded by Tom Rowan
456  * Department of Computer Sciences
457  * University of Texas at Austin
458  *
459  * subopt sets options for optimize_subplex.
460  *
461  * input
462  *
463  * n - problem dimension
464  *
465  * common
466  * */
467 
468 
469 
470  /* subroutines and functions
471  *
472  * fortran
473  *
474  *-----------------------------------------------------------
475  *
476  ************************************************************
477  * simplex method strategy parameters
478  ************************************************************
479  *
480  * alpha - reflection coefficient
481  * alpha .gt. 0
482  * */
483  usubc.alpha = 1.0f;
484 
485  /* beta - contraction coefficient
486  * 0 .lt. beta .lt. 1
487  * */
488  usubc.beta = .50f;
489 
490  /* gamma - expansion coefficient
491  * gamma .gt. 1
492  * */
493  usubc.gamma = 2.0f;
494 
495  /* delta - shrinkage (massive contraction) coefficient
496  * 0 .lt. delta .lt. 1
497  * */
498  usubc.delta = .50f;
499 
500  /************************************************************
501  * subplex method strategy parameters
502  ************************************************************
503  *
504  * psi - simplex reduction coefficient
505  * 0 .lt. psi .lt. 1
506  * */
507  usubc.psi = .250f;
508 
509  /* omega - step reduction coefficient
510  * 0 .lt. omega .lt. 1
511  * */
512  usubc.omega = .10f;
513 
514  /* nsmin and nsmax specify a range of subspace dimensions.
515  * In addition to satisfying 1 .le. nsmin .le. nsmax .le. n,
516  * nsmin and nsmax must be chosen so that n can be expressed
517  * as a sum of positive integers where each of these integers
518  * ns(i) satisfies nsmin .le. ns(i) .ge. nsmax.
519  * Specifically,
520  * nsmin*ceil(n/nsmax) .le. n must be true.
521  *
522  * nsmin - subspace dimension minimum
523  * */
524  usubc.nsmin = MIN2(2,n);
525 
526  /* nsmax - subspace dimension maximum
527  * */
528  usubc.nsmax = MIN2(5,n);
529 
530  /************************************************************
531  * subplex method special cases
532  ************************************************************
533  * nelder-mead simplex method with periodic restarts
534  * nsmin = nsmax = n
535  ************************************************************
536  * nelder-mead simplex method
537  * nsmin = nsmax = n, psi = small positive
538  ************************************************************
539  *
540  * irepl, ifxsw, and bonus deal with measurement replication.
541  * Objective functions subject to large amounts of noise can
542  * cause an optimization method to halt at a false optimum.
543  * An expensive solution to this problem is to evaluate f
544  * several times at each point and return the average (or max
545  * or min) of these trials as the function value. optimize_subplex
546  * performs measurement replication only at the current best
547  * point. The longer a point is retained as best, the more
548  * accurate its function value becomes.
549  *
550  * The common variable nfxe contains the number of function
551  * evaluations at the current best point. fxstat contains the
552  * mean, max, min, and standard deviation of these trials.
553  *
554  * irepl - measurement replication switch
555  * irepl = 0, 1, or 2
556  * = 0 : no measurement replication
557  * = 1 : optimize_subplex performs measurement replication
558  * = 2 : user performs measurement replication
559  * (This is useful when optimizing on the mean,
560  * max, or min of trials is insufficient. Common
561  * variable initx is true for first function
562  * evaluation. newx is true for first trial at
563  * this point. The user uses subroutine fstats
564  * within his objective function to maintain
565  * fxstat. By monitoring newx, the user can tell
566  * whether to return the function evaluation
567  * (newx = .true.) or to use the new function
568  * evaluation to refine the function evaluation
569  * of the current best point (newx = .false.).
570  * The common variable ftest gives the function
571  * value that a new point must beat to be
572  * considered the new best point.)
573  * */
574  usubc.irepl = 0;
575 
576  /* ifxsw - measurement replication optimization switch
577  * ifxsw = 1, 2, or 3
578  * = 1 : retain mean of trials as best function value
579  * = 2 : retain max
580  * = 3 : retain min
581  * */
582  usubc.ifxsw = 1;
583 
584  /* Since the current best point will also be the most
585  * accurately evaluated point whenever irepl .gt. 0, a bonus
586  * should be added to the function value of the best point
587  * so that the best point is not replaced by a new point
588  * that only appears better because of noise.
589  * optimize_subplex uses bonus to determine how many multiples of
590  * fxstat(4) should be added as a bonus to the function
591  * evaluation. (The bonus is adjusted automatically by
592  * optimize_subplex when ifxsw or minf is changed.)
593  *
594  * bonus - measurement replication bonus coefficient
595  * bonus .ge. 0 (normally, bonus = 0 or 1)
596  * = 0 : bonus not used
597  * = 1 : bonus used
598  * */
599  usubc.bonus = 1.0f;
600 
601  /* nfstop = 0 : f(x) is not tested against fstop
602  * = 1 : if f(x) has reached fstop, optimize_subplex returns
603  * iflag = 2
604  * = 2 : (only valid when irepl .gt. 0)
605  * if f(x) has reached fstop and
606  * nfxe .gt. nfstop, optimize_subplex returns iflag = 2
607  * */
608  usubc.nfstop = 0;
609 
610  /* fstop - f target value
611  * Its usage is determined by the value of nfstop.
612  *
613  * minf - logical switch
614  * = .true. : optimize_subplex performs minimization
615  * = .false. : optimize_subplex performs maximization
616  * */
617  usubc.minf = true;
618  return;
619 }
620 /**********************************************************************/
621 /* >>chng 01 jan 03, cleaned up -1 and formatting in this routine */
622 STATIC void cdcopy(long int n,
623  realnum dx[],
624  long int incx,
625  realnum dy[],
626  long int incy)
627 {
628  long int i,
629  ix,
630  iy,
631  m;
632 
633  DEBUG_ENTRY( "cdcopy()" );
634 
635  /*
636  * copies a vector, x, to a vector, y.
637  * uses unrolled loops for increments equal to one.
638  * Jack Dongarra, lapack, 3/11/78.
639  */
640 
641  if( n > 0 )
642  {
643  if( incx == 1 && incy == 1 )
644  {
645 
646  /* code for both increments equal to 1 */
647 
648  /* first the clean-up loop */
649  m = n%7;
650  if( m != 0 )
651  {
652  for( i=0; i < m; i++ )
653  {
654  dy[i] = dx[i];
655  }
656  if( n < 7 )
657  {
658  return;
659  }
660  }
661 
662  for( i=m; i < n; i += 7 )
663  {
664  dy[i] = dx[i];
665  dy[i+1] = dx[i+1];
666  dy[i+2] = dx[i+2];
667  dy[i+3] = dx[i+3];
668  dy[i+4] = dx[i+4];
669  dy[i+5] = dx[i+5];
670  dy[i+6] = dx[i+6];
671  }
672  }
673  else
674  {
675 
676  /* code for unequal increments or equal increments
677  * not equal to 1 */
678 
679  ix = 1;
680  iy = 1;
681  if( incx < 0 )
682  ix = (-n + 1)*incx + 1;
683  if( incy < 0 )
684  iy = (-n + 1)*incy + 1;
685  for( i=0; i < n; i++ )
686  {
687  dy[iy-1] = dx[ix-1];
688  ix += incx;
689  iy += incy;
690  }
691  }
692  }
693  return;
694 }
695 /********************************************************************* */
696 STATIC void evalf(long int ns,
697  long int ips[],
698  realnum xs[],
699  long int n,
700  realnum x[],
701  realnum *sfx,
702  long int *nfe)
703 {
704  static int newbst;
705  static long int i,
706  i_;
707  static realnum fx;
708  /* gary delete since in header */
709  /*double optimize_func();*/
710 
711  DEBUG_ENTRY( "evalf()" );
712  /* gary add, not used, so trick compiler notes */
713  i_ = n;
714 
715  /*>>chng 97 dec 07, implicit nonte, rid of first function argument
716  *
717  * Coded by Tom Rowan
718  * Department of Computer Sciences
719  * University of Texas at Austin
720  *
721  * evalf evaluates the function f at a point defined by x
722  * with ns of its components replaced by those in xs.
723  *
724  * input
725  *
726  * f - user supplied function f(n,x) to be optimized
727  * changed to optimize_func - cannot specify arbitrary function now
728  *
729  * ns - subspace dimension
730  *
731  * ips - permutation vector
732  *
733  * xs - real ns-vector to be mapped to x
734  *
735  * n - problem dimension
736  *
737  * x - real n-vector
738  *
739  * nfe - number of function evaluations
740  *
741  * output
742  *
743  * sfx - signed value of f evaluated at x
744  *
745  * nfe - incremented number of function evaluations
746  *
747  * common
748  * */
749 
750 
751 
752 
753  /* local variables
754  * */
755 
756 
757  /* subroutines and functions
758  * */
759 
760  /*-----------------------------------------------------------
761  * */
762  for( i=1; i <= ns; i++ )
763  {
764  i_ = i - 1;
765  x[ips[i_]-1] = xs[i_];
766  }
767  usubc.newx = isubc.IntNew || usubc.irepl != 2;
768  fx = (realnum)optimize_func(x);
769  /* fx = f(n,x) */
770  if( usubc.irepl == 0 )
771  {
772  if( usubc.minf )
773  {
774  *sfx = fx;
775  }
776  else
777  {
778  *sfx = -fx;
779  }
780  }
781  else if( isubc.IntNew )
782  {
783  if( usubc.minf )
784  {
785  *sfx = fx;
786  newbst = fx < usubc.ftest;
787  }
788  else
789  {
790  *sfx = -fx;
791  newbst = fx > usubc.ftest;
792  }
793  if( usubc.initx || newbst )
794  {
795  if( usubc.irepl == 1 )
796  fstats(fx,1,true);
797  usubc.ftest = fx;
798  isubc.sfbest = *sfx;
799  }
800  }
801  else
802  {
803  if( usubc.irepl == 1 )
804  {
805  fstats(fx,1,false);
806  fx = usubc.fxstat[usubc.ifxsw-1];
807  }
808  usubc.ftest = fx + isubc.fbonus*usubc.fxstat[3];
809  if( usubc.minf )
810  {
811  *sfx = usubc.ftest;
812  isubc.sfbest = fx;
813  }
814  else
815  {
816  *sfx = -usubc.ftest;
817  isubc.sfbest = -fx;
818  }
819  }
820  *nfe += 1;
821  return;
822 }
823 
824 /********************************************************************* */
825 STATIC void setstp(long int nsubs,
826  long int n,
827  realnum deltax[],
828  realnum step[])
829 {
830  static long int i,
831  i_;
832  static realnum stpfac;
833 
834  DEBUG_ENTRY( "setstp()" );
835 
836 
837  /* Coded by Tom Rowan
838  * Department of Computer Sciences
839  * University of Texas at Austin
840  *
841  * setstp sets the stepsizes for the corresponding components
842  * of the solution vector.
843  *
844  * input
845  *
846  * nsubs - number of subspaces
847  *
848  * n - number of components (problem dimension)
849  *
850  * deltax - vector of change in solution vector
851  *
852  * step - stepsizes for corresponding components of
853  * solution vector
854  *
855  * output
856  *
857  * step - new stepsizes
858  *
859  * common
860  * */
861 
862 
863  /* local variables
864  * */
865 
866 
867  /* subroutines and functions
868  *
869  * blas */
870  /* fortran
871  *
872  *-----------------------------------------------------------
873  *
874  * set new step
875  * */
876  if( nsubs > 1 )
877  {
878  double a , b , c;
879  a = cdasum(n,deltax,1);
880  b = cdasum(n,step,1);
881  c = MAX2(a/b ,usubc.omega);
882  stpfac = (realnum)MIN2(c , 1.0f/usubc.omega);
883  }
884  else
885  {
886  stpfac = usubc.psi;
887  }
888  csscal(n,stpfac,step,1);
889 
890  /* reorient simplex
891  * */
892  for( i=1; i <= n; i++ )
893  {
894  i_ = i - 1;
895  if( deltax[i_] == 0.f )
896  {
897  step[i_] = -step[i_];
898  }
899  else
900  {
901  step[i_] = (realnum)sign(step[i_],deltax[i_]);
902  }
903  }
904  return;
905 }
906 
907 /********************************************************************* */
908 STATIC void sortd(long int n,
909  realnum xkey[],
910  long int ix[])
911 {
912  long int i,
913  i_,
914  ifirst,
915  ilast,
916  iswap,
917  ixi,
918  ixip1;
919 
920  DEBUG_ENTRY( "sortd()" );
921 
922 
923  /* Coded by Tom Rowan
924  * Department of Computer Sciences
925  * University of Texas at Austin
926  *
927  * sortd uses the shakersort method to sort an array of keys
928  * in decreasing order. The sort is performed implicitly by
929  * modifying a vector of indices.
930  *
931  * For nearly sorted arrays, sortd requires O(n) comparisons.
932  * for completely unsorted arrays, sortd requires O(n**2)
933  * comparisons and will be inefficient unless n is small.
934  *
935  * input
936  *
937  * n - number of components
938  *
939  * xkey - real vector of keys
940  *
941  * ix - integer vector of indices
942  *
943  * output
944  *
945  * ix - indices satisfy xkey(ix(i)) .ge. xkey(ix(i+1))
946  * for i = 1,...,n-1
947  *
948  * local variables
949  * */
950 
951  /*-----------------------------------------------------------
952  * */
953  ifirst = 1;
954  iswap = 1;
955  ilast = n - 1;
956  while( ifirst <= ilast )
957  {
958  for( i=ifirst; i <= ilast; i++ )
959  {
960  i_ = i - 1;
961  ixi = ix[i_];
962  ixip1 = ix[i_+1];
963  if( xkey[ixi-1] < xkey[ixip1-1] )
964  {
965  ix[i_] = ixip1;
966  ix[i_+1] = ixi;
967  iswap = i;
968  }
969  }
970  ilast = iswap - 1;
971  for( i=ilast; i >= ifirst; i-- )
972  {
973  i_ = i - 1;
974  ixi = ix[i_];
975  ixip1 = ix[i_+1];
976  if( xkey[ixi-1] < xkey[ixip1-1] )
977  {
978  ix[i_] = ixip1;
979  ix[i_+1] = ixi;
980  iswap = i;
981  }
982  }
983  ifirst = iswap + 1;
984  }
985  return;
986 }
987 
988 /********************************************************************* */
989 STATIC void partx(long int n,
990  long int ip[],
991  realnum absdx[],
992  long int *nsubs,
993  long int nsvals[])
994 {
995  static long int i,
996  limit,
997  nleft,
998  ns1,
999  ns1_,
1000  ns2,
1001  nused;
1002  static realnum as1,
1003  as1max,
1004  as2,
1005  asleft,
1006  gap,
1007  gapmax;
1008 
1009  DEBUG_ENTRY( "partx()" );
1010 
1011 
1012  /* Coded by Tom Rowan
1013  * Department of Computer Sciences
1014  * University of Texas at Austin
1015  *
1016  * partx partitions the vector x by grouping components of
1017  * similar magnitude of change.
1018  *
1019  * input
1020  *
1021  * n - number of components (problem dimension)
1022  *
1023  * ip - permutation vector
1024  *
1025  * absdx - vector of magnitude of change in x
1026  *
1027  * nsvals - integer array dimensioned .ge. INT(n/nsmin)
1028  *
1029  * output
1030  *
1031  * nsubs - number of subspaces
1032  *
1033  * nsvals - integer array of subspace dimensions
1034  *
1035  * common
1036  * */
1037 
1038 
1039  /* local variables
1040  * */
1041 
1042 
1043  /* subroutines and functions
1044  *
1045  *
1046  *-----------------------------------------------------------
1047  * */
1048  *nsubs = 0;
1049  nused = 0;
1050  nleft = n;
1051  asleft = absdx[0];
1052  for( i=1; i < n; i++ )
1053  {
1054  asleft += absdx[i];
1055  }
1056 
1057  while( nused < n )
1058  {
1059  *nsubs += 1;
1060  as1 = 0.0f;
1061  for( i=0; i < (usubc.nsmin - 1); i++ )
1062  {
1063  as1 += absdx[ip[nused+i]-1];
1064  }
1065 
1066  gapmax = -1.0f;
1067  limit = MIN2(usubc.nsmax,nleft);
1068  for( ns1=usubc.nsmin; ns1 <= limit; ns1++ )
1069  {
1070  ns1_ = ns1 - 1;
1071  as1 += absdx[ip[nused+ns1_]-1];
1072  ns2 = nleft - ns1;
1073  if( ns2 > 0 )
1074  {
1075  if( ns2 >= ((ns2 - 1)/usubc.nsmax + 1)*usubc.nsmin )
1076  {
1077  as2 = asleft - as1;
1078  gap = as1/ns1 - as2/ns2;
1079  if( gap > gapmax )
1080  {
1081  gapmax = gap;
1082  nsvals[*nsubs-1] = ns1;
1083  as1max = as1;
1084  }
1085  }
1086  }
1087  else if( as1/ns1 > gapmax )
1088  {
1089  goto L_21;
1090  }
1091  }
1092  nused += nsvals[*nsubs-1];
1093  nleft = n - nused;
1094  asleft -= as1max;
1095  }
1096  return;
1097 L_21:
1098  nsvals[*nsubs-1] = ns1;
1099  return;
1100 }
1101 
1102 /********************************************************************* */
1103 
1104 #undef S
1105 #define S(I_,J_) (*(s+(I_)*(ns)+(J_)))
1106 
1107 STATIC void simplx(long int n,
1108  realnum step[],
1109  long int ns,
1110  long int ips[],
1111  long int maxnfe,
1112  int cmode,
1113  realnum x[],
1114  realnum *fx,
1115  long int *nfe,
1116  realnum *s,
1117  realnum fs[],
1118  long int *iflag)
1119 {
1120  static int small,
1121  updatc;
1122 
1123  static long int i,
1124  icent,
1125  ih,
1126  il,
1127  inew = 0,
1128  is,
1129  itemp,
1130  j,
1131  j_,
1132  npts;
1133 
1134  static realnum dum,
1135  fc,
1136  fe,
1137  fr,
1138  tol;
1139 
1140  DEBUG_ENTRY( "simplx()" );
1141 
1142 
1143  /* Coded by Tom Rowan
1144  * Department of Computer Sciences
1145  * University of Texas at Austin
1146  *
1147  * simplx uses the Nelder-Mead simplex method to minimize the
1148  * function f on a subspace.
1149  *
1150  * input
1151  *
1152  * ffff - function to be minimized, declared external in
1153  * calling routine
1154  * removed - always calls optimize_func
1155  *
1156  * n - problem dimension
1157  *
1158  * step - stepsizes for corresponding components of x
1159  *
1160  * ns - subspace dimension
1161  *
1162  * ips - permutation vector
1163  *
1164  * maxnfe - maximum number of function evaluations
1165  *
1166  * cmode - logical switch
1167  * = .true. : continuation of previous call
1168  * = .false. : first call
1169  *
1170  * x - starting guess for minimum
1171  *
1172  * fx - value of f at x
1173  *
1174  * nfe - number of function evaluations
1175  *
1176  * s - real work array of dimension .ge.
1177  * ns*(ns+3) used to store simplex
1178  *
1179  * fs - real work array of dimension .ge.
1180  * ns+1 used to store function values of simplex
1181  * vertices
1182  *
1183  * output
1184  *
1185  * x - computed minimum
1186  *
1187  * fx - value of f at x
1188  *
1189  * nfe - incremented number of function evaluations
1190  *
1191  * iflag - error flag
1192  * = -1 : maxnfe exceeded
1193  * = 0 : simplex reduced by factor of psi
1194  * = 1 : limit of machine precision
1195  * = 2 : reached fstop
1196  *
1197  * common
1198  * */
1199 
1200 
1201 
1202 
1203  /* local variables
1204  * */
1205 
1206 
1207  /* subroutines and functions
1208  *
1209  * external optimize_func,calcc,dist,evalf,newpt,order,start */
1210  /* blas */
1211 
1212  /*-----------------------------------------------------------
1213  * */
1214  if( cmode )
1215  goto L_50;
1216  npts = ns + 1;
1217  icent = ns + 2;
1218  itemp = ns + 3;
1219  updatc = false;
1220  start(n,x,step,ns,ips,s,&small);
1221  if( small )
1222  {
1223  *iflag = 1;
1224  return;
1225  }
1226  else
1227  {
1228  if( usubc.irepl > 0 )
1229  {
1230  isubc.IntNew = false;
1231  evalf(ns,ips,&S(0,0),n,x,&fs[0],nfe);
1232  }
1233  else
1234  {
1235  fs[0] = *fx;
1236  }
1237  isubc.IntNew = true;
1238  for( j=2; j <= npts; j++ )
1239  {
1240  j_ = j - 1;
1241  evalf(ns,ips,&S(j_,0),n,x,&fs[j_],nfe);
1242  }
1243  il = 1;
1244  order(npts,fs,&il,&is,&ih);
1245  tol = (realnum)(usubc.psi*dist(ns,&S(ih-1,0),&S(il-1,0)));
1246  }
1247 
1248  /* main loop
1249  * */
1250 L_20:
1251  calcc(ns,s,ih,inew,updatc,&S(icent-1,0));
1252  updatc = true;
1253  inew = ih;
1254 
1255  /* reflect
1256  * */
1257  newpt(ns,usubc.alpha,&S(icent-1,0),&S(ih-1,0),true,&S(itemp-1,0),
1258  &small);
1259  if( !small )
1260  {
1261  evalf(ns,ips,&S(itemp-1,0),n,x,&fr,nfe);
1262  if( fr < fs[il-1] )
1263  {
1264 
1265  /* expand
1266  * */
1267  newpt(ns,-usubc.gamma,&S(icent-1,0),&S(itemp-1,0),true,
1268  &S(ih-1,0),&small);
1269  if( small )
1270  goto L_40;
1271  evalf(ns,ips,&S(ih-1,0),n,x,&fe,nfe);
1272  if( fe < fr )
1273  {
1274  fs[ih-1] = fe;
1275  }
1276  else
1277  {
1278  cdcopy(ns,&S(itemp-1,0),1,&S(ih-1,0),1);
1279  fs[ih-1] = fr;
1280  }
1281  }
1282  else if( fr < fs[is-1] )
1283  {
1284 
1285  /* accept reflected point
1286  * */
1287  cdcopy(ns,&S(itemp-1,0),1,&S(ih-1,0),1);
1288  fs[ih-1] = fr;
1289  }
1290  else
1291  {
1292 
1293  /* contract
1294  * */
1295  if( fr > fs[ih-1] )
1296  {
1297  newpt(ns,-usubc.beta,&S(icent-1,0),&S(ih-1,0),true,
1298  &S(itemp-1,0),&small);
1299  }
1300  else
1301  {
1302  newpt(ns,-usubc.beta,&S(icent-1,0),&S(itemp-1,0),false,
1303  (realnum*)&dum,&small);
1304  }
1305  if( small )
1306  goto L_40;
1307  evalf(ns,ips,&S(itemp-1,0),n,x,&fc,nfe);
1308  if( fc < (realnum)MIN2(fr,fs[ih-1]) )
1309  {
1310  cdcopy(ns,&S(itemp-1,0),1,&S(ih-1,0),1);
1311  fs[ih-1] = fc;
1312  }
1313  else
1314  {
1315 
1316  /* shrink simplex
1317  * */
1318  for( j=1; j <= npts; j++ )
1319  {
1320  j_ = j - 1;
1321  if( j != il )
1322  {
1323  newpt(ns,-usubc.delta,&S(il-1,0),&S(j_,0),
1324  false,(realnum*)&dum,&small);
1325  if( small )
1326  goto L_40;
1327  evalf(ns,ips,&S(j_,0),n,x,&fs[j_],nfe);
1328  }
1329  }
1330  }
1331  updatc = false;
1332  }
1333  order(npts,fs,&il,&is,&ih);
1334  }
1335 
1336  /* check termination
1337  * */
1338 
1339 L_40:
1340  if( usubc.irepl == 0 )
1341  {
1342  *fx = fs[il-1];
1343  }
1344  else
1345  {
1346  *fx = isubc.sfbest;
1347  }
1348 
1349 L_50:
1350  if( (usubc.nfstop > 0 && *fx <= isubc.sfstop) && usubc.nfxe >=
1351  usubc.nfstop )
1352  goto L_51;
1353  if( *nfe >= maxnfe )
1354  goto L_52;
1355  if( !(dist(ns,&S(ih-1,0),&S(il-1,0)) <= tol || small) )
1356  goto L_20;
1357  *iflag = 0;
1358  goto L_53;
1359 
1360 L_52:
1361  *iflag = -1;
1362  goto L_53;
1363 
1364 L_51:
1365  *iflag = 2;
1366 
1367  /* end main loop, return best point
1368  * */
1369 
1370 L_53:
1371  for( i=0; i < ns; i++ )
1372  {
1373  x[ips[i]-1] = S(il-1,i);
1374  }
1375  return;
1376 }
1377 
1378 /********************************************************************* */
1379 STATIC void fstats(double fx,
1380  long int ifxwt,
1381  int reset)
1382 {
1383  static long int nsv;
1384  static realnum f1sv,
1385  fscale;
1386 
1387  DEBUG_ENTRY( "fstats()" );
1388 
1389 
1390  /* Coded by Tom Rowan
1391  * Department of Computer Sciences
1392  * University of Texas at Austin
1393  *
1394  * fstats modifies the common /usubc/ variables nfxe,fxstat.
1395  *
1396  * input
1397  *
1398  * fx - most recent evaluation of f at best x
1399  *
1400  * ifxwt - integer weight for fx
1401  *
1402  * reset - logical switch
1403  * = .true. : initialize nfxe,fxstat
1404  * = .false. : update nfxe,fxstat
1405  *
1406  * common
1407  * */
1408 
1409 
1410  /* local variables
1411  * */
1412 
1413 
1414  /* subroutines and functions
1415  *
1416  *
1417  *-----------------------------------------------------------
1418  * */
1419  if( reset )
1420  {
1421  usubc.nfxe = ifxwt;
1422  usubc.fxstat[0] = (realnum)fx;
1423  usubc.fxstat[1] = (realnum)fx;
1424  usubc.fxstat[2] = (realnum)fx;
1425  usubc.fxstat[3] = 0.0f;
1426  }
1427  else
1428  {
1429  nsv = usubc.nfxe;
1430  f1sv = usubc.fxstat[0];
1431  usubc.nfxe += ifxwt;
1432  usubc.fxstat[0] += (realnum)(ifxwt*(fx - usubc.fxstat[0])/usubc.nfxe);
1433  usubc.fxstat[1] = MAX2(usubc.fxstat[1],(realnum)fx);
1434  usubc.fxstat[2] = MIN2(usubc.fxstat[2],(realnum)fx);
1435  fscale = (realnum)MAX3(fabs(usubc.fxstat[1]),fabs(usubc.fxstat[2]),1.);
1436  usubc.fxstat[3] = (realnum)(fscale*sqrt(((nsv-1)*POW2(usubc.fxstat[3]/
1437  fscale)+nsv*POW2((usubc.fxstat[0]-f1sv)/fscale)+ifxwt*
1438  POW2((fx-usubc.fxstat[0])/fscale))/(usubc.nfxe-1)));
1439  }
1440  return;
1441 }
1442 
1443 STATIC double cdasum(long int n,
1444  realnum dx[],
1445  long int incx)
1446 {
1447  /*
1448  *
1449  * takes the sum of the absolute values.
1450  * uses unrolled loops for increment equal to one.
1451  * jack dongarra, lapack, 3/11/78.
1452  * modified to correct problem with negative increment, 8/21/90.
1453  *
1454  */
1455  long int i,
1456  ix,
1457  m;
1458  realnum cdasum_v,
1459  dtemp;
1460 
1461  DEBUG_ENTRY( "cdasum()" );
1462 
1463  cdasum_v = 0.00f;
1464  dtemp = 0.00f;
1465  if( n > 0 )
1466  {
1467  if( incx == 1 )
1468  {
1469 
1470  /* code for increment equal to 1
1471  *
1472  *
1473  * clean-up loop
1474  * */
1475  m = n%6;
1476  if( m != 0 )
1477  {
1478  for( i=0; i < m; i++ )
1479  {
1480  dtemp += (realnum)fabs(dx[i]);
1481  }
1482  if( n < 6 )
1483  goto L_60;
1484  }
1485 
1486  for( i=m; i < n; i += 6 )
1487  {
1488  dtemp += (realnum)(fabs(dx[i]) + fabs(dx[i+1]) + fabs(dx[i+2]) +
1489  fabs(dx[i+3]) + fabs(dx[i+4]) + fabs(dx[i+5]));
1490  }
1491 L_60:
1492  cdasum_v = dtemp;
1493  }
1494  else
1495  {
1496 
1497  /* code for increment not equal to 1
1498  * */
1499  ix = 1;
1500 
1501  if( incx < 0 )
1502  ix = (-n + 1)*incx + 1;
1503 
1504  for( i=0; i < n; i++ )
1505  {
1506  dtemp += (realnum)fabs(dx[ix-1]);
1507  ix += incx;
1508  }
1509  cdasum_v = dtemp;
1510  }
1511  }
1512  return( cdasum_v );
1513 }
1514 
1515 STATIC void csscal(long int n,
1516  double da,
1517  realnum dx[],
1518  long int incx)
1519 {
1520  long int
1521  i,
1522  i_,
1523  m,
1524  mp1,
1525  nincx;
1526 
1527  DEBUG_ENTRY( "csscal()" );
1528 
1529  /* scales a vector by a constant.
1530  * uses unrolled loops for increment equal to one.
1531  * jack dongarra, lapack, 3/11/78.
1532  * modified 3/93 to return if incx .le. 0.
1533  * modified 12/3/93, array(1) declarations changed to array(*)
1534  * changed to single precisions
1535  * */
1536 
1537  if( !(n <= 0 || incx <= 0) )
1538  {
1539  if( incx == 1 )
1540  {
1541 
1542  /* code for increment equal to 1
1543  *
1544  *
1545  * clean-up loop
1546  * */
1547  m = n%5;
1548  if( m != 0 )
1549  {
1550  for( i=1; i <= m; i++ )
1551  {
1552  i_ = i - 1;
1553  dx[i_] *= (realnum)da;
1554  }
1555  if( n < 5 )
1556  {
1557  return;
1558  }
1559  }
1560  mp1 = m + 1;
1561  for( i=mp1; i <= n; i += 5 )
1562  {
1563  i_ = i - 1;
1564  dx[i_] *= (realnum)da;
1565  dx[i_+1] *= (realnum)da;
1566  dx[i_+2] *= (realnum)da;
1567  dx[i_+3] *= (realnum)da;
1568  dx[i_+4] *= (realnum)da;
1569  }
1570  }
1571  else
1572  {
1573 
1574  /* code for increment not equal to 1
1575  * */
1576  nincx = n*incx;
1577  /*for( i=1, _do0=DOCNT(i,nincx,_do1 = incx); _do0 > 0; i += _do1, _do0-- )*/
1578  /* gary change forc */
1579  for( i=0; i<nincx; i=i+incx)
1580  /*for( i=1, _do0=DOCNT(i,nincx,_do1 = incx); _do0 > 0; i += _do1, _do0-- )*/
1581  {
1582  dx[i] *= (realnum)da;
1583  }
1584  }
1585  }
1586  return;
1587 }
1588 
1589 /********************************************************************* */
1590 
1591 #undef S
1592 #define S(I_,J_) (*(s+(I_)*(ns)+(J_)))
1593 
1594 STATIC void start(long int,
1595  realnum x[],
1596  realnum step[],
1597  long int ns,
1598  long int ips[],
1599  realnum *s,
1600  int *small)
1601 {
1602  DEBUG_ENTRY( "start()" );
1603 
1604  /* Coded by Tom Rowan
1605  * Department of Computer Sciences
1606  * University of Texas at Austin
1607  *
1608  * start creates the initial simplex for simplx minimization.
1609  *
1610  * input
1611  *
1612  * n - problem dimension (not used)
1613  *
1614  * x - current best point
1615  *
1616  * step - stepsizes for corresponding components of x
1617  *
1618  * ns - subspace dimension
1619  *
1620  * ips - permutation vector
1621  *
1622  *
1623  * output
1624  *
1625  * s - first ns+1 columns contain initial simplex
1626  *
1627  * small - logical flag
1628  * = .true. : coincident points
1629  * = .false. : otherwise
1630  *
1631  * local variables
1632  * */
1633 
1634  /* subroutines and functions
1635  *
1636  * blas */
1637  /* fortran */
1638 
1639  /*-----------------------------------------------------------
1640  * */
1641  for( long i=1; i <= ns; i++ )
1642  {
1643  long i_ = i - 1;
1644  S(0,i_) = x[ips[i_]-1];
1645  }
1646 
1647  for( long j=2; j <= (ns + 1); j++ )
1648  {
1649  long j_ = j - 1;
1650  cdcopy(ns,&S(0,0),1,&S(j_,0),1);
1651  S(j_,j_-1) = S(0,j_-1) + step[ips[j_-1]-1];
1652  }
1653 
1654  /* check for coincident points
1655  * */
1656  for( long j=2; j <= (ns + 1); j++ )
1657  {
1658  long j_ = j - 1;
1659  if( (double)(S(j_,j_-1)) == (double)(S(0,j_-1)) )
1660  goto L_40;
1661  }
1662  *small = false;
1663  return;
1664 
1665  /* coincident points
1666  * */
1667 L_40:
1668  *small = true;
1669  return;
1670 }
1671 
1672 /********************************************************************* */
1673 STATIC void order(long int npts,
1674  realnum fs[],
1675  long int *il,
1676  long int *is,
1677  long int *ih)
1678 {
1679  long int i,
1680  il0,
1681  j;
1682 
1683  DEBUG_ENTRY( "order()" );
1684 
1685 
1686  /* Coded by Tom Rowan
1687  * Department of Computer Sciences
1688  * University of Texas at Austin
1689  *
1690  * order determines the indices of the vertices with the
1691  * lowest, second highest, and highest function values.
1692  *
1693  * input
1694  *
1695  * npts - number of points in simplex
1696  *
1697  * fs - real vector of function values of
1698  * simplex
1699  *
1700  * il - index to vertex with lowest function value
1701  *
1702  * output
1703  *
1704  * il - new index to vertex with lowest function value
1705  *
1706  * is - new index to vertex with second highest
1707  * function value
1708  *
1709  * ih - new index to vertex with highest function value
1710  *
1711  * local variables
1712  * */
1713 
1714  /* subroutines and functions
1715  *
1716  *
1717  *-----------------------------------------------------------
1718  * */
1719  il0 = *il;
1720  j = (il0%npts) + 1;
1721 
1722  if( fs[j-1] >= fs[*il-1] )
1723  {
1724  *ih = j;
1725  *is = il0;
1726  }
1727  else
1728  {
1729  *ih = il0;
1730  *is = j;
1731  *il = j;
1732  }
1733 
1734  for( i=il0 + 1; i <= (il0 + npts - 2); i++ )
1735  {
1736  j = (i%npts) + 1;
1737  if( fs[j-1] >= fs[*ih-1] )
1738  {
1739  *is = *ih;
1740  *ih = j;
1741  }
1742  else if( fs[j-1] > fs[*is-1] )
1743  {
1744  *is = j;
1745  }
1746  else if( fs[j-1] < fs[*il-1] )
1747  {
1748  *il = j;
1749  }
1750  }
1751  return;
1752 }
1753 
1754 STATIC double dist(long int n,
1755  realnum x[],
1756  realnum y[])
1757 {
1758  /*
1759  *
1760  */
1761  long int i,
1762  i_;
1763  realnum absxmy,
1764  dist_v,
1765  scale,
1766  sum;
1767 
1768  DEBUG_ENTRY( "dist()" );
1769 
1770  /* Coded by Tom Rowan
1771  * Department of Computer Sciences
1772  * University of Texas at Austin
1773  *
1774  * dist calculates the distance between the points x,y.
1775  *
1776  * input
1777  *
1778  * n - number of components
1779  *
1780  * x - point in n-space
1781  *
1782  * y - point in n-space
1783  *
1784  * local variables
1785  * */
1786 
1787  /* subroutines and functions
1788  *
1789  * fortran
1790  *
1791  *-----------------------------------------------------------
1792  * */
1793  absxmy = (realnum)fabs(x[0]-y[0]);
1794  if( absxmy <= 1.0f )
1795  {
1796  sum = absxmy*absxmy;
1797  scale = 1.0f;
1798  }
1799  else
1800  {
1801  sum = 1.0f;
1802  scale = absxmy;
1803  }
1804 
1805  for( i=2; i <= n; i++ )
1806  {
1807  i_ = i - 1;
1808  absxmy = (realnum)fabs(x[i_]-y[i_]);
1809  if( absxmy <= scale )
1810  {
1811  sum += POW2(absxmy/scale);
1812  }
1813  else
1814  {
1815  sum = 1.0f + sum*POW2(scale/absxmy);
1816  scale = absxmy;
1817  }
1818  }
1819  dist_v = (realnum)(scale*sqrt(sum));
1820  return( dist_v );
1821 }
1822 
1823 /********************************************************************* */
1824 
1825 #undef S
1826 #define S(I_,J_) (*(s+(I_)*(ns)+(J_)))
1827 
1828 STATIC void calcc(long int ns,
1829  realnum *s,
1830  long int ih,
1831  long int inew,
1832  int updatc,
1833  realnum c[])
1834 {
1835  long int i,
1836  j,
1837  j_;
1838  /* >>chng 99 apr 21, following was not initialized, caught by Peter van Hoof */
1839  realnum xNothing[1] = { 0.0f };
1840 
1841  DEBUG_ENTRY( "calcc()" );
1842 
1843 
1844  /* Coded by Tom Rowan
1845  * Department of Computer Sciences
1846  * University of Texas at Austin
1847  *
1848  * calcc calculates the centroid of the simplex without the
1849  * vertex with highest function value.
1850  *
1851  * input
1852  *
1853  * ns - subspace dimension
1854  *
1855  * s - real work space of dimension .ge.
1856  * ns*(ns+3) used to store simplex
1857  *
1858  * ih - index to vertex with highest function value
1859  *
1860  * inew - index to new point
1861  *
1862  * updatc - logical switch
1863  * = .true. : update centroid
1864  * = .false. : calculate centroid from scratch
1865  *
1866  * c - centroid of the simplex without vertex with
1867  * highest function value
1868  *
1869  * output
1870  *
1871  * c - new centroid
1872  *
1873  * local variables
1874  * */
1875  /* added to get prototypes to work */
1876 
1877  /* subroutines and functions
1878  *
1879  * blas */
1880 
1881  /*-----------------------------------------------------------
1882  * */
1883  if( !updatc )
1884  {
1885  /* call cdcopy (ns,0.0,0,c,1)
1886  * xNothing will not be used since 0 is increment */
1887  cdcopy(ns,xNothing,0,c,1);
1888  for( j=1; j <= (ns + 1); j++ )
1889  {
1890  j_ = j - 1;
1891  if( j != ih )
1892  cdaxpy(ns,1.0f,&S(j_,0),1,c,1);
1893  }
1894  csscal(ns,1.0f/ns,c,1);
1895  }
1896  else if( ih != inew )
1897  {
1898  for( i=0; i < ns; i++ )
1899  {
1900  c[i] += (S(inew-1,i) - S(ih-1,i))/ns;
1901  }
1902  }
1903  return;
1904 }
1905 
1906 /********************************************************************* */
1907 STATIC void newpt(long int ns,
1908  double coef,
1909  realnum xbase[],
1910  realnum xold[],
1911  int IntNew,
1912  realnum xnew[],
1913  int *small)
1914 {
1915  int eqbase,
1916  eqold;
1917  long int i;
1918  realnum xoldi;
1919 
1920  DEBUG_ENTRY( "newpt()" );
1921 
1922 
1923  /* Coded by Tom Rowan
1924  * Department of Computer Sciences
1925  * University of Texas at Austin
1926  *
1927  * newpt performs reflections, expansions, contractions, and
1928  * shrinkages (massive contractions) by computing:
1929  *
1930  * xbase + coef * (xbase - xold)
1931  *
1932  * The result is stored in xnew if IntNew .eq. .true.,
1933  * in xold otherwise.
1934  *
1935  * use : coef .gt. 0 to reflect
1936  * coef .lt. 0 to expand, contract, or shrink
1937  *
1938  * input
1939  *
1940  * ns - number of components (subspace dimension)
1941  *
1942  * coef - one of four simplex method coefficients
1943  *
1944  * xbase - real ns-vector representing base
1945  * point
1946  *
1947  * xold - real ns-vector representing old
1948  * point
1949  *
1950  * IntNew - logical switch
1951  * = .true. : store result in xnew
1952  * = .false. : store result in xold, xnew is not
1953  * referenced
1954  *
1955  * output
1956  *
1957  * xold - unchanged if IntNew .eq. .true., contains new
1958  * point otherwise
1959  *
1960  * xnew - real ns-vector representing new
1961  * point if new .eq. .true., not referenced
1962  * otherwise
1963  *
1964  * small - logical flag
1965  * = .true. : coincident points
1966  * = .false. : otherwise
1967  *
1968  * local variables
1969  * */
1970 
1971  /* subroutines and functions
1972  *
1973  * fortran */
1974 
1975  /*-----------------------------------------------------------
1976  * */
1977  eqbase = true;
1978  eqold = true;
1979  if( IntNew )
1980  {
1981  for( i=0; i < ns; i++ )
1982  {
1983  xnew[i] = (realnum)(xbase[i] + coef*(xbase[i] - xold[i]));
1984  eqbase = eqbase && ((double)(xnew[i]) == (double)(xbase[i]));
1985  eqold = eqold && ((double)(xnew[i]) == (double)(xold[i]));
1986  }
1987  }
1988  else
1989  {
1990  for( i=0; i < ns; i++ )
1991  {
1992  xoldi = xold[i];
1993  xold[i] = (realnum)(xbase[i] + coef*(xbase[i] - xold[i]));
1994  eqbase = eqbase && ((double)(xold[i]) == (double)(xbase[i]));
1995  eqold = eqold && ((double)(xold[i]) == (double)(xoldi));
1996  }
1997  }
1998  *small = eqbase || eqold;
1999  return;
2000 }
2001 /********************************************************************* */
2002 STATIC void cdaxpy(long int n,
2003  double da,
2004  realnum dx[],
2005  long int incx,
2006  realnum dy[],
2007  long int incy)
2008 {
2009  long int i,
2010  i_,
2011  ix,
2012  iy,
2013  m;
2014 
2015  DEBUG_ENTRY( "cdaxpy()" );
2016 
2017  /* constant times a vector plus a vector.
2018  * uses unrolled loops for increments equal to one.
2019  * jack dongarra, lapack, 3/11/78.
2020  * */
2021 
2022  if( n > 0 )
2023  {
2024  if( da != 0.00f )
2025  {
2026  if( incx == 1 && incy == 1 )
2027  {
2028 
2029  /* code for both increments equal to 1
2030  *
2031  *
2032  * clean-up loop
2033  * */
2034  m = n%4;
2035  if( m != 0 )
2036  {
2037  for( i=1; i <= m; i++ )
2038  {
2039  i_ = i - 1;
2040  dy[i_] += (realnum)(da*dx[i_]);
2041  }
2042  if( n < 4 )
2043  {
2044  return;
2045  }
2046  }
2047 
2048  for( i=m; i < n; i += 4 )
2049  {
2050  dy[i] += (realnum)(da*dx[i]);
2051  dy[i+1] += (realnum)(da*dx[i+1]);
2052  dy[i+2] += (realnum)(da*dx[i+2]);
2053  dy[i+3] += (realnum)(da*dx[i+3]);
2054  }
2055  }
2056  else
2057  {
2058 
2059  /* code for unequal increments or equal increments
2060  * not equal to 1
2061  * */
2062  ix = 1;
2063  iy = 1;
2064  if( incx < 0 )
2065  ix = (-n + 1)*incx + 1;
2066  if( incy < 0 )
2067  iy = (-n + 1)*incy + 1;
2068  for( i=0; i < n; i++ )
2069  {
2070  dy[iy-1] += (realnum)(da*dx[ix-1]);
2071  ix += incx;
2072  iy += incy;
2073  }
2074  }
2075  }
2076  }
2077  return;
2078 }
2079 /*lint +e801 goto is deprecated */
2080 /*lint +e64 type mismatch */
STATIC void partx(long, long[], realnum[], long *, long[])
realnum gamma
STATIC void setstp(long, long, realnum[], realnum[])
realnum ftest
STATIC double cdasum(long, realnum[], long)
realnum delta
long int nsmax
t_fe fe
Definition: fe.cpp:5
long int nfstop
STATIC void subopt(long)
realnum psi
chi2_type optimize_func(const realnum param[], int grid_index=-1)
T sign(T x, T y)
Definition: cddefines.h:842
STATIC void fstats(double, long, int)
STATIC void order(long, realnum[], long *, long *, long *)
realnum alpha
#define MIN2(a, b)
Definition: cddefines.h:803
void optimize_subplex(long int n, double tol, long int maxnfe, long int mode, realnum scale[], realnum x[], realnum *fx, long int *nfe, realnum work[], long int iwork[], long int *iflag)
realnum sfbest
STATIC void newpt(long, double, realnum[], realnum[], int, realnum[], int *)
STATIC void calcc(long, realnum *, long, long, int, realnum[])
long int irepl
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:854
struct t_isubc isubc
#define POW2
Definition: cddefines.h:979
long int nsmin
#define STATIC
Definition: cddefines.h:118
STATIC void cdaxpy(long, double, realnum[], long, realnum[], long)
float realnum
Definition: cddefines.h:124
STATIC double da(double z, double temp, double eden)
STATIC void cdcopy(long, realnum[], long, realnum[], long)
STATIC void sortd(long, realnum[], long[])
long int nfxe
STATIC void csscal(long, double, realnum[], long)
STATIC double dist(long, realnum[], realnum[])
realnum sfstop
realnum fstop
struct t_usubc usubc
STATIC void simplx(long, realnum[], long, long[], long, int, realnum[], realnum *, long *, realnum *, realnum[], long *)
STATIC void start(long, realnum[], realnum[], long, long[], realnum *, int *)
long int ifxsw
realnum bonus
static int nleft
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
realnum fxstat[4]
#define MAX2(a, b)
Definition: cddefines.h:824
#define S(I_, J_)
realnum omega
STATIC void evalf(long, long[], realnum[], long, realnum[], realnum *, long *)
#define MAX3(a, b, c)
Definition: cddefines.h:829
realnum fbonus
realnum beta