cg.c

DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.

        goto L530;
    }
/*< 520   B = A >*/
L520:
    b = a;
/*<       FB = FA >*/
    fb = fa;
/*<       DB = DA >*/
    db = da;
/*< 530   A = C >*/
L530:
    a = c__;
/*<       FA = F >*/
    fa = f;
/*<       DA = D >*/
    da = d__;
/*<       IF ( D .GT. A6*G ) GOTO 450 >*/
    if (d__ > a6 * g) {
        goto L450;
    }
/*<       GOTO 560 >*/
    goto L560;
/*< 540   IF ( D .LT. 0. ) GOTO 520 >*/
L540:
    if (d__ < (float)0.) {
        goto L520;
    }
/*<       GOTO 530 >*/
    goto L530;
/*< 550   C = .5*(A+B) >*/
L550:
    c__ = (a + b) * (float).5;
/*<       NB = NB + 1 >*/
    ++nb;
/*<       W = DABS(B-A) >*/
    w = (d__1 = b - a, abs(d__1));
/*<       GOTO 500 >*/
    goto L500;
/*< 560   E = 0. >*/
L560:
    *e = (float)0.;
/*<       DO 570 I = 1,N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<            IF ( DABS(H(I,3)) .GT. E ) E = DABS(H(I,3)) >*/
        if ((d__1 = h__[i__ + h_dim1 * 3], abs(d__1)) > *e) {
            *e = (d__2 = h__[i__ + h_dim1 * 3], abs(d__2));
        }
/*< 570        X(I) = H(I,2) >*/
/* L570: */
        x[i__] = h__[i__ + (h_dim1 << 1)];
    }
/*<       IT = IT + 1 >*/
    ++(*it);
/*<       IF ( E .LE. T ) GOTO 660 >*/
    if (*e <= *t) {
        goto L660;
    }
/*<       IF ( IT .GE. LIMIT ) GOTO 660 >*/
    if (*it >= *limit) {
        goto L660;
    }
/*<       F = FA >*/
    f = fa;
/*<       D = DA >*/
    d__ = da;
/*<       A = A7*A >*/
    a = a7 * a;
/*<       CALL PRE(H(1,2),H(1,3)) >*/
    (*pre)(&h__[(h_dim1 << 1) + 1], &h__[h_dim1 * 3 + 1], userdata);
/*<       R = 0. >*/
    r__ = (float)0.;
/*<       DO 580 I = 1,N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*< 580        R = R + H(I,2)*H(I,3) >*/
/* L580: */
        r__ += h__[i__ + (h_dim1 << 1)] * h__[i__ + h_dim1 * 3];
    }
/*<       IF ( R .LT. 0. ) GOTO 620 >*/
    if (r__ < (float)0.) {
        goto L620;
    }
/*<       S = R/G >*/
    s = r__ / g;
/*<       G = R >*/
    g = r__;
/*<       L = L + 1 >*/
    ++l;
/*<       IF ( L .GE. M ) GOTO 50 >*/
    if (l >= *m) {
        goto L50;
    }
/*<       D = 0. >*/
    d__ = (float)0.;
/*<       DO 590 I = 1,N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<            H(I,1) = -H(I,2) + S*H(I,1) >*/
        h__[i__ + h_dim1] = -h__[i__ + (h_dim1 << 1)] + s * h__[i__ + h_dim1];
/*< 590        D = D + H(I,1)*H(I,3) >*/
/* L590: */
        d__ += h__[i__ + h_dim1] * h__[i__ + h_dim1 * 3];
    }
/*<       GOTO 70 >*/
    goto L70;
/*< 600   IF ( D .LT. G ) GOTO 560 >*/
L600:
    if (d__ < g) {
        goto L560;
    }
/*<       WRITE(6,*) 'UNABLE TO OBTAIN DESCENT DIRECTION' >*/
    if(error_code)
      {
      *error_code = 1;
      return 0;
      }
    else
      {
      printf("UNABLE TO OBTAIN DESCENT DIRECTION\n");
      }
/*<       STOP >*/
    /*assert(0);*/
    return 0;
/*< 610   WRITE(6,*) 'THE FUNCTION DECREASES WITH NO MINIMUM' >*/
L610:
    if(error_code)
      {
      *error_code = 2;
      }
    else
      {
      printf("THE FUNCTION DECREASES WITH NO MINIMUM\n");
      }
/*<       STOP >*/
    /*assert(0);*/
    return 0;
/*< 620   WRITE(6,*) 'PRECONDITIONER NOT POSITIVE DEFINITE' >*/
L620:
    if(error_code)
      {
      *error_code = 3;
      return 0;
      }
    else
      {
      printf("PRECONDITIONER NOT POSITIVE DEFINITE\n");
      }
/*<       STOP >*/
    /*assert(0);*/
    return 0;
/*< 630   Q = Q*A3**25 >*/
L630:
/* Computing 25th power */
    d__1 = a3, d__2 = d__1, d__1 *= d__1, d__1 *= d__1, d__1 *= d__1, d__2 *= 
            d__1;
    q *= d__2 * (d__1 * d__1);
/*<       ND = 0 >*/
    nd = 0;
/*< 640   ND = ND + 1 >*/
L640:
    ++nd;
/*<       IF ( ND .GT. 25 ) GOTO 650 >*/
    if (nd > 25) {
        goto L650;
    }
/*<       Q = A3*Q >*/
    q = a3 * q;
/*<       P = FV(Q,X,H,N,VALUE) >*/
    p = fv_(&q, &x[1], &h__[h_offset], n, value, userdata);
/*<       CALL INS(Q,P,A,B,C,FA,FB,FC,J,Y,Z) >*/
    ins_(&q, &p, &a, &b, &c__, &fa, &fb, &fc, &j, y, z__);
/*<       IF ( P-F .GT. V*Q ) GOTO 640 >*/
    if (p - f > v * q) {
        goto L640;
    }
/*<       GOTO 135 >*/
    goto L135;
/*< 650   WRITE(6,*) 'UNABLE TO SATISFY ARMIJO CONDITION' >*/
L650:
    printf("UNABLE TO SATISFY ARMIJO CONDITION\n");
/*<       RETURN >*/
    return 0;
/*< 660   STEP = A >*/
L660:
    *step = a;
/*<       RETURN >*/
    return 0;
/*<       END >*/
} /* cg_ */

/*<       DOUBLE PRECISION FUNCTION FV(A,X,H,N,VALUE) >*/
doublereal fv_(doublereal *a, doublereal *x, doublereal *h__, integer *n, 
               double (*value)(double*,void*), void* userdata)
{
    /* System generated locals */
    integer h_dim1, h_offset, i__1;
    doublereal ret_val;

    /* Local variables */
    integer i__;

/*<       REAL*8 H(N,1),X(1),A,VALUE >*/
/*<       EXTERNAL VALUE >*/
/*<       DO 10 I = 1 , N >*/
    /* Parameter adjustments */
    --x;
    h_dim1 = *n;
    h_offset = 1 + h_dim1;
    h__ -= h_offset;

    /* Function Body */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*< 10         H(I,2) = X(I) + A*H(I,1) >*/
/* L10: */
        h__[i__ + (h_dim1 << 1)] = x[i__] + *a * h__[i__ + h_dim1];
    }
/*<       FV = VALUE(H(1,2)) >*/
    ret_val = (*value)(&h__[(h_dim1 << 1) + 1], userdata);
/*<       RETURN >*/
    return ret_val;
/*<       END >*/
} /* fv_ */

/*<       DOUBLE PRECISION FUNCTION FD(A,X,H,N,GRAD) >*/
doublereal fd_(doublereal *a, doublereal *x, doublereal *h__, integer *n, 
        void (*grad)(double*,double*,void*), void* userdata)
{
    /* System generated locals */
    integer h_dim1, h_offset, i__1;
    doublereal ret_val;

    /* Local variables */
    doublereal d__;
    integer i__;

/*<       REAL*8 H(N,1),X(1),A,D >*/
/*<       EXTERNAL GRAD >*/
/*<       DO 10 I = 1 , N >*/
    /* Parameter adjustments */
    --x;
    h_dim1 = *n;
    h_offset = 1 + h_dim1;
    h__ -= h_offset;

    /* Function Body */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*< 10         H(I,2) = X(I) + A*H(I,1) >*/
/* L10: */
        h__[i__ + (h_dim1 << 1)] = x[i__] + *a * h__[i__ + h_dim1];
    }
/*<       CALL GRAD(H(1,3),H(1,2)) >*/
    (*grad)(&h__[h_dim1 * 3 + 1], &h__[(h_dim1 << 1) + 1], userdata);
/*<       D = 0. >*/
    d__ = (float)0.;
/*<       DO 20 I = 1,N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*< 20         D = D + H(I,1)*H(I,3) >*/
/* L20: */
        d__ += h__[i__ + h_dim1] * h__[i__ + h_dim1 * 3];
    }
/*<       FD = D >*/
    ret_val = d__;
/*<       RETURN >*/
    return ret_val;
/*<       END >*/
} /* fd_ */

/*<       SUBROUTINE FVD(V,D,A,X,H,N,BOTH) >*/
/* Subroutine */ int fvd_(doublereal *v, doublereal *d__, doublereal *a, 
                          doublereal *x, doublereal *h__, integer *n,
                          void (*both)(double*,double*,double*,void*),
                          void* userdata)
{
    /* System generated locals */
    integer h_dim1, h_offset, i__1;

    /* Local variables */
    integer i__;

/*<       REAL*8 H(N,1),X(1),A,D,V >*/
/*<       EXTERNAL BOTH >*/
/*<       DO 10 I = 1 , N >*/
    /* Parameter adjustments */
    --x;
    h_dim1 = *n;
    h_offset = 1 + h_dim1;
    h__ -= h_offset;

    /* Function Body */
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*< 10         H(I,2) = X(I) + A*H(I,1) >*/
/* L10: */
        h__[i__ + (h_dim1 << 1)] = x[i__] + *a * h__[i__ + h_dim1];
    }
/*<       CALL BOTH(V,H(1,3),H(1,2)) >*/
    (*both)(v, &h__[h_dim1 * 3 + 1], &h__[(h_dim1 << 1) + 1], userdata);
/*<       D = 0. >*/
    *d__ = (float)0.;
/*<       DO 20 I = 1,N >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*< 20         D = D + H(I,1)*H(I,3) >*/
/* L20: */
        *d__ += h__[i__ + h_dim1] * h__[i__ + h_dim1 * 3];
    }
/*<       RETURN >*/
    return 0;
/*<       END >*/
} /* fvd_ */

/*<       SUBROUTINE CUB(X,A,B,C,D,E,F) >*/
/* Subroutine */ int cub_(doublereal *x, doublereal *a, doublereal *b, 
        doublereal *c__, doublereal *d__, doublereal *e, doublereal *f)
{
    /* System generated locals */
    doublereal d__1;

    /* Builtin functions */
    double sqrt(doublereal), d_sign(doublereal *, doublereal *);

    /* Local variables */
    doublereal g, v, w, y, z__;

/*<       REAL*8 A,B,C,D,E,F,G,V,W,X,Y,Z >*/
/*<       G = B - A >*/
    g = *b - *a;
/*<       IF ( G .EQ. 0. ) GOTO 50 >*/
    if (g == (float)0.) {
        goto L50;
    }
/*<       V = E + F - 3*(D-C)/G >*/
    v = *e + *f - (*d__ - *c__) * 3 / g;
/*<       W = V*V-E*F >*/
    w = v * v - *e * *f;
/*<       IF ( W .LT. 0. ) W = 0. >*/
    if (w < (float)0.) {
        w = (float)0.;
    }
/*<       W = DSIGN(DSQRT(W),G) >*/
    d__1 = sqrt(w);
    w = d_sign(&d__1, &g);
/*<       Y = E + V >*/
    y = *e + v;
/*<       Z = F + V >*/
    z__ = *f + v;
/*<       IF ( DSIGN(Y,G) .NE. Y ) GOTO 30 >*/
    if (d_sign(&y, &g) != y) {
        goto L30;
    }
/*<       IF ( DSIGN(Z,G) .NE. Z ) GOTO 20 >*/
    if (d_sign(&z__, &g) != z__) {
        goto L20;
    }
/*<       IF ( Z .EQ. 0. ) GOTO 20 >*/
    if (z__ == (float)0.) {
        goto L20;
    }
/*< 10    X = B - G*F/(Z+W) >*/
L10:
    *x = *b - g * *f / (z__ + w);
/*<       RETURN >*/
    return 0;
/*< 20    IF ( C .LT. D ) X = A >*/
L20:
    if (*c__ < *d__) {
        *x = *a;
    }
/*<       IF ( C .GE. D ) X = B >*/
    if (*c__ >= *d__) {
        *x = *b;
    }
/*<       RETURN >*/
    return 0;
/*< 30    IF ( DSIGN(Z,G) .NE. Z ) GOTO 40 >*/
L30:
    if (d_sign(&z__, &g) != z__) {
        goto L40;
    }
/*<       IF ( DABS(E) .GT. DABS(F) ) GOTO 10 >*/
    if (abs(*e) > abs(*f)) {
        goto L10;
    }
/*< 40    X = A + G*E/(Y-W) >*/
L40:
    *x = *a + g * *e / (y - w);
/*<       RETURN >*/
    return 0;
/*< 50    X = A >*/
L50:
    *x = *a;
/*<       RETURN >*/
    return 0;
/*<       END >*/
} /* cub_ */

/*<       SUBROUTINE INS(S,F,A,B,C,FA,FB,FC,J,Y,Z) >*/
/* Subroutine */ int ins_(doublereal *s, doublereal *f, doublereal *a, 
        doublereal *b, doublereal *c__, doublereal *fa, doublereal *fb, 
        doublereal *fc, integer *j, doublereal *y, doublereal *z__)
{
/*<       REAL*8 A,B,C,F,FA,FB,FC,S,Y(1),Z(1) >*/
/*<       INTEGER J >*/
/*<       J = J + 1 >*/
    /* Parameter adjustments */
    --z__;
    --y;

    /* Function Body */
    ++(*j);
/*<       Y(J) = S >*/
    y[*j] = *s;
/*<       Z(J) = F >*/
    z__[*j] = *f;
/*<       IF ( F .LE. FA ) GOTO 20 >*/
    if (*f <= *fa) {
        goto L20;
    }
/*<       IF ( F .LE. FB ) GOTO 10 >*/
    if (*f <= *fb) {
        goto L10;
    }
/*<       IF ( F .GT. FC ) RETURN >*/
    if (*f > *fc) {
        return 0;
    }
/*<       C = S >*/
    *c__ = *s;
/*<       FC = F >*/
    *fc = *f;
/*<       RETURN >*/
    return 0;
/*< 10    C = B >*/
L10:
    *c__ = *b;
/*<       B = S >*/
    *b = *s;
/*<       FC = FB >*/
    *fc = *fb;
/*<       FB = F >*/
    *fb = *f;
/*<       RETURN >*/
    return 0;
/*< 20    C = B >*/
L20:
    *c__ = *b;
/*<       B = A >*/
    *b = *a;
/*<       A = S >*/
    *a = *s;
/*<       FC = FB >*/
    *fc = *fb;
/*<       FB = FA >*/
    *fb = *fa;
/*<       FA = F >*/
    *fa = *f;
/*<       RETURN >*/
    return 0;
/*<       END >*/
} /* ins_ */

#ifdef __cplusplus
        }
#endif