1 /* DIRsubrout.f -- translated by f2c (version 20050501).
3 f2c output hand-cleaned by SGJ (August 2007).
7 #include "direct-internal.h"
9 /* Table of constant values */
11 static integer c__1 = 1;
12 static integer c__32 = 32;
13 static integer c__0 = 0;
15 /* +-----------------------------------------------------------------------+ */
16 /* | INTEGER Function DIRGetlevel | */
17 /* | Returns the level of the hyperrectangle. Depending on the value of the| */
18 /* | global variable JONES. IF JONES equals 0, the level is given by | */
19 /* | kN + p, where the rectangle has p sides with a length of| */
20 /* | 1/3^(k+1), and N-p sides with a length of 1/3^k. | */
21 /* | If JONES equals 1, the level is the power of 1/3 of the length of the | */
22 /* | longest side hyperrectangle. | */
25 /* | the maximal length | */
27 /* | pos -- the position of the midpoint in the array length | */
28 /* | length -- the array with the dimensions | */
29 /* | maxfunc -- the leading dimension of length | */
30 /* | n -- the dimension of the problem | */
32 /* +-----------------------------------------------------------------------+ */
33 integer direct_dirgetlevel_(integer *pos, integer *length, integer *maxfunc, integer
36 /* System generated locals */
37 integer length_dim1, length_offset, ret_val, i__1;
40 integer i__, k, p, help;
42 /* JG 09/15/00 Added variable JONES (see above) */
43 /* Parameter adjustments */
45 length_offset = 1 + length_dim1;
46 length -= length_offset;
50 help = length[*pos * length_dim1 + 1];
54 for (i__ = 2; i__ <= i__1; ++i__) {
55 if (length[i__ + *pos * length_dim1] < k) {
56 k = length[i__ + *pos * length_dim1];
58 if (length[i__ + *pos * length_dim1] == help) {
64 ret_val = k * *n + *n - p;
69 help = length[*pos * length_dim1 + 1];
71 for (i__ = 2; i__ <= i__1; ++i__) {
72 if (length[i__ + *pos * length_dim1] < help) {
73 help = length[i__ + *pos * length_dim1];
82 /* +-----------------------------------------------------------------------+ */
83 /* | Program : Direct.f (subfile DIRsubrout.f) | */
84 /* | Last modified : 07-16-2001 | */
85 /* | Written by : Joerg Gablonsky | */
86 /* | Subroutines used by the algorithm DIRECT. | */
87 /* +-----------------------------------------------------------------------+ */
88 /* +-----------------------------------------------------------------------+ */
89 /* | SUBROUTINE DIRChoose | */
90 /* | Decide, which is the next sampling point. | */
91 /* | Changed 09/25/00 JG | */
92 /* | Added maxdiv to call and changed S to size maxdiv. | */
93 /* | Changed 01/22/01 JG | */
94 /* | Added Ifeasiblef to call to keep track if a feasible point has| */
96 /* | Changed 07/16/01 JG | */
97 /* | Changed if statement to prevent run-time errors. |
99 /* +-----------------------------------------------------------------------+ */
100 /* Subroutine */ void direct_dirchoose_(integer *anchor, integer *s, integer *actdeep,
101 doublereal *f, doublereal *minf, doublereal epsrel, doublereal epsabs, doublereal *thirds,
102 integer *maxpos, integer *length, integer *maxfunc, const integer *maxdeep,
103 const integer *maxdiv, integer *n, FILE *logfile,
104 integer *cheat, doublereal *kmax, integer *ifeasiblef, integer jones)
106 /* System generated locals */
107 integer s_dim1, s_offset, length_dim1, length_offset, i__1;
109 /* Local variables */
111 doublereal helplower;
113 doublereal helpgreater;
114 integer novaluedeep = 0;
118 /* Parameter adjustments */
122 s_offset = 1 + s_dim1;
125 length_offset = 1 + length_dim1;
126 length -= length_offset;
132 if (*ifeasiblef >= 1) {
134 for (j = 0; j <= i__1; ++j) {
136 s[k + s_dim1] = anchor[j];
137 s[k + (s_dim1 << 1)] = direct_dirgetlevel_(&s[k + s_dim1], &length[
138 length_offset], maxfunc, n, jones);
149 for (j = 0; j <= i__1; ++j) {
151 s[k + s_dim1] = anchor[j];
152 s[k + (s_dim1 << 1)] = direct_dirgetlevel_(&s[k + s_dim1], &length[
153 length_offset], maxfunc, n, jones);
160 if (anchor[-1] > 0) {
161 novalue = anchor[-1];
162 novaluedeep = direct_dirgetlevel_(&novalue, &length[length_offset], maxfunc,
167 for (j = k - 1; j <= i__1; ++j) {
171 for (j = *maxpos; j >= 1; --j) {
174 j___ = s[j + s_dim1];
176 for (i__ = 1; i__ <= i__1; ++i__) {
177 i___ = s[i__ + s_dim1];
178 /* +-----------------------------------------------------------------------+ */
179 /* | JG 07/16/01 Changed IF statement into two to prevent run-time errors | */
180 /* | which could occur if the compiler checks the second | */
181 /* | expression in an .AND. statement although the first | */
182 /* | statement is already not true. | */
183 /* +-----------------------------------------------------------------------+ */
184 if (i___ > 0 && ! (i__ == j)) {
185 if (f[(i___ << 1) + 2] <= 1.) {
186 help2 = thirds[s[i__ + (s_dim1 << 1)]] - thirds[s[j + (
188 help2 = (f[(i___ << 1) + 1] - f[(j___ << 1) + 1]) / help2;
191 fprintf(logfile, "thirds > 0, help2 <= 0\n");
194 if (help2 < helplower) {
196 fprintf(logfile, "helplower = %g\n", help2);
204 for (i__ = j + 1; i__ <= i__1; ++i__) {
205 i___ = s[i__ + s_dim1];
206 /* +-----------------------------------------------------------------------+ */
207 /* | JG 07/16/01 Changed IF statement into two to prevent run-time errors | */
208 /* | which could occur if the compiler checks the second | */
209 /* | expression in an .AND. statement although the first | */
210 /* | statement is already not true. | */
211 /* +-----------------------------------------------------------------------+ */
212 if (i___ > 0 && ! (i__ == j)) {
213 if (f[(i___ << 1) + 2] <= 1.) {
214 help2 = thirds[s[i__ + (s_dim1 << 1)]] - thirds[s[j + (
216 help2 = (f[(i___ << 1) + 1] - f[(j___ << 1) + 1]) / help2;
219 fprintf(logfile, "thirds < 0, help2 <= 0\n");
222 if (help2 > helpgreater) {
224 fprintf(logfile, "helpgreater = %g\n", help2);
231 if (helplower > 1e20 && helpgreater > 0.) {
232 helplower = helpgreater;
235 if (helpgreater <= helplower) {
236 if (*cheat == 1 && helplower > *kmax) {
239 if (f[(j___ << 1) + 1] - helplower * thirds[s[j + (s_dim1 << 1)]] >
240 MIN(*minf - epsrel * fabs(*minf),
243 fprintf(logfile, "> minf - epslminfl\n");
248 fprintf(logfile, "helpgreater > helplower: %g %g %g\n",
249 helpgreater, helplower, helpgreater - helplower);
260 s[*maxpos + s_dim1] = novalue;
261 s[*maxpos + (s_dim1 << 1)] = novaluedeep;
265 /* +-----------------------------------------------------------------------+ */
266 /* | SUBROUTINE DIRDoubleInsert | */
267 /* | Routine to make sure that if there are several potential optimal | */
268 /* | hyperrectangles of the same level (i.e. hyperrectangles that have| */
269 /* | the same level and the same function value at the center), all of| */
270 /* | them are divided. This is the way as originally described in | */
271 /* | Jones et.al. | */
272 /* | JG 07/16/01 Added errorflag to calling sequence. We check if more | */
273 /* | we reach the capacity of the array S. If this happens, we | */
274 /* | return to the main program with an error. | */
275 /* +-----------------------------------------------------------------------+ */
276 /* Subroutine */ void direct_dirdoubleinsert_(integer *anchor, integer *s, integer *
277 maxpos, integer *point, doublereal *f, const integer *maxdeep, integer *
278 maxfunc, const integer *maxdiv, integer *ierror)
280 /* System generated locals */
281 integer s_dim1, s_offset, i__1;
283 /* Local variables */
284 integer i__, oldmaxpos, pos, help, iflag, actdeep;
286 /* +-----------------------------------------------------------------------+ */
287 /* | JG 07/16/01 Added flag to prevent run time-errors on some systems. | */
288 /* +-----------------------------------------------------------------------+ */
289 /* Parameter adjustments */
294 s_offset = 1 + s_dim1;
300 for (i__ = 1; i__ <= i__1; ++i__) {
301 if (s[i__ + s_dim1] > 0) {
302 actdeep = s[i__ + (s_dim1 << 1)];
303 help = anchor[actdeep];
306 /* +-----------------------------------------------------------------------+ */
307 /* | JG 07/16/01 Added flag to prevent run time-errors on some systems. On | */
308 /* | some systems the second conditions in an AND statement is | */
309 /* | evaluated even if the first one is already not true. | */
310 /* +-----------------------------------------------------------------------+ */
311 while(pos > 0 && iflag == 0) {
312 if (f[(pos << 1) + 1] - f[(help << 1) + 1] <= 1e-13) {
313 if (*maxpos < *maxdiv) {
315 s[*maxpos + s_dim1] = pos;
316 s[*maxpos + (s_dim1 << 1)] = actdeep;
319 /* +-----------------------------------------------------------------------+ */
320 /* | JG 07/16/01 Maximum number of elements possible in S has been reached!| */
321 /* +-----------------------------------------------------------------------+ */
332 } /* dirdoubleinsert_ */
334 /* +-----------------------------------------------------------------------+ */
335 /* | INTEGER Function GetmaxDeep | */
336 /* | function to get the maximal length (1/length) of the n-dimensional | */
337 /* | rectangle with midpoint pos. | */
339 /* | On Return : | */
340 /* | the maximal length | */
342 /* | pos -- the position of the midpoint in the array length | */
343 /* | length -- the array with the dimensions | */
344 /* | maxfunc -- the leading dimension of length | */
345 /* | n -- the dimension of the problem | */
347 /* +-----------------------------------------------------------------------+ */
348 integer direct_dirgetmaxdeep_(integer *pos, integer *length, integer *maxfunc,
351 /* System generated locals */
352 integer length_dim1, length_offset, i__1, i__2, i__3;
354 /* Local variables */
357 /* Parameter adjustments */
359 length_offset = 1 + length_dim1;
360 length -= length_offset;
363 help = length[*pos * length_dim1 + 1];
365 for (i__ = 2; i__ <= i__1; ++i__) {
367 i__2 = help, i__3 = length[i__ + *pos * length_dim1];
368 help = MIN(i__2,i__3);
372 } /* dirgetmaxdeep_ */
374 static integer isinbox_(doublereal *x, doublereal *a, doublereal *b, integer *n,
377 /* System generated locals */
378 integer ret_val, i__1;
380 /* Local variables */
381 integer outofbox, i__;
383 /* Parameter adjustments */
391 for (i__ = 1; i__ <= i__1; ++i__) {
392 if (a[i__] > x[i__] || b[i__] < x[i__]) {
403 /* +-----------------------------------------------------------------------+ */
404 /* | JG Added 09/25/00 | */
406 /* | SUBROUTINE DIRResortlist | */
408 /* | Resort the list so that the infeasible point is in the list with the | */
409 /* | replaced value. | */
410 /* +-----------------------------------------------------------------------+ */
411 /* Subroutine */ static void dirresortlist_(integer *replace, integer *anchor,
412 doublereal *f, integer *point, integer *length, integer *n, integer *
413 maxfunc, integer *maxdim, const integer *maxdeep, FILE *logfile,
416 /* System generated locals */
417 integer length_dim1, length_offset, i__1;
419 /* Local variables */
423 /* +-----------------------------------------------------------------------+ */
424 /* | Get the length of the hyper rectangle with infeasible mid point and | */
425 /* | Index of the corresponding list. | */
426 /* +-----------------------------------------------------------------------+ */
427 /* JG 09/25/00 Replaced with DIRgetlevel */
428 /* l = DIRgetmaxDeep(replace,length,maxfunc,n) */
429 /* Parameter adjustments */
433 length_offset = 1 + length_dim1;
434 length -= length_offset;
438 l = direct_dirgetlevel_(replace, &length[length_offset], maxfunc, n, jones);
440 /* +-----------------------------------------------------------------------+ */
441 /* | If the hyper rectangle with infeasibel midpoint is already the start | */
442 /* | of the list, give out message, nothing to do. | */
443 /* +-----------------------------------------------------------------------+ */
444 if (*replace == start) {
445 /* write(logfile,*) 'No resorting of list necessarry, since new ', */
446 /* + 'point is already anchor of list .',l */
448 /* +-----------------------------------------------------------------------+ */
449 /* | Take the hyper rectangle with infeasible midpoint out of the list. | */
450 /* +-----------------------------------------------------------------------+ */
453 for (i__ = 1; i__ <= i__1; ++i__) {
454 if (point[pos] == *replace) {
455 point[pos] = point[*replace];
462 fprintf(logfile, "Error in DIRREsortlist: "
463 "We went through the whole list\n"
464 "and could not find the point to replace!!\n");
469 /* +-----------------------------------------------------------------------+ */
470 /* | If the anchor of the list has a higher value than the value of a | */
471 /* | nearby point, put the infeasible point at the beginning of the list. | */
472 /* +-----------------------------------------------------------------------+ */
474 if (f[(start << 1) + 1] > f[(*replace << 1) + 1]) {
475 anchor[l] = *replace;
476 point[*replace] = start;
477 /* write(logfile,*) 'Point is replacing current anchor for ' */
478 /* + , 'this list ',l,replace,start */
480 /* +-----------------------------------------------------------------------+ */
481 /* | Insert the point into the list according to its (replaced) function | */
483 /* +-----------------------------------------------------------------------+ */
486 for (i__ = 1; i__ <= i__1; ++i__) {
487 /* +-----------------------------------------------------------------------+ */
488 /* | The point has to be added at the end of the list. | */
489 /* +-----------------------------------------------------------------------+ */
490 if (point[pos] == 0) {
491 point[*replace] = point[pos];
492 point[pos] = *replace;
493 /* write(logfile,*) 'Point is added at the end of the ' */
494 /* + , 'list ',l, replace */
497 if (f[(point[pos] << 1) + 1] > f[(*replace << 1) + 1]) {
498 point[*replace] = point[pos];
499 point[pos] = *replace;
500 /* write(logfile,*) 'There are points with a higher ' */
501 /* + ,'f-value in the list ',l,replace, pos */
512 } /* dirresortlist_ */
514 /* +-----------------------------------------------------------------------+ */
515 /* | JG Added 09/25/00 | */
516 /* | SUBROUTINE DIRreplaceInf | */
518 /* | Find out if there are infeasible points which are near feasible ones. | */
519 /* | If this is the case, replace the function value at the center of the | */
520 /* | hyper rectangle by the lowest function value of a nearby function. | */
521 /* +-----------------------------------------------------------------------+ */
522 /* Subroutine */ void direct_dirreplaceinf_(integer *free, integer *freeold,
523 doublereal *f, doublereal *c__, doublereal *thirds, integer *length,
524 integer *anchor, integer *point, doublereal *c1, doublereal *c2,
525 integer *maxfunc, const integer *maxdeep, integer *maxdim, integer *n,
526 FILE *logfile, doublereal *fmax, integer jones)
528 /* System generated locals */
529 integer c_dim1, c_offset, length_dim1, length_offset, i__1, i__2, i__3;
530 doublereal d__1, d__2;
532 /* Local variables */
533 doublereal a[32], b[32];
534 integer i__, j, k, l;
535 doublereal x[32], sidelength;
538 /* +-----------------------------------------------------------------------+ */
539 /* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
540 /* +-----------------------------------------------------------------------+ */
541 /* Parameter adjustments */
545 length_dim1 = *maxdim;
546 length_offset = 1 + length_dim1;
547 length -= length_offset;
549 c_offset = 1 + c_dim1;
556 for (i__ = 1; i__ <= i__1; ++i__) {
557 if (f[(i__ << 1) + 2] > 0.) {
558 /* +-----------------------------------------------------------------------+ */
559 /* | Get the maximum side length of the hyper rectangle and then set the | */
560 /* | new side length to this lengths times the growth factor. | */
561 /* +-----------------------------------------------------------------------+ */
562 help = direct_dirgetmaxdeep_(&i__, &length[length_offset], maxfunc, n);
563 sidelength = thirds[help] * 2.;
564 /* +-----------------------------------------------------------------------+ */
565 /* | Set the Center and the upper and lower bounds of the rectangles. | */
566 /* +-----------------------------------------------------------------------+ */
568 for (j = 1; j <= i__2; ++j) {
569 sidelength = thirds[length[i__ + j * length_dim1]];
570 a[j - 1] = c__[j + i__ * c_dim1] - sidelength;
571 b[j - 1] = c__[j + i__ * c_dim1] + sidelength;
574 /* +-----------------------------------------------------------------------+ */
575 /* | The function value is reset to 'Inf', since it may have been changed | */
576 /* | in an earlier iteration and now the feasible point which was close | */
577 /* | is not close anymore (since the hyper rectangle surrounding the | */
578 /* | current point may have shrunk). | */
579 /* +-----------------------------------------------------------------------+ */
580 f[(i__ << 1) + 1] = HUGE_VAL;
581 f[(i__ << 1) + 2] = 2.;
582 /* +-----------------------------------------------------------------------+ */
583 /* | Check if any feasible point is near this infeasible point. | */
584 /* +-----------------------------------------------------------------------+ */
586 for (k = 1; k <= i__2; ++k) {
587 /* +-----------------------------------------------------------------------+ */
588 /* | If the point k is feasible, check if it is near. | */
589 /* +-----------------------------------------------------------------------+ */
590 if (f[(k << 1) + 2] == 0.) {
591 /* +-----------------------------------------------------------------------+ */
592 /* | Copy the coordinates of the point k into x. | */
593 /* +-----------------------------------------------------------------------+ */
595 for (l = 1; l <= i__3; ++l) {
596 x[l - 1] = c__[l + k * c_dim1];
599 /* +-----------------------------------------------------------------------+ */
600 /* | Check if the point k is near the infeasible point, if so, replace the | */
602 /* +-----------------------------------------------------------------------+ */
603 if (isinbox_(x, a, b, n, &c__32) == 1) {
605 d__1 = f[(i__ << 1) + 1], d__2 = f[(k << 1) + 1];
606 f[(i__ << 1) + 1] = MIN(d__1,d__2);
607 f[(i__ << 1) + 2] = 1.;
612 if (f[(i__ << 1) + 2] == 1.) {
613 f[(i__ << 1) + 1] += (d__1 = f[(i__ << 1) + 1], fabs(d__1)) *
616 for (l = 1; l <= i__2; ++l) {
617 x[l - 1] = c__[l + i__ * c_dim1] * c1[l] + c__[l + i__ *
621 dirresortlist_(&i__, &anchor[-1], &f[3], &point[1],
622 &length[length_offset], n,
623 maxfunc, maxdim, maxdeep, logfile, jones);
625 /* +-----------------------------------------------------------------------+ */
626 /* | JG 01/22/01 | */
627 /* | Replaced fixed value for infeasible points with maximum value found, | */
628 /* | increased by 1. | */
629 /* +-----------------------------------------------------------------------+ */
630 if (! (*fmax == f[(i__ << 1) + 1])) {
632 d__1 = *fmax + 1., d__2 = f[(i__ << 1) + 1];
633 f[(i__ << 1) + 1] = MAX(d__1,d__2);
640 } /* dirreplaceinf_ */
642 /* +-----------------------------------------------------------------------+ */
643 /* | SUBROUTINE DIRInsert | */
644 /* +-----------------------------------------------------------------------+ */
645 /* Subroutine */ static void dirinsert_(integer *start, integer *ins, integer *point,
646 doublereal *f, integer *maxfunc)
648 /* System generated locals */
651 /* Local variables */
654 /* JG 09/17/00 Rewrote this routine. */
655 /* DO 10,i = 1,maxfunc */
656 /* IF (f(ins,1) .LT. f(point(start),1)) THEN */
657 /* help = point(start) */
658 /* point(start) = ins */
659 /* point(ins) = help */
662 /* IF (point(start) .EQ. 0) THEN */
663 /* point(start) = ins */
667 /* start = point(start) */
670 /* Parameter adjustments */
676 for (i__ = 1; i__ <= i__1; ++i__) {
677 if (point[*start] == 0) {
678 point[*start] = *ins;
681 } else if (f[(*ins << 1) + 1] < f[(point[*start] << 1) + 1]) {
682 help = point[*start];
683 point[*start] = *ins;
687 *start = point[*start];
692 /* +-----------------------------------------------------------------------+ */
693 /* | SUBROUTINE DIRInsertList | */
694 /* | Changed 02-24-2000 | */
695 /* | Got rid of the distinction between feasible and infeasible points| */
696 /* | I could do this since infeasible points get set to a high | */
697 /* | function value, which may be replaced by a function value of a | */
698 /* | nearby function at the end of the main loop. | */
699 /* +-----------------------------------------------------------------------+ */
700 /* Subroutine */ void direct_dirinsertlist_(integer *new__, integer *anchor, integer *
701 point, doublereal *f, integer *maxi, integer *length, integer *
702 maxfunc, const integer *maxdeep, integer *n, integer *samp,
705 /* System generated locals */
706 integer length_dim1, length_offset, i__1;
708 /* Local variables */
711 integer pos1, pos2, deep;
713 /* JG 09/24/00 Changed this to Getlevel */
714 /* Parameter adjustments */
719 length_offset = 1 + length_dim1;
720 length -= length_offset;
724 for (j = 1; j <= i__1; ++j) {
727 *new__ = point[pos2];
728 /* JG 09/24/00 Changed this to Getlevel */
729 /* deep = DIRGetMaxdeep(pos1,length,maxfunc,n) */
730 deep = direct_dirgetlevel_(&pos1, &length[length_offset], maxfunc, n, jones);
731 if (anchor[deep] == 0) {
732 if (f[(pos2 << 1) + 1] < f[(pos1 << 1) + 1]) {
742 if (f[(pos2 << 1) + 1] < f[(pos1 << 1) + 1]) {
743 if (f[(pos2 << 1) + 1] < f[(pos << 1) + 1]) {
745 /* JG 08/30/00 Fixed bug. Sorting was not correct when */
746 /* f(1,pos2) < f(1,pos1) < f(1,pos) */
747 if (f[(pos1 << 1) + 1] < f[(pos << 1) + 1]) {
752 dirinsert_(&pos, &pos1, &point[1], &f[3], maxfunc);
755 dirinsert_(&pos, &pos2, &point[1], &f[3], maxfunc);
756 dirinsert_(&pos, &pos1, &point[1], &f[3], maxfunc);
759 if (f[(pos1 << 1) + 1] < f[(pos << 1) + 1]) {
760 /* JG 08/30/00 Fixed bug. Sorting was not correct when */
761 /* f(pos1,1) < f(pos2,1) < f(pos,1) */
763 if (f[(pos << 1) + 1] < f[(pos2 << 1) + 1]) {
765 dirinsert_(&pos, &pos2, &point[1], &f[3], maxfunc);
771 dirinsert_(&pos, &pos1, &point[1], &f[3], maxfunc);
772 dirinsert_(&pos, &pos2, &point[1], &f[3], maxfunc);
778 /* JG 09/24/00 Changed this to Getlevel */
779 /* deep = DIRGetMaxdeep(samp,length,maxfunc,n) */
780 deep = direct_dirgetlevel_(samp, &length[length_offset], maxfunc, n, jones);
782 if (f[(*samp << 1) + 1] < f[(pos << 1) + 1]) {
783 anchor[deep] = *samp;
786 dirinsert_(&pos, samp, &point[1], &f[3], maxfunc);
788 } /* dirinsertlist_ */
790 /* +-----------------------------------------------------------------------+ */
791 /* | SUBROUTINE DIRInsertList2 (Old way to do it.) | */
792 /* +-----------------------------------------------------------------------+ */
793 /* Subroutine */ static void dirinsertlist_2__(integer *start, integer *j, integer *k,
794 integer *list2, doublereal *w, integer *maxi, integer *n)
796 /* System generated locals */
797 integer list2_dim1, list2_offset, i__1;
799 /* Local variables */
802 /* Parameter adjustments */
805 list2_offset = 1 + list2_dim1;
806 list2 -= list2_offset;
811 list2[*j + list2_dim1] = 0;
815 if (w[*start] > w[*j]) {
816 list2[*j + list2_dim1] = *start;
820 for (i__ = 1; i__ <= i__1; ++i__) {
821 if (list2[pos + list2_dim1] == 0) {
822 list2[*j + list2_dim1] = 0;
823 list2[pos + list2_dim1] = *j;
826 if (w[*j] < w[list2[pos + list2_dim1]]) {
827 list2[*j + list2_dim1] = list2[pos + list2_dim1];
828 list2[pos + list2_dim1] = *j;
832 pos = list2[pos + list2_dim1];
837 list2[*j + (list2_dim1 << 1)] = *k;
838 } /* dirinsertlist_2__ */
840 /* +-----------------------------------------------------------------------+ */
841 /* | SUBROUTINE DIRSearchmin | */
842 /* | Search for the minimum in the list. ! */
843 /* +-----------------------------------------------------------------------+ */
844 /* Subroutine */ static void dirsearchmin_(integer *start, integer *list2, integer *
845 pos, integer *k, integer *n)
847 /* System generated locals */
848 integer list2_dim1, list2_offset;
850 /* Parameter adjustments */
852 list2_offset = 1 + list2_dim1;
853 list2 -= list2_offset;
857 *pos = list2[*start + (list2_dim1 << 1)];
858 *start = list2[*start + list2_dim1];
859 } /* dirsearchmin_ */
861 /* +-----------------------------------------------------------------------+ */
862 /* | SUBROUTINE DIRSamplepoints | */
863 /* | Subroutine to sample the new points. | */
864 /* +-----------------------------------------------------------------------+ */
865 /* Subroutine */ void direct_dirsamplepoints_(doublereal *c__, integer *arrayi,
866 doublereal *delta, integer *sample, integer *start, integer *length,
867 FILE *logfile, doublereal *f, integer *free,
868 integer *maxi, integer *point, doublereal *x, doublereal *l,
869 doublereal *minf, integer *minpos, doublereal *u, integer *n,
870 integer *maxfunc, const integer *maxdeep, integer *oops)
872 /* System generated locals */
873 integer length_dim1, length_offset, c_dim1, c_offset, i__1, i__2;
875 /* Local variables */
879 /* Parameter adjustments */
887 length_offset = 1 + length_dim1;
888 length -= length_offset;
890 c_offset = 1 + c_dim1;
897 i__1 = *maxi + *maxi;
898 for (k = 1; k <= i__1; ++k) {
900 for (j = 1; j <= i__2; ++j) {
901 length[j + *free * length_dim1] = length[j + *sample *
903 c__[j + *free * c_dim1] = c__[j + *sample * c_dim1];
907 *free = point[*free];
910 fprintf(logfile, "Error, no more free positions! "
911 "Increase maxfunc!\n");
920 for (j = 1; j <= i__1; ++j) {
921 c__[arrayi[j] + pos * c_dim1] = c__[arrayi[j] + *sample * c_dim1] + *
924 c__[arrayi[j] + pos * c_dim1] = c__[arrayi[j] + *sample * c_dim1] - *
930 } /* dirsamplepoints_ */
932 /* +-----------------------------------------------------------------------+ */
933 /* | SUBROUTINE DIRDivide | */
934 /* | Subroutine to divide the hyper rectangles according to the rules. | */
935 /* | Changed 02-24-2000 | */
936 /* | Replaced if statement by min (line 367) | */
937 /* +-----------------------------------------------------------------------+ */
938 /* Subroutine */ void direct_dirdivide_(integer *new__, integer *currentlength,
939 integer *length, integer *point, integer *arrayi, integer *sample,
940 integer *list2, doublereal *w, integer *maxi, doublereal *f, integer *
941 maxfunc, const integer *maxdeep, integer *n)
943 /* System generated locals */
944 integer length_dim1, length_offset, list2_dim1, list2_offset, i__1, i__2;
945 doublereal d__1, d__2;
947 /* Local variables */
948 integer i__, j, k, pos, pos2;
952 /* Parameter adjustments */
957 list2_offset = 1 + list2_dim1;
958 list2 -= list2_offset;
961 length_offset = 1 + length_dim1;
962 length -= length_offset;
968 for (i__ = 1; i__ <= i__1; ++i__) {
970 w[j] = f[(pos << 1) + 1];
974 d__1 = f[(pos << 1) + 1], d__2 = w[j];
975 w[j] = MIN(d__1,d__2);
977 dirinsertlist_2__(&start, &j, &k, &list2[list2_offset], &w[1], maxi,
983 for (j = 1; j <= i__1; ++j) {
984 dirsearchmin_(&start, &list2[list2_offset], &pos, &k, n);
986 length[k + *sample * length_dim1] = *currentlength + 1;
987 i__2 = *maxi - j + 1;
988 for (i__ = 1; i__ <= i__2; ++i__) {
989 length[k + pos * length_dim1] = *currentlength + 1;
991 length[k + pos * length_dim1] = *currentlength + 1;
992 /* JG 07/10/01 pos2 = 0 at the end of the 30-loop. Since we end */
993 /* the loop now, we do not need to reassign pos and pos2. */
995 pos = list2[pos2 + (list2_dim1 << 1)];
996 pos2 = list2[pos2 + list2_dim1];
1004 /* +-----------------------------------------------------------------------+ */
1006 /* | SUBROUTINE DIRINFCN | */
1008 /* | Subroutine DIRinfcn unscales the variable x for use in the | */
1009 /* | user-supplied function evaluation subroutine fcn. After fcn returns | */
1010 /* | to DIRinfcn, DIRinfcn then rescales x for use by DIRECT. | */
1014 /* | fcn -- The argument containing the name of the user-supplied | */
1015 /* | subroutine that returns values for the function to be | */
1016 /* | minimized. | */
1018 /* | x -- A double-precision vector of length n. The point at | */
1019 /* | which the derivative is to be evaluated. | */
1021 /* | xs1 -- A double-precision vector of length n. Used for | */
1022 /* | scaling and unscaling the vector x by DIRinfcn. | */
1024 /* | xs2 -- A double-precision vector of length n. Used for | */
1025 /* | scaling and unscaling the vector x by DIRinfcn. | */
1027 /* | n -- An integer. The dimension of the problem. | */
1028 /* | kret -- An Integer. If kret = 1, the point is infeasible, | */
1029 /* | kret = -1, bad problem set up, | */
1030 /* | kret = 0, feasible. | */
1034 /* | f -- A double-precision scalar. | */
1036 /* | Subroutines and Functions | */
1038 /* | The subroutine whose name is passed through the argument fcn. | */
1040 /* +-----------------------------------------------------------------------+ */
1041 /* Subroutine */ void direct_dirinfcn_(fp fcn, doublereal *x, doublereal *c1,
1042 doublereal *c2, integer *n, doublereal *f, integer *flag__,
1045 /* System generated locals */
1048 /* Local variables */
1051 /* +-----------------------------------------------------------------------+ */
1052 /* | Variables to pass user defined data to the function to be optimized. | */
1053 /* +-----------------------------------------------------------------------+ */
1054 /* +-----------------------------------------------------------------------+ */
1055 /* | Unscale the variable x. | */
1056 /* +-----------------------------------------------------------------------+ */
1057 /* Parameter adjustments */
1064 for (i__ = 1; i__ <= i__1; ++i__) {
1065 x[i__] = (x[i__] + c2[i__]) * c1[i__];
1068 /* +-----------------------------------------------------------------------+ */
1069 /* | Call the function-evaluation subroutine fcn. | */
1070 /* +-----------------------------------------------------------------------+ */
1072 *f = fcn(*n, &x[1], flag__, fcn_data);
1073 /* +-----------------------------------------------------------------------+ */
1074 /* | Rescale the variable x. | */
1075 /* +-----------------------------------------------------------------------+ */
1077 for (i__ = 1; i__ <= i__1; ++i__) {
1078 x[i__] = x[i__] / c1[i__] - c2[i__];
1083 /* +-----------------------------------------------------------------------+ */
1084 /* | SUBROUTINE DIRGet_I | */
1085 /* +-----------------------------------------------------------------------+ */
1086 /* Subroutine */ void direct_dirget_i__(integer *length, integer *pos, integer *
1087 arrayi, integer *maxi, integer *n, integer *maxfunc)
1089 /* System generated locals */
1090 integer length_dim1, length_offset, i__1;
1092 /* Local variables */
1093 integer i__, j, help;
1095 /* Parameter adjustments */
1098 length_offset = 1 + length_dim1;
1099 length -= length_offset;
1103 help = length[*pos * length_dim1 + 1];
1105 for (i__ = 2; i__ <= i__1; ++i__) {
1106 if (length[i__ + *pos * length_dim1] < help) {
1107 help = length[i__ + *pos * length_dim1];
1112 for (i__ = 1; i__ <= i__1; ++i__) {
1113 if (length[i__ + *pos * length_dim1] == help) {
1122 /* +-----------------------------------------------------------------------+ */
1123 /* | SUBROUTINE DIRInit | */
1124 /* | Initialise all needed variables and do the first run of the | */
1125 /* | algorithm. | */
1126 /* | Changed 02/24/2000 | */
1127 /* | Changed fcn Double precision to fcn external! | */
1128 /* | Changed 09/15/2000 | */
1129 /* | Added distinction between Jones way to characterize rectangles | */
1130 /* | and our way. Common variable JONES controls which way we use. | */
1131 /* | JONES = 0 Jones way (Distance from midpoint to corner) | */
1132 /* | JONES = 1 Our way (Length of longest side) | */
1133 /* | Changed 09/24/00 | */
1134 /* | Added array levels. Levels contain the values to characterize | */
1135 /* | the hyperrectangles. | */
1136 /* | Changed 01/22/01 | */
1137 /* | Added variable fmax to keep track of maximum value found. | */
1138 /* | Added variable Ifeasiblef to keep track if feasibel point has | */
1139 /* | been found. | */
1140 /* | Changed 01/23/01 | */
1141 /* | Added variable Ierror to keep track of errors. | */
1142 /* +-----------------------------------------------------------------------+ */
1143 /* Subroutine */ void direct_dirinit_(doublereal *f, fp fcn, doublereal *c__,
1144 integer *length, integer *actdeep, integer *point, integer *anchor,
1145 integer *free, FILE *logfile, integer *arrayi,
1146 integer *maxi, integer *list2, doublereal *w, doublereal *x,
1147 doublereal *l, doublereal *u, doublereal *minf, integer *minpos,
1148 doublereal *thirds, doublereal *levels, integer *maxfunc, const integer *
1149 maxdeep, integer *n, integer *maxor, doublereal *fmax, integer *
1150 ifeasiblef, integer *iinfeasible, integer *ierror, void *fcndata,
1151 integer jones, int *force_stop)
1153 /* System generated locals */
1154 integer c_dim1, c_offset, length_dim1, length_offset, list2_dim1,
1155 list2_offset, i__1, i__2;
1157 /* Local variables */
1159 integer new__, help, oops;
1160 doublereal help2, delta;
1163 /* +-----------------------------------------------------------------------+ */
1164 /* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
1165 /* +-----------------------------------------------------------------------+ */
1166 /* +-----------------------------------------------------------------------+ */
1167 /* | JG 01/22/01 Added variable Ifeasiblef to keep track if feasibel point | */
1168 /* | has been found. | */
1169 /* | JG 01/23/01 Added variable Ierror to keep track of errors. | */
1170 /* | JG 03/09/01 Added IInfeasible to keep track if an infeasible point has| */
1171 /* | been found. | */
1172 /* +-----------------------------------------------------------------------+ */
1173 /* JG 09/15/00 Added variable JONES (see above) */
1174 /* +-----------------------------------------------------------------------+ */
1175 /* | Variables to pass user defined data to the function to be optimized. | */
1176 /* +-----------------------------------------------------------------------+ */
1177 /* Parameter adjustments */
1185 list2_dim1 = *maxor;
1186 list2_offset = 1 + list2_dim1;
1187 list2 -= list2_offset;
1190 length_offset = 1 + length_dim1;
1191 length -= length_offset;
1193 c_offset = 1 + c_dim1;
1199 /* JG 09/15/00 If Jones way of characterising rectangles is used, */
1200 /* initialise thirds to reflect this. */
1203 for (j = 0; j <= i__1; ++j) {
1204 w[j + 1] = sqrt(*n - j + j / 9.) * .5;
1208 i__1 = *maxdeep / *n;
1209 for (i__ = 1; i__ <= i__1; ++i__) {
1211 for (j = 0; j <= i__2; ++j) {
1212 levels[(i__ - 1) * *n + j] = w[j + 1] / help2;
1219 /* JG 09/15/00 Initialiase levels to contain 1/j */
1222 for (i__ = 1; i__ <= i__1; ++i__) {
1223 levels[i__] = 1. / help2;
1231 for (i__ = 1; i__ <= i__1; ++i__) {
1232 thirds[i__] = 1. / help2;
1238 for (i__ = 1; i__ <= i__1; ++i__) {
1239 c__[i__ + c_dim1] = .5;
1241 length[i__ + length_dim1] = 0;
1244 direct_dirinfcn_(fcn, &x[1], &l[1], &u[1], n, &f[3], &help, fcndata);
1245 if (force_stop && *force_stop) {
1249 f[4] = (doublereal) help;
1250 *iinfeasible = help;
1252 /* 09/25/00 Added this */
1253 /* if (f(1,1) .ge. 1.E+6) then */
1261 /* JG 09/25/00 Remove IF */
1269 direct_dirget_i__(&length[length_offset], &c__1, &arrayi[1], maxi, n, maxfunc);
1271 direct_dirsamplepoints_(&c__[c_offset], &arrayi[1], &delta, &c__1, &new__, &
1272 length[length_offset], logfile, &f[3], free, maxi, &
1273 point[1], &x[1], &l[1], minf, minpos, &u[1], n,
1274 maxfunc, maxdeep, &oops);
1275 /* +-----------------------------------------------------------------------+ */
1276 /* | JG 01/23/01 Added error checking. | */
1277 /* +-----------------------------------------------------------------------+ */
1282 /* +-----------------------------------------------------------------------+ */
1283 /* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
1284 /* | Added variable to keep track if feasible point was found. | */
1285 /* +-----------------------------------------------------------------------+ */
1286 direct_dirsamplef_(&c__[c_offset], &arrayi[1], &delta, &c__1, &new__, &length[
1287 length_offset], logfile, &f[3], free, maxi, &point[
1288 1], fcn, &x[1], &l[1], minf, minpos, &u[1], n, maxfunc,
1289 maxdeep, &oops, fmax, ifeasiblef, iinfeasible, fcndata,
1291 if (force_stop && *force_stop) {
1295 /* +-----------------------------------------------------------------------+ */
1296 /* | JG 01/23/01 Added error checking. | */
1297 /* +-----------------------------------------------------------------------+ */
1302 direct_dirdivide_(&new__, &c__0, &length[length_offset], &point[1], &arrayi[1], &
1303 c__1, &list2[list2_offset], &w[1], maxi, &f[3], maxfunc,
1305 direct_dirinsertlist_(&new__, &anchor[-1], &point[1], &f[3], maxi, &
1306 length[length_offset], maxfunc, maxdeep, n, &c__1, jones);
1309 /* +-----------------------------------------------------------------------+ */
1310 /* | SUBROUTINE DIRInitList | */
1311 /* | Initialise the list. | */
1312 /* +-----------------------------------------------------------------------+ */
1313 /* Subroutine */ void direct_dirinitlist_(integer *anchor, integer *free, integer *
1314 point, doublereal *f, integer *maxfunc, const integer *maxdeep)
1316 /* System generated locals */
1319 /* Local variables */
1322 /* f -- values of functions. */
1323 /* anchor -- anchors of lists with deep i */
1324 /* point -- lists */
1325 /* free -- first free position */
1326 /* Parameter adjustments */
1333 for (i__ = -1; i__ <= i__1; ++i__) {
1338 for (i__ = 1; i__ <= i__1; ++i__) {
1339 f[(i__ << 1) + 1] = 0.;
1340 f[(i__ << 1) + 2] = 0.;
1341 point[i__] = i__ + 1;
1345 point[*maxfunc] = 0;
1347 } /* dirinitlist_ */
1349 /* +-----------------------------------------------------------------------+ */
1351 /* | SUBROUTINE DIRPREPRC | */
1353 /* | Subroutine DIRpreprc uses an afine mapping to map the hyper-box given | */
1354 /* | by the constraints on the variable x onto the n-dimensional unit cube.| */
1355 /* | This mapping is done using the following equation: | */
1357 /* | x(i)=x(i)/(u(i)-l(i))-l(i)/(u(i)-l(i)). | */
1359 /* | DIRpreprc checks if the bounds l and u are well-defined. That is, if | */
1361 /* | l(i) < u(i) forevery i. | */
1365 /* | u -- A double-precision vector of length n. The vector | */
1366 /* | containing the upper bounds for the n independent | */
1367 /* | variables. | */
1369 /* | l -- A double-precision vector of length n. The vector | */
1370 /* | containing the lower bounds for the n independent | */
1371 /* | variables. | */
1373 /* | n -- An integer. The dimension of the problem. | */
1377 /* | xs1 -- A double-precision vector of length n, used for scaling | */
1378 /* | and unscaling the vector x. | */
1380 /* | xs2 -- A double-precision vector of length n, used for scaling | */
1381 /* | and unscaling the vector x. | */
1384 /* | oops -- An integer. If an upper bound is less than a lower | */
1385 /* | bound or if the initial point is not in the | */
1386 /* | hyper-box oops is set to 1 and iffco terminates. | */
1388 /* +-----------------------------------------------------------------------+ */
1389 /* Subroutine */ void direct_dirpreprc_(doublereal *u, doublereal *l, integer *n,
1390 doublereal *xs1, doublereal *xs2, integer *oops)
1392 /* System generated locals */
1395 /* Local variables */
1399 /* Parameter adjustments */
1408 for (i__ = 1; i__ <= i__1; ++i__) {
1409 /* +-----------------------------------------------------------------------+ */
1410 /* | Check if the hyper-box is well-defined. | */
1411 /* +-----------------------------------------------------------------------+ */
1412 if (u[i__] <= l[i__]) {
1418 /* +-----------------------------------------------------------------------+ */
1419 /* | Scale the initial iterate so that it is in the unit cube. | */
1420 /* +-----------------------------------------------------------------------+ */
1422 for (i__ = 1; i__ <= i__1; ++i__) {
1423 help = u[i__] - l[i__];
1424 xs2[i__] = l[i__] / help;
1430 /* Subroutine */ void direct_dirheader_(FILE *logfile, integer *version,
1431 doublereal *x, integer *n, doublereal *eps, integer *maxf, integer *
1432 maxt, doublereal *l, doublereal *u, integer *algmethod, integer *
1433 maxfunc, const integer *maxdeep, doublereal *fglobal, doublereal *fglper,
1434 integer *ierror, doublereal *epsfix, integer *iepschange, doublereal *
1435 volper, doublereal *sigmaper)
1437 /* System generated locals */
1440 /* Local variables */
1441 integer imainver, i__, numerrors, isubsubver, ihelp, isubver;
1444 /* +-----------------------------------------------------------------------+ */
1445 /* | Variables to pass user defined data to the function to be optimized. | */
1446 /* +-----------------------------------------------------------------------+ */
1447 /* Parameter adjustments */
1454 fprintf(logfile, "------------------- Log file ------------------\n");
1458 imainver = *version / 100;
1459 ihelp = *version - imainver * 100;
1460 isubver = ihelp / 10;
1461 ihelp -= isubver * 10;
1463 /* +-----------------------------------------------------------------------+ */
1464 /* | JG 01/13/01 Added check for epsilon. If epsilon is smaller 0, we use | */
1465 /* | the update formula from Jones. We then set the flag | */
1466 /* | iepschange to 1, and store the absolute value of eps in | */
1467 /* | epsfix. epsilon is then changed after each iteration. | */
1468 /* +-----------------------------------------------------------------------+ */
1478 /* +-----------------------------------------------------------------------+ */
1479 /* | JG 07/16/01 Removed printout of contents in cdata(1). | */
1480 /* +-----------------------------------------------------------------------+ */
1481 /* write(logfile,*) cdata(1) */
1484 fprintf(logfile, "DIRECT Version %d.%d.%d\n"
1485 " Problem dimension n: %d\n"
1487 " Maximum number of f-evaluations (maxf): %d\n"
1488 " Maximum number of iterations (MaxT): %d\n"
1489 " Value of f_global: %e\n"
1490 " Global percentage wanted: %e\n"
1491 " Volume percentage wanted: %e\n"
1492 " Measure percentage wanted: %e\n",
1493 imainver, isubver, isubsubver, *n, *eps, *maxf, *maxt,
1494 *fglobal, *fglper, *volper, *sigmaper);
1495 fprintf(logfile, *iepschange == 1
1496 ? "Epsilon is changed using the Jones formula.\n"
1497 : "Epsilon is constant.\n");
1498 fprintf(logfile, *algmethod == 0
1499 ? "Jones original DIRECT algorithm is used.\n"
1500 : "Our modification of the DIRECT algorithm is used.\n");
1504 for (i__ = 1; i__ <= i__1; ++i__) {
1505 if (u[i__] <= l[i__]) {
1508 fprintf(logfile, "WARNING: bounds on variable x%d: "
1509 "%g <= xi <= %g\n", i__, l[i__], u[i__]);
1513 fprintf(logfile, "Bounds on variable x%d: "
1514 "%g <= xi <= %g\n", i__, l[i__], u[i__]);
1518 /* +-----------------------------------------------------------------------+ */
1519 /* | If there are to many function evaluations or to many iteration, note | */
1520 /* | this and set the error flag accordingly. Note: If more than one error | */
1521 /* | occurred, we give out an extra message. | */
1522 /* +-----------------------------------------------------------------------+ */
1523 if (*maxf + 20 > *maxfunc) {
1526 "WARNING: The maximum number of function evaluations (%d) is higher than\n"
1527 " the constant maxfunc (%d). Increase maxfunc in subroutine DIRECT\n"
1528 " or decrease the maximum number of function evaluations.\n",
1534 if (logfile) fprintf(logfile, "----------------------------------\n");
1535 if (numerrors == 1) {
1537 fprintf(logfile, "WARNING: One error in the input!\n");
1540 fprintf(logfile, "WARNING: %d errors in the input!\n",
1544 if (logfile) fprintf(logfile, "----------------------------------\n");
1547 fprintf(logfile, "Iteration # of f-eval. minf\n");
1552 /* Subroutine */ void direct_dirsummary_(FILE *logfile, doublereal *x, doublereal *
1553 l, doublereal *u, integer *n, doublereal *minf, doublereal *fglobal,
1554 integer *numfunc, integer *ierror)
1556 /* Local variables */
1559 /* Parameter adjustments */
1566 fprintf(logfile, "-----------------------Summary------------------\n"
1567 "Final function value: %g\n"
1568 "Number of function evaluations: %d\n", *minf, *numfunc);
1569 if (*fglobal > -1e99)
1570 fprintf(logfile, "Final function value is within %g%% of global optimum\n", 100*(*minf - *fglobal) / MAX(1.0, fabs(*fglobal)));
1571 fprintf(logfile, "Index, final solution, x(i)-l(i), u(i)-x(i)\n");
1572 for (i__ = 1; i__ <= *n; ++i__)
1573 fprintf(logfile, "%d, %g, %g, %g\n", i__, x[i__],
1574 x[i__]-l[i__], u[i__] - x[i__]);
1575 fprintf(logfile, "-----------------------------------------------\n");