« back to all changes in this revision
Viewing changes to base/subplex.c
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Bazaar Package Importer
- Author(s): Josselin Mouette
- Date: 2002-04-17 10:36:45 UTC
- Revision ID: james.westby@ubuntu.com-20020417103645-29vomjspk4yf4olw
Tags: upstream-2.1
Import upstream version 2.1
- files added:
- autom4te.cache
- base
- doc
- examples
- utils
added
removed
base/subplex.c
1
/*
2
Downloaded from http://www.netlib.org/opt/subplex.tgz
3
4
README file for SUBPLEX
5
6
NAME
7
subplex - subspace-searching simplex method for unconstrained
8
optimization
9
10
DESCRIPTION
11
Subplex is a subspace-searching simplex method for the
12
unconstrained optimization of general multivariate functions.
13
Like the Nelder-Mead simplex method it generalizes, the subplex
14
method is well suited for optimizing noisy objective functions.
15
The number of function evaluations required for convergence
16
typically increases only linearly with the problem size, so for
17
most applications the subplex method is much more efficient than
18
the simplex method.
19
20
INSTALLATION
21
To build subplex on UNIX systems, edit the Makefile as necessary
22
and type:
23
24
make
25
26
This will create a linkable library named subplex.a and a
27
demonstration executable named demo.
28
29
EXAMPLE
30
To run subplex on a simple objective function type:
31
32
demo < demo.in
33
34
To run subplex on other problems, edit a copy of the sample driver
35
demo.f as necessary.
36
37
AUTHOR
38
Tom Rowan
39
Oak Ridge National Laboratory
40
Mathematical Sciences Section
41
P.O. Box 2008, Bldg. 6012
42
Oak Ridge, TN 37831-6367
43
44
Phone: (423) 574-3131
45
Fax : (423) 574-0680
46
Email: na.rowan@na-net.ornl.gov
47
48
REFERENCE
49
T. Rowan, "Functional Stability Analysis of Numerical Algorithms",
50
Ph.D. thesis, Department of Computer Sciences, University of Texas
51
at Austin, 1990.
52
53
COMMENTS
54
Please send comments, suggestions, or bug reports to
55
na.rowan@na-net.ornl.gov.
56
*/
57
58
#include <stdio.h>
59
#include <stdlib.h>
60
#include <math.h>
61
62
#include "ctl.h"
63
64
typedef number doublereal;
65
typedef boolean logical;
66
67
#define TRUE_ 1
68
#define FALSE_ 0
69
70
typedef subplex_func D_fp;
71
72
#define max(a,b) ((a) > (b) ? (a) : (b))
73
#define min(a,b) ((a) < (b) ? (a) : (b))
74
#define abs(x) fabs(x)
75
76
/****************************************************************************/
77
/****************************************************************************/
78
79
/* dasum.f -- translated by f2c (version 19991025).
80
You must link the resulting object file with the libraries:
81
-lf2c -lm (in that order)
82
*/
83
84
static doublereal dasum_(integer *n, doublereal *dx, integer *incx)
85
{
86
/* System generated locals */
87
integer i__1;
88
doublereal ret_val, d__1, d__2, d__3, d__4, d__5, d__6;
89
90
/* Local variables */
91
integer i__, m;
92
doublereal dtemp;
93
integer ix, mp1;
94
95
96
/* takes the sum of the absolute values. */
97
/* uses unrolled loops for increment equal to one. */
98
/* jack dongarra, linpack, 3/11/78. */
99
/* modified to correct problem with negative increment, 8/21/90. */
100
101
102
/* Parameter adjustments */
103
--dx;
104
105
/* Function Body */
106
ret_val = 0.;
107
dtemp = 0.;
108
if (*n <= 0) {
109
return ret_val;
110
}
111
if (*incx == 1) {
112
goto L20;
113
}
114
115
/* code for increment not equal to 1 */
116
117
ix = 1;
118
if (*incx < 0) {
119
ix = (-(*n) + 1) * *incx + 1;
120
}
121
i__1 = *n;
122
for (i__ = 1; i__ <= i__1; ++i__) {
123
dtemp += (d__1 = dx[ix], abs(d__1));
124
ix += *incx;
125
/* L10: */
126
}
127
ret_val = dtemp;
128
return ret_val;
129
130
/* code for increment equal to 1 */
131
132
133
/* clean-up loop */
134
135
L20:
136
m = *n % 6;
137
if (m == 0) {
138
goto L40;
139
}
140
i__1 = m;
141
for (i__ = 1; i__ <= i__1; ++i__) {
142
dtemp += (d__1 = dx[i__], abs(d__1));
143
/* L30: */
144
}
145
if (*n < 6) {
146
goto L60;
147
}
148
L40:
149
mp1 = m + 1;
150
i__1 = *n;
151
for (i__ = mp1; i__ <= i__1; i__ += 6) {
152
dtemp = dtemp + (d__1 = dx[i__], abs(d__1)) + (d__2 = dx[i__ + 1],
153
abs(d__2)) + (d__3 = dx[i__ + 2], abs(d__3)) + (d__4 = dx[i__
154
+ 3], abs(d__4)) + (d__5 = dx[i__ + 4], abs(d__5)) + (d__6 =
155
dx[i__ + 5], abs(d__6));
156
/* L50: */
157
}
158
L60:
159
ret_val = dtemp;
160
return ret_val;
161
} /* dasum_ */
162
163
/* daxpy.f -- translated by f2c (version 19991025).
164
You must link the resulting object file with the libraries:
165
-lf2c -lm (in that order)
166
*/
167
168
static int daxpy_(integer *n, doublereal *da, doublereal *dx,
169
integer *incx, doublereal *dy, integer *incy)
170
{
171
/* System generated locals */
172
integer i__1;
173
174
/* Local variables */
175
integer i__, m, ix, iy, mp1;
176
177
178
/* constant times a vector plus a vector. */
179
/* uses unrolled loops for increments equal to one. */
180
/* jack dongarra, linpack, 3/11/78. */
181
182
183
/* Parameter adjustments */
184
--dy;
185
--dx;
186
187
/* Function Body */
188
if (*n <= 0) {
189
return 0;
190
}
191
if (*da == 0.) {
192
return 0;
193
}
194
if (*incx == 1 && *incy == 1) {
195
goto L20;
196
}
197
198
/* code for unequal increments or equal increments */
199
/* not equal to 1 */
200
201
ix = 1;
202
iy = 1;
203
if (*incx < 0) {
204
ix = (-(*n) + 1) * *incx + 1;
205
}
206
if (*incy < 0) {
207
iy = (-(*n) + 1) * *incy + 1;
208
}
209
i__1 = *n;
210
for (i__ = 1; i__ <= i__1; ++i__) {
211
dy[iy] += *da * dx[ix];
212
ix += *incx;
213
iy += *incy;
214
/* L10: */
215
}
216
return 0;
217
218
/* code for both increments equal to 1 */
219
220
221
/* clean-up loop */
222
223
L20:
224
m = *n % 4;
225
if (m == 0) {
226
goto L40;
227
}
228
i__1 = m;
229
for (i__ = 1; i__ <= i__1; ++i__) {
230
dy[i__] += *da * dx[i__];
231
/* L30: */
232
}
233
if (*n < 4) {
234
return 0;
235
}
236
L40:
237
mp1 = m + 1;
238
i__1 = *n;
239
for (i__ = mp1; i__ <= i__1; i__ += 4) {
240
dy[i__] += *da * dx[i__];
241
dy[i__ + 1] += *da * dx[i__ + 1];
242
dy[i__ + 2] += *da * dx[i__ + 2];
243
dy[i__ + 3] += *da * dx[i__ + 3];
244
/* L50: */
245
}
246
return 0;
247
} /* daxpy_ */
248
249
/* dcopy.f -- translated by f2c (version 19991025).
250
You must link the resulting object file with the libraries:
251
-lf2c -lm (in that order)
252
*/
253
254
static int dcopy_(integer *n, doublereal *dx, integer *incx,
255
doublereal *dy, integer *incy)
256
{
257
/* System generated locals */
258
integer i__1;
259
260
/* Local variables */
261
integer i__, m, ix, iy, mp1;
262
263
264
/* copies a vector, x, to a vector, y. */
265
/* uses unrolled loops for increments equal to one. */
266
/* jack dongarra, linpack, 3/11/78. */
267
268
269
/* Parameter adjustments */
270
--dy;
271
--dx;
272
273
/* Function Body */
274
if (*n <= 0) {
275
return 0;
276
}
277
if (*incx == 1 && *incy == 1) {
278
goto L20;
279
}
280
281
/* code for unequal increments or equal increments */
282
/* not equal to 1 */
283
284
ix = 1;
285
iy = 1;
286
if (*incx < 0) {
287
ix = (-(*n) + 1) * *incx + 1;
288
}
289
if (*incy < 0) {
290
iy = (-(*n) + 1) * *incy + 1;
291
}
292
i__1 = *n;
293
for (i__ = 1; i__ <= i__1; ++i__) {
294
dy[iy] = dx[ix];
295
ix += *incx;
296
iy += *incy;
297
/* L10: */
298
}
299
return 0;
300
301
/* code for both increments equal to 1 */
302
303
304
/* clean-up loop */
305
306
L20:
307
m = *n % 7;
308
if (m == 0) {
309
goto L40;
310
}
311
i__1 = m;
312
for (i__ = 1; i__ <= i__1; ++i__) {
313
dy[i__] = dx[i__];
314
/* L30: */
315
}
316
if (*n < 7) {
317
return 0;
318
}
319
L40:
320
mp1 = m + 1;
321
i__1 = *n;
322
for (i__ = mp1; i__ <= i__1; i__ += 7) {
323
dy[i__] = dx[i__];
324
dy[i__ + 1] = dx[i__ + 1];
325
dy[i__ + 2] = dx[i__ + 2];
326
dy[i__ + 3] = dx[i__ + 3];
327
dy[i__ + 4] = dx[i__ + 4];
328
dy[i__ + 5] = dx[i__ + 5];
329
dy[i__ + 6] = dx[i__ + 6];
330
/* L50: */
331
}
332
return 0;
333
} /* dcopy_ */
334
335
/* dscal.f -- translated by f2c (version 19991025).
336
You must link the resulting object file with the libraries:
337
-lf2c -lm (in that order)
338
*/
339
340
static int dscal_(integer *n, doublereal *da, doublereal *dx,
341
integer *incx)
342
{
343
/* System generated locals */
344
integer i__1;
345
346
/* Local variables */
347
integer i__, m, ix, mp1;
348
349
350
/* scales a vector by a constant. */
351
/* uses unrolled loops for increment equal to one. */
352
/* jack dongarra, linpack, 3/11/78. */
353
/* modified to correct problem with negative increment, 8/21/90. */
354
355
356
/* Parameter adjustments */
357
--dx;
358
359
/* Function Body */
360
if (*n <= 0) {
361
return 0;
362
}
363
if (*incx == 1) {
364
goto L20;
365
}
366
367
/* code for increment not equal to 1 */
368
369
ix = 1;
370
if (*incx < 0) {
371
ix = (-(*n) + 1) * *incx + 1;
372
}
373
i__1 = *n;
374
for (i__ = 1; i__ <= i__1; ++i__) {
375
dx[ix] = *da * dx[ix];
376
ix += *incx;
377
/* L10: */
378
}
379
return 0;
380
381
/* code for increment equal to 1 */
382
383
384
/* clean-up loop */
385
386
L20:
387
m = *n % 5;
388
if (m == 0) {
389
goto L40;
390
}
391
i__1 = m;
392
for (i__ = 1; i__ <= i__1; ++i__) {
393
dx[i__] = *da * dx[i__];
394
/* L30: */
395
}
396
if (*n < 5) {
397
return 0;
398
}
399
L40:
400
mp1 = m + 1;
401
i__1 = *n;
402
for (i__ = mp1; i__ <= i__1; i__ += 5) {
403
dx[i__] = *da * dx[i__];
404
dx[i__ + 1] = *da * dx[i__ + 1];
405
dx[i__ + 2] = *da * dx[i__ + 2];
406
dx[i__ + 3] = *da * dx[i__ + 3];
407
dx[i__ + 4] = *da * dx[i__ + 4];
408
/* L50: */
409
}
410
return 0;
411
} /* dscal_ */
412
413
/* dist.f -- translated by f2c (version 19991025).
414
You must link the resulting object file with the libraries:
415
-lf2c -lm (in that order)
416
*/
417
418
static doublereal dist_(integer *n, doublereal *x, doublereal *y)
419
{
420
/* System generated locals */
421
integer i__1;
422
doublereal ret_val, d__1;
423
424
/* Builtin functions */
425
double sqrt(doublereal);
426
427
/* Local variables */
428
integer i__;
429
doublereal scale, absxmy, sum;
430
431
432
433
/* Coded by Tom Rowan */
434
/* Department of Computer Sciences */
435
/* University of Texas at Austin */
436
437
/* dist calculates the distance between the points x,y. */
438
439
/* input */
440
441
/* n - number of components */
442
443
/* x - point in n-space */
444
445
/* y - point in n-space */
446
447
/* local variables */
448
449
450
/* subroutines and functions */
451
452
/* fortran */
453
454
/* ----------------------------------------------------------- */
455
456
/* Parameter adjustments */
457
--y;
458
--x;
459
460
/* Function Body */
461
absxmy = (d__1 = x[1] - y[1], abs(d__1));
462
if (absxmy <= 1.) {
463
sum = absxmy * absxmy;
464
scale = 1.;
465
} else {
466
sum = 1.;
467
scale = absxmy;
468
}
469
i__1 = *n;
470
for (i__ = 2; i__ <= i__1; ++i__) {
471
absxmy = (d__1 = x[i__] - y[i__], abs(d__1));
472
if (absxmy <= scale) {
473
/* Computing 2nd power */
474
d__1 = absxmy / scale;
475
sum += d__1 * d__1;
476
} else {
477
/* Computing 2nd power */
478
d__1 = scale / absxmy;
479
sum = sum * (d__1 * d__1) + 1.;
480
scale = absxmy;
481
}
482
/* L10: */
483
}
484
ret_val = scale * sqrt(sum);
485
return ret_val;
486
} /* dist_ */
487
488
/* calcc.f -- translated by f2c (version 19991025).
489
You must link the resulting object file with the libraries:
490
-lf2c -lm (in that order)
491
*/
492
493
/* Table of constant values */
494
495
static doublereal c_b3 = 0.;
496
static integer c__0 = 0;
497
static integer c__1 = 1;
498
static doublereal c_b7 = 1.;
499
500
static int calcc_(integer *ns, doublereal *s, integer *ih, integer *
501
inew, logical *updatc, doublereal *c__)
502
{
503
/* System generated locals */
504
integer s_dim1, s_offset, i__1;
505
doublereal d__1;
506
507
/* Local variables */
508
integer i__, j;
509
510
/* Coded by Tom Rowan */
511
/* Department of Computer Sciences */
512
/* University of Texas at Austin */
513
514
/* calcc calculates the centroid of the simplex without the */
515
/* vertex with highest function value. */
516
517
/* input */
518
519
/* ns - subspace dimension */
520
521
/* s - double precision work space of dimension .ge. */
522
/* ns*(ns+3) used to store simplex */
523
524
/* ih - index to vertex with highest function value */
525
526
/* inew - index to new point */
527
528
/* updatc - logical switch */
529
/* = .true. : update centroid */
530
/* = .false. : calculate centroid from scratch */
531
532
/* c - centroid of the simplex without vertex with */
533
/* highest function value */
534
535
/* output */
536
537
/* c - new centroid */
538
539
/* local variables */
540
541
542
/* subroutines and functions */
543
544
/* blas */
545
546
/* ----------------------------------------------------------- */
547
548
/* Parameter adjustments */
549
--c__;
550
s_dim1 = *ns;
551
s_offset = 1 + s_dim1 * 1;
552
s -= s_offset;
553
554
/* Function Body */
555
if (*updatc) {
556
if (*ih == *inew) {
557
return 0;
558
}
559
i__1 = *ns;
560
for (i__ = 1; i__ <= i__1; ++i__) {
561
c__[i__] += (s[i__ + *inew * s_dim1] - s[i__ + *ih * s_dim1]) / *
562
ns;
563
/* L10: */
564
}
565
} else {
566
dcopy_(ns, &c_b3, &c__0, &c__[1], &c__1);
567
i__1 = *ns + 1;
568
for (j = 1; j <= i__1; ++j) {
569
if (j != *ih) {
570
daxpy_(ns, &c_b7, &s[j * s_dim1 + 1], &c__1, &c__[1], &c__1);
571
}
572
/* L20: */
573
}
574
d__1 = 1. / *ns;
575
dscal_(ns, &d__1, &c__[1], &c__1);
576
}
577
return 0;
578
} /* calcc_ */
579
580
/* order.f -- translated by f2c (version 19991025).
581
You must link the resulting object file with the libraries:
582
-lf2c -lm (in that order)
583
*/
584
585
static int order_(integer *npts, doublereal *fs, integer *il,
586
integer *is, integer *ih)
587
{
588
/* System generated locals */
589
integer i__1;
590
591
/* Local variables */
592
integer i__, j, il0;
593
594
595
596
/* Coded by Tom Rowan */
597
/* Department of Computer Sciences */
598
/* University of Texas at Austin */
599
600
/* order determines the indices of the vertices with the */
601
/* lowest, second highest, and highest function values. */
602
603
/* input */
604
605
/* npts - number of points in simplex */
606
607
/* fs - double precision vector of function values of */
608
/* simplex */
609
610
/* il - index to vertex with lowest function value */
611
612
/* output */
613
614
/* il - new index to vertex with lowest function value */
615
616
/* is - new index to vertex with second highest */
617
/* function value */
618
619
/* ih - new index to vertex with highest function value */
620
621
/* local variables */
622
623
624
/* subroutines and functions */
625
626
/* fortran */
627
628
/* ----------------------------------------------------------- */
629
630
/* Parameter adjustments */
631
--fs;
632
633
/* Function Body */
634
il0 = *il;
635
j = il0 % *npts + 1;
636
if (fs[j] >= fs[*il]) {
637
*ih = j;
638
*is = il0;
639
} else {
640
*ih = il0;
641
*is = j;
642
*il = j;
643
}
644
i__1 = il0 + *npts - 2;
645
for (i__ = il0 + 1; i__ <= i__1; ++i__) {
646
j = i__ % *npts + 1;
647
if (fs[j] >= fs[*ih]) {
648
*is = *ih;
649
*ih = j;
650
} else if (fs[j] > fs[*is]) {
651
*is = j;
652
} else if (fs[j] < fs[*il]) {
653
*il = j;
654
}
655
/* L10: */
656
}
657
return 0;
658
} /* order_ */
659
660
/* partx.f -- translated by f2c (version 19991025).
661
You must link the resulting object file with the libraries:
662
-lf2c -lm (in that order)
663
*/
664
665
/* Common Block Declarations */
666
667
static struct {
668
doublereal alpha, beta, gamma, delta, psi, omega;
669
integer nsmin, nsmax, irepl, ifxsw;
670
doublereal bonus, fstop;
671
integer nfstop, nfxe;
672
doublereal fxstat[4], ftest;
673
logical minf, initx, newx;
674
} usubc_;
675
676
#define usubc_1 usubc_
677
678
static int partx_(integer *n, integer *ip, doublereal *absdx,
679
integer *nsubs, integer *nsvals)
680
{
681
/* System generated locals */
682
integer i__1;
683
684
/* Local variables */
685
static integer i__, nleft, nused;
686
static doublereal as1max, gapmax, asleft, as1, as2;
687
static integer ns1, ns2;
688
static doublereal gap;
689
690
691
692
/* Coded by Tom Rowan */
693
/* Department of Computer Sciences */
694
/* University of Texas at Austin */
695
696
/* partx partitions the vector x by grouping components of */
697
/* similar magnitude of change. */
698
699
/* input */
700
701
/* n - number of components (problem dimension) */
702
703
/* ip - permutation vector */
704
705
/* absdx - vector of magnitude of change in x */
706
707
/* nsvals - integer array dimensioned .ge. int(n/nsmin) */
708
709
/* output */
710
711
/* nsubs - number of subspaces */
712
713
/* nsvals - integer array of subspace dimensions */
714
715
/* common */
716
717
718
719
/* local variables */
720
721
722
723
/* subroutines and functions */
724
725
/* fortran */
726
727
/* ----------------------------------------------------------- */
728
729
/* Parameter adjustments */
730
--absdx;
731
--ip;
732
--nsvals;
733
734
/* Function Body */
735
*nsubs = 0;
736
nused = 0;
737
nleft = *n;
738
asleft = absdx[1];
739
i__1 = *n;
740
for (i__ = 2; i__ <= i__1; ++i__) {
741
asleft += absdx[i__];
742
/* L10: */
743
}
744
L20:
745
if (nused < *n) {
746
++(*nsubs);
747
as1 = 0.;
748
i__1 = usubc_1.nsmin - 1;
749
for (i__ = 1; i__ <= i__1; ++i__) {
750
as1 += absdx[ip[nused + i__]];
751
/* L30: */
752
}
753
gapmax = -1.;
754
i__1 = min(usubc_1.nsmax,nleft);
755
for (ns1 = usubc_1.nsmin; ns1 <= i__1; ++ns1) {
756
as1 += absdx[ip[nused + ns1]];
757
ns2 = nleft - ns1;
758
if (ns2 > 0) {
759
if (ns2 >= ((ns2 - 1) / usubc_1.nsmax + 1) * usubc_1.nsmin) {
760
as2 = asleft - as1;
761
gap = as1 / ns1 - as2 / ns2;
762
if (gap > gapmax) {
763
gapmax = gap;
764
nsvals[*nsubs] = ns1;
765
as1max = as1;
766
}
767
}
768
} else {
769
if (as1 / ns1 > gapmax) {
770
nsvals[*nsubs] = ns1;
771
return 0;
772
}
773
}
774
/* L40: */
775
}
776
nused += nsvals[*nsubs];
777
nleft = *n - nused;
778
asleft -= as1max;
779
goto L20;
780
}
781
return 0;
782
} /* partx_ */
783
784
/* sortd.f -- translated by f2c (version 19991025).
785
You must link the resulting object file with the libraries:
786
-lf2c -lm (in that order)
787
*/
788
789
static int sortd_(integer *n, doublereal *xkey, integer *ix)
790
{
791
/* System generated locals */
792
integer i__1;
793
794
/* Local variables */
795
integer ixip1, i__, ilast, iswap, ifirst, ixi;
796
797
798
799
/* Coded by Tom Rowan */
800
/* Department of Computer Sciences */
801
/* University of Texas at Austin */
802
803
/* sortd uses the shakersort method to sort an array of keys */
804
/* in decreasing order. The sort is performed implicitly by */
805
/* modifying a vector of indices. */
806
807
/* For nearly sorted arrays, sortd requires O(n) comparisons. */
808
/* for completely unsorted arrays, sortd requires O(n**2) */
809
/* comparisons and will be inefficient unless n is small. */
810
811
/* input */
812
813
/* n - number of components */
814
815
/* xkey - double precision vector of keys */
816
817
/* ix - integer vector of indices */
818
819
/* output */
820
821
/* ix - indices satisfy xkey(ix(i)) .ge. xkey(ix(i+1)) */
822
/* for i = 1,...,n-1 */
823
824
/* local variables */
825
826
827
/* ----------------------------------------------------------- */
828
829
/* Parameter adjustments */
830
--ix;
831
--xkey;
832
833
/* Function Body */
834
ifirst = 1;
835
iswap = 1;
836
ilast = *n - 1;
837
L10:
838
if (ifirst <= ilast) {
839
i__1 = ilast;
840
for (i__ = ifirst; i__ <= i__1; ++i__) {
841
ixi = ix[i__];
842
ixip1 = ix[i__ + 1];
843
if (xkey[ixi] < xkey[ixip1]) {
844
ix[i__] = ixip1;
845
ix[i__ + 1] = ixi;
846
iswap = i__;
847
}
848
/* L20: */
849
}
850
ilast = iswap - 1;
851
i__1 = ifirst;
852
for (i__ = ilast; i__ >= i__1; --i__) {
853
ixi = ix[i__];
854
ixip1 = ix[i__ + 1];
855
if (xkey[ixi] < xkey[ixip1]) {
856
ix[i__] = ixip1;
857
ix[i__ + 1] = ixi;
858
iswap = i__;
859
}
860
/* L30: */
861
}
862
ifirst = iswap + 1;
863
goto L10;
864
}
865
return 0;
866
} /* sortd_ */
867
868
/* newpt.f -- translated by f2c (version 19991025).
869
You must link the resulting object file with the libraries:
870
-lf2c -lm (in that order)
871
*/
872
873
static int newpt_(integer *ns, doublereal *coef, doublereal *xbase,
874
doublereal *xold, logical *new__, doublereal *xnew, logical *small)
875
{
876
/* System generated locals */
877
integer i__1;
878
879
/* Local variables */
880
integer i__;
881
logical eqold;
882
doublereal xoldi;
883
logical eqbase;
884
885
886
887
/* Coded by Tom Rowan */
888
/* Department of Computer Sciences */
889
/* University of Texas at Austin */
890
891
/* newpt performs reflections, expansions, contractions, and */
892
/* shrinkages (massive contractions) by computing: */
893
894
/* xbase + coef * (xbase - xold) */
895
896
/* The result is stored in xnew if new .eq. .true., */
897
/* in xold otherwise. */
898
899
/* use : coef .gt. 0 to reflect */
900
/* coef .lt. 0 to expand, contract, or shrink */
901
902
/* input */
903
904
/* ns - number of components (subspace dimension) */
905
906
/* coef - one of four simplex method coefficients */
907
908
/* xbase - double precision ns-vector representing base */
909
/* point */
910
911
/* xold - double precision ns-vector representing old */
912
/* point */
913
914
/* new - logical switch */
915
/* = .true. : store result in xnew */
916
/* = .false. : store result in xold, xnew is not */
917
/* referenced */
918
919
/* output */
920
921
/* xold - unchanged if new .eq. .true., contains new */
922
/* point otherwise */
923
924
/* xnew - double precision ns-vector representing new */
925
/* point if new .eq. .true., not referenced */
926
/* otherwise */
927
928
/* small - logical flag */
929
/* = .true. : coincident points */
930
/* = .false. : otherwise */
931
932
/* local variables */
933
934
935
/* subroutines and functions */
936
937
/* fortran */
938
939
/* ----------------------------------------------------------- */
940
941
/* Parameter adjustments */
942
--xold;
943
--xbase;
944
--xnew;
945
946
/* Function Body */
947
eqbase = TRUE_;
948
eqold = TRUE_;
949
if (*new__) {
950
i__1 = *ns;
951
for (i__ = 1; i__ <= i__1; ++i__) {
952
xnew[i__] = xbase[i__] + *coef * (xbase[i__] - xold[i__]);
953
eqbase = eqbase && xnew[i__] == xbase[i__];
954
eqold = eqold && xnew[i__] == xold[i__];
955
/* L10: */
956
}
957
} else {
958
i__1 = *ns;
959
for (i__ = 1; i__ <= i__1; ++i__) {
960
xoldi = xold[i__];
961
xold[i__] = xbase[i__] + *coef * (xbase[i__] - xold[i__]);
962
eqbase = eqbase && xold[i__] == xbase[i__];
963
eqold = eqold && xold[i__] == xoldi;
964
/* L20: */
965
}
966
}
967
*small = eqbase || eqold;
968
return 0;
969
} /* newpt_ */
970
971
/* start.f -- translated by f2c (version 19991025).
972
You must link the resulting object file with the libraries:
973
-lf2c -lm (in that order)
974
*/
975
976
static int start_(integer *n, doublereal *x, doublereal *step,
977
integer *ns, integer *ips, doublereal *s, logical *small)
978
{
979
/* System generated locals */
980
integer s_dim1, s_offset, i__1;
981
982
/* Local variables */
983
integer i__, j;
984
985
986
/* Coded by Tom Rowan */
987
/* Department of Computer Sciences */
988
/* University of Texas at Austin */
989
990
/* start creates the initial simplex for simplx minimization. */
991
992
/* input */
993
994
/* n - problem dimension */
995
996
/* x - current best point */
997
998
/* step - stepsizes for corresponding components of x */
999
1000
/* ns - subspace dimension */
1001
1002
/* ips - permutation vector */
1003
1004
1005
/* output */
1006
1007
/* s - first ns+1 columns contain initial simplex */
1008
1009
/* small - logical flag */
1010
/* = .true. : coincident points */
1011
/* = .false. : otherwise */
1012
1013
/* local variables */
1014
1015
1016
/* subroutines and functions */
1017
1018
/* blas */
1019
/* fortran */
1020
1021
/* ----------------------------------------------------------- */
1022
1023
/* Parameter adjustments */
1024
--ips;
1025
--step;
1026
--x;
1027
s_dim1 = *ns;
1028
s_offset = 1 + s_dim1 * 1;
1029
s -= s_offset;
1030
1031
/* Function Body */
1032
i__1 = *ns;
1033
for (i__ = 1; i__ <= i__1; ++i__) {
1034
s[i__ + s_dim1] = x[ips[i__]];
1035
/* L10: */
1036
}
1037
i__1 = *ns + 1;
1038
for (j = 2; j <= i__1; ++j) {
1039
dcopy_(ns, &s[s_dim1 + 1], &c__1, &s[j * s_dim1 + 1], &c__1);
1040
s[j - 1 + j * s_dim1] = s[j - 1 + s_dim1] + step[ips[j - 1]];
1041
/* L20: */
1042
}
1043
1044
/* check for coincident points */
1045
1046
i__1 = *ns + 1;
1047
for (j = 2; j <= i__1; ++j) {
1048
if (s[j - 1 + j * s_dim1] == s[j - 1 + s_dim1]) {
1049
goto L40;
1050
}
1051
/* L30: */
1052
}
1053
*small = FALSE_;
1054
return 0;
1055
1056
/* coincident points */
1057
1058
L40:
1059
*small = TRUE_;
1060
return 0;
1061
} /* start_ */
1062
1063
/* fstats.f -- translated by f2c (version 19991025).
1064
You must link the resulting object file with the libraries:
1065
-lf2c -lm (in that order)
1066
*/
1067
1068
static int fstats_(doublereal *fx, integer *ifxwt, logical *reset)
1069
{
1070
/* System generated locals */
1071
doublereal d__1, d__2, d__3;
1072
1073
/* Builtin functions */
1074
double sqrt(doublereal);
1075
1076
/* Local variables */
1077
static doublereal fscale;
1078
static integer nsv;
1079
static doublereal f1sv;
1080
1081
1082
1083
/* Coded by Tom Rowan */
1084
/* Department of Computer Sciences */
1085
/* University of Texas at Austin */
1086
1087
/* fstats modifies the common /usubc/ variables nfxe,fxstat. */
1088
1089
/* input */
1090
1091
/* fx - most recent evaluation of f at best x */
1092
1093
/* ifxwt - integer weight for fx */
1094
1095
/* reset - logical switch */
1096
/* = .true. : initialize nfxe,fxstat */
1097
/* = .false. : update nfxe,fxstat */
1098
1099
/* common */
1100
1101
1102
1103
/* local variables */
1104
1105
1106
1107
/* subroutines and functions */
1108
1109
/* fortran */
1110
1111
/* ----------------------------------------------------------- */
1112
1113
if (*reset) {
1114
usubc_1.nfxe = *ifxwt;
1115
usubc_1.fxstat[0] = *fx;
1116
usubc_1.fxstat[1] = *fx;
1117
usubc_1.fxstat[2] = *fx;
1118
usubc_1.fxstat[3] = 0.;
1119
} else {
1120
nsv = usubc_1.nfxe;
1121
f1sv = usubc_1.fxstat[0];
1122
usubc_1.nfxe += *ifxwt;
1123
usubc_1.fxstat[0] += *ifxwt * (*fx - usubc_1.fxstat[0]) /
1124
usubc_1.nfxe;
1125
usubc_1.fxstat[1] = max(usubc_1.fxstat[1],*fx);
1126
usubc_1.fxstat[2] = min(usubc_1.fxstat[2],*fx);
1127
/* Computing MAX */
1128
d__1 = abs(usubc_1.fxstat[1]), d__2 = abs(usubc_1.fxstat[2]), d__1 =
1129
max(d__1,d__2);
1130
fscale = max(d__1,1.);
1131
/* Computing 2nd power */
1132
d__1 = usubc_1.fxstat[3] / fscale;
1133
/* Computing 2nd power */
1134
d__2 = (usubc_1.fxstat[0] - f1sv) / fscale;
1135
/* Computing 2nd power */
1136
d__3 = (*fx - usubc_1.fxstat[0]) / fscale;
1137
usubc_1.fxstat[3] = fscale * sqrt(((nsv - 1) * (d__1 * d__1) + nsv * (
1138
d__2 * d__2) + *ifxwt * (d__3 * d__3)) / (usubc_1.nfxe - 1));
1139
}
1140
return 0;
1141
} /* fstats_ */
1142
1143
/* evalf.f -- translated by f2c (version 19991025).
1144
You must link the resulting object file with the libraries:
1145
-lf2c -lm (in that order)
1146
*/
1147
1148
/* Common Block Declarations */
1149
1150
static struct {
1151
doublereal fbonus, sfstop, sfbest;
1152
logical new__;
1153
} isubc_;
1154
1155
#define isubc_1 isubc_
1156
1157
static logical c_true = TRUE_;
1158
static logical c_false = FALSE_;
1159
1160
static int evalf_(D_fp f,void*fdata, integer *ns, integer *ips, doublereal *xs,
1161
integer *n, doublereal *x, doublereal *sfx, integer *nfe)
1162
{
1163
/* System generated locals */
1164
integer i__1;
1165
1166
/* Local variables */
1167
static integer i__;
1168
static doublereal fx;
1169
static logical newbst;
1170
1171
/* Coded by Tom Rowan */
1172
/* Department of Computer Sciences */
1173
/* University of Texas at Austin */
1174
1175
/* evalf evaluates the function f at a point defined by x */
1176
/* with ns of its components replaced by those in xs. */
1177
1178
/* input */
1179
1180
/* f - user supplied function f(n,x) to be optimized */
1181
1182
/* ns - subspace dimension */
1183
1184
/* ips - permutation vector */
1185
1186
/* xs - double precision ns-vector to be mapped to x */
1187
1188
/* n - problem dimension */
1189
1190
/* x - double precision n-vector */
1191
1192
/* nfe - number of function evaluations */
1193
1194
/* output */
1195
1196
/* sfx - signed value of f evaluated at x */
1197
1198
/* nfe - incremented number of function evaluations */
1199
1200
/* common */
1201
1202
1203
1204
1205
1206
/* local variables */
1207
1208
1209
1210
/* subroutines and functions */
1211
1212
1213
/* ----------------------------------------------------------- */
1214
1215
/* Parameter adjustments */
1216
--ips;
1217
--xs;
1218
--x;
1219
1220
/* Function Body */
1221
i__1 = *ns;
1222
for (i__ = 1; i__ <= i__1; ++i__) {
1223
x[ips[i__]] = xs[i__];
1224
/* L10: */
1225
}
1226
usubc_1.newx = isubc_1.new__ || usubc_1.irepl != 2;
1227
fx = (*f)(*n, &x[1], fdata);
1228
if (usubc_1.irepl == 0) {
1229
if (usubc_1.minf) {
1230
*sfx = fx;
1231
} else {
1232
*sfx = -fx;
1233
}
1234
} else if (isubc_1.new__) {
1235
if (usubc_1.minf) {
1236
*sfx = fx;
1237
newbst = fx < usubc_1.ftest;
1238
} else {
1239
*sfx = -fx;
1240
newbst = fx > usubc_1.ftest;
1241
}
1242
if (usubc_1.initx || newbst) {
1243
if (usubc_1.irepl == 1) {
1244
fstats_(&fx, &c__1, &c_true);
1245
}
1246
usubc_1.ftest = fx;
1247
isubc_1.sfbest = *sfx;
1248
}
1249
} else {
1250
if (usubc_1.irepl == 1) {
1251
fstats_(&fx, &c__1, &c_false);
1252
fx = usubc_1.fxstat[usubc_1.ifxsw - 1];
1253
}
1254
usubc_1.ftest = fx + isubc_1.fbonus * usubc_1.fxstat[3];
1255
if (usubc_1.minf) {
1256
*sfx = usubc_1.ftest;
1257
isubc_1.sfbest = fx;
1258
} else {
1259
*sfx = -usubc_1.ftest;
1260
isubc_1.sfbest = -fx;
1261
}
1262
}
1263
++(*nfe);
1264
return 0;
1265
} /* evalf_ */
1266
1267
/* simplx.f -- translated by f2c (version 19991025).
1268
You must link the resulting object file with the libraries:
1269
-lf2c -lm (in that order)
1270
*/
1271
1272
static int simplx_(D_fp f, void *fdata, integer *n, doublereal *step, integer *
1273
ns, integer *ips, integer *maxnfe, logical *cmode, doublereal *x,
1274
doublereal *fx, integer *nfe, doublereal *s, doublereal *fs, integer *
1275
iflag)
1276
{
1277
/* System generated locals */
1278
integer s_dim1, s_offset, i__1;
1279
doublereal d__1, d__2;
1280
1281
/* Local variables */
1282
static integer inew;
1283
static integer npts;
1284
static integer i__, j;
1285
static integer icent;
1286
static logical small;
1287
static integer itemp;
1288
static doublereal fc, fe;
1289
static integer ih, il;
1290
static doublereal fr;
1291
static integer is;
1292
static logical updatc;
1293
static doublereal dum, tol;
1294
1295
1296
1297
/* Coded by Tom Rowan */
1298
/* Department of Computer Sciences */
1299
/* University of Texas at Austin */
1300
1301
/* simplx uses the Nelder-Mead simplex method to minimize the */
1302
/* function f on a subspace. */
1303
1304
/* input */
1305
1306
/* f - function to be minimized, declared external in */
1307
/* calling routine */
1308
1309
/* n - problem dimension */
1310
1311
/* step - stepsizes for corresponding components of x */
1312
1313
/* ns - subspace dimension */
1314
1315
/* ips - permutation vector */
1316
1317
/* maxnfe - maximum number of function evaluations */
1318
1319
/* cmode - logical switch */
1320
/* = .true. : continuation of previous call */
1321
/* = .false. : first call */
1322
1323
/* x - starting guess for minimum */
1324
1325
/* fx - value of f at x */
1326
1327
/* nfe - number of function evaluations */
1328
1329
/* s - double precision work array of dimension .ge. */
1330
/* ns*(ns+3) used to store simplex */
1331
1332
/* fs - double precision work array of dimension .ge. */
1333
/* ns+1 used to store function values of simplex */
1334
/* vertices */
1335
1336
/* output */
1337
1338
/* x - computed minimum */
1339
1340
/* fx - value of f at x */
1341
1342
/* nfe - incremented number of function evaluations */
1343
1344
/* iflag - error flag */
1345
/* = -1 : maxnfe exceeded */
1346
/* = 0 : simplex reduced by factor of psi */
1347
/* = 1 : limit of machine precision */
1348
/* = 2 : reached fstop */
1349
1350
/* common */
1351
1352
1353
1354
1355
1356
/* local variables */
1357
1358
1359
1360
/* subroutines and functions */
1361
1362
/* blas */
1363
/* fortran */
1364
1365
/* ----------------------------------------------------------- */
1366
1367
/* Parameter adjustments */
1368
--x;
1369
--step;
1370
--fs;
1371
s_dim1 = *ns;
1372
s_offset = 1 + s_dim1 * 1;
1373
s -= s_offset;
1374
--ips;
1375
1376
/* Function Body */
1377
if (*cmode) {
1378
goto L50;
1379
}
1380
npts = *ns + 1;
1381
icent = *ns + 2;
1382
itemp = *ns + 3;
1383
updatc = FALSE_;
1384
start_(n, &x[1], &step[1], ns, &ips[1], &s[s_offset], &small);
1385
if (small) {
1386
*iflag = 1;
1387
return 0;
1388
}
1389
if (usubc_1.irepl > 0) {
1390
isubc_1.new__ = FALSE_;
1391
evalf_((D_fp)f,fdata, ns, &ips[1], &s[s_dim1 + 1], n, &x[1], &fs[1], nfe);
1392
} else {
1393
fs[1] = *fx;
1394
}
1395
isubc_1.new__ = TRUE_;
1396
i__1 = npts;
1397
for (j = 2; j <= i__1; ++j) {
1398
evalf_((D_fp)f, fdata,ns, &ips[1], &s[j * s_dim1 + 1], n, &x[1], &fs[j],
1399
nfe);
1400
/* L10: */
1401
}
1402
il = 1;
1403
order_(&npts, &fs[1], &il, &is, &ih);
1404
tol = usubc_1.psi * dist_(ns, &s[ih * s_dim1 + 1], &s[il * s_dim1 + 1]);
1405
1406
/* main loop */
1407
1408
L20:
1409
calcc_(ns, &s[s_offset], &ih, &inew, &updatc, &s[icent * s_dim1 + 1]);
1410
updatc = TRUE_;
1411
inew = ih;
1412
1413
/* reflect */
1414
1415
newpt_(ns, &usubc_1.alpha, &s[icent * s_dim1 + 1], &s[ih * s_dim1 + 1], &
1416
c_true, &s[itemp * s_dim1 + 1], &small);
1417
if (small) {
1418
goto L40;
1419
}
1420
evalf_((D_fp)f,fdata, ns, &ips[1], &s[itemp * s_dim1 + 1], n, &x[1], &fr, nfe);
1421
if (fr < fs[il]) {
1422
1423
/* expand */
1424
1425
d__1 = -usubc_1.gamma;
1426
newpt_(ns, &d__1, &s[icent * s_dim1 + 1], &s[itemp * s_dim1 + 1], &
1427
c_true, &s[ih * s_dim1 + 1], &small);
1428
if (small) {
1429
goto L40;
1430
}
1431
evalf_((D_fp)f,fdata, ns, &ips[1], &s[ih * s_dim1 + 1], n, &x[1], &fe, nfe);
1432
if (fe < fr) {
1433
fs[ih] = fe;
1434
} else {
1435
dcopy_(ns, &s[itemp * s_dim1 + 1], &c__1, &s[ih * s_dim1 + 1], &
1436
c__1);
1437
fs[ih] = fr;
1438
}
1439
} else if (fr < fs[is]) {
1440
1441
/* accept reflected point */
1442
1443
dcopy_(ns, &s[itemp * s_dim1 + 1], &c__1, &s[ih * s_dim1 + 1], &c__1);
1444
fs[ih] = fr;
1445
} else {
1446
1447
/* contract */
1448
1449
if (fr > fs[ih]) {
1450
d__1 = -usubc_1.beta;
1451
newpt_(ns, &d__1, &s[icent * s_dim1 + 1], &s[ih * s_dim1 + 1], &
1452
c_true, &s[itemp * s_dim1 + 1], &small);
1453
} else {
1454
d__1 = -usubc_1.beta;
1455
newpt_(ns, &d__1, &s[icent * s_dim1 + 1], &s[itemp * s_dim1 + 1],
1456
&c_false, &dum, &small);
1457
}
1458
if (small) {
1459
goto L40;
1460
}
1461
evalf_((D_fp)f,fdata, ns, &ips[1], &s[itemp * s_dim1 + 1], n, &x[1], &fc,
1462
nfe);
1463
/* Computing MIN */
1464
d__1 = fr, d__2 = fs[ih];
1465
if (fc < min(d__1,d__2)) {
1466
dcopy_(ns, &s[itemp * s_dim1 + 1], &c__1, &s[ih * s_dim1 + 1], &
1467
c__1);
1468
fs[ih] = fc;
1469
} else {
1470
1471
/* shrink simplex */
1472
1473
i__1 = npts;
1474
for (j = 1; j <= i__1; ++j) {
1475
if (j != il) {
1476
d__1 = -usubc_1.delta;
1477
newpt_(ns, &d__1, &s[il * s_dim1 + 1], &s[j * s_dim1 + 1],
1478
&c_false, &dum, &small);
1479
if (small) {
1480
goto L40;
1481
}
1482
evalf_((D_fp)f,fdata, ns, &ips[1], &s[j * s_dim1 + 1], n, &x[1],
1483
&fs[j], nfe);
1484
}
1485
/* L30: */
1486
}
1487
}
1488
updatc = FALSE_;
1489
}
1490
order_(&npts, &fs[1], &il, &is, &ih);
1491
1492
/* check termination */
1493
1494
L40:
1495
if (usubc_1.irepl == 0) {
1496
*fx = fs[il];
1497
} else {
1498
*fx = isubc_1.sfbest;
1499
}
1500
L50:
1501
if (usubc_1.nfstop > 0 && *fx <= isubc_1.sfstop && usubc_1.nfxe >=
1502
usubc_1.nfstop) {
1503
*iflag = 2;
1504
} else if (*nfe >= *maxnfe) {
1505
*iflag = -1;
1506
} else if (dist_(ns, &s[ih * s_dim1 + 1], &s[il * s_dim1 + 1]) <= tol ||
1507
small) {
1508
*iflag = 0;
1509
} else {
1510
goto L20;
1511
}
1512
1513
/* end main loop, return best point */
1514
1515
i__1 = *ns;
1516
for (i__ = 1; i__ <= i__1; ++i__) {
1517
x[ips[i__]] = s[i__ + il * s_dim1];
1518
/* L60: */
1519
}
1520
return 0;
1521
} /* simplx_ */
1522
1523
/* subopt.f -- translated by f2c (version 19991025).
1524
You must link the resulting object file with the libraries:
1525
-lf2c -lm (in that order)
1526
*/
1527
1528
static int subopt_(integer *n)
1529
{
1530
1531
1532
/* Coded by Tom Rowan */
1533
/* Department of Computer Sciences */
1534
/* University of Texas at Austin */
1535
1536
/* subopt sets options for subplx. */
1537
1538
/* input */
1539
1540
/* n - problem dimension */
1541
1542
/* common */
1543
1544
1545
1546
1547
/* subroutines and functions */
1548
1549
/* fortran */
1550
1551
/* ----------------------------------------------------------- */
1552
1553
/* *********************************************************** */
1554
/* simplex method strategy parameters */
1555
/* *********************************************************** */
1556
1557
/* alpha - reflection coefficient */
1558
/* alpha .gt. 0 */
1559
1560
usubc_1.alpha = 1.;
1561
1562
/* beta - contraction coefficient */
1563
/* 0 .lt. beta .lt. 1 */
1564
1565
usubc_1.beta = .5;
1566
1567
/* gamma - expansion coefficient */
1568
/* gamma .gt. 1 */
1569
1570
usubc_1.gamma = 2.;
1571
1572
/* delta - shrinkage (massive contraction) coefficient */
1573
/* 0 .lt. delta .lt. 1 */
1574
1575
usubc_1.delta = .5;
1576
1577
/* *********************************************************** */
1578
/* subplex method strategy parameters */
1579
/* *********************************************************** */
1580
1581
/* psi - simplex reduction coefficient */
1582
/* 0 .lt. psi .lt. 1 */
1583
1584
usubc_1.psi = .25;
1585
1586
/* omega - step reduction coefficient */
1587
/* 0 .lt. omega .lt. 1 */
1588
1589
usubc_1.omega = .1;
1590
1591
/* nsmin and nsmax specify a range of subspace dimensions. */
1592
/* In addition to satisfying 1 .le. nsmin .le. nsmax .le. n, */
1593
/* nsmin and nsmax must be chosen so that n can be expressed */
1594
/* as a sum of positive integers where each of these integers */
1595
/* ns(i) satisfies nsmin .le. ns(i) .ge. nsmax. */
1596
/* Specifically, */
1597
/* nsmin*ceil(n/nsmax) .le. n must be true. */
1598
1599
/* nsmin - subspace dimension minimum */
1600
1601
usubc_1.nsmin = min(2,*n);
1602
1603
/* nsmax - subspace dimension maximum */
1604
1605
usubc_1.nsmax = min(5,*n);
1606
1607
/* *********************************************************** */
1608
/* subplex method special cases */
1609
/* *********************************************************** */
1610
/* nelder-mead simplex method with periodic restarts */
1611
/* nsmin = nsmax = n */
1612
/* *********************************************************** */
1613
/* nelder-mead simplex method */
1614
/* nsmin = nsmax = n, psi = small positive */
1615
/* *********************************************************** */
1616
1617
/* irepl, ifxsw, and bonus deal with measurement replication. */
1618
/* Objective functions subject to large amounts of noise can */
1619
/* cause an optimization method to halt at a false optimum. */
1620
/* An expensive solution to this problem is to evaluate f */
1621
/* several times at each point and return the average (or max */
1622
/* or min) of these trials as the function value. subplx */
1623
/* performs measurement replication only at the current best */
1624
/* point. The longer a point is retained as best, the more */
1625
/* accurate its function value becomes. */
1626
1627
/* The common variable nfxe contains the number of function */
1628
/* evaluations at the current best point. fxstat contains the */
1629
/* mean, max, min, and standard deviation of these trials. */
1630
1631
/* irepl - measurement replication switch */
1632
/* irepl = 0, 1, or 2 */
1633
/* = 0 : no measurement replication */
1634
/* = 1 : subplx performs measurement replication */
1635
/* = 2 : user performs measurement replication */
1636
/* (This is useful when optimizing on the mean, */
1637
/* max, or min of trials is insufficient. Common */
1638
/* variable initx is true for first function */
1639
/* evaluation. newx is true for first trial at */
1640
/* this point. The user uses subroutine fstats */
1641
/* within his objective function to maintain */
1642
/* fxstat. By monitoring newx, the user can tell */
1643
/* whether to return the function evaluation */
1644
/* (newx = .true.) or to use the new function */
1645
/* evaluation to refine the function evaluation */
1646
/* of the current best point (newx = .false.). */
1647
/* The common variable ftest gives the function */
1648
/* value that a new point must beat to be */
1649
/* considered the new best point.) */
1650
1651
usubc_1.irepl = 0;
1652
1653
/* ifxsw - measurement replication optimization switch */
1654
/* ifxsw = 1, 2, or 3 */
1655
/* = 1 : retain mean of trials as best function value */
1656
/* = 2 : retain max */
1657
/* = 3 : retain min */
1658
1659
usubc_1.ifxsw = 1;
1660
1661
/* Since the current best point will also be the most */
1662
/* accurately evaluated point whenever irepl .gt. 0, a bonus */
1663
/* should be added to the function value of the best point */
1664
/* so that the best point is not replaced by a new point */
1665
/* that only appears better because of noise. */
1666
/* subplx uses bonus to determine how many multiples of */
1667
/* fxstat(4) should be added as a bonus to the function */
1668
/* evaluation. (The bonus is adjusted automatically by */
1669
/* subplx when ifxsw or minf is changed.) */
1670
1671
/* bonus - measurement replication bonus coefficient */
1672
/* bonus .ge. 0 (normally, bonus = 0 or 1) */
1673
/* = 0 : bonus not used */
1674
/* = 1 : bonus used */
1675
1676
usubc_1.bonus = 1.;
1677
1678
/* nfstop = 0 : f(x) is not tested against fstop */
1679
/* = 1 : if f(x) has reached fstop, subplx returns */
1680
/* iflag = 2 */
1681
/* = 2 : (only valid when irepl .gt. 0) */
1682
/* if f(x) has reached fstop and */
1683
/* nfxe .gt. nfstop, subplx returns iflag = 2 */
1684
1685
usubc_1.nfstop = 0;
1686
1687
/* fstop - f target value */
1688
/* Its usage is determined by the value of nfstop. */
1689
1690
/* minf - logical switch */
1691
/* = .true. : subplx performs minimization */
1692
/* = .false. : subplx performs maximization */
1693
1694
usubc_1.minf = TRUE_;
1695
return 0;
1696
} /* subopt_ */
1697
1698
/* setstp.f -- translated by f2c (version 19991025).
1699
You must link the resulting object file with the libraries:
1700
-lf2c -lm (in that order)
1701
*/
1702
1703
static double d_sign(doublereal *x, doublereal *y)
1704
{
1705
return copysign(*x, *y);
1706
}
1707
1708
static int setstp_(integer *nsubs, integer *n, doublereal *deltax,
1709
doublereal *step)
1710
{
1711
/* System generated locals */
1712
integer i__1;
1713
doublereal d__1, d__2, d__3;
1714
1715
/* Builtin functions */
1716
/* double d_sign(doublereal *, doublereal *); */
1717
1718
/* Local variables */
1719
static integer i__;
1720
static doublereal stpfac;
1721
1722
1723
1724
/* Coded by Tom Rowan */
1725
/* Department of Computer Sciences */
1726
/* University of Texas at Austin */
1727
1728
/* setstp sets the stepsizes for the corresponding components */
1729
/* of the solution vector. */
1730
1731
/* input */
1732
1733
/* nsubs - number of subspaces */
1734
1735
/* n - number of components (problem dimension) */
1736
1737
/* deltax - vector of change in solution vector */
1738
1739
/* step - stepsizes for corresponding components of */
1740
/* solution vector */
1741
1742
/* output */
1743
1744
/* step - new stepsizes */
1745
1746
/* common */
1747
1748
1749
1750
/* local variables */
1751
1752
1753
1754
/* subroutines and functions */
1755
1756
/* blas */
1757
/* fortran */
1758
1759
/* ----------------------------------------------------------- */
1760
1761
/* set new step */
1762
1763
/* Parameter adjustments */
1764
--step;
1765
--deltax;
1766
1767
/* Function Body */
1768
if (*nsubs > 1) {
1769
/* Computing MIN */
1770
/* Computing MAX */
1771
d__3 = dasum_(n, &deltax[1], &c__1) / dasum_(n, &step[1], &c__1);
1772
d__1 = max(d__3,usubc_1.omega), d__2 = 1. / usubc_1.omega;
1773
stpfac = min(d__1,d__2);
1774
} else {
1775
stpfac = usubc_1.psi;
1776
}
1777
dscal_(n, &stpfac, &step[1], &c__1);
1778
1779
/* reorient simplex */
1780
1781
i__1 = *n;
1782
for (i__ = 1; i__ <= i__1; ++i__) {
1783
if (deltax[i__] != 0.) {
1784
step[i__] = d_sign(&step[i__], &deltax[i__]);
1785
} else {
1786
step[i__] = -step[i__];
1787
}
1788
/* L10: */
1789
}
1790
return 0;
1791
} /* setstp_ */
1792
1793
/* subplx.f -- translated by f2c (version 19991025).
1794
You must link the resulting object file with the libraries:
1795
-lf2c -lm (in that order)
1796
*/
1797
1798
static int subplx_(D_fp f, void *fdata, integer *n, doublereal *tol, integer *
1799
maxnfe, integer *mode, doublereal *scale, doublereal *x, doublereal *
1800
fx, integer *nfe, doublereal *work, integer *iwork, integer *iflag)
1801
{
1802
/* Initialized data */
1803
1804
static doublereal bnsfac[6] /* was [3][2] */ = { -1.,-2.,0.,1.,0.,2. };
1805
1806
/* System generated locals */
1807
integer i__1;
1808
doublereal d__1, d__2, d__3, d__4, d__5, d__6;
1809
1810
/* Local variables */
1811
static integer i__;
1812
static logical cmode;
1813
static integer istep;
1814
static doublereal xpscl;
1815
static integer nsubs, ipptr;
1816
static integer isptr;
1817
static integer ns, insfnl, ifsptr;
1818
static integer insptr;
1819
static integer istptr;
1820
static doublereal scl, dum;
1821
static integer ins;
1822
static doublereal sfx;
1823
1824
1825
1826
/* Coded by Tom Rowan */
1827
/* Department of Computer Sciences */
1828
/* University of Texas at Austin */
1829
1830
/* subplx uses the subplex method to solve unconstrained */
1831
/* optimization problems. The method is well suited for */
1832
/* optimizing objective functions that are noisy or are */
1833
/* discontinuous at the solution. */
1834
1835
/* subplx sets default optimization options by calling the */
1836
/* subroutine subopt. The user can override these defaults */
1837
/* by calling subopt prior to calling subplx, changing the */
1838
/* appropriate common variables, and setting the value of */
1839
/* mode as indicated below. */
1840
1841
/* By default, subplx performs minimization. */
1842
1843
/* input */
1844
1845
/* f - user supplied function f(n,x) to be optimized, */
1846
/* declared external in calling routine */
1847
1848
/* n - problem dimension */
1849
1850
/* tol - relative error tolerance for x (tol .ge. 0.) */
1851
1852
/* maxnfe - maximum number of function evaluations */
1853
1854
/* mode - integer mode switch with binary expansion */
1855
/* (bit 1) (bit 0) : */
1856
/* bit 0 = 0 : first call to subplx */
1857
/* = 1 : continuation of previous call */
1858
/* bit 1 = 0 : use default options */
1859
/* = 1 : user set options */
1860
1861
/* scale - scale and initial stepsizes for corresponding */
1862
/* components of x */
1863
/* (If scale(1) .lt. 0., */
1864
/* abs(scale(1)) is used for all components of x, */
1865
/* and scale(2),...,scale(n) are not referenced.) */
1866
1867
/* x - starting guess for optimum */
1868
1869
/* work - double precision work array of dimension .ge. */
1870
/* 2*n + nsmax*(nsmax+4) + 1 */
1871
/* (nsmax is set in subroutine subopt. */
1872
/* default: nsmax = min(5,n)) */
1873
1874
/* iwork - integer work array of dimension .ge. */
1875
/* n + int(n/nsmin) */
1876
/* (nsmin is set in subroutine subopt. */
1877
/* default: nsmin = min(2,n)) */
1878
1879
/* output */
1880
1881
/* x - computed optimum */
1882
1883
/* fx - value of f at x */
1884
1885
/* nfe - number of function evaluations */
1886
1887
/* iflag - error flag */
1888
/* = -2 : invalid input */
1889
/* = -1 : maxnfe exceeded */
1890
/* = 0 : tol satisfied */
1891
/* = 1 : limit of machine precision */
1892
/* = 2 : fstop reached (fstop usage is determined */
1893
/* by values of options minf, nfstop, and */
1894
/* irepl. default: f(x) not tested against */
1895
/* fstop) */
1896
/* iflag should not be reset between calls to */
1897
/* subplx. */
1898
1899
/* common */
1900
1901
1902
1903
1904
1905
/* local variables */
1906
1907
1908
1909
/* subroutines and functions */
1910
1911
/* blas */
1912
/* fortran */
1913
1914
/* data */
1915
1916
/* Parameter adjustments */
1917
--x;
1918
--scale;
1919
--work;
1920
--iwork;
1921
1922
/* Function Body */
1923
/* ----------------------------------------------------------- */
1924
1925
if (*mode % 2 == 0) {
1926
1927
/* first call, check input */
1928
1929
if (*n < 1) {
1930
goto L120;
1931
}
1932
if (*tol < 0.) {
1933
goto L120;
1934
}
1935
if (*maxnfe < 1) {
1936
goto L120;
1937
}
1938
if (scale[1] >= 0.) {
1939
i__1 = *n;
1940
for (i__ = 1; i__ <= i__1; ++i__) {
1941
xpscl = x[i__] + scale[i__];
1942
if (xpscl == x[i__]) {
1943
goto L120;
1944
}
1945
/* L10: */
1946
}
1947
} else {
1948
scl = abs(scale[1]);
1949
i__1 = *n;
1950
for (i__ = 1; i__ <= i__1; ++i__) {
1951
xpscl = x[i__] + scl;
1952
if (xpscl == x[i__]) {
1953
goto L120;
1954
}
1955
/* L20: */
1956
}
1957
}
1958
if (*mode / 2 % 2 == 0) {
1959
subopt_(n);
1960
} else {
1961
if (usubc_1.alpha <= 0.) {
1962
goto L120;
1963
}
1964
if (usubc_1.beta <= 0. || usubc_1.beta >= 1.) {
1965
goto L120;
1966
}
1967
if (usubc_1.gamma <= 1.) {
1968
goto L120;
1969
}
1970
if (usubc_1.delta <= 0. || usubc_1.delta >= 1.) {
1971
goto L120;
1972
}
1973
if (usubc_1.psi <= 0. || usubc_1.psi >= 1.) {
1974
goto L120;
1975
}
1976
if (usubc_1.omega <= 0. || usubc_1.omega >= 1.) {
1977
goto L120;
1978
}
1979
if (usubc_1.nsmin < 1 || usubc_1.nsmax < usubc_1.nsmin || *n <
1980
usubc_1.nsmax) {
1981
goto L120;
1982
}
1983
if (*n < ((*n - 1) / usubc_1.nsmax + 1) * usubc_1.nsmin) {
1984
goto L120;
1985
}
1986
if (usubc_1.irepl < 0 || usubc_1.irepl > 2) {
1987
goto L120;
1988
}
1989
if (usubc_1.ifxsw < 1 || usubc_1.ifxsw > 3) {
1990
goto L120;
1991
}
1992
if (usubc_1.bonus < 0.) {
1993
goto L120;
1994
}
1995
if (usubc_1.nfstop < 0) {
1996
goto L120;
1997
}
1998
}
1999
2000
/* initialization */
2001
2002
istptr = *n + 1;
2003
isptr = istptr + *n;
2004
ifsptr = isptr + usubc_1.nsmax * (usubc_1.nsmax + 3);
2005
insptr = *n + 1;
2006
if (scale[1] > 0.) {
2007
dcopy_(n, &scale[1], &c__1, &work[1], &c__1);
2008
dcopy_(n, &scale[1], &c__1, &work[istptr], &c__1);
2009
} else {
2010
dcopy_(n, &scl, &c__0, &work[1], &c__1);
2011
dcopy_(n, &scl, &c__0, &work[istptr], &c__1);
2012
}
2013
i__1 = *n;
2014
for (i__ = 1; i__ <= i__1; ++i__) {
2015
iwork[i__] = i__;
2016
/* L30: */
2017
}
2018
*nfe = 0;
2019
usubc_1.nfxe = 1;
2020
if (usubc_1.irepl == 0) {
2021
isubc_1.fbonus = 0.;
2022
} else if (usubc_1.minf) {
2023
isubc_1.fbonus = bnsfac[usubc_1.ifxsw - 1] * usubc_1.bonus;
2024
} else {
2025
isubc_1.fbonus = bnsfac[usubc_1.ifxsw + 2] * usubc_1.bonus;
2026
}
2027
if (usubc_1.nfstop == 0) {
2028
isubc_1.sfstop = 0.;
2029
} else if (usubc_1.minf) {
2030
isubc_1.sfstop = usubc_1.fstop;
2031
} else {
2032
isubc_1.sfstop = -usubc_1.fstop;
2033
}
2034
usubc_1.ftest = 0.;
2035
cmode = FALSE_;
2036
isubc_1.new__ = TRUE_;
2037
usubc_1.initx = TRUE_;
2038
evalf_((D_fp)f, fdata, &c__0, &iwork[1], &dum, n, &x[1], &sfx, nfe);
2039
usubc_1.initx = FALSE_;
2040
} else {
2041
2042
/* continuation of previous call */
2043
2044
if (*iflag == 2) {
2045
if (usubc_1.minf) {
2046
isubc_1.sfstop = usubc_1.fstop;
2047
} else {
2048
isubc_1.sfstop = -usubc_1.fstop;
2049
}
2050
cmode = TRUE_;
2051
goto L70;
2052
} else if (*iflag == -1) {
2053
cmode = TRUE_;
2054
goto L70;
2055
} else if (*iflag == 0) {
2056
cmode = FALSE_;
2057
goto L90;
2058
} else {
2059
return 0;
2060
}
2061
}
2062
2063
/* subplex loop */
2064
2065
L40:
2066
i__1 = *n;
2067
for (i__ = 1; i__ <= i__1; ++i__) {
2068
work[i__] = (d__1 = work[i__], abs(d__1));
2069
/* L50: */
2070
}
2071
sortd_(n, &work[1], &iwork[1]);
2072
partx_(n, &iwork[1], &work[1], &nsubs, &iwork[insptr]);
2073
dcopy_(n, &x[1], &c__1, &work[1], &c__1);
2074
ins = insptr;
2075
insfnl = insptr + nsubs - 1;
2076
ipptr = 1;
2077
2078
/* simplex loop */
2079
2080
L60:
2081
ns = iwork[ins];
2082
L70:
2083
simplx_((D_fp)f, fdata, n, &work[istptr], &ns, &iwork[ipptr], maxnfe, &cmode, &x[
2084
1], &sfx, nfe, &work[isptr], &work[ifsptr], iflag);
2085
cmode = FALSE_;
2086
if (*iflag != 0) {
2087
goto L110;
2088
}
2089
if (ins < insfnl) {
2090
++ins;
2091
ipptr += ns;
2092
goto L60;
2093
}
2094
2095
/* end simplex loop */
2096
2097
i__1 = *n;
2098
for (i__ = 1; i__ <= i__1; ++i__) {
2099
work[i__] = x[i__] - work[i__];
2100
/* L80: */
2101
}
2102
2103
/* check termination */
2104
2105
L90:
2106
istep = istptr;
2107
i__1 = *n;
2108
for (i__ = 1; i__ <= i__1; ++i__) {
2109
/* Computing MAX */
2110
d__4 = (d__2 = work[i__], abs(d__2)), d__5 = (d__1 = work[istep], abs(
2111
d__1)) * usubc_1.psi;
2112
/* Computing MAX */
2113
d__6 = (d__3 = x[i__], abs(d__3));
2114
if (max(d__4,d__5) / max(d__6,1.) > *tol) {
2115
setstp_(&nsubs, n, &work[1], &work[istptr]);
2116
goto L40;
2117
}
2118
++istep;
2119
/* L100: */
2120
}
2121
2122
/* end subplex loop */
2123
2124
*iflag = 0;
2125
L110:
2126
if (usubc_1.minf) {
2127
*fx = sfx;
2128
} else {
2129
*fx = -sfx;
2130
}
2131
return 0;
2132
2133
/* invalid input */
2134
2135
L120:
2136
*iflag = -2;
2137
return 0;
2138
} /* subplx_ */
2139
2140
/****************************************************************************/
2141
/****************************************************************************/
2142
2143
/* front-end for subplex routines */
2144
2145
/* Wrapper around f2c'ed subplx_ routine, for multidimensinal
2146
unconstrained optimization:
2147
2148
Parameters:
2149
2150
f: function f(n,x,fdata) to be optimized
2151
n: problem dimension
2152
x[n]: (input) starting guess position, (output) computed minimum
2153
fdata: data pointer passed to f
2154
tol: relative error tolerance for x
2155
maxnfe: maximum number of function evaluations
2156
scale[n]: (input) scale & initial stepsizes for components of x
2157
(if *scale < 0, |*scale| is used for all components)
2158
nfe: (output) number of function evaluations
2159
errflag: (output)
2160
= -2 : invalid input
2161
= -1 : maxnfe exceeded
2162
= 0 : tol satisfied
2163
= 1 : limit of machine precision
2164
= 2 : fstop reached (fstop usage is determined by values of
2165
options minf, nfstop, and irepl. default: f(x) not
2166
tested against fstop)
2167
2168
Return value: value of f at minimum.
2169
*/
2170
number subplex(subplex_func f, number *x, integer n, void *fdata,
2171
number tol, integer maxnfe,
2172
number fmin, boolean use_fmin,
2173
number *scale,
2174
integer *nfe, integer *errflag)
2175
{
2176
integer mode = 0, *iwork, nsmax, nsmin;
2177
number *work, fx;
2178
2179
nsmax = min(5,n);
2180
nsmin = min(2,n);
2181
work = (number*) malloc(sizeof(number) * (2*n + nsmax*(nsmax+4) + 1));
2182
iwork = (integer*) malloc(sizeof(integer) * (n + n/nsmin + 1));
2183
if (!work || !iwork) {
2184
fprintf(stderr, "subplex: error, out of memory!\n");
2185
exit(EXIT_FAILURE);
2186
}
2187
2188
if (use_fmin) { /* stop when fmin is reached */
2189
subopt_(&n);
2190
usubc_1.nfstop = 1;
2191
usubc_1.fstop = fmin;
2192
mode = 2;
2193
}
2194
2195
subplx_(f,fdata, &n,
2196
&tol, &maxnfe, &mode,
2197
scale, x,
2198
&fx, nfe,
2199
work, iwork, errflag);
2200
2201
free(iwork);
2202
free(work);
2203
2204
return fx;
2205
}
2206
2207
number f_scm_wrapper(integer n, number *x, void *f_scm_p)
2208
{
2209
SCM *f_scm = (SCM *) f_scm_p;
2210
return gh_scm2double(gh_call1(*f_scm, make_number_list(n, x)));
2211
}
2212
2213
/* Scheme-callable wrapper for subplex() function, above. */
2214
SCM subplex_scm(SCM f_scm, SCM x_scm,
2215
SCM tol_scm, SCM maxnfe_scm,
2216
SCM fmin_scm, SCM use_fmin_scm,
2217
SCM scale_scm)
2218
{
2219
number *x, tol, *scale, fx, fmin;
2220
integer i, n, maxnfe, nfe, errflag, scale_len;
2221
boolean use_fmin;
2222
SCM retval;
2223
2224
n = list_length(x_scm);
2225
tol = fabs(gh_scm2double(tol_scm));
2226
maxnfe = gh_scm2int(maxnfe_scm);
2227
fmin = gh_scm2double(fmin_scm);
2228
use_fmin = gh_scm2bool(use_fmin_scm);
2229
2230
scale_len = list_length(scale_scm);
2231
if (scale_len != 1 && scale_len != n) {
2232
fprintf(stderr, "subplex: invalid scale argument length %d\n",
2233
scale_len);
2234
return SCM_UNDEFINED;
2235
}
2236
2237
x = (number*) malloc(sizeof(number) * n);
2238
scale = (number*) malloc(sizeof(number) * scale_len);
2239
if (!x || !scale) {
2240
fprintf(stderr, "subplex: error, out of memory!\n");
2241
exit(EXIT_FAILURE);
2242
}
2243
2244
for (i = 0; i < n; ++i)
2245
x[i] = number_list_ref(x_scm, i);
2246
for (i = 0; i < scale_len; ++i)
2247
scale[i] = fabs(number_list_ref(scale_scm, i));
2248
if (scale_len == 1 && scale_len != n)
2249
*scale *= -1;
2250
2251
fx = subplex(f_scm_wrapper, x, n, &f_scm,
2252
tol, maxnfe,
2253
fmin, use_fmin,
2254
scale,
2255
&nfe, &errflag);
2256
2257
switch (errflag) {
2258
case -2:
2259
fprintf(stderr, "subplex error: invalid inputs\n");
2260
return SCM_UNDEFINED;
2261
case -1:
2262
fprintf(stderr, "subplex warning: max # iterations exceeded\n");
2263
break;
2264
case 1:
2265
fprintf(stderr, "subplex warning: machine precision reached\n");
2266
break;
2267
case 2:
2268
fprintf(stderr, "subplex warning: fstop reached\n");
2269
break;
2270
}
2271
2272
retval = gh_cons(make_number_list(n, x), gh_double2scm(fx));
2273
2274
free(scale);
2275
free(x);
2276
2277
return retval;
2278
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1