76 #define min2(a, b) ((a) <= (b) ? (a) : (b))
77 #define max2(a, b) ((a) >= (b) ? (a) : (b))
78 #define max3(a, b, c) max2(max2((a), (b)), (c));
99 1e-20, 1e20, 1e-4, 0.9, 0.9, 1.0e-16,
202 memcpy(param, &
_defparam,
sizeof(*param));
216 int32_t i, j, k, ls, end, bound;
221 const int32_t m = param.
m;
224 float64_t *g = NULL, *gp = NULL, *pg = NULL;
225 float64_t *d = NULL, *w = NULL, *pf = NULL;
247 if (param.
past < 0) {
250 if (param.
delta < 0.) {
259 if (param.
ftol < 0.) {
268 if (param.
gtol < 0.) {
271 if (param.
xtol < 0.) {
319 if (xp == NULL || g == NULL || gp == NULL || d == NULL || w == NULL) {
341 for (i = 0;i < m;++i) {
347 if (it->s == NULL || it->y == NULL) {
354 if (0 < param.
past) {
383 std::copy(pg,pg+n,d);
396 if (xnorm < 1.0) xnorm = 1.0;
397 if (gnorm / xnorm <= param.
epsilon) {
416 ls = linesearch(n, x, &fx, g, d, &step, xp, gp, w, &cd, ¶m);
418 ls = linesearch(n, x, &fx, g, d, &step, xp, pg, w, &cd, ¶m);
426 std::copy(xp,xp+n,x);
427 std::copy(gp,gp+n,g);
452 if (xnorm < 1.0) xnorm = 1.0;
453 if (gnorm / xnorm <= param.
epsilon) {
466 if (param.
past <= k) {
468 rate = (pf[k % param.
past] - fx) / fx;
471 if (rate < param.
delta) {
478 pf[k % param.
past] = fx;
514 bound = (m <= k) ? m : k;
521 std::copy(g, g+n, d);
524 std::copy(pg, pg+n, d);
529 for (i = 0;i < bound;++i) {
541 for (i = 0;i < bound;++i) {
556 if (d[i] * pg[i] >= 0) {
570 if (ptr_fx != NULL) {
578 for (i = 0;i < m;++i) {
630 dgtest = param->
ftol * dginit;
633 std::copy(xp,xp+n,x);
641 if (*f > finit + *stp * dgtest) {
652 if (dg < param->wolfe * dginit) {
661 if(dg > -param->
wolfe * dginit) {
670 if (*stp < param->min_step) {
703 int32_t i, count = 0;
713 for (i = 0;i < n;++i) {
714 wp[i] = (xp[i] == 0.) ? -gp[i] : xp[i];
719 std::copy(xp,xp+n,x);
735 for (i = 0;i < n;++i) {
736 dgtest += (x[i] - xp[i]) * gp[i];
739 if (*f <= finit + param->ftol * dgtest) {
744 if (*stp < param->min_step) {
778 int32_t brackt, stage1, uinfo = 0;
804 dgtest = param->
ftol * dginit;
806 prev_width = 2.0 * width;
827 stmin =
min2(stx, sty);
828 stmax =
max2(stx, sty);
831 stmax = *stp + 4.0 * (*stp - stx);
835 if (*stp < param->min_step) *stp = param->
min_step;
842 if ((brackt && ((*stp <= stmin || stmax <= *stp) || param->
max_linesearch <= count + 1 || uinfo != 0)) || (brackt && (stmax - stmin <= param->
xtol * stmax))) {
850 std::copy(xp,xp+n,x);
857 ftest1 = finit + *stp * dgtest;
861 if (brackt && ((*stp <= stmin || stmax <= *stp) || uinfo != 0)) {
865 if (*stp == param->
max_step && *f <= ftest1 && dg <= dgtest) {
869 if (*stp == param->
min_step && (ftest1 < *f || dgtest <= dg)) {
873 if (brackt && (stmax - stmin) <= param->
xtol * stmax) {
881 if (*f <= ftest1 && fabs(dg) <= param->
gtol * (-dginit)) {
890 if (stage1 && *f <= ftest1 &&
min2(param->
ftol, param->
gtol) * dginit <= dg) {
901 if (stage1 && ftest1 < *f && *f <= fx) {
903 fm = *f - *stp * dgtest;
904 fxm = fx - stx * dgtest;
905 fym = fy - sty * dgtest;
918 stmin, stmax, &brackt
922 fx = fxm + stx * dgtest;
923 fy = fym + sty * dgtest;
935 stmin, stmax, &brackt
943 if (0.66 * prev_width <= fabs(sty - stx)) {
944 *stp = stx + 0.5 * (sty - stx);
947 width = fabs(sty - stx);
959 #define USES_MINIMIZER \
960 float64_t a, d, gamma, theta, p, q, r, s;
972 #define CUBIC_MINIMIZER(cm, u, fu, du, v, fv, dv) \
974 theta = ((fu) - (fv)) * 3 / d + (du) + (dv); \
981 gamma = s * sqrt(a * a - ((du) / s) * ((dv) / s)); \
982 if ((v) < (u)) gamma = -gamma; \
983 p = gamma - (du) + theta; \
984 q = gamma - (du) + gamma + (dv); \
1000 #define CUBIC_MINIMIZER2(cm, u, fu, du, v, fv, dv, xmin, xmax) \
1002 theta = ((fu) - (fv)) * 3 / d + (du) + (dv); \
1006 s = max3(p, q, r); \
1009 gamma = s * sqrt(max2(0, a * a - ((du) / s) * ((dv) / s))); \
1010 if ((u) < (v)) gamma = -gamma; \
1011 p = gamma - (dv) + theta; \
1012 q = gamma - (dv) + gamma + (du); \
1014 if (r < 0. && gamma != 0.) { \
1015 (cm) = (v) - r * d; \
1016 } else if (a < 0) { \
1031 #define QUARD_MINIMIZER(qm, u, fu, du, v, fv) \
1033 (qm) = (u) + (du) / (((fu) - (fv)) / a + (du)) / 2 * a;
1043 #define QUARD_MINIMIZER2(qm, u, du, v, dv) \
1045 (qm) = (v) + (dv) / ((dv) - (du)) * a;
1047 #define fsigndiff(x, y) (*(x) * (*(y) / fabs(*(y))) < 0.)
1102 if (*t <=
min2(*x, *y) ||
max2(*x, *y) <= *t) {
1106 if (0. <= *dx * (*t - *x)) {
1130 if (fabs(mc - *x) < fabs(mq - *x)) {
1133 newt = mc + 0.5 * (mq - mc);
1146 if (fabs(mc - *t) > fabs(mq - *t)) {
1151 }
else if (fabs(*dt) < fabs(*dx)) {
1167 if (fabs(*t - mc) < fabs(*t - mq)) {
1173 if (fabs(*t - mc) > fabs(*t - mq)) {
1189 }
else if (*x < *t) {
1226 if (tmax < newt) newt = tmax;
1227 if (newt < tmin) newt = tmin;
1233 if (*brackt && bound) {
1234 mq = *x + 0.66 * (*y - *x);
1236 if (mq < newt) newt = mq;
1238 if (newt < mq) newt = mq;
1253 const int32_t start,
1260 for (i = start;i < n;++i) {
1273 const int32_t start,
1280 for (i = 0;i < start;++i) {
1285 for (i = start;i < end;++i) {
1289 }
else if (0. < x[i]) {
1296 }
else if (c < g[i]) {
1305 for (i = end;i < n;++i) {
1313 const int32_t start,
1319 for (i = start;i < end;++i) {
1320 if (d[i] * sign[i] <= 0) {