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;
129 helplower = HUGE_VAL;
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) {
172 helplower = HUGE_VAL;
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 (helpgreater <= helplower) {
232 if (*cheat == 1 && helplower > *kmax) {
235 if (f[(j___ << 1) + 1] - helplower * thirds[s[j + (s_dim1 << 1)]] >
236 MIN(*minf - epsrel * fabs(*minf),
239 fprintf(logfile, "> minf - epslminfl\n");
244 fprintf(logfile, "helpgreater > helplower: %g %g %g\n",
245 helpgreater, helplower, helpgreater - helplower);
256 s[*maxpos + s_dim1] = novalue;
257 s[*maxpos + (s_dim1 << 1)] = novaluedeep;
261 /* +-----------------------------------------------------------------------+ */
262 /* | SUBROUTINE DIRDoubleInsert | */
263 /* | Routine to make sure that if there are several potential optimal | */
264 /* | hyperrectangles of the same level (i.e. hyperrectangles that have| */
265 /* | the same level and the same function value at the center), all of| */
266 /* | them are divided. This is the way as originally described in | */
267 /* | Jones et.al. | */
268 /* | JG 07/16/01 Added errorflag to calling sequence. We check if more | */
269 /* | we reach the capacity of the array S. If this happens, we | */
270 /* | return to the main program with an error. | */
271 /* +-----------------------------------------------------------------------+ */
272 /* Subroutine */ void direct_dirdoubleinsert_(integer *anchor, integer *s, integer *
273 maxpos, integer *point, doublereal *f, const integer *maxdeep, integer *
274 maxfunc, const integer *maxdiv, integer *ierror)
276 /* System generated locals */
277 integer s_dim1, s_offset, i__1;
279 /* Local variables */
280 integer i__, oldmaxpos, pos, help, iflag, actdeep;
282 /* +-----------------------------------------------------------------------+ */
283 /* | JG 07/16/01 Added flag to prevent run time-errors on some systems. | */
284 /* +-----------------------------------------------------------------------+ */
285 /* Parameter adjustments */
290 s_offset = 1 + s_dim1;
296 for (i__ = 1; i__ <= i__1; ++i__) {
297 if (s[i__ + s_dim1] > 0) {
298 actdeep = s[i__ + (s_dim1 << 1)];
299 help = anchor[actdeep];
302 /* +-----------------------------------------------------------------------+ */
303 /* | JG 07/16/01 Added flag to prevent run time-errors on some systems. On | */
304 /* | some systems the second conditions in an AND statement is | */
305 /* | evaluated even if the first one is already not true. | */
306 /* +-----------------------------------------------------------------------+ */
307 while(pos > 0 && iflag == 0) {
308 if (f[(pos << 1) + 1] - f[(help << 1) + 1] <= 1e-13) {
309 if (*maxpos < *maxdiv) {
311 s[*maxpos + s_dim1] = pos;
312 s[*maxpos + (s_dim1 << 1)] = actdeep;
315 /* +-----------------------------------------------------------------------+ */
316 /* | JG 07/16/01 Maximum number of elements possible in S has been reached!| */
317 /* +-----------------------------------------------------------------------+ */
328 } /* dirdoubleinsert_ */
330 /* +-----------------------------------------------------------------------+ */
331 /* | INTEGER Function GetmaxDeep | */
332 /* | function to get the maximal length (1/length) of the n-dimensional | */
333 /* | rectangle with midpoint pos. | */
335 /* | On Return : | */
336 /* | the maximal length | */
338 /* | pos -- the position of the midpoint in the array length | */
339 /* | length -- the array with the dimensions | */
340 /* | maxfunc -- the leading dimension of length | */
341 /* | n -- the dimension of the problem | */
343 /* +-----------------------------------------------------------------------+ */
344 integer direct_dirgetmaxdeep_(integer *pos, integer *length, integer *maxfunc,
347 /* System generated locals */
348 integer length_dim1, length_offset, i__1, i__2, i__3;
350 /* Local variables */
353 /* Parameter adjustments */
355 length_offset = 1 + length_dim1;
356 length -= length_offset;
359 help = length[*pos * length_dim1 + 1];
361 for (i__ = 2; i__ <= i__1; ++i__) {
363 i__2 = help, i__3 = length[i__ + *pos * length_dim1];
364 help = MIN(i__2,i__3);
368 } /* dirgetmaxdeep_ */
370 static integer isinbox_(doublereal *x, doublereal *a, doublereal *b, integer *n,
373 /* System generated locals */
374 integer ret_val, i__1;
376 /* Local variables */
377 integer outofbox, i__;
379 /* Parameter adjustments */
387 for (i__ = 1; i__ <= i__1; ++i__) {
388 if (a[i__] > x[i__] || b[i__] < x[i__]) {
399 /* +-----------------------------------------------------------------------+ */
400 /* | JG Added 09/25/00 | */
402 /* | SUBROUTINE DIRResortlist | */
404 /* | Resort the list so that the infeasible point is in the list with the | */
405 /* | replaced value. | */
406 /* +-----------------------------------------------------------------------+ */
407 /* Subroutine */ static void dirresortlist_(integer *replace, integer *anchor,
408 doublereal *f, integer *point, integer *length, integer *n, integer *
409 maxfunc, integer *maxdim, const integer *maxdeep, FILE *logfile,
412 /* System generated locals */
413 integer length_dim1, length_offset, i__1;
415 /* Local variables */
419 /* +-----------------------------------------------------------------------+ */
420 /* | Get the length of the hyper rectangle with infeasible mid point and | */
421 /* | Index of the corresponding list. | */
422 /* +-----------------------------------------------------------------------+ */
423 /* JG 09/25/00 Replaced with DIRgetlevel */
424 /* l = DIRgetmaxDeep(replace,length,maxfunc,n) */
425 /* Parameter adjustments */
429 length_offset = 1 + length_dim1;
430 length -= length_offset;
434 l = direct_dirgetlevel_(replace, &length[length_offset], maxfunc, n, jones);
436 /* +-----------------------------------------------------------------------+ */
437 /* | If the hyper rectangle with infeasibel midpoint is already the start | */
438 /* | of the list, give out message, nothing to do. | */
439 /* +-----------------------------------------------------------------------+ */
440 if (*replace == start) {
441 /* write(logfile,*) 'No resorting of list necessarry, since new ', */
442 /* + 'point is already anchor of list .',l */
444 /* +-----------------------------------------------------------------------+ */
445 /* | Take the hyper rectangle with infeasible midpoint out of the list. | */
446 /* +-----------------------------------------------------------------------+ */
449 for (i__ = 1; i__ <= i__1; ++i__) {
450 if (point[pos] == *replace) {
451 point[pos] = point[*replace];
458 fprintf(logfile, "Error in DIRREsortlist: "
459 "We went through the whole list\n"
460 "and could not find the point to replace!!\n");
465 /* +-----------------------------------------------------------------------+ */
466 /* | If the anchor of the list has a higher value than the value of a | */
467 /* | nearby point, put the infeasible point at the beginning of the list. | */
468 /* +-----------------------------------------------------------------------+ */
470 if (f[(start << 1) + 1] > f[(*replace << 1) + 1]) {
471 anchor[l] = *replace;
472 point[*replace] = start;
473 /* write(logfile,*) 'Point is replacing current anchor for ' */
474 /* + , 'this list ',l,replace,start */
476 /* +-----------------------------------------------------------------------+ */
477 /* | Insert the point into the list according to its (replaced) function | */
479 /* +-----------------------------------------------------------------------+ */
482 for (i__ = 1; i__ <= i__1; ++i__) {
483 /* +-----------------------------------------------------------------------+ */
484 /* | The point has to be added at the end of the list. | */
485 /* +-----------------------------------------------------------------------+ */
486 if (point[pos] == 0) {
487 point[*replace] = point[pos];
488 point[pos] = *replace;
489 /* write(logfile,*) 'Point is added at the end of the ' */
490 /* + , 'list ',l, replace */
493 if (f[(point[pos] << 1) + 1] > f[(*replace << 1) + 1]) {
494 point[*replace] = point[pos];
495 point[pos] = *replace;
496 /* write(logfile,*) 'There are points with a higher ' */
497 /* + ,'f-value in the list ',l,replace, pos */
508 } /* dirresortlist_ */
510 /* +-----------------------------------------------------------------------+ */
511 /* | JG Added 09/25/00 | */
512 /* | SUBROUTINE DIRreplaceInf | */
514 /* | Find out if there are infeasible points which are near feasible ones. | */
515 /* | If this is the case, replace the function value at the center of the | */
516 /* | hyper rectangle by the lowest function value of a nearby function. | */
517 /* +-----------------------------------------------------------------------+ */
518 /* Subroutine */ void direct_dirreplaceinf_(integer *free, integer *freeold,
519 doublereal *f, doublereal *c__, doublereal *thirds, integer *length,
520 integer *anchor, integer *point, doublereal *c1, doublereal *c2,
521 integer *maxfunc, const integer *maxdeep, integer *maxdim, integer *n,
522 FILE *logfile, doublereal *fmax, integer jones)
524 /* System generated locals */
525 integer c_dim1, c_offset, length_dim1, length_offset, i__1, i__2, i__3;
526 doublereal d__1, d__2;
528 /* Local variables */
529 doublereal a[32], b[32];
530 integer i__, j, k, l;
531 doublereal x[32], sidelength;
534 /* +-----------------------------------------------------------------------+ */
535 /* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
536 /* +-----------------------------------------------------------------------+ */
537 /* Parameter adjustments */
541 length_dim1 = *maxdim;
542 length_offset = 1 + length_dim1;
543 length -= length_offset;
545 c_offset = 1 + c_dim1;
552 for (i__ = 1; i__ <= i__1; ++i__) {
553 if (f[(i__ << 1) + 2] > 0.) {
554 /* +-----------------------------------------------------------------------+ */
555 /* | Get the maximum side length of the hyper rectangle and then set the | */
556 /* | new side length to this lengths times the growth factor. | */
557 /* +-----------------------------------------------------------------------+ */
558 help = direct_dirgetmaxdeep_(&i__, &length[length_offset], maxfunc, n);
559 sidelength = thirds[help] * 2.;
560 /* +-----------------------------------------------------------------------+ */
561 /* | Set the Center and the upper and lower bounds of the rectangles. | */
562 /* +-----------------------------------------------------------------------+ */
564 for (j = 1; j <= i__2; ++j) {
565 sidelength = thirds[length[i__ + j * length_dim1]];
566 a[j - 1] = c__[j + i__ * c_dim1] - sidelength;
567 b[j - 1] = c__[j + i__ * c_dim1] + sidelength;
570 /* +-----------------------------------------------------------------------+ */
571 /* | The function value is reset to 'Inf', since it may have been changed | */
572 /* | in an earlier iteration and now the feasible point which was close | */
573 /* | is not close anymore (since the hyper rectangle surrounding the | */
574 /* | current point may have shrunk). | */
575 /* +-----------------------------------------------------------------------+ */
576 f[(i__ << 1) + 1] = HUGE_VAL;
577 f[(i__ << 1) + 2] = 2.;
578 /* +-----------------------------------------------------------------------+ */
579 /* | Check if any feasible point is near this infeasible point. | */
580 /* +-----------------------------------------------------------------------+ */
582 for (k = 1; k <= i__2; ++k) {
583 /* +-----------------------------------------------------------------------+ */
584 /* | If the point k is feasible, check if it is near. | */
585 /* +-----------------------------------------------------------------------+ */
586 if (f[(k << 1) + 2] == 0.) {
587 /* +-----------------------------------------------------------------------+ */
588 /* | Copy the coordinates of the point k into x. | */
589 /* +-----------------------------------------------------------------------+ */
591 for (l = 1; l <= i__3; ++l) {
592 x[l - 1] = c__[l + k * c_dim1];
595 /* +-----------------------------------------------------------------------+ */
596 /* | Check if the point k is near the infeasible point, if so, replace the | */
598 /* +-----------------------------------------------------------------------+ */
599 if (isinbox_(x, a, b, n, &c__32) == 1) {
601 d__1 = f[(i__ << 1) + 1], d__2 = f[(k << 1) + 1];
602 f[(i__ << 1) + 1] = MIN(d__1,d__2);
603 f[(i__ << 1) + 2] = 1.;
608 if (f[(i__ << 1) + 2] == 1.) {
609 f[(i__ << 1) + 1] += (d__1 = f[(i__ << 1) + 1], fabs(d__1)) *
612 for (l = 1; l <= i__2; ++l) {
613 x[l - 1] = c__[l + i__ * c_dim1] * c1[l] + c__[l + i__ *
617 dirresortlist_(&i__, &anchor[-1], &f[3], &point[1],
618 &length[length_offset], n,
619 maxfunc, maxdim, maxdeep, logfile, jones);
621 /* +-----------------------------------------------------------------------+ */
622 /* | JG 01/22/01 | */
623 /* | Replaced fixed value for infeasible points with maximum value found, | */
624 /* | increased by 1. | */
625 /* +-----------------------------------------------------------------------+ */
626 if (! (*fmax == f[(i__ << 1) + 1])) {
628 d__1 = *fmax + 1., d__2 = f[(i__ << 1) + 1];
629 f[(i__ << 1) + 1] = MAX(d__1,d__2);
636 } /* dirreplaceinf_ */
638 /* +-----------------------------------------------------------------------+ */
639 /* | SUBROUTINE DIRInsert | */
640 /* +-----------------------------------------------------------------------+ */
641 /* Subroutine */ static void dirinsert_(integer *start, integer *ins, integer *point,
642 doublereal *f, integer *maxfunc)
644 /* System generated locals */
647 /* Local variables */
650 /* JG 09/17/00 Rewrote this routine. */
651 /* DO 10,i = 1,maxfunc */
652 /* IF (f(ins,1) .LT. f(point(start),1)) THEN */
653 /* help = point(start) */
654 /* point(start) = ins */
655 /* point(ins) = help */
658 /* IF (point(start) .EQ. 0) THEN */
659 /* point(start) = ins */
663 /* start = point(start) */
666 /* Parameter adjustments */
672 for (i__ = 1; i__ <= i__1; ++i__) {
673 if (point[*start] == 0) {
674 point[*start] = *ins;
677 } else if (f[(*ins << 1) + 1] < f[(point[*start] << 1) + 1]) {
678 help = point[*start];
679 point[*start] = *ins;
683 *start = point[*start];
688 /* +-----------------------------------------------------------------------+ */
689 /* | SUBROUTINE DIRInsertList | */
690 /* | Changed 02-24-2000 | */
691 /* | Got rid of the distinction between feasible and infeasible points| */
692 /* | I could do this since infeasible points get set to a high | */
693 /* | function value, which may be replaced by a function value of a | */
694 /* | nearby function at the end of the main loop. | */
695 /* +-----------------------------------------------------------------------+ */
696 /* Subroutine */ void direct_dirinsertlist_(integer *new__, integer *anchor, integer *
697 point, doublereal *f, integer *maxi, integer *length, integer *
698 maxfunc, const integer *maxdeep, integer *n, integer *samp,
701 /* System generated locals */
702 integer length_dim1, length_offset, i__1;
704 /* Local variables */
707 integer pos1, pos2, deep;
709 /* JG 09/24/00 Changed this to Getlevel */
710 /* Parameter adjustments */
715 length_offset = 1 + length_dim1;
716 length -= length_offset;
720 for (j = 1; j <= i__1; ++j) {
723 *new__ = point[pos2];
724 /* JG 09/24/00 Changed this to Getlevel */
725 /* deep = DIRGetMaxdeep(pos1,length,maxfunc,n) */
726 deep = direct_dirgetlevel_(&pos1, &length[length_offset], maxfunc, n, jones);
727 if (anchor[deep] == 0) {
728 if (f[(pos2 << 1) + 1] < f[(pos1 << 1) + 1]) {
738 if (f[(pos2 << 1) + 1] < f[(pos1 << 1) + 1]) {
739 if (f[(pos2 << 1) + 1] < f[(pos << 1) + 1]) {
741 /* JG 08/30/00 Fixed bug. Sorting was not correct when */
742 /* f(1,pos2) < f(1,pos1) < f(1,pos) */
743 if (f[(pos1 << 1) + 1] < f[(pos << 1) + 1]) {
748 dirinsert_(&pos, &pos1, &point[1], &f[3], maxfunc);
751 dirinsert_(&pos, &pos2, &point[1], &f[3], maxfunc);
752 dirinsert_(&pos, &pos1, &point[1], &f[3], maxfunc);
755 if (f[(pos1 << 1) + 1] < f[(pos << 1) + 1]) {
756 /* JG 08/30/00 Fixed bug. Sorting was not correct when */
757 /* f(pos1,1) < f(pos2,1) < f(pos,1) */
759 if (f[(pos << 1) + 1] < f[(pos2 << 1) + 1]) {
761 dirinsert_(&pos, &pos2, &point[1], &f[3], maxfunc);
767 dirinsert_(&pos, &pos1, &point[1], &f[3], maxfunc);
768 dirinsert_(&pos, &pos2, &point[1], &f[3], maxfunc);
774 /* JG 09/24/00 Changed this to Getlevel */
775 /* deep = DIRGetMaxdeep(samp,length,maxfunc,n) */
776 deep = direct_dirgetlevel_(samp, &length[length_offset], maxfunc, n, jones);
778 if (f[(*samp << 1) + 1] < f[(pos << 1) + 1]) {
779 anchor[deep] = *samp;
782 dirinsert_(&pos, samp, &point[1], &f[3], maxfunc);
784 } /* dirinsertlist_ */
786 /* +-----------------------------------------------------------------------+ */
787 /* | SUBROUTINE DIRInsertList2 (Old way to do it.) | */
788 /* +-----------------------------------------------------------------------+ */
789 /* Subroutine */ static void dirinsertlist_2__(integer *start, integer *j, integer *k,
790 integer *list2, doublereal *w, integer *maxi, integer *n)
792 /* System generated locals */
793 integer list2_dim1, list2_offset, i__1;
795 /* Local variables */
798 /* Parameter adjustments */
801 list2_offset = 1 + list2_dim1;
802 list2 -= list2_offset;
807 list2[*j + list2_dim1] = 0;
811 if (w[*start] > w[*j]) {
812 list2[*j + list2_dim1] = *start;
816 for (i__ = 1; i__ <= i__1; ++i__) {
817 if (list2[pos + list2_dim1] == 0) {
818 list2[*j + list2_dim1] = 0;
819 list2[pos + list2_dim1] = *j;
822 if (w[*j] < w[list2[pos + list2_dim1]]) {
823 list2[*j + list2_dim1] = list2[pos + list2_dim1];
824 list2[pos + list2_dim1] = *j;
828 pos = list2[pos + list2_dim1];
833 list2[*j + (list2_dim1 << 1)] = *k;
834 } /* dirinsertlist_2__ */
836 /* +-----------------------------------------------------------------------+ */
837 /* | SUBROUTINE DIRSearchmin | */
838 /* | Search for the minimum in the list. ! */
839 /* +-----------------------------------------------------------------------+ */
840 /* Subroutine */ static void dirsearchmin_(integer *start, integer *list2, integer *
841 pos, integer *k, integer *n)
843 /* System generated locals */
844 integer list2_dim1, list2_offset;
846 /* Parameter adjustments */
848 list2_offset = 1 + list2_dim1;
849 list2 -= list2_offset;
853 *pos = list2[*start + (list2_dim1 << 1)];
854 *start = list2[*start + list2_dim1];
855 } /* dirsearchmin_ */
857 /* +-----------------------------------------------------------------------+ */
858 /* | SUBROUTINE DIRSamplepoints | */
859 /* | Subroutine to sample the new points. | */
860 /* +-----------------------------------------------------------------------+ */
861 /* Subroutine */ void direct_dirsamplepoints_(doublereal *c__, integer *arrayi,
862 doublereal *delta, integer *sample, integer *start, integer *length,
863 FILE *logfile, doublereal *f, integer *free,
864 integer *maxi, integer *point, doublereal *x, doublereal *l,
865 doublereal *minf, integer *minpos, doublereal *u, integer *n,
866 integer *maxfunc, const integer *maxdeep, integer *oops)
868 /* System generated locals */
869 integer length_dim1, length_offset, c_dim1, c_offset, i__1, i__2;
871 /* Local variables */
875 /* Parameter adjustments */
883 length_offset = 1 + length_dim1;
884 length -= length_offset;
886 c_offset = 1 + c_dim1;
893 i__1 = *maxi + *maxi;
894 for (k = 1; k <= i__1; ++k) {
896 for (j = 1; j <= i__2; ++j) {
897 length[j + *free * length_dim1] = length[j + *sample *
899 c__[j + *free * c_dim1] = c__[j + *sample * c_dim1];
903 *free = point[*free];
906 fprintf(logfile, "Error, no more free positions! "
907 "Increase maxfunc!\n");
916 for (j = 1; j <= i__1; ++j) {
917 c__[arrayi[j] + pos * c_dim1] = c__[arrayi[j] + *sample * c_dim1] + *
920 c__[arrayi[j] + pos * c_dim1] = c__[arrayi[j] + *sample * c_dim1] - *
926 } /* dirsamplepoints_ */
928 /* +-----------------------------------------------------------------------+ */
929 /* | SUBROUTINE DIRDivide | */
930 /* | Subroutine to divide the hyper rectangles according to the rules. | */
931 /* | Changed 02-24-2000 | */
932 /* | Replaced if statement by min (line 367) | */
933 /* +-----------------------------------------------------------------------+ */
934 /* Subroutine */ void direct_dirdivide_(integer *new__, integer *currentlength,
935 integer *length, integer *point, integer *arrayi, integer *sample,
936 integer *list2, doublereal *w, integer *maxi, doublereal *f, integer *
937 maxfunc, const integer *maxdeep, integer *n)
939 /* System generated locals */
940 integer length_dim1, length_offset, list2_dim1, list2_offset, i__1, i__2;
941 doublereal d__1, d__2;
943 /* Local variables */
944 integer i__, j, k, pos, pos2;
948 /* Parameter adjustments */
953 list2_offset = 1 + list2_dim1;
954 list2 -= list2_offset;
957 length_offset = 1 + length_dim1;
958 length -= length_offset;
964 for (i__ = 1; i__ <= i__1; ++i__) {
966 w[j] = f[(pos << 1) + 1];
970 d__1 = f[(pos << 1) + 1], d__2 = w[j];
971 w[j] = MIN(d__1,d__2);
973 dirinsertlist_2__(&start, &j, &k, &list2[list2_offset], &w[1], maxi,
979 for (j = 1; j <= i__1; ++j) {
980 dirsearchmin_(&start, &list2[list2_offset], &pos, &k, n);
982 length[k + *sample * length_dim1] = *currentlength + 1;
983 i__2 = *maxi - j + 1;
984 for (i__ = 1; i__ <= i__2; ++i__) {
985 length[k + pos * length_dim1] = *currentlength + 1;
987 length[k + pos * length_dim1] = *currentlength + 1;
988 /* JG 07/10/01 pos2 = 0 at the end of the 30-loop. Since we end */
989 /* the loop now, we do not need to reassign pos and pos2. */
991 pos = list2[pos2 + (list2_dim1 << 1)];
992 pos2 = list2[pos2 + list2_dim1];
1000 /* +-----------------------------------------------------------------------+ */
1002 /* | SUBROUTINE DIRINFCN | */
1004 /* | Subroutine DIRinfcn unscales the variable x for use in the | */
1005 /* | user-supplied function evaluation subroutine fcn. After fcn returns | */
1006 /* | to DIRinfcn, DIRinfcn then rescales x for use by DIRECT. | */
1010 /* | fcn -- The argument containing the name of the user-supplied | */
1011 /* | subroutine that returns values for the function to be | */
1012 /* | minimized. | */
1014 /* | x -- A double-precision vector of length n. The point at | */
1015 /* | which the derivative is to be evaluated. | */
1017 /* | xs1 -- A double-precision vector of length n. Used for | */
1018 /* | scaling and unscaling the vector x by DIRinfcn. | */
1020 /* | xs2 -- A double-precision vector of length n. Used for | */
1021 /* | scaling and unscaling the vector x by DIRinfcn. | */
1023 /* | n -- An integer. The dimension of the problem. | */
1024 /* | kret -- An Integer. If kret = 1, the point is infeasible, | */
1025 /* | kret = -1, bad problem set up, | */
1026 /* | kret = 0, feasible. | */
1030 /* | f -- A double-precision scalar. | */
1032 /* | Subroutines and Functions | */
1034 /* | The subroutine whose name is passed through the argument fcn. | */
1036 /* +-----------------------------------------------------------------------+ */
1037 /* Subroutine */ void direct_dirinfcn_(fp fcn, doublereal *x, doublereal *c1,
1038 doublereal *c2, integer *n, doublereal *f, integer *flag__,
1041 /* System generated locals */
1044 /* Local variables */
1047 /* +-----------------------------------------------------------------------+ */
1048 /* | Variables to pass user defined data to the function to be optimized. | */
1049 /* +-----------------------------------------------------------------------+ */
1050 /* +-----------------------------------------------------------------------+ */
1051 /* | Unscale the variable x. | */
1052 /* +-----------------------------------------------------------------------+ */
1053 /* Parameter adjustments */
1060 for (i__ = 1; i__ <= i__1; ++i__) {
1061 x[i__] = (x[i__] + c2[i__]) * c1[i__];
1064 /* +-----------------------------------------------------------------------+ */
1065 /* | Call the function-evaluation subroutine fcn. | */
1066 /* +-----------------------------------------------------------------------+ */
1068 *f = fcn(*n, &x[1], flag__, fcn_data);
1069 /* +-----------------------------------------------------------------------+ */
1070 /* | Rescale the variable x. | */
1071 /* +-----------------------------------------------------------------------+ */
1073 for (i__ = 1; i__ <= i__1; ++i__) {
1074 x[i__] = x[i__] / c1[i__] - c2[i__];
1079 /* +-----------------------------------------------------------------------+ */
1080 /* | SUBROUTINE DIRGet_I | */
1081 /* +-----------------------------------------------------------------------+ */
1082 /* Subroutine */ void direct_dirget_i__(integer *length, integer *pos, integer *
1083 arrayi, integer *maxi, integer *n, integer *maxfunc)
1085 /* System generated locals */
1086 integer length_dim1, length_offset, i__1;
1088 /* Local variables */
1089 integer i__, j, help;
1091 /* Parameter adjustments */
1094 length_offset = 1 + length_dim1;
1095 length -= length_offset;
1099 help = length[*pos * length_dim1 + 1];
1101 for (i__ = 2; i__ <= i__1; ++i__) {
1102 if (length[i__ + *pos * length_dim1] < help) {
1103 help = length[i__ + *pos * length_dim1];
1108 for (i__ = 1; i__ <= i__1; ++i__) {
1109 if (length[i__ + *pos * length_dim1] == help) {
1118 /* +-----------------------------------------------------------------------+ */
1119 /* | SUBROUTINE DIRInit | */
1120 /* | Initialise all needed variables and do the first run of the | */
1121 /* | algorithm. | */
1122 /* | Changed 02/24/2000 | */
1123 /* | Changed fcn Double precision to fcn external! | */
1124 /* | Changed 09/15/2000 | */
1125 /* | Added distinction between Jones way to characterize rectangles | */
1126 /* | and our way. Common variable JONES controls which way we use. | */
1127 /* | JONES = 0 Jones way (Distance from midpoint to corner) | */
1128 /* | JONES = 1 Our way (Length of longest side) | */
1129 /* | Changed 09/24/00 | */
1130 /* | Added array levels. Levels contain the values to characterize | */
1131 /* | the hyperrectangles. | */
1132 /* | Changed 01/22/01 | */
1133 /* | Added variable fmax to keep track of maximum value found. | */
1134 /* | Added variable Ifeasiblef to keep track if feasibel point has | */
1135 /* | been found. | */
1136 /* | Changed 01/23/01 | */
1137 /* | Added variable Ierror to keep track of errors. | */
1138 /* +-----------------------------------------------------------------------+ */
1139 /* Subroutine */ void direct_dirinit_(doublereal *f, fp fcn, doublereal *c__,
1140 integer *length, integer *actdeep, integer *point, integer *anchor,
1141 integer *free, FILE *logfile, integer *arrayi,
1142 integer *maxi, integer *list2, doublereal *w, doublereal *x,
1143 doublereal *l, doublereal *u, doublereal *minf, integer *minpos,
1144 doublereal *thirds, doublereal *levels, integer *maxfunc, const integer *
1145 maxdeep, integer *n, integer *maxor, doublereal *fmax, integer *
1146 ifeasiblef, integer *iinfeasible, integer *ierror, void *fcndata,
1147 integer jones, int *force_stop)
1149 /* System generated locals */
1150 integer c_dim1, c_offset, length_dim1, length_offset, list2_dim1,
1151 list2_offset, i__1, i__2;
1153 /* Local variables */
1155 integer new__, help, oops;
1156 doublereal help2, delta;
1159 /* +-----------------------------------------------------------------------+ */
1160 /* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
1161 /* +-----------------------------------------------------------------------+ */
1162 /* +-----------------------------------------------------------------------+ */
1163 /* | JG 01/22/01 Added variable Ifeasiblef to keep track if feasibel point | */
1164 /* | has been found. | */
1165 /* | JG 01/23/01 Added variable Ierror to keep track of errors. | */
1166 /* | JG 03/09/01 Added IInfeasible to keep track if an infeasible point has| */
1167 /* | been found. | */
1168 /* +-----------------------------------------------------------------------+ */
1169 /* JG 09/15/00 Added variable JONES (see above) */
1170 /* +-----------------------------------------------------------------------+ */
1171 /* | Variables to pass user defined data to the function to be optimized. | */
1172 /* +-----------------------------------------------------------------------+ */
1173 /* Parameter adjustments */
1181 list2_dim1 = *maxor;
1182 list2_offset = 1 + list2_dim1;
1183 list2 -= list2_offset;
1186 length_offset = 1 + length_dim1;
1187 length -= length_offset;
1189 c_offset = 1 + c_dim1;
1195 /* JG 09/15/00 If Jones way of characterising rectangles is used, */
1196 /* initialise thirds to reflect this. */
1199 for (j = 0; j <= i__1; ++j) {
1200 w[j + 1] = sqrt(*n - j + j / 9.) * .5;
1204 i__1 = *maxdeep / *n;
1205 for (i__ = 1; i__ <= i__1; ++i__) {
1207 for (j = 0; j <= i__2; ++j) {
1208 levels[(i__ - 1) * *n + j] = w[j + 1] / help2;
1215 /* JG 09/15/00 Initialiase levels to contain 1/j */
1218 for (i__ = 1; i__ <= i__1; ++i__) {
1219 levels[i__] = 1. / help2;
1227 for (i__ = 1; i__ <= i__1; ++i__) {
1228 thirds[i__] = 1. / help2;
1234 for (i__ = 1; i__ <= i__1; ++i__) {
1235 c__[i__ + c_dim1] = .5;
1237 length[i__ + length_dim1] = 0;
1240 direct_dirinfcn_(fcn, &x[1], &l[1], &u[1], n, &f[3], &help, fcndata);
1241 if (force_stop && *force_stop) {
1245 f[4] = (doublereal) help;
1246 *iinfeasible = help;
1248 /* 09/25/00 Added this */
1249 /* if (f(1,1) .ge. 1.E+6) then */
1257 /* JG 09/25/00 Remove IF */
1265 direct_dirget_i__(&length[length_offset], &c__1, &arrayi[1], maxi, n, maxfunc);
1267 direct_dirsamplepoints_(&c__[c_offset], &arrayi[1], &delta, &c__1, &new__, &
1268 length[length_offset], logfile, &f[3], free, maxi, &
1269 point[1], &x[1], &l[1], minf, minpos, &u[1], n,
1270 maxfunc, maxdeep, &oops);
1271 /* +-----------------------------------------------------------------------+ */
1272 /* | JG 01/23/01 Added error checking. | */
1273 /* +-----------------------------------------------------------------------+ */
1278 /* +-----------------------------------------------------------------------+ */
1279 /* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
1280 /* | Added variable to keep track if feasible point was found. | */
1281 /* +-----------------------------------------------------------------------+ */
1282 direct_dirsamplef_(&c__[c_offset], &arrayi[1], &delta, &c__1, &new__, &length[
1283 length_offset], logfile, &f[3], free, maxi, &point[
1284 1], fcn, &x[1], &l[1], minf, minpos, &u[1], n, maxfunc,
1285 maxdeep, &oops, fmax, ifeasiblef, iinfeasible, fcndata,
1287 if (force_stop && *force_stop) {
1291 /* +-----------------------------------------------------------------------+ */
1292 /* | JG 01/23/01 Added error checking. | */
1293 /* +-----------------------------------------------------------------------+ */
1298 direct_dirdivide_(&new__, &c__0, &length[length_offset], &point[1], &arrayi[1], &
1299 c__1, &list2[list2_offset], &w[1], maxi, &f[3], maxfunc,
1301 direct_dirinsertlist_(&new__, &anchor[-1], &point[1], &f[3], maxi, &
1302 length[length_offset], maxfunc, maxdeep, n, &c__1, jones);
1305 /* +-----------------------------------------------------------------------+ */
1306 /* | SUBROUTINE DIRInitList | */
1307 /* | Initialise the list. | */
1308 /* +-----------------------------------------------------------------------+ */
1309 /* Subroutine */ void direct_dirinitlist_(integer *anchor, integer *free, integer *
1310 point, doublereal *f, integer *maxfunc, const integer *maxdeep)
1312 /* System generated locals */
1315 /* Local variables */
1318 /* f -- values of functions. */
1319 /* anchor -- anchors of lists with deep i */
1320 /* point -- lists */
1321 /* free -- first free position */
1322 /* Parameter adjustments */
1329 for (i__ = -1; i__ <= i__1; ++i__) {
1334 for (i__ = 1; i__ <= i__1; ++i__) {
1335 f[(i__ << 1) + 1] = 0.;
1336 f[(i__ << 1) + 2] = 0.;
1337 point[i__] = i__ + 1;
1341 point[*maxfunc] = 0;
1343 } /* dirinitlist_ */
1345 /* +-----------------------------------------------------------------------+ */
1347 /* | SUBROUTINE DIRPREPRC | */
1349 /* | Subroutine DIRpreprc uses an afine mapping to map the hyper-box given | */
1350 /* | by the constraints on the variable x onto the n-dimensional unit cube.| */
1351 /* | This mapping is done using the following equation: | */
1353 /* | x(i)=x(i)/(u(i)-l(i))-l(i)/(u(i)-l(i)). | */
1355 /* | DIRpreprc checks if the bounds l and u are well-defined. That is, if | */
1357 /* | l(i) < u(i) forevery i. | */
1361 /* | u -- A double-precision vector of length n. The vector | */
1362 /* | containing the upper bounds for the n independent | */
1363 /* | variables. | */
1365 /* | l -- A double-precision vector of length n. The vector | */
1366 /* | containing the lower bounds for the n independent | */
1367 /* | variables. | */
1369 /* | n -- An integer. The dimension of the problem. | */
1373 /* | xs1 -- A double-precision vector of length n, used for scaling | */
1374 /* | and unscaling the vector x. | */
1376 /* | xs2 -- A double-precision vector of length n, used for scaling | */
1377 /* | and unscaling the vector x. | */
1380 /* | oops -- An integer. If an upper bound is less than a lower | */
1381 /* | bound or if the initial point is not in the | */
1382 /* | hyper-box oops is set to 1 and iffco terminates. | */
1384 /* +-----------------------------------------------------------------------+ */
1385 /* Subroutine */ void direct_dirpreprc_(doublereal *u, doublereal *l, integer *n,
1386 doublereal *xs1, doublereal *xs2, integer *oops)
1388 /* System generated locals */
1391 /* Local variables */
1395 /* Parameter adjustments */
1404 for (i__ = 1; i__ <= i__1; ++i__) {
1405 /* +-----------------------------------------------------------------------+ */
1406 /* | Check if the hyper-box is well-defined. | */
1407 /* +-----------------------------------------------------------------------+ */
1408 if (u[i__] <= l[i__]) {
1414 /* +-----------------------------------------------------------------------+ */
1415 /* | Scale the initial iterate so that it is in the unit cube. | */
1416 /* +-----------------------------------------------------------------------+ */
1418 for (i__ = 1; i__ <= i__1; ++i__) {
1419 help = u[i__] - l[i__];
1420 xs2[i__] = l[i__] / help;
1426 /* Subroutine */ void direct_dirheader_(FILE *logfile, integer *version,
1427 doublereal *x, integer *n, doublereal *eps, integer *maxf, integer *
1428 maxt, doublereal *l, doublereal *u, integer *algmethod, integer *
1429 maxfunc, const integer *maxdeep, doublereal *fglobal, doublereal *fglper,
1430 integer *ierror, doublereal *epsfix, integer *iepschange, doublereal *
1431 volper, doublereal *sigmaper)
1433 /* System generated locals */
1436 /* Local variables */
1437 integer imainver, i__, numerrors, isubsubver, ihelp, isubver;
1440 /* +-----------------------------------------------------------------------+ */
1441 /* | Variables to pass user defined data to the function to be optimized. | */
1442 /* +-----------------------------------------------------------------------+ */
1443 /* Parameter adjustments */
1450 fprintf(logfile, "------------------- Log file ------------------\n");
1454 imainver = *version / 100;
1455 ihelp = *version - imainver * 100;
1456 isubver = ihelp / 10;
1457 ihelp -= isubver * 10;
1459 /* +-----------------------------------------------------------------------+ */
1460 /* | JG 01/13/01 Added check for epsilon. If epsilon is smaller 0, we use | */
1461 /* | the update formula from Jones. We then set the flag | */
1462 /* | iepschange to 1, and store the absolute value of eps in | */
1463 /* | epsfix. epsilon is then changed after each iteration. | */
1464 /* +-----------------------------------------------------------------------+ */
1474 /* +-----------------------------------------------------------------------+ */
1475 /* | JG 07/16/01 Removed printout of contents in cdata(1). | */
1476 /* +-----------------------------------------------------------------------+ */
1477 /* write(logfile,*) cdata(1) */
1480 fprintf(logfile, "DIRECT Version %d.%d.%d\n"
1481 " Problem dimension n: %d\n"
1483 " Maximum number of f-evaluations (maxf): %d\n"
1484 " Maximum number of iterations (MaxT): %d\n"
1485 " Value of f_global: %e\n"
1486 " Global percentage wanted: %e\n"
1487 " Volume percentage wanted: %e\n"
1488 " Measure percentage wanted: %e\n",
1489 imainver, isubver, isubsubver, *n, *eps, *maxf, *maxt,
1490 *fglobal, *fglper, *volper, *sigmaper);
1491 fprintf(logfile, *iepschange == 1
1492 ? "Epsilon is changed using the Jones formula.\n"
1493 : "Epsilon is constant.\n");
1494 fprintf(logfile, *algmethod == 0
1495 ? "Jones original DIRECT algorithm is used.\n"
1496 : "Our modification of the DIRECT algorithm is used.\n");
1500 for (i__ = 1; i__ <= i__1; ++i__) {
1501 if (u[i__] <= l[i__]) {
1504 fprintf(logfile, "WARNING: bounds on variable x%d: "
1505 "%g <= xi <= %g\n", i__, l[i__], u[i__]);
1509 fprintf(logfile, "Bounds on variable x%d: "
1510 "%g <= xi <= %g\n", i__, l[i__], u[i__]);
1514 /* +-----------------------------------------------------------------------+ */
1515 /* | If there are to many function evaluations or to many iteration, note | */
1516 /* | this and set the error flag accordingly. Note: If more than one error | */
1517 /* | occurred, we give out an extra message. | */
1518 /* +-----------------------------------------------------------------------+ */
1519 if (*maxf + 20 > *maxfunc) {
1522 "WARNING: The maximum number of function evaluations (%d) is higher than\n"
1523 " the constant maxfunc (%d). Increase maxfunc in subroutine DIRECT\n"
1524 " or decrease the maximum number of function evaluations.\n",
1530 if (logfile) fprintf(logfile, "----------------------------------\n");
1531 if (numerrors == 1) {
1533 fprintf(logfile, "WARNING: One error in the input!\n");
1536 fprintf(logfile, "WARNING: %d errors in the input!\n",
1540 if (logfile) fprintf(logfile, "----------------------------------\n");
1543 fprintf(logfile, "Iteration # of f-eval. minf\n");
1548 /* Subroutine */ void direct_dirsummary_(FILE *logfile, doublereal *x, doublereal *
1549 l, doublereal *u, integer *n, doublereal *minf, doublereal *fglobal,
1550 integer *numfunc, integer *ierror)
1552 /* Local variables */
1555 /* Parameter adjustments */
1562 fprintf(logfile, "-----------------------Summary------------------\n"
1563 "Final function value: %g\n"
1564 "Number of function evaluations: %d\n", *minf, *numfunc);
1565 if (*fglobal > -1e99)
1566 fprintf(logfile, "Final function value is within %g%% of global optimum\n", 100*(*minf - *fglobal) / MAX(1.0, fabs(*fglobal)));
1567 fprintf(logfile, "Index, final solution, x(i)-l(i), u(i)-x(i)\n");
1568 for (i__ = 1; i__ <= *n; ++i__)
1569 fprintf(logfile, "%d, %g, %g, %g\n", i__, x[i__],
1570 x[i__]-l[i__], u[i__] - x[i__]);
1571 fprintf(logfile, "-----------------------------------------------\n");