cg.c

InsightToolkit-1.4.0(有大量的优化算法程序)

    if (da > a6 * g) {
        goto L410;
    }
    if (da >= 0.) {
        goto L560;
    }
    r = a;
    q = 0.f;
    for (i = 0; i < j; ++i) {
        if (y[i] > a) {
            goto L370;
        }
        if (y[i] <= q || y[i] == a) {
            continue; /* next i */
        }
        q = y[i];
    }
    if (a <= a8 * q) {
        goto L560;
    }
    q = a;
L340:
    ++nd;
    if (nd > 25) {
        goto L610;
    }
    q *= a3;
    p = fv_(&q, x, h, n, value);
    f1 = fa;
    ins_(&q, &p, &a, &b, &c, &fa, &fb, &fc, &j, y, z);
    if (p < f1) {
        goto L340;
    }
    if (a > r) {
        goto L360;
    }
    for (i = 0; i < *n; ++i) {
        h[i+*n] = x[i] + a * h[i];
    }
    goto L560;
L360:
    da = fd_(&a, x, h, n, grad);
    if (da > a6 * g) {
        goto L410;
    }
    goto L560;
L370:
    q = y[i];
    for (k = i; k < j; ++k) {
        if (y[k] <= a) {
            continue; /* next k */
        }
        if (y[k] < q) {
            q = y[k];
        }
    }
    if (q <= a5 * a) {
        goto L560;
    }
    f0 = log(q / a);
    s = (doublereal) ((integer) (f0 * l3 + .999f));
    f0 = exp(f0 / s) * 1.001f;
    s = a;
L390:
    s *= f0;
    if (s >= q) {
        goto L320;
    }
    p = fv_(&s, x, h, n, value);
    f1 = fa;
    ins_(&s, &p, &a, &b, &c, &fa, &fb, &fc, &j, y, z);
    if (p < f1) {
        goto L390;
    }
    if (a > r) {
        goto L320;
    }
    for (i = 0; i < *n; ++i) {
        h[i+*n] = x[i] + a * h[i];
    }
    goto L560;
L410:
    b = 0.f;
    k = 0;
    i = k;
L420:
    ++i;
    if (i+1 > j) {
        goto L430;
    }
    if (y[i] >= a) {
        goto L420;
    }
    if (y[i] < b) {
        goto L420;
    }
    b = y[i];
    k = i;
    goto L420;
L430:
    fb = z[k];
    db = d;
    if (b != 0.) {
        db = fd_(&b, x, h, n, grad);
    }
    w = abs(b - a) * 2.f;
    cub_(&c, &a, &b, &fa, &fb, &da, &db);
    nc = 1;
    goto L480;
L450:
    w *= .5f;
    if (w < abs(c0 - c)) {
        goto L550;
    }
    if (c0 < c) {
        goto L460;
    }
    if (d0 >= d) {
        goto L470;
    }
    goto L550;
L460:
    if (d0 > d) {
        goto L550;
    }
L470:
    cub_(&c, &c, &c0, &f, &f0, &d, &d0);
    ++nc;
    if (nc > 30) {
        goto L600;
    }
L480:
    r = max(a,b);
    s = min(a,b);
    if (c > r) {
        goto L490;
    }
    if (c > s) {
        goto L500;
    }
    c = s + (s - c);
    s = (a + b) * .5f;
    if (c > s) {
        c = s;
    }
    goto L500;
L490:
    c = r - (c - r);
    s = (a + b) * .5f;
    if (c < s) {
        c = s;
    }
L500:
    c0 = a;
    f0 = fa;
    d0 = da;
    fvd_(&f, &d, &c, x, h, n, both);
    if (f < fa) {
        goto L510;
    }
    b = c;
    fb = f;
    db = d;
    goto L450;
L510:
    if (c < a) {
        goto L540;
    }
    if (d < 0.) {
        goto L530;
    }
L520:
    b = a;
    fb = fa;
    db = da;
L530:
    a = c;
    fa = f;
    da = d;
    if (d > a6 * g) {
        goto L450;
    }
    goto L560;
L540:
    if (d < 0.) {
        goto L520;
    }
    goto L530;
L550:
    c = (a + b) * .5f;
    ++nb;
    w = abs(b - a);
    goto L500;
L560:
    *e = 0.f;
    for (i = 0; i < *n; ++i) {
        if (abs(h[i+*n*2]) > *e) {
            *e = abs(h[i+*n*2]);
        }
        x[i] = h[i+*n];
    }
    ++(*it);
    if (*e <= *t) {
        goto L660;
    }
    if (*it >= *limit) {
        goto L660;
    }
    f = fa;
    d = da;
    a *= a7;
    (*pre)(&h[*n], &h[*n*2]);
    r = 0.f;
    for (i = 0; i < *n; ++i) {
        r += h[i+*n] * h[i+*n*2];
    }
    if (r < 0.) {
        goto L620;
    }
    s = r / g;
    g = r;
    ++l;
    if (l >= *m) {
        goto L50;
    }
    d = 0.f;
    for (i = 0; i < *n; ++i) {
        h[i] = -h[i+*n] + s * h[i];
        d += h[i] * h[i+*n*2];
    }
    goto L70;
L600:
    if (d < g) {
        goto L560;
    }
#ifdef DEBUG
    s_wsle(&io___43);
    do_lio(&c__9, &c__1, "UNABLE TO OBTAIN DESCENT DIRECTION", 34L);
    e_wsle();
#endif
    printf("UNABLE TO OBTAIN DESCENT DIRECTION\n"); assert(0);
/*     s_stop("", 0L); */
L610:
#ifdef DEBUG
    s_wsle(&io___44);
    do_lio(&c__9, &c__1, "THE FUNCTION DECREASES WITH NO MINIMUM", 38L);
    e_wsle();
#endif
    printf("THE FUNCTION DECREASES WITH NO MINIMUM\n"); assert(0);
/*     s_stop("", 0L); */
L620:
#ifdef DEBUG
    s_wsle(&io___45);
    do_lio(&c__9, &c__1, "PRECONDITIONER NOT POSITIVE DEFINITE", 36L);
    e_wsle();
#endif
    printf("PRECONDITIONER NOT POSITIVE DEFINITE\n"); assert(0);
/*     s_stop("", 0L); */
L630:
    /* Computing 25th power */
    d__1 = a3, d__1 *= d__1, d__1 *= d__1, d__1 *= d__1,
    q *= a3 * d__1 * d__1 * d__1;
    nd = 0;
L640:
    ++nd;
    if (nd > 25) {
        goto L650;
    }
    q *= a3;
    p = fv_(&q, x, h, n, value);
    ins_(&q, &p, &a, &b, &c, &fa, &fb, &fc, &j, y, z);
    if (p - f > v * q) {
        goto L640;
    }
    goto L135;
L650:
#ifdef DEBUG
    s_wsle(&io___46);
    do_lio(&c__9, &c__1, "UNABLE TO SATISFY ARMIJO CONDITION", 34L);
    e_wsle();
#endif
    printf("UNABLE TO SATISFY ARMIJO CONDITION\n");
    return;
L660:
    *step = a;
} /* cg_ */

doublereal fv_(a, x, h, n, value)
doublereal *a, *x, *h;
const integer *n;
doublereal (*value) (doublereal*);
{
    /* Local variables */
    static integer i;

    for (i = 0; i < *n; ++i) {
        h[i+*n] = x[i] + *a * h[i];
    }
    return (*value)(&h[*n]);
} /* fv_ */

doublereal fd_(a, x, h, n, grad)
doublereal *a, *x, *h;
const integer *n;
void (*grad) (doublereal*,doublereal*);
{
    /* Local variables */
    static doublereal d;
    static integer i;

    for (i = 0; i < *n; ++i) {
        h[i+*n] = x[i] + *a * h[i];
    }
    (*grad)(&h[*n*2], &h[*n]);
    d = 0.f;
    for (i = 0; i < *n; ++i) {
        d += h[i] * h[i+*n*2];
    }
    return d;
} /* fd_ */

/* Subroutine */ void fvd_(v, d, a, x, h, n, both)
doublereal *v, *d, *a, *x, *h;
const integer *n;
/* Subroutine */ void (*both) (doublereal*,doublereal*,doublereal*);
{
    /* Local variables */
    static integer i;

    for (i = 0; i < *n; ++i) {
        h[i+*n] = x[i] + *a * h[i];
    }
    (*both)(v, &h[*n*2], &h[*n]);
    *d = 0.f;
    for (i = 0; i < *n; ++i) {
        *d += h[i] * h[i+*n*2];
    }
    return;
} /* fvd_ */

/* Subroutine */ void cub_(x, a, b, c, d, e, f)
doublereal *x, *a, *b, *c, *d, *e, *f;
{
    /* Local variables */
    static doublereal g, v, w, y, z;

    g = *b - *a;
    if (g == 0.) {
        goto L50;
    }
    v = *e + *f - (*d - *c) * 3 / g;
    w = v * v - *e * *f;
    if (w < 0.) {
        w = 0.f;
    }
    w = sqrt(w);
    w = d_sign(&w, &g);
    y = *e + v;
    z = *f + v;
    if (d_sign(&y, &g) != y) {
        goto L30;
    }
    if (d_sign(&z, &g) != z) {
        goto L20;
    }
    if (z == 0.) {
        goto L20;
    }
L10:
    *x = *b - g * *f / (z + w);
    return;
L20:
    if (*c < *d) {
        *x = *a;
    }
    if (*c >= *d) {
        *x = *b;
    }
    return;
L30:
    if (d_sign(&z, &g) != z) {
        goto L40;
    }
    if (abs(*e) > abs(*f)) {
        goto L10;
    }
L40:
    *x = *a + g * *e / (y - w);
    return;
L50:
    *x = *a;
    return;
} /* cub_ */

/* Subroutine */ void ins_(s, f, a, b, c, fa, fb, fc, j, y, z)
doublereal *s, *f, *a, *b, *c, *fa, *fb, *fc;
integer *j;
doublereal *y, *z;
{
    y[*j] = *s;
    z[*j] = *f;
    ++(*j);
    if (*f <= *fa) {
        *c = *b; *b = *a; *a = *s;
        *fc = *fb; *fb = *fa; *fa = *f;
        return;
    }
    if (*f <= *fb) {
        *c = *b; *b = *s;
        *fc = *fb; *fb = *f;
        return;
    }
    if (*f > *fc) {
        return;
    }
    *c = *s;
    *fc = *f;
} /* ins_ */