lbfgsb.c

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

/*<       updatd = .true. >*/
    updatd = TRUE_;
/*<       iupdat = iupdat + 1 >*/
    ++iupdat;
/*     Update matrices WS and WY and form the middle matrix in B. */
/*<    >*/
    matupd_(n, m, &ws[ws_offset], &wy[wy_offset], &sy[sy_offset], &ss[
            ss_offset], &d__[1], &r__[1], &itail, &iupdat, &col, &head, &
            theta, &rr, &dr, &stp, &dtd);
/*     Form the upper half of the pds T = theta*SS + L*D^(-1)*L'; */
/*        Store T in the upper triangular of the array wt; */
/*        Cholesky factorize T to J*J' with */
/*           J' stored in the upper triangular of wt. */
/*<       call formt(m,wt,sy,ss,col,theta,info) >*/
    formt_(m, &wt[wt_offset], &sy[sy_offset], &ss[ss_offset], &col, &theta, &
            info);
/*<       if (info .ne. 0) then  >*/
    if (info != 0) {
/*          nonpositive definiteness in Cholesky factorization; */
/*          refresh the lbfgs memory and restart the iteration. */
/*<          if(iprint .ge. 1) write (6, 1007) >*/
/*
 1007 format (/,
     +' Nonpositive definiteness in Cholesky factorization in formt;',/,
     +'   refresh the lbfgs memory and restart the iteration.')
*/
        if (*iprint >= 1) {
            printf(" Nonpositive definiteness in Cholesky factorization in formt;\n"
                   "   refresh the lbfgs memory and restart the iteration.\n");
        }
/*<          info = 0 >*/
        info = 0;
/*<          col = 0 >*/
        col = 0;
/*<          head = 1 >*/
        head = 1;
/*<          theta = one >*/
        theta = 1.;
/*<          iupdat = 0 >*/
        iupdat = 0;
/*<          updatd = .false. >*/
        updatd = FALSE_;
/*<          goto 222 >*/
        goto L222;
/*<       endif >*/
    }
/*     Now the inverse of the middle matrix in B is */
/*       [  D^(1/2)      O ] [ -D^(1/2)  D^(-1/2)*L' ] */
/*       [ -L*D^(-1/2)   J ] [  0        J'          ] */
/*<  888  continue >*/
L888:
/* -------------------- the end of the loop ----------------------------- */
/*<       goto 222 >*/
    goto L222;
/*<  999  continue >*/
L999:
/*<       call timer(time2) >*/
    timer_(&time2);
/*<       time = time2 - time1 >*/
    time = time2 - time1;
/*<    >*/
    prn3lb_(n, &x[1], f, task, iprint, &info, &itfile, &iter, &nfgv, &nintol, 
            &nskip, &nact, &sbgnrm, &time, &nint, word, &iback, &stp, &xstep, 
            &k, &cachyt, &sbtime, &lnscht, (ftnlen)60, (ftnlen)3);
/*<  1000 continue >*/
L1000:
/*     Save local variables. */
/*<       lsave(1)  = prjctd >*/
    lsave[1] = prjctd;
/*<       lsave(2)  = cnstnd >*/
    lsave[2] = cnstnd;
/*<       lsave(3)  = boxed >*/
    lsave[3] = boxed;
/*<       lsave(4)  = updatd >*/
    lsave[4] = updatd;
/*<       isave(1)  = nintol  >*/
    isave[1] = nintol;
/*<       isave(3)  = itfile  >*/
    isave[3] = itfile;
/*<       isave(4)  = iback  >*/
    isave[4] = iback;
/*<       isave(5)  = nskip  >*/
    isave[5] = nskip;
/*<       isave(6)  = head  >*/
    isave[6] = head;
/*<       isave(7)  = col  >*/
    isave[7] = col;
/*<       isave(8)  = itail  >*/
    isave[8] = itail;
/*<       isave(9)  = iter  >*/
    isave[9] = iter;
/*<       isave(10) = iupdat  >*/
    isave[10] = iupdat;
/*<       isave(12) = nint  >*/
    isave[12] = nint;
/*<       isave(13) = nfgv  >*/
    isave[13] = nfgv;
/*<       isave(14) = info  >*/
    isave[14] = info;
/*<       isave(15) = ifun  >*/
    isave[15] = ifun;
/*<       isave(16) = iword  >*/
    isave[16] = iword;
/*<       isave(17) = nfree  >*/
    isave[17] = nfree;
/*<       isave(18) = nact  >*/
    isave[18] = nact;
/*<       isave(19) = ileave  >*/
    isave[19] = ileave;
/*<       isave(20) = nenter  >*/
    isave[20] = nenter;
/*<       dsave(1)  = theta  >*/
    dsave[1] = theta;
/*<       dsave(2)  = fold  >*/
    dsave[2] = fold;
/*<       dsave(3)  = tol  >*/
    dsave[3] = tol;
/*<       dsave(4)  = dnorm  >*/
    dsave[4] = dnorm;
/*<       dsave(5)  = epsmch  >*/
    dsave[5] = epsmch;
/*<       dsave(6)  = cpu1  >*/
    dsave[6] = cpu1;
/*<       dsave(7)  = cachyt  >*/
    dsave[7] = cachyt;
/*<       dsave(8)  = sbtime  >*/
    dsave[8] = sbtime;
/*<       dsave(9)  = lnscht  >*/
    dsave[9] = lnscht;
/*<       dsave(10) = time1  >*/
    dsave[10] = time1;
/*<       dsave(11) = gd  >*/
    dsave[11] = gd;
/*<       dsave(12) = stpmx  >*/
    dsave[12] = stpmx;
/*<       dsave(13) = sbgnrm >*/
    dsave[13] = sbgnrm;
/*<       dsave(14) = stp >*/
    dsave[14] = stp;
/*<       dsave(15) = gdold >*/
    dsave[15] = gdold;
/*<       dsave(16) = dtd   >*/
    dsave[16] = dtd;
/*<  1001 format (//,'ITERATION ',i5) >*/
/*<  1 >*/
/*<  1 >*/
/*<  1004 format ('  ys=',1p,e10.3,'  -gs=',1p,e10.3,' BFGS update SKIPPED') >*/
/*<  1 >*/
/*<  1 >*/
/*<  1 >*/
/*<  1 >*/
/*<       return    >*/
    return 0;
/*<       end >*/
} /* mainlb_ */

/* ======================= The end of mainlb ============================= */
/*<    >*/
/* Subroutine */ int active_(integer *n, doublereal *l, doublereal *u, 
        integer *nbd, doublereal *x, integer *iwhere, integer *iprint, 
        logical *prjctd, logical *cnstnd, logical *boxed)
{
    /* System generated locals */
    integer i__1;

    /* Local variables */
    integer i__, nbdd;

/*<       logical          prjctd, cnstnd, boxed >*/
/*<       integer          n, iprint, nbd(n), iwhere(n) >*/
/*<       double precision x(n), l(n), u(n) >*/
/*     ************ */

/*     Subroutine active */

/*     This subroutine initializes iwhere and projects the initial x to */
/*       the feasible set if necessary. */

/*     iwhere is an integer array of dimension n. */
/*       On entry iwhere is unspecified. */
/*       On exit iwhere(i)=-1  if x(i) has no bounds */
/*                         3   if l(i)=u(i) */
/*                         0   otherwise. */
/*       In cauchy, iwhere is given finer gradations. */


/*                           *  *  * */

/*     NEOS, November 1994. (Latest revision June 1996.) */
/*     Optimization Technology Center. */
/*     Argonne National Laboratory and Northwestern University. */
/*     Written by */
/*                        Ciyou Zhu */
/*     in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal. */


/*     ************ */
/*<       integer          nbdd,i >*/
/*<       double precision zero >*/
/*<       parameter        (zero=0.0d0) >*/
/*     Initialize nbdd, prjctd, cnstnd and boxed. */
/*<       nbdd = 0 >*/
    /* Parameter adjustments */
    --iwhere;
    --x;
    --nbd;
    --u;
    --l;

    /* Function Body */
    nbdd = 0;
/*<       prjctd = .false. >*/
    *prjctd = FALSE_;
/*<       cnstnd = .false. >*/
    *cnstnd = FALSE_;
/*<       boxed = .true. >*/
    *boxed = TRUE_;
/*     Project the initial x to the easible set if necessary. */
/*<       do 10 i = 1, n >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          if (nbd(i) .gt. 0) then >*/
        if (nbd[i__] > 0) {
/*<             if (nbd(i) .le. 2 .and. x(i) .le. l(i)) then >*/
            if (nbd[i__] <= 2 && x[i__] <= l[i__]) {
/*<            if (x(i) .lt. l(i)) then >*/
                if (x[i__] < l[i__]) {
/*<                   prjctd = .true. >*/
                    *prjctd = TRUE_;
/*<               x(i) = l(i) >*/
                    x[i__] = l[i__];
/*<                endif >*/
                }
/*<                nbdd = nbdd + 1 >*/
                ++nbdd;
/*<             else if (nbd(i) .ge. 2 .and. x(i) .ge. u(i)) then >*/
            } else if (nbd[i__] >= 2 && x[i__] >= u[i__]) {
/*<            if (x(i) .gt. u(i)) then >*/
                if (x[i__] > u[i__]) {
/*<                   prjctd = .true. >*/
                    *prjctd = TRUE_;
/*<               x(i) = u(i) >*/
                    x[i__] = u[i__];
/*<                endif >*/
                }
/*<                nbdd = nbdd + 1 >*/
                ++nbdd;
/*<             endif >*/
            }
/*<          endif >*/
        }
/*<   10  continue >*/
/* L10: */
    }
/*     Initialize iwhere and assign values to cnstnd and boxed. */
/*<       do 20 i = 1, n >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          if (nbd(i) .ne. 2) boxed = .false. >*/
        if (nbd[i__] != 2) {
            *boxed = FALSE_;
        }
/*<          if (nbd(i) .eq. 0) then >*/
        if (nbd[i__] == 0) {
/*                                this variable is always free */
/*<         iwhere(i) = -1 >*/
            iwhere[i__] = -1;
/*           otherwise set x(i)=mid(x(i), u(i), l(i)). */
/*<          else >*/
        } else {
/*<         cnstnd = .true. >*/
            *cnstnd = TRUE_;
/*<             if (nbd(i) .eq. 2 .and. u(i) - l(i) .le. zero) then >*/
            if (nbd[i__] == 2 && u[i__] - l[i__] <= 0.) {
/*                   this variable is always fixed */
/*<            iwhere(i) = 3 >*/
                iwhere[i__] = 3;
/*<             else  >*/
            } else {
/*<                iwhere(i) = 0 >*/
                iwhere[i__] = 0;
/*<             endif >*/
            }
/*<          endif >*/
        }
/*<   20  continue >*/
/* L20: */
    }
/*<       if (iprint .ge. 0) then >*/
    if (*iprint >= 0) {
/*<    >*/
        if (*prjctd) {
            printf("The initial X is infeasible.  Restart with its projection.\n");
        }
/*<    >*/
        if (! (*cnstnd)) {
            printf("This problem is unconstrained.\n");
        }
/*<       endif >*/
    }
/*<       if (iprint .gt. 0) write (6,1001) nbdd >*/
    if (*iprint > 0) {
        printf("At X0 %9ld variables are exactly at the bounds\n", nbdd);
    }
/*<  1001 format (/,'At X0 ',i9,' variables are exactly at the bounds')  >*/
/*<       return >*/
    return 0;
/*<       end >*/
} /* active_ */

/* ======================= The end of active ============================= */
/*<       subroutine bmv(m, sy, wt, col, v, p, info) >*/
/* Subroutine */ int bmv_(integer *m, doublereal *sy, doublereal *wt, integer 
        *col, doublereal *v, doublereal *p, integer *info)
{
    /* System generated locals */
    integer sy_dim1, sy_offset, wt_dim1, wt_offset, i__1, i__2;

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

    /* Local variables */
    integer i__, k, i2;
    doublereal sum;
    extern /* Subroutine */ int dtrsl_(doublereal *, integer *, integer *, 
            doublereal *, integer *, integer *);

/*<       integer m, col, info >*/
/*<       double precision sy(m, m), wt(m, m), v(2*col), p(2*col) >*/
/*     ************ */

/*     Subroutine bmv */

/*     This subroutine computes the product of the 2m x 2m middle matrix */
/*      in the compact L-BFGS formula of B and a 2m vector v; */
/*      it returns the product in p. */

/*     m is an integer variable. */
/*       On entry m is the maximum number of variable metric corrections */
/*         used to define the limited memory matrix. */
/*       On exit m is unchanged. */

/*     sy is a double precision array of dimension m x m. */
/*       On entry sy specifies the matrix S'Y. */
/*       On exit sy is unchanged. */

/*     wt is a double precision array of dimension m x m. */
/*       On entry wt specifies the u