cloudy trunk
Loading...
Searching...
No Matches
optimize_subplx.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 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
6STATIC void calcc(long,realnum*,long,long,int,realnum[]);
7STATIC double cdasum(long,realnum[],long);
8STATIC void cdaxpy(long,double,realnum[],long,realnum[],long);
9STATIC void cdcopy(long,realnum[],long,realnum[],long);
10STATIC void csscal(long,double,realnum[],long);
11STATIC double dist(long,realnum[],realnum[]);
12STATIC void evalf(long,long[],realnum[],long,realnum[],realnum*,long*);
13STATIC void fstats(double,long,int);
14STATIC void newpt(long,double,realnum[],realnum[],int,realnum[],int*);
15STATIC void order(long,realnum[],long*,long*,long*);
16STATIC void partx(long,long[],realnum[],long*,long[]);
17STATIC void setstp(long,long,realnum[],realnum[]);
18STATIC void simplx(long,realnum[],long,long[],long,int,realnum[],
19 realnum*,long*,realnum*,realnum[],long*);
20STATIC void sortd(long,realnum[],long[]);
21STATIC void start(long,realnum[],realnum[],long,long[],realnum*,int*);
22STATIC 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 *********************************************************************** */
63
64void 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 {
311 isubc.sfstop = usubc.fstop;
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 {
332 isubc.sfstop = usubc.fstop;
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 * */
358L_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
375L_60:
376 ns = iwork[ins-1];
377
378L_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 * */
402L_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;
430L_121:
431 if( usubc.minf )
432 {
433 *fx = sfx;
434 }
435 else
436 {
437 *fx = -sfx;
438 }
439 return;
440
441 /* invalid input
442 * */
443L_120:
444
445 *iflag = -2;
446 return;
447}
448/********************************************************************* */
449STATIC 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 */
622STATIC 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/********************************************************************* */
696STATIC 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/********************************************************************* */
825STATIC 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/********************************************************************* */
908STATIC 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/********************************************************************* */
989STATIC 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;
1097L_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
1107STATIC 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 * */
1250L_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
1339L_40:
1340 if( usubc.irepl == 0 )
1341 {
1342 *fx = fs[il-1];
1343 }
1344 else
1345 {
1346 *fx = isubc.sfbest;
1347 }
1348
1349L_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
1360L_52:
1361 *iflag = -1;
1362 goto L_53;
1363
1364L_51:
1365 *iflag = 2;
1366
1367 /* end main loop, return best point
1368 * */
1369
1370L_53:
1371 for( i=0; i < ns; i++ )
1372 {
1373 x[ips[i]-1] = S(il-1,i);
1374 }
1375 return;
1376}
1377
1378/********************************************************************* */
1379STATIC 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
1443STATIC 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 }
1491L_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
1515STATIC 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
1594STATIC void start(long int n,
1595 realnum x[],
1596 realnum step[],
1597 long int ns,
1598 long int ips[],
1599 realnum *s,
1600 int *small)
1601{
1602 long int i,
1603 i_,
1604 j,
1605 j_;
1606
1607 DEBUG_ENTRY( "start()" );
1608
1609 /* gary add, not used so set to i_ to mute compiler note */
1610 i_ = n;
1611
1612
1613 /* Coded by Tom Rowan
1614 * Department of Computer Sciences
1615 * University of Texas at Austin
1616 *
1617 * start creates the initial simplex for simplx minimization.
1618 *
1619 * input
1620 *
1621 * n - problem dimension
1622 *
1623 * x - current best point
1624 *
1625 * step - stepsizes for corresponding components of x
1626 *
1627 * ns - subspace dimension
1628 *
1629 * ips - permutation vector
1630 *
1631 *
1632 * output
1633 *
1634 * s - first ns+1 columns contain initial simplex
1635 *
1636 * small - logical flag
1637 * = .true. : coincident points
1638 * = .false. : otherwise
1639 *
1640 * local variables
1641 * */
1642
1643 /* subroutines and functions
1644 *
1645 * blas */
1646 /* fortran */
1647
1648 /*-----------------------------------------------------------
1649 * */
1650 for( i=1; i <= ns; i++ )
1651 {
1652 i_ = i - 1;
1653 S(0,i_) = x[ips[i_]-1];
1654 }
1655
1656 for( j=2; j <= (ns + 1); j++ )
1657 {
1658 j_ = j - 1;
1659 cdcopy(ns,&S(0,0),1,&S(j_,0),1);
1660 S(j_,j_-1) = S(0,j_-1) + step[ips[j_-1]-1];
1661 }
1662
1663 /* check for coincident points
1664 * */
1665 for( j=2; j <= (ns + 1); j++ )
1666 {
1667 j_ = j - 1;
1668 if( (double)(S(j_,j_-1)) == (double)(S(0,j_-1)) )
1669 goto L_40;
1670 }
1671 *small = false;
1672 return;
1673
1674 /* coincident points
1675 * */
1676L_40:
1677 *small = true;
1678 return;
1679}
1680
1681/********************************************************************* */
1682STATIC void order(long int npts,
1683 realnum fs[],
1684 long int *il,
1685 long int *is,
1686 long int *ih)
1687{
1688 long int i,
1689 il0,
1690 j;
1691
1692 DEBUG_ENTRY( "order()" );
1693
1694
1695 /* Coded by Tom Rowan
1696 * Department of Computer Sciences
1697 * University of Texas at Austin
1698 *
1699 * order determines the indices of the vertices with the
1700 * lowest, second highest, and highest function values.
1701 *
1702 * input
1703 *
1704 * npts - number of points in simplex
1705 *
1706 * fs - real vector of function values of
1707 * simplex
1708 *
1709 * il - index to vertex with lowest function value
1710 *
1711 * output
1712 *
1713 * il - new index to vertex with lowest function value
1714 *
1715 * is - new index to vertex with second highest
1716 * function value
1717 *
1718 * ih - new index to vertex with highest function value
1719 *
1720 * local variables
1721 * */
1722
1723 /* subroutines and functions
1724 *
1725 *
1726 *-----------------------------------------------------------
1727 * */
1728 il0 = *il;
1729 j = (il0%npts) + 1;
1730
1731 if( fs[j-1] >= fs[*il-1] )
1732 {
1733 *ih = j;
1734 *is = il0;
1735 }
1736 else
1737 {
1738 *ih = il0;
1739 *is = j;
1740 *il = j;
1741 }
1742
1743 for( i=il0 + 1; i <= (il0 + npts - 2); i++ )
1744 {
1745 j = (i%npts) + 1;
1746 if( fs[j-1] >= fs[*ih-1] )
1747 {
1748 *is = *ih;
1749 *ih = j;
1750 }
1751 else if( fs[j-1] > fs[*is-1] )
1752 {
1753 *is = j;
1754 }
1755 else if( fs[j-1] < fs[*il-1] )
1756 {
1757 *il = j;
1758 }
1759 }
1760 return;
1761}
1762
1763STATIC double dist(long int n,
1764 realnum x[],
1765 realnum y[])
1766{
1767 /*
1768 *
1769 */
1770 long int i,
1771 i_;
1772 realnum absxmy,
1773 dist_v,
1774 scale,
1775 sum;
1776
1777 DEBUG_ENTRY( "dist()" );
1778
1779 /* Coded by Tom Rowan
1780 * Department of Computer Sciences
1781 * University of Texas at Austin
1782 *
1783 * dist calculates the distance between the points x,y.
1784 *
1785 * input
1786 *
1787 * n - number of components
1788 *
1789 * x - point in n-space
1790 *
1791 * y - point in n-space
1792 *
1793 * local variables
1794 * */
1795
1796 /* subroutines and functions
1797 *
1798 * fortran
1799 *
1800 *-----------------------------------------------------------
1801 * */
1802 absxmy = (realnum)fabs(x[0]-y[0]);
1803 if( absxmy <= 1.0f )
1804 {
1805 sum = absxmy*absxmy;
1806 scale = 1.0f;
1807 }
1808 else
1809 {
1810 sum = 1.0f;
1811 scale = absxmy;
1812 }
1813
1814 for( i=2; i <= n; i++ )
1815 {
1816 i_ = i - 1;
1817 absxmy = (realnum)fabs(x[i_]-y[i_]);
1818 if( absxmy <= scale )
1819 {
1820 sum += POW2(absxmy/scale);
1821 }
1822 else
1823 {
1824 sum = 1.0f + sum*POW2(scale/absxmy);
1825 scale = absxmy;
1826 }
1827 }
1828 dist_v = (realnum)(scale*sqrt(sum));
1829 return( dist_v );
1830}
1831
1832/********************************************************************* */
1833
1834#undef S
1835#define S(I_,J_) (*(s+(I_)*(ns)+(J_)))
1836
1837STATIC void calcc(long int ns,
1838 realnum *s,
1839 long int ih,
1840 long int inew,
1841 int updatc,
1842 realnum c[])
1843{
1844 long int i,
1845 j,
1846 j_;
1847 /* >>chng 99 apr 21, following was not initialized, caught by Peter van Hoof */
1848 realnum xNothing[1] = { 0.0f };
1849
1850 DEBUG_ENTRY( "calcc()" );
1851
1852
1853 /* Coded by Tom Rowan
1854 * Department of Computer Sciences
1855 * University of Texas at Austin
1856 *
1857 * calcc calculates the centroid of the simplex without the
1858 * vertex with highest function value.
1859 *
1860 * input
1861 *
1862 * ns - subspace dimension
1863 *
1864 * s - real work space of dimension .ge.
1865 * ns*(ns+3) used to store simplex
1866 *
1867 * ih - index to vertex with highest function value
1868 *
1869 * inew - index to new point
1870 *
1871 * updatc - logical switch
1872 * = .true. : update centroid
1873 * = .false. : calculate centroid from scratch
1874 *
1875 * c - centroid of the simplex without vertex with
1876 * highest function value
1877 *
1878 * output
1879 *
1880 * c - new centroid
1881 *
1882 * local variables
1883 * */
1884 /* added to get prototypes to work */
1885
1886 /* subroutines and functions
1887 *
1888 * blas */
1889
1890 /*-----------------------------------------------------------
1891 * */
1892 if( !updatc )
1893 {
1894 /* call cdcopy (ns,0.0,0,c,1)
1895 * xNothing will not be used since 0 is increment */
1896 cdcopy(ns,xNothing,0,c,1);
1897 for( j=1; j <= (ns + 1); j++ )
1898 {
1899 j_ = j - 1;
1900 if( j != ih )
1901 cdaxpy(ns,1.0f,&S(j_,0),1,c,1);
1902 }
1903 csscal(ns,1.0f/ns,c,1);
1904 }
1905 else if( ih != inew )
1906 {
1907 for( i=0; i < ns; i++ )
1908 {
1909 c[i] += (S(inew-1,i) - S(ih-1,i))/ns;
1910 }
1911 }
1912 return;
1913}
1914
1915/********************************************************************* */
1916STATIC void newpt(long int ns,
1917 double coef,
1918 realnum xbase[],
1919 realnum xold[],
1920 int IntNew,
1921 realnum xnew[],
1922 int *small)
1923{
1924 int eqbase,
1925 eqold;
1926 long int i;
1927 realnum xoldi;
1928
1929 DEBUG_ENTRY( "newpt()" );
1930
1931
1932 /* Coded by Tom Rowan
1933 * Department of Computer Sciences
1934 * University of Texas at Austin
1935 *
1936 * newpt performs reflections, expansions, contractions, and
1937 * shrinkages (massive contractions) by computing:
1938 *
1939 * xbase + coef * (xbase - xold)
1940 *
1941 * The result is stored in xnew if IntNew .eq. .true.,
1942 * in xold otherwise.
1943 *
1944 * use : coef .gt. 0 to reflect
1945 * coef .lt. 0 to expand, contract, or shrink
1946 *
1947 * input
1948 *
1949 * ns - number of components (subspace dimension)
1950 *
1951 * coef - one of four simplex method coefficients
1952 *
1953 * xbase - real ns-vector representing base
1954 * point
1955 *
1956 * xold - real ns-vector representing old
1957 * point
1958 *
1959 * IntNew - logical switch
1960 * = .true. : store result in xnew
1961 * = .false. : store result in xold, xnew is not
1962 * referenced
1963 *
1964 * output
1965 *
1966 * xold - unchanged if IntNew .eq. .true., contains new
1967 * point otherwise
1968 *
1969 * xnew - real ns-vector representing new
1970 * point if new .eq. .true., not referenced
1971 * otherwise
1972 *
1973 * small - logical flag
1974 * = .true. : coincident points
1975 * = .false. : otherwise
1976 *
1977 * local variables
1978 * */
1979
1980 /* subroutines and functions
1981 *
1982 * fortran */
1983
1984 /*-----------------------------------------------------------
1985 * */
1986 eqbase = true;
1987 eqold = true;
1988 if( IntNew )
1989 {
1990 for( i=0; i < ns; i++ )
1991 {
1992 xnew[i] = (realnum)(xbase[i] + coef*(xbase[i] - xold[i]));
1993 eqbase = eqbase && ((double)(xnew[i]) == (double)(xbase[i]));
1994 eqold = eqold && ((double)(xnew[i]) == (double)(xold[i]));
1995 }
1996 }
1997 else
1998 {
1999 for( i=0; i < ns; i++ )
2000 {
2001 xoldi = xold[i];
2002 xold[i] = (realnum)(xbase[i] + coef*(xbase[i] - xold[i]));
2003 eqbase = eqbase && ((double)(xold[i]) == (double)(xbase[i]));
2004 eqold = eqold && ((double)(xold[i]) == (double)(xoldi));
2005 }
2006 }
2007 *small = eqbase || eqold;
2008 return;
2009}
2010/********************************************************************* */
2011STATIC void cdaxpy(long int n,
2012 double da,
2013 realnum dx[],
2014 long int incx,
2015 realnum dy[],
2016 long int incy)
2017{
2018 long int i,
2019 i_,
2020 ix,
2021 iy,
2022 m;
2023
2024 DEBUG_ENTRY( "cdaxpy()" );
2025
2026 /* constant times a vector plus a vector.
2027 * uses unrolled loops for increments equal to one.
2028 * jack dongarra, lapack, 3/11/78.
2029 * */
2030
2031 if( n > 0 )
2032 {
2033 if( da != 0.00f )
2034 {
2035 if( incx == 1 && incy == 1 )
2036 {
2037
2038 /* code for both increments equal to 1
2039 *
2040 *
2041 * clean-up loop
2042 * */
2043 m = n%4;
2044 if( m != 0 )
2045 {
2046 for( i=1; i <= m; i++ )
2047 {
2048 i_ = i - 1;
2049 dy[i_] += (realnum)(da*dx[i_]);
2050 }
2051 if( n < 4 )
2052 {
2053 return;
2054 }
2055 }
2056
2057 for( i=m; i < n; i += 4 )
2058 {
2059 dy[i] += (realnum)(da*dx[i]);
2060 dy[i+1] += (realnum)(da*dx[i+1]);
2061 dy[i+2] += (realnum)(da*dx[i+2]);
2062 dy[i+3] += (realnum)(da*dx[i+3]);
2063 }
2064 }
2065 else
2066 {
2067
2068 /* code for unequal increments or equal increments
2069 * not equal to 1
2070 * */
2071 ix = 1;
2072 iy = 1;
2073 if( incx < 0 )
2074 ix = (-n + 1)*incx + 1;
2075 if( incy < 0 )
2076 iy = (-n + 1)*incy + 1;
2077 for( i=0; i < n; i++ )
2078 {
2079 dy[iy-1] += (realnum)(da*dx[ix-1]);
2080 ix += incx;
2081 iy += incy;
2082 }
2083 }
2084 }
2085 }
2086 return;
2087}
2088/*lint +e801 goto is deprecated */
2089/*lint +e64 type mismatch */
STATIC double da(double z, double temp, double eden)
#define STATIC
Definition cddefines.h:97
#define MAX3(a, b, c)
Definition cddefines.h:787
#define MIN2
Definition cddefines.h:761
#define POW2
Definition cddefines.h:929
T sign(T x, T y)
Definition cddefines.h:800
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
t_fe fe
Definition fe.cpp:5
chi2_type optimize_func(const realnum param[], int grid_index=-1)
STATIC double dist(long, realnum[], realnum[])
STATIC void sortd(long, realnum[], long[])
STATIC void start(long, realnum[], realnum[], long, long[], realnum *, int *)
STATIC void cdcopy(long, realnum[], long, realnum[], long)
struct t_usubc usubc
STATIC void cdaxpy(long, double, realnum[], long, realnum[], long)
STATIC void newpt(long, double, realnum[], realnum[], int, realnum[], int *)
#define S(I_, J_)
STATIC void csscal(long, double, realnum[], long)
STATIC void evalf(long, long[], realnum[], long, realnum[], realnum *, long *)
STATIC void fstats(double, long, int)
STATIC void simplx(long, realnum[], long, long[], long, int, realnum[], realnum *, long *, realnum *, realnum[], long *)
STATIC void calcc(long, realnum *, long, long, int, realnum[])
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)
struct t_isubc isubc
STATIC void setstp(long, long, realnum[], realnum[])
STATIC double cdasum(long, realnum[], long)
STATIC void partx(long, long[], realnum[], long *, long[])
STATIC void subopt(long)
STATIC void order(long, realnum[], long *, long *, long *)
realnum fbonus
realnum sfstop
realnum sfbest
realnum omega
realnum bonus
long int ifxsw
long int nfstop
realnum gamma
long int nsmin
realnum alpha
long int irepl
long int nsmax
realnum ftest
realnum fstop
realnum delta
long int nfxe
realnum fxstat[4]
static int nleft