75 #define min2(a, b) ((a) <= (b) ? (a) : (b))
76 #define max2(a, b) ((a) >= (b) ? (a) : (b))
77 #define max3(a, b, c) max2(max2((a), (b)), (c));
98 1e-20, 1e20, 1e-4, 0.9, 0.9, 1.0e-16,
201 memcpy(param, &
_defparam,
sizeof(*param));
215 int32_t i, j, k, ls, end, bound;
220 const int32_t m = param.
m;
223 float64_t *g = NULL, *gp = NULL, *pg = NULL;
224 float64_t *d = NULL, *w = NULL, *pf = NULL;
246 if (param.
past < 0) {
249 if (param.
delta < 0.) {
258 if (param.
ftol < 0.) {
267 if (param.
gtol < 0.) {
270 if (param.
xtol < 0.) {
318 if (xp == NULL || g == NULL || gp == NULL || d == NULL || w == NULL) {
340 for (i = 0;i < m;++i) {
346 if (it->s == NULL || it->y == NULL) {
353 if (0 < param.
past) {
382 std::copy(pg,pg+n,d);
395 if (xnorm < 1.0) xnorm = 1.0;
396 if (gnorm / xnorm <= param.
epsilon) {
415 ls = linesearch(n, x, &fx, g, d, &step, xp, gp, w, &cd, ¶m);
417 ls = linesearch(n, x, &fx, g, d, &step, xp, pg, w, &cd, ¶m);
425 std::copy(xp,xp+n,x);
426 std::copy(gp,gp+n,g);
451 if (xnorm < 1.0) xnorm = 1.0;
452 if (gnorm / xnorm <= param.
epsilon) {
465 if (param.
past <= k) {
467 rate = (pf[k % param.
past] - fx) / fx;
470 if (rate < param.
delta) {
477 pf[k % param.
past] = fx;
513 bound = (m <= k) ? m : k;
520 std::copy(g, g+n, d);
523 std::copy(pg, pg+n, d);
528 for (i = 0;i < bound;++i) {
540 for (i = 0;i < bound;++i) {
555 if (d[i] * pg[i] >= 0) {
569 if (ptr_fx != NULL) {
577 for (i = 0;i < m;++i) {
629 dgtest = param->
ftol * dginit;
632 std::copy(xp,xp+n,x);
640 if (*f > finit + *stp * dgtest) {
651 if (dg < param->wolfe * dginit) {
660 if(dg > -param->
wolfe * dginit) {
669 if (*stp < param->min_step) {
702 int32_t i, count = 0;
712 for (i = 0;i < n;++i) {
713 wp[i] = (xp[i] == 0.) ? -gp[i] : xp[i];
718 std::copy(xp,xp+n,x);
734 for (i = 0;i < n;++i) {
735 dgtest += (x[i] - xp[i]) * gp[i];
738 if (*f <= finit + param->ftol * dgtest) {
743 if (*stp < param->min_step) {
777 int32_t brackt, stage1, uinfo = 0;
803 dgtest = param->
ftol * dginit;
805 prev_width = 2.0 * width;
826 stmin =
min2(stx, sty);
827 stmax =
max2(stx, sty);
830 stmax = *stp + 4.0 * (*stp - stx);
834 if (*stp < param->min_step) *stp = param->
min_step;
841 if ((brackt && ((*stp <= stmin || stmax <= *stp) || param->
max_linesearch <= count + 1 || uinfo != 0)) || (brackt && (stmax - stmin <= param->
xtol * stmax))) {
849 std::copy(xp,xp+n,x);
856 ftest1 = finit + *stp * dgtest;
860 if (brackt && ((*stp <= stmin || stmax <= *stp) || uinfo != 0)) {
864 if (*stp == param->
max_step && *f <= ftest1 && dg <= dgtest) {
868 if (*stp == param->
min_step && (ftest1 < *f || dgtest <= dg)) {
872 if (brackt && (stmax - stmin) <= param->
xtol * stmax) {
880 if (*f <= ftest1 && fabs(dg) <= param->
gtol * (-dginit)) {
889 if (stage1 && *f <= ftest1 &&
min2(param->
ftol, param->
gtol) * dginit <= dg) {
900 if (stage1 && ftest1 < *f && *f <= fx) {
902 fm = *f - *stp * dgtest;
903 fxm = fx - stx * dgtest;
904 fym = fy - sty * dgtest;
917 stmin, stmax, &brackt
921 fx = fxm + stx * dgtest;
922 fy = fym + sty * dgtest;
934 stmin, stmax, &brackt
942 if (0.66 * prev_width <= fabs(sty - stx)) {
943 *stp = stx + 0.5 * (sty - stx);
946 width = fabs(sty - stx);
958 #define USES_MINIMIZER \
959 float64_t a, d, gamma, theta, p, q, r, s;
971 #define CUBIC_MINIMIZER(cm, u, fu, du, v, fv, dv) \
973 theta = ((fu) - (fv)) * 3 / d + (du) + (dv); \
980 gamma = s * sqrt(a * a - ((du) / s) * ((dv) / s)); \
981 if ((v) < (u)) gamma = -gamma; \
982 p = gamma - (du) + theta; \
983 q = gamma - (du) + gamma + (dv); \
999 #define CUBIC_MINIMIZER2(cm, u, fu, du, v, fv, dv, xmin, xmax) \
1001 theta = ((fu) - (fv)) * 3 / d + (du) + (dv); \
1005 s = max3(p, q, r); \
1008 gamma = s * sqrt(max2(0, a * a - ((du) / s) * ((dv) / s))); \
1009 if ((u) < (v)) gamma = -gamma; \
1010 p = gamma - (dv) + theta; \
1011 q = gamma - (dv) + gamma + (du); \
1013 if (r < 0. && gamma != 0.) { \
1014 (cm) = (v) - r * d; \
1015 } else if (a < 0) { \
1030 #define QUARD_MINIMIZER(qm, u, fu, du, v, fv) \
1032 (qm) = (u) + (du) / (((fu) - (fv)) / a + (du)) / 2 * a;
1042 #define QUARD_MINIMIZER2(qm, u, du, v, dv) \
1044 (qm) = (v) + (dv) / ((dv) - (du)) * a;
1046 #define fsigndiff(x, y) (*(x) * (*(y) / fabs(*(y))) < 0.)
1101 if (*t <=
min2(*x, *y) ||
max2(*x, *y) <= *t) {
1105 if (0. <= *dx * (*t - *x)) {
1129 if (fabs(mc - *x) < fabs(mq - *x)) {
1132 newt = mc + 0.5 * (mq - mc);
1145 if (fabs(mc - *t) > fabs(mq - *t)) {
1150 }
else if (fabs(*dt) < fabs(*dx)) {
1166 if (fabs(*t - mc) < fabs(*t - mq)) {
1172 if (fabs(*t - mc) > fabs(*t - mq)) {
1188 }
else if (*x < *t) {
1225 if (tmax < newt) newt = tmax;
1226 if (newt < tmin) newt = tmin;
1232 if (*brackt && bound) {
1233 mq = *x + 0.66 * (*y - *x);
1235 if (mq < newt) newt = mq;
1237 if (newt < mq) newt = mq;
1252 const int32_t start,
1259 for (i = start;i < n;++i) {
1272 const int32_t start,
1279 for (i = 0;i < start;++i) {
1284 for (i = start;i < end;++i) {
1288 }
else if (0. < x[i]) {
1295 }
else if (c < g[i]) {
1304 for (i = end;i < n;++i) {
1312 const int32_t start,
1318 for (i = start;i < end;++i) {
1319 if (d[i] * sign[i] <= 0) {