45 #define dmax(a1,a2) ((a1 > a2) ? a1 : a2)
46 #define dmin(a1,a2) ((a1 < a2) ? a1 : a2)
73 C_INT i__1, i__2, i__3;
82 static C_INT i__, j, k;
83 static C_FLOAT64 s, t, y[100], large, z__[100], small, value, f1;
91 static C_INT km1, im1;
95 C_FLOAT64 lastValue = std::numeric_limits<C_FLOAT64>::infinity();
161 small = *machep * *machep;
162 vsmall = small * small;
164 vlarge = 1. / vsmall;
193 t = small + fabs(*t0);
208 for (i__ = 1; i__ <= i__1; ++i__)
212 for (j = 1; j <= i__2; ++j)
215 q_1.
v[i__ + j * 100 - 101] = 0.;
219 q_1.
v[i__ + i__ * 100 - 101] = 1.;
226 for (i__ = 1; i__ <= i__1; ++i__)
228 q_1.
q0[i__ - 1] = x[i__];
230 q_1.
q1[i__ - 1] = x[i__];
235 print_(n, &x[1], prin, fmin);
247 min_(n, &
c__1, &
c__2, d__, &s, &value, &
c_false, f, &x[1], &t,
257 for (i__ = 1; i__ <= i__1; ++i__)
265 if (sf > d__[0] * .9 && sf * .9 < d__[0])
272 for (i__ = 2; i__ <= i__1; ++i__)
282 for (k = 2; k <= i__1; ++k)
286 for (i__ = 1; i__ <= i__2; ++i__)
314 for (i__ = 1; i__ <= i__2; ++i__)
320 for (j = 1; j <= i__3; ++j)
323 x[j] += s *
q_1.
v[j + i__ * 100 - 101];
338 for (k2 = k; k2 <= i__2; ++k2)
344 x[1], &t, machep, &h__);
355 d__1 = s + z__[k2 - 1];
356 s = d__[k2 - 1] * (d__1 * d__1);
370 if (illc || df >= (d__1 = *machep * 100 *
global_1.
fx, fabs(d__1)))
382 if (k == 2 && *prin > 1)
393 for (k2 = 1; k2 <= i__2; ++k2)
398 x[1], &t, machep, &h__);
407 for (i__ = 1; i__ <= i__2; ++i__)
437 for (ii = 1; ii <= i__2; ++ii)
442 for (j = 1; j <= i__3; ++j)
445 q_1.
v[j + (i__ + 1) * 100 - 101] =
q_1.
v[j + i__ * 100 - 101];
449 d__[i__] = d__[i__ - 1];
456 for (i__ = 1; i__ <= i__2; ++i__)
459 q_1.
v[i__ + k * 100 - 101] = y[i__ - 1] / lds;
477 for (i__ = 1; i__ <= i__2; ++i__)
480 q_1.
v[i__ + k * 100 - 101] = -
q_1.
v[i__ + k * 100 - 101];
493 print_(n, &x[1], prin, fmin);
499 for (i__ = 1; i__ <= i__2; ++i__)
507 t2 = m2 * sqrt(t2) + t;
534 quad_(n, f, &x[1], &t, machep, &h__);
538 for (i__ = 1; i__ <= i__1; ++i__)
540 d__[i__ - 1] = 1. / sqrt(d__[i__ - 1]);
542 if (dn < d__[i__ - 1])
557 for (j = 1; j <= i__1; ++j)
562 for (i__ = 1; i__ <= i__2; ++i__)
565 q_1.
v[i__ + j * 100 - 101] = s *
q_1.
v[i__ + j * 100 - 101];
579 for (i__ = 1; i__ <= i__2; ++i__)
584 for (j = 1; j <= i__1; ++j)
587 sl +=
q_1.
v[i__ + j * 100 - 101] *
q_1.
v[i__ + j * 100 - 101];
590 z__[i__ - 1] = sqrt(sl);
592 if (z__[i__ - 1] < m4)
597 if (s > z__[i__ - 1])
607 for (i__ = 1; i__ <= i__2; ++i__)
609 sl = s / z__[i__ - 1];
610 z__[i__ - 1] = 1. / sl;
612 if (z__[i__ - 1] <= scbd)
622 for (j = 1; j <= i__1; ++j)
625 q_1.
v[i__ + j * 100 - 101] = sl *
q_1.
v[i__ + j * 100 - 101];
638 for (i__ = 2; i__ <= i__2; ++i__)
643 for (j = 1; j <= i__1; ++j)
645 s =
q_1.
v[i__ + j * 100 - 101];
646 q_1.
v[i__ + j * 100 - 101] =
q_1.
v[j + i__ * 100 - 101];
648 q_1.
v[j + i__ * 100 - 101] = s;
670 for (i__ = 1; i__ <= i__2; ++i__)
675 for (j = 1; j <= i__1; ++j)
678 q_1.
v[i__ + j * 100 - 101] = s *
q_1.
v[i__ + j * 100 - 101];
686 for (i__ = 1; i__ <= i__2; ++i__)
691 for (j = 1; j <= i__1; ++j)
695 d__1 =
q_1.
v[j + i__ * 100 - 101];
700 d__[i__ - 1] = s * d__[i__ - 1];
704 for (j = 1; j <= i__1; ++j)
707 q_1.
v[j + i__ * 100 - 101] = s *
q_1.
v[j + i__ * 100 - 101];
716 for (i__ = 1; i__ <= i__2; ++i__)
718 dni = dn * d__[i__ - 1];
730 d__[i__ - 1] = 1 / (dni * dni);
733 d__[i__ - 1] = vlarge;
736 d__[i__ - 1] = vsmall;
758 if (*prin > 1 && scbd > 1.)
804 C_INT ab_dim1, ab_offset, i__1, i__2, i__3;
808 static C_FLOAT64 temp, c__, e[100], f, g, h__;
809 static C_INT i__, j, k, l;
811 static C_INT l2, ii, kk, kt, ll2, lp1;
824 ab_offset = ab_dim1 + 1;
838 for (i__ = 1; i__ <= i__1; ++i__)
845 for (j = i__; j <= i__2; ++j)
849 d__1 = ab[j + i__ * ab_dim1];
860 f = ab[i__ + i__ * ab_dim1];
869 ab[i__ + i__ * ab_dim1] = f - g;
878 for (j = l; j <= i__2; ++j)
883 for (k = i__; k <= i__3; ++k)
886 f += ab[k + i__ * ab_dim1] * ab[k + j * ab_dim1];
892 for (k = i__; k <= i__3; ++k)
895 ab[k + j * ab_dim1] += f * ab[k + i__ * ab_dim1];
910 for (j = l; j <= i__3; ++j)
913 s += ab[i__ + j * ab_dim1] * ab[i__ + j * ab_dim1];
929 f = ab[i__ + (i__ + 1) * ab_dim1];
945 ab[i__ + (i__ + 1) * ab_dim1] = f - g;
948 for (j = l; j <= i__3; ++j)
951 e[j - 1] = ab[i__ + j * ab_dim1] / h__;
956 for (j = l; j <= i__3; ++j)
961 for (k = l; k <= i__2; ++k)
964 s += ab[j + k * ab_dim1] * ab[i__ + k * ab_dim1];
969 for (k = l; k <= i__2; ++k)
972 ab[j + k * ab_dim1] += s * e[k - 1];
977 y = (d__1 = q[i__], fabs(d__1)) + (d__2 = e[i__ - 1], fabs(d__2));
987 ab[*n + *n * ab_dim1] = 1.;
992 for (ii = 2; ii <= i__1; ++ii)
1001 h__ = ab[i__ + (i__ + 1) * ab_dim1] * g;
1004 for (j = l; j <= i__2; ++j)
1007 ab[j + i__ * ab_dim1] = ab[i__ + j * ab_dim1] / h__;
1012 for (j = l; j <= i__2; ++j)
1017 for (k = l; k <= i__3; ++k)
1020 s += ab[i__ + k * ab_dim1] * ab[k + j * ab_dim1];
1025 for (k = l; k <= i__3; ++k)
1028 ab[k + j * ab_dim1] += s * ab[k + i__ * ab_dim1];
1035 for (j = l; j <= i__3; ++j)
1037 ab[i__ + j * ab_dim1] = 0.;
1039 ab[j + i__ * ab_dim1] = 0.;
1042 ab[i__ + i__ * ab_dim1] = 1.;
1053 for (kk = 1; kk <= i__1; ++kk)
1066 printf(
"QR FAILED\n");
1070 for (ll2 = 1; ll2 <= i__3; ++ll2)
1075 if ((d__1 = e[l - 1], fabs(d__1)) <= eps)
1085 if ((d__1 = q[l - 1], fabs(d__1)) <= eps)
1100 for (i__ = l; i__ <= i__3; ++i__)
1103 e[i__ - 1] = c__ * e[i__ - 1];
1113 if (fabs(f) < fabs(g))
1133 h__ = fabs(f) * sqrt(d__1 * d__1 + 1);
1138 h__ = fabs(g) * sqrt(d__1 * d__1 + 1);
1169 f = ((y - z__) * (y + z__) + (g - h__) * (g + h__)) / (h__ * 2 * y);
1170 g = sqrt(f * f + 1.);
1178 f = ((x - z__) * (x + z__) + h__ * (y / temp - h__)) / x;
1191 for (i__ = lp1; i__ <= i__3; ++i__)
1198 if (fabs(f) < fabs(h__))
1218 z__ = fabs(f) * sqrt(d__1 * d__1 + 1);
1223 z__ = fabs(h__) * sqrt(d__1 * d__1 + 1);
1237 f = x * c__ + g * s;
1238 g = -x * s + g * c__;
1243 for (j = 1; j <= i__2; ++j)
1245 x = ab[j + (i__ - 1) * ab_dim1];
1246 z__ = ab[j + i__ * ab_dim1];
1247 ab[j + (i__ - 1) * ab_dim1] = x * c__ + z__ * s;
1249 ab[j + i__ * ab_dim1] = -x * s + z__ * c__;
1252 if (fabs(f) < fabs(h__))
1272 z__ = fabs(f) * sqrt(d__1 * d__1 + 1);
1277 z__ = fabs(h__) * sqrt(d__1 * d__1 + 1);
1291 f = c__ * g + s * y;
1293 x = -s * g + c__ * y;
1312 for (j = 1; j <= i__3; ++j)
1315 ab[j + k * ab_dim1] = -ab[j + k * ab_dim1];
1324 q[1] = ab[ab_dim1 + 1];
1325 ab[ab_dim1 + 1] = 1.;
1339 static C_INT i__, k;
1340 static C_FLOAT64 s, small, d1, f0, f2, m2, m4, t2, x2, fm;
1361 small = d__1 * d__1;
1375 for (i__ = 1; i__ <= i__1; ++i__)
1400 t2 =
dmax(t2, small);
1402 d__1 = t2, d__2 = *h__ * .01;
1403 t2 =
dmin(d__1, d__2);
1405 if (!(*fk) || *f1 > fm)
1414 if (*fk && fabs(*x1) >= t2)
1463 *d2 = (x2 * (*f1 - f0) - *x1 * (f2 - f0)) / (*x1 * x2 * (*x1 - x2));
1466 d1 = (*f1 - f0) / *x1 - *x1 * *d2;
1484 x2 = d1 * -.5 / *d2;
1487 if (fabs(x2) <= *h__)
1510 if (k >= *nits || f2 <= f0)
1518 if (f0 < *f1 && *x1 * x2 > 0.)
1541 if ((d__1 = x2 * (x2 - *x1), fabs(d__1)) <= small)
1546 *d2 = (x2 * (*f1 - f0) - *x1 * (fm - f0)) / (*x1 * x2 * (*x1 - x2));
1582 for (i__ = 1; i__ <= i__1; ++i__)
1585 x[i__] += *x1 *
q_1.
v[i__ + *j * 100 - 101];
1616 for (i__ = 1; i__ <= i__1; ++i__)
1619 t[i__ - 1] = x[i__] + *l *
q_1.
v[i__ + *j * 100 - 101];
1630 for (i__ = 1; i__ <= i__1; ++i__)
1640 ret_val = (*f)(t, n);
1648 C_INT v_dim1, v_offset, i__1, i__2;
1651 static C_INT i__, j, k;
1653 static C_INT ip1, nm1;
1660 v_offset = v_dim1 + 1;
1673 for (i__ = 1; i__ <= i__1; ++i__)
1680 for (j = ip1; j <= i__2; ++j)
1702 for (j = 1; j <= i__2; ++j)
1704 s = v[j + i__ * v_dim1];
1705 v[j + i__ * v_dim1] = v[j + k * v_dim1];
1707 v[j + k * v_dim1] = s;
1739 for (i__ = 1; i__ <= i__1; ++i__)
1742 l =
q_1.
q1[i__ - 1];
1744 q_1.
q1[i__ - 1] = s;
1761 min_(n, &
c__0, &
c__2, &s, &l, &value, &
c_true, f, &x[1], t, machep,
1776 for (i__ = 1; i__ <= i__1; ++i__)
1778 s =
q_1.
q0[i__ - 1];
1779 q_1.
q0[i__ - 1] = x[i__];
1799 switch ((
int)*option)
1808 printf(
"THE SECOND DIFFERENCE ARRAY D[*] IS :\n");
1811 for (i = 1; i <= i__1; ++i)
1813 printf(
"%g\n", v[i]);
1818 printf(
"THE SCALE FACTORS ARE:\n");
1821 for (i = 1; i <= i__1; ++i)
1823 printf(
"%g\n", v[i]);
1828 printf(
"THE APPROXIMATING QUADRATIC FORM HAS THE PRINCEPAL VALUES:\n");
1831 for (i = 1; i <= i__1; ++i)
1833 printf(
"%g\n", v[i]);
1841 for (i = 1; i <= i__1; ++i)
1843 printf(
"%g\n", v[i]);
1869 printf(
" LINEAR SEARCHES, THE FUNCTION HAS BEEN EVALUATED ");
1871 printf(
"THE SMALLEST VALUE FOUND IS F(X) = ");
1881 printf(
"log (f(x)) - ");
1882 printf(
"%g", *fmin);
1887 printf(
"LOG (F(x)) -- ");
1888 printf(
"%g", *fmin);
1889 printf(
" IS UNDERFINED\n");
1892 if (*n > 4 && *prin <= 2)
1899 for (i = 1; i <= i__1; ++i)
1902 printf(
"%g\n", x[i]);
1915 C_INT v_dim1, v_offset, i__1, i__2;
1918 C_INT i, j, low, upp;
1925 v_offset = v_dim1 + 1;
1932 switch ((
int)*option)
1939 printf(
"HE NEW DIRECTIONS ARE:\n");
1942 printf(
"AND THE PRINCIPAL AXES:\n");
1952 for (i = 1; i <= i__1; ++i)
1958 for (j = low; j <= i__2; ++j)
1960 printf(
" %12g", v[i*v_dim1 + j]);
C_FLOAT64 praxis_(C_FLOAT64 *t0, C_FLOAT64 *machep, C_FLOAT64 *h0, C_INT *n, C_INT *prin, C_FLOAT64 *x, FPraxis *f, C_FLOAT64 *fmin)
int maprnt_(C_INT *, C_FLOAT64 *, C_INT *, C_INT *)
int sort_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *)
int minfit_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
int vcprnt_(C_INT *, C_FLOAT64 *, C_INT *)
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
virtual C_FLOAT64 getRandomCC()
int quad_(C_INT *, FPraxis *f, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
int min_(C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, bool *, FPraxis *f, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
int print_(C_INT *n, C_FLOAT64 *x, C_INT *prin, C_FLOAT64 *fmin)
C_FLOAT64 flin_(C_INT *, C_INT *, C_FLOAT64 *, FPraxis *, C_FLOAT64 *, C_INT *)