lbfgsb.c

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

        }
/*<          info   = 0 >*/
        info = 0;
/*<          col    = 0 >*/
        col = 0;
/*<          head   = 1 >*/
        head = 1;
/*<          theta  = one >*/
        theta = 1.;
/*<          iupdat = 0 >*/
        iupdat = 0;
/*<          updatd = .false. >*/
        updatd = FALSE_;
/*<          call timer(cpu2)  >*/
        timer_(&cpu2);
/*<          cachyt = cachyt + cpu2 - cpu1 >*/
        cachyt = cachyt + cpu2 - cpu1;
/*<          goto 222 >*/
        goto L222;
/*<       endif >*/
    }
/*<       call timer(cpu2)  >*/
    timer_(&cpu2);
/*<       cachyt = cachyt + cpu2 - cpu1 >*/
    cachyt = cachyt + cpu2 - cpu1;
/*<       nintol = nintol + nint >*/
    nintol += nint;
/*     Count the entering and leaving variables for iter > 0; */
/*     find the index set of free and active variables at the GCP. */
/*<    >*/
    freev_(n, &nfree, &index[1], &nenter, &ileave, &indx2[1], &iwhere[1], &
            wrk, &updatd, &cnstnd, iprint, &iter);
/*<       nact = n - nfree >*/
    nact = *n - nfree;
/*<  333  continue >*/
L333:
/*     If there are no free variables or B=theta*I, then */
/*                                        skip the subspace minimization. */
/*<       if (nfree .eq. 0 .or. col .eq. 0) goto 555 >*/
    if (nfree == 0 || col == 0) {
        goto L555;
    }
/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */

/*     Subspace minimization. */

/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
/*<       call timer(cpu1)  >*/
    timer_(&cpu1);
/*     Form  the LEL^T factorization of the indefinite */
/*       matrix    K = [-D -Y'ZZ'Y/theta     L_a'-R_z'  ] */
/*                     [L_a -R_z           theta*S'AA'S ] */
/*       where     E = [-I  0] */
/*                     [ 0  I] */
/*<    >*/
    if (wrk) {
        formk_(n, &nfree, &index[1], &nenter, &ileave, &indx2[1], &iupdat, &
                updatd, &wn[wn_offset], &snd[snd_offset], m, &ws[ws_offset], &
                wy[wy_offset], &sy[sy_offset], &theta, &col, &head, &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, 1006) >*/
/*
 1006 format (/, 
     +' Nonpositive definiteness in Cholesky factorization in formk;',/,
     +'   refresh the lbfgs memory and restart the iteration.')
 */
        if (*iprint >= 1) {
            printf(" Nonpositive definiteness in Cholesky factorization in formk;\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_;
/*<          call timer(cpu2)  >*/
        timer_(&cpu2);
/*<          sbtime = sbtime + cpu2 - cpu1  >*/
        sbtime = sbtime + cpu2 - cpu1;
/*<          goto 222 >*/
        goto L222;
/*<       endif  >*/
    }
/*        compute r=-Z'B(xcp-xk)-Z'g (using wa(2m+1)=W'(xcp-x) */
/*                                                   from 'cauchy'). */
/*<    >*/
    cmprlb_(n, m, &x[1], &g[1], &ws[ws_offset], &wy[wy_offset], &sy[sy_offset]
            , &wt[wt_offset], &z__[1], &r__[1], &wa[1], &index[1], &theta, &
            col, &head, &nfree, &cnstnd, &info);
/*<       if (info .ne. 0) goto 444 >*/
    if (info != 0) {
        goto L444;
    }
/*       call the direct method. */
/*<    >*/
    subsm_(n, m, &nfree, &index[1], &l[1], &u[1], &nbd[1], &z__[1], &r__[1], &
            ws[ws_offset], &wy[wy_offset], &theta, &col, &head, &iword, &wa[1]
            , &wn[wn_offset], iprint, &info);
/*<  444  continue >*/
L444:
/*<       if (info .ne. 0) then  >*/
    if (info != 0) {
/*          singular triangular system detected; */
/*          refresh the lbfgs memory and restart the iteration. */
/*<          if(iprint .ge. 1) write (6, 1005) >*/
        if (*iprint >= 1) {
            printf(" Singular triangular system detected;\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_;
/*<          call timer(cpu2)  >*/
        timer_(&cpu2);
/*<          sbtime = sbtime + cpu2 - cpu1  >*/
        sbtime = sbtime + cpu2 - cpu1;
/*<          goto 222 >*/
        goto L222;
/*<       endif >*/
    }
/*<       call timer(cpu2)  >*/
    timer_(&cpu2);
/*<       sbtime = sbtime + cpu2 - cpu1  >*/
    sbtime = sbtime + cpu2 - cpu1;
/*<  555  continue >*/
L555:
/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */

/*     Line search and optimality tests. */

/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
/*     Generate the search direction d:=z-x. */
/*<       do 40 i = 1, n >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<      d(i) = z(i) - x(i) >*/
        d__[i__] = z__[i__] - x[i__];
/*<   40  continue >*/
/* L40: */
    }
/*<       call timer(cpu1)  >*/
    timer_(&cpu1);
/*<  666  continue >*/
L666:
/*<    >*/
    lnsrlb_(n, &l[1], &u[1], &nbd[1], &x[1], f, &fold, &gd, &gdold, &g[1], &
            d__[1], &r__[1], &t[1], &z__[1], &stp, &dnorm, &dtd, &xstep, &
            stpmx, &iter, &ifun, &iback, &nfgv, &info, task, &boxed, &cnstnd, 
            csave, &isave[22], &dsave[17], (ftnlen)60, (ftnlen)60);
/*<       if (info .ne. 0 .or. iback .ge. 20) then >*/
    if (info != 0 || iback >= 20) {
/*          restore the previous iterate. */
/*<          call dcopy(n,t,1,x,1) >*/
        dcopy_(n, &t[1], &c__1, &x[1], &c__1);
/*<          call dcopy(n,r,1,g,1) >*/
        dcopy_(n, &r__[1], &c__1, &g[1], &c__1);
/*<          f = fold >*/
        *f = fold;
/*<          if (col .eq. 0) then >*/
        if (col == 0) {
/*             abnormal termination. */
/*<             if (info .eq. 0) then >*/
            if (info == 0) {
/*<                info = -9 >*/
                info = -9;
/*                restore the actual number of f and g evaluations etc. */
/*<                nfgv = nfgv - 1 >*/
                --nfgv;
/*<                ifun = ifun - 1 >*/
                --ifun;
/*<                iback = iback - 1 >*/
                --iback;
/*<             endif >*/
            }
/*<             task = 'ABNORMAL_TERMINATION_IN_LNSRCH' >*/
            s_copy(task, "ABNORMAL_TERMINATION_IN_LNSRCH", (ftnlen)60, (
                    ftnlen)30);
/*<             iter = iter + 1 >*/
            ++iter;
/*<             goto 999 >*/
            goto L999;
/*<          else >*/
        } else {
/*             refresh the lbfgs memory and restart the iteration. */
/*<             if(iprint .ge. 1) write (6, 1008) >*/
/*
 1008 format (/,
     +' Bad direction in the line search;',/,
     +'   refresh the lbfgs memory and restart the iteration.')
*/
            if (*iprint >= 1) {
                printf(" Bad direction in the line search;\n"
                       "   refresh the lbfgs memory and restart the iteration.\n");
            }
/*<             if (info .eq. 0) nfgv = nfgv - 1 >*/
            if (info == 0) {
                --nfgv;
            }
/*<             info   = 0 >*/
            info = 0;
/*<             col    = 0 >*/
            col = 0;
/*<             head   = 1 >*/
            head = 1;
/*<             theta  = one >*/
            theta = 1.;
/*<             iupdat = 0 >*/
            iupdat = 0;
/*<             updatd = .false. >*/
            updatd = FALSE_;
/*<             task   = 'RESTART_FROM_LNSRCH' >*/
            s_copy(task, "RESTART_FROM_LNSRCH", (ftnlen)60, (ftnlen)19);
/*<             call timer(cpu2) >*/
            timer_(&cpu2);
/*<             lnscht = lnscht + cpu2 - cpu1 >*/
            lnscht = lnscht + cpu2 - cpu1;
/*<             goto 222 >*/
            goto L222;
/*<          endif >*/
        }
/*<       else if (task(1:5) .eq. 'FG_LN') then >*/
    } else if (s_cmp(task, "FG_LN", (ftnlen)5, (ftnlen)5) == 0) {
/*          return to the driver for calculating f and g; reenter at 666. */
/*<      goto 1000 >*/
        goto L1000;
/*<       else  >*/
    } else {
/*          calculate and print out the quantities related to the new X. */
/*<          call timer(cpu2)  >*/
        timer_(&cpu2);
/*<          lnscht = lnscht + cpu2 - cpu1 >*/
        lnscht = lnscht + cpu2 - cpu1;
/*<          iter = iter + 1 >*/
        ++iter;
/*        Compute the infinity norm of the projected (-)gradient. */
/*<          call projgr(n,l,u,nbd,x,g,sbgnrm) >*/
        projgr_(n, &l[1], &u[1], &nbd[1], &x[1], &g[1], &sbgnrm);
/*        Print iteration information. */
/*<    >*/
        prn2lb_(n, &x[1], f, &g[1], iprint, &itfile, &iter, &nfgv, &nact, &
                sbgnrm, &nint, word, &iword, &iback, &stp, &xstep, (ftnlen)3);
/*<          goto 1000 >*/
        goto L1000;
/*<       endif >*/
    }
/*<  777  continue >*/
L777:
/*     Test for termination. */
/*<       if (sbgnrm .le. pgtol) then >*/
    if (sbgnrm <= *pgtol) {
/*                                terminate the algorithm. */
/*<          task = 'CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL' >*/
        s_copy(task, "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL", (
                ftnlen)60, (ftnlen)48);
/*<          goto 999 >*/
        goto L999;
/*<       endif  >*/
    }
/*<       ddum = max(abs(fold), abs(f), one) >*/
/* Computing MAX */
    d__1 = abs(fold), d__2 = abs(*f), d__1 = max(d__1,d__2);
    ddum = max(d__1,1.);
/*<       if ((fold - f) .le. tol*ddum) then >*/
    if (fold - *f <= tol * ddum) {
/*                                        terminate the algorithm. */
/*<          task = 'CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH' >*/
        s_copy(task, "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH", (
                ftnlen)60, (ftnlen)47);
/*<          if (iback .ge. 10) info = -5 >*/
        if (iback >= 10) {
            info = -5;
        }
/*           i.e., to issue a warning if iback>10 in the line search. */
/*<          goto 999 >*/
        goto L999;
/*<       endif  >*/
    }
/*     Compute d=newx-oldx, r=newg-oldg, rr=y'y and dr=y's. */
/*<       do 42 i = 1, n >*/
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<          r(i) = g(i) - r(i) >*/
        r__[i__] = g[i__] - r__[i__];
/*<   42  continue >*/
/* L42: */
    }
/*<       rr = ddot(n,r,1,r,1) >*/
    rr = ddot_(n, &r__[1], &c__1, &r__[1], &c__1);
/*<       if (stp .eq. one) then   >*/
    if (stp == 1.) {
/*<          dr = gd - gdold >*/
        dr = gd - gdold;
/*<          ddum = -gdold >*/
        ddum = -gdold;
/*<       else >*/
    } else {
/*<          dr = (gd - gdold)*stp >*/
        dr = (gd - gdold) * stp;
/*<      call dscal(n,stp,d,1) >*/
        dscal_(n, &stp, &d__[1], &c__1);
/*<          ddum = -gdold*stp >*/
        ddum = -gdold * stp;
/*<       endif >*/
    }
/*<       if (dr .le. epsmch*ddum) then >*/
    if (dr <= epsmch * ddum) {
/*                            skip the L-BFGS update. */
/*<          nskip = nskip + 1 >*/
        ++nskip;
/*<          updatd = .false. >*/
        updatd = FALSE_;
/*<          if (iprint .ge. 1) write (6,1004) dr, ddum >*/
/*
 1004 format ('  ys=',1p,e10.3,'  -gs=',1p,e10.3,' BFGS update SKIPPED')
 */
        if (*iprint >= 1) {
            printf("  ys=%10.3g  -gs=%10.3g BFGS update SKIPPED\n",
                   dr, ddum);
        }
/*<          goto 888 >*/
        goto L888;
/*<       endif  >*/
    }
/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */

/*     Update the L-BFGS matrix. */

/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */