lbfgsb.c

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

/*<    >*/
/*<       double precision one,zero >*/
/*<       parameter        (one=1.0d0,zero=0.0d0) >*/
/*<       if (task .eq. 'START') then >*/
    /* Parameter adjustments */
    --indx2;
    --iwhere;
    --index;
    --t;
    --d__;
    --r__;
    --z__;
    --g;
    --nbd;
    --u;
    --l;
    --x;
    --ygo;
    --yg;
    --sgo;
    --sg;
    --wa;
    snd_dim1 = 2 * *m;
    snd_offset = 1 + snd_dim1;
    snd -= snd_offset;
    wn_dim1 = 2 * *m;
    wn_offset = 1 + wn_dim1;
    wn -= wn_offset;
    wt_dim1 = *m;
    wt_offset = 1 + wt_dim1;
    wt -= wt_offset;
    yy_dim1 = *m;
    yy_offset = 1 + yy_dim1;
    yy -= yy_offset;
    ss_dim1 = *m;
    ss_offset = 1 + ss_dim1;
    ss -= ss_offset;
    sy_dim1 = *m;
    sy_offset = 1 + sy_dim1;
    sy -= sy_offset;
    wy_dim1 = *n;
    wy_offset = 1 + wy_dim1;
    wy -= wy_offset;
    ws_dim1 = *n;
    ws_offset = 1 + ws_dim1;
    ws -= ws_offset;
    --lsave;
    --isave;
    --dsave;

    /* Function Body */
    if (s_cmp(task, "START", (ftnlen)60, (ftnlen)5) == 0) {
/*<          call timer(time1) >*/
        timer_(&time1);
/*        Generate the current machine precision. */
/*<          epsmch = dpmeps() >*/
        epsmch = dpmeps_();
/*        Initialize counters and scalars when task='START'. */
/*           for the limited memory BFGS matrices: */
/*<          col    = 0 >*/
        col = 0;
/*<          head   = 1 >*/
        head = 1;
/*<          theta  = one >*/
        theta = 1.;
/*<          iupdat = 0 >*/
        iupdat = 0;
/*<          updatd = .false. >*/
        updatd = FALSE_;
/*           for operation counts: */
/*<          iter   = 0 >*/
        iter = 0;
/*<          nfgv   = 0 >*/
        nfgv = 0;
/*<          nint   = 0 >*/
        nint = 0;
/*<          nintol = 0 >*/
        nintol = 0;
/*<          nskip  = 0 >*/
        nskip = 0;
/*<          nfree  = n >*/
        nfree = *n;
/*           for stopping tolerance: */
/*<          tol = factr*epsmch >*/
        tol = *factr * epsmch;
/*           for measuring running time: */
/*<          cachyt = 0 >*/
        cachyt = 0.;
/*<          sbtime = 0 >*/
        sbtime = 0.;
/*<          lnscht = 0 >*/
        lnscht = 0.;
/*           'word' records the status of subspace solutions. */
/*<          word = '---' >*/
        s_copy(word, "---", (ftnlen)3, (ftnlen)3);
/*           'info' records the termination information. */
/*<          info = 0 >*/
        info = 0;
/*<          if (iprint .ge. 1) then >*/
        if (*iprint >= 1) {
/*                                open a summary file 'iterate.dat' */
/*<             open (8, file = 'iterate.dat', status = 'unknown') >*/
/*
            o__1.oerr = 0;
            o__1.ounit = 8;
            o__1.ofnmlen = 11;
            o__1.ofnm = "iterate.dat";
            o__1.orl = 0;
            o__1.osta = "unknown";
            o__1.oacc = 0;
            o__1.ofm = 0;
            o__1.oblnk = 0;
            f_open(&o__1);
*/
/*<             itfile = 8 >*/
            itfile = 8;
/*<          endif             >*/
        }
/*        Check the input arguments for errors. */
/*<      call errclb(n,m,factr,l,u,nbd,task,info,k) >*/
        errclb_(n, m, factr, &l[1], &u[1], &nbd[1], task, &info, &k, (ftnlen)
                60);
/*<          if (task(1:5) .eq. 'ERROR') then >*/
        if (s_cmp(task, "ERROR", (ftnlen)5, (ftnlen)5) == 0) {
/*<    >*/
            prn3lb_(n, &x[1], f, task, iprint, &info, &itfile, &iter, &nfgv, &
                    nintol, &nskip, &nact, &sbgnrm, &c_b9, &nint, word, &
                    iback, &stp, &xstep, &k, &cachyt, &sbtime, &lnscht, (
                    ftnlen)60, (ftnlen)3);
/*<             return >*/
            return 0;
/*<          endif >*/
        }
/*<          call prn1lb(n,m,l,u,x,iprint,itfile,epsmch) >*/
        prn1lb_(n, m, &l[1], &u[1], &x[1], iprint, &itfile, &epsmch);
/*        Initialize iwhere & project x onto the feasible set. */
/*<          call active(n,l,u,nbd,x,iwhere,iprint,prjctd,cnstnd,boxed)  >*/
        active_(n, &l[1], &u[1], &nbd[1], &x[1], &iwhere[1], iprint, &prjctd, 
                &cnstnd, &boxed);
/*        The end of the initialization. */
/*<       else >*/
    } else {
/*          restore local variables. */
/*<          prjctd = lsave(1) >*/
        prjctd = lsave[1];
/*<          cnstnd = lsave(2) >*/
        cnstnd = lsave[2];
/*<          boxed  = lsave(3) >*/
        boxed = lsave[3];
/*<          updatd = lsave(4) >*/
        updatd = lsave[4];
/*<          nintol = isave(1) >*/
        nintol = isave[1];
/*<          itfile = isave(3) >*/
        itfile = isave[3];
/*<          iback  = isave(4) >*/
        iback = isave[4];
/*<          nskip  = isave(5) >*/
        nskip = isave[5];
/*<          head   = isave(6) >*/
        head = isave[6];
/*<          col    = isave(7) >*/
        col = isave[7];
/*<          itail  = isave(8) >*/
        itail = isave[8];
/*<          iter   = isave(9) >*/
        iter = isave[9];
/*<          iupdat = isave(10) >*/
        iupdat = isave[10];
/*<          nint   = isave(12) >*/
        nint = isave[12];
/*<          nfgv   = isave(13) >*/
        nfgv = isave[13];
/*<          info   = isave(14) >*/
        info = isave[14];
/*<          ifun   = isave(15) >*/
        ifun = isave[15];
/*<          iword  = isave(16) >*/
        iword = isave[16];
/*<          nfree  = isave(17) >*/
        nfree = isave[17];
/*<          nact   = isave(18) >*/
        nact = isave[18];
/*<          ileave = isave(19) >*/
        ileave = isave[19];
/*<          nenter = isave(20) >*/
        nenter = isave[20];
/*<          theta  = dsave(1) >*/
        theta = dsave[1];
/*<          fold   = dsave(2) >*/
        fold = dsave[2];
/*<          tol    = dsave(3) >*/
        tol = dsave[3];
/*<          dnorm  = dsave(4) >*/
        dnorm = dsave[4];
/*<          epsmch = dsave(5) >*/
        epsmch = dsave[5];
/*<          cpu1   = dsave(6) >*/
        cpu1 = dsave[6];
/*<          cachyt = dsave(7) >*/
        cachyt = dsave[7];
/*<          sbtime = dsave(8) >*/
        sbtime = dsave[8];
/*<          lnscht = dsave(9) >*/
        lnscht = dsave[9];
/*<          time1  = dsave(10) >*/
        time1 = dsave[10];
/*<          gd     = dsave(11) >*/
        gd = dsave[11];
/*<          stpmx  = dsave(12) >*/
        stpmx = dsave[12];
/*<          sbgnrm = dsave(13) >*/
        sbgnrm = dsave[13];
/*<          stp    = dsave(14) >*/
        stp = dsave[14];
/*<          gdold  = dsave(15) >*/
        gdold = dsave[15];
/*<          dtd    = dsave(16) >*/
        dtd = dsave[16];
/*        After returning from the driver go to the point where execution */
/*        is to resume. */
/*<          if (task(1:5) .eq. 'FG_LN') goto 666 >*/
        if (s_cmp(task, "FG_LN", (ftnlen)5, (ftnlen)5) == 0) {
            goto L666;
        }
/*<          if (task(1:5) .eq. 'NEW_X') goto 777 >*/
        if (s_cmp(task, "NEW_X", (ftnlen)5, (ftnlen)5) == 0) {
            goto L777;
        }
/*<          if (task(1:5) .eq. 'FG_ST') goto 111 >*/
        if (s_cmp(task, "FG_ST", (ftnlen)5, (ftnlen)5) == 0) {
            goto L111;
        }
/*<          if (task(1:4) .eq. 'STOP') then >*/
        if (s_cmp(task, "STOP", (ftnlen)4, (ftnlen)4) == 0) {
/*<             if (task(7:9) .eq. 'CPU') then >*/
            if (s_cmp(task + 6, "CPU", (ftnlen)3, (ftnlen)3) == 0) {
/*                                          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;
/*<             endif >*/
            }
/*<             goto 999 >*/
            goto L999;
/*<          endif >*/
        }
/*<       endif  >*/
    }
/*     Compute f0 and g0. */
/*<       task = 'FG_START'  >*/
    s_copy(task, "FG_START", (ftnlen)60, (ftnlen)8);
/*          return to the driver to calculate f and g; reenter at 111. */
/*<       goto 1000 >*/
    goto L1000;
/*<  111  continue >*/
L111:
/*<       nfgv = 1 >*/
    nfgv = 1;
/*     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);
/*<       if (iprint .ge. 1) then >*/
    if (*iprint >= 1) {
/*<          write (6,1002) iter,f,sbgnrm >*/
/*
 1002 format
     +  (/,'At iterate',i5,4x,'f= ',1p,d12.5,4x,'|proj g|= ',1p,d12.5)
*/
        printf("At iterate %5ld    f= %12.5g    |proj g|= %12.5g\n",
               iter, *f, sbgnrm);
/*<          write (itfile,1003) iter,nfgv,sbgnrm,f >*/
/*
 1003 format (2(1x,i4),5x,'-',5x,'-',3x,'-',5x,'-',5x,'-',8x,'-',3x,
     +        1p,2(1x,d10.3))
*/
#ifdef LBFGSB_ENABLE_ITERATE_FILE
        fprintf(0,
                " %4ld %4ld     -     -   -     -     -        -    %10.3g %10.3g\n",
                iter, nfgv, sbgnrm, *f);
#endif
/*<       endif >*/
    }
/*<       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  >*/
    }
/* ----------------- the beginning of the loop -------------------------- */
/*<  222  continue >*/
L222:
/*<       if (iprint .ge. 99) write (6,1001) iter + 1 >*/
/*
 1001 format (//,'ITERATION ',i5)
*/
    if (*iprint >= 99) {
        i__1 = iter + 1;
        printf("ITERATION %5ld\n", i__1);
    }
/*<       iword = -1 >*/
    iword = -1;

/*<       if (.not. cnstnd .and. col .gt. 0) then  >*/
    if (! cnstnd && col > 0) {
/*                                            skip the search for GCP. */
/*<          call dcopy(n,x,1,z,1) >*/
        dcopy_(n, &x[1], &c__1, &z__[1], &c__1);
/*<      wrk = updatd >*/
        wrk = updatd;
/*<          nint = 0 >*/
        nint = 0;
/*<          goto 333 >*/
        goto L333;
/*<       endif >*/
    }
/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */

/*     Compute the Generalized Cauchy Point (GCP). */

/* ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
/*<       call timer(cpu1)  >*/
    timer_(&cpu1);
/*<    >*/
    cauchy_(n, &x[1], &l[1], &u[1], &nbd[1], &g[1], &indx2[1], &iwhere[1], &t[
            1], &d__[1], &z__[1], m, &wy[wy_offset], &ws[ws_offset], &sy[
            sy_offset], &wt[wt_offset], &theta, &col, &head, &wa[1], &wa[(*m 
            << 1) + 1], &wa[(*m << 2) + 1], &wa[*m * 6 + 1], &nint, &sg[1], &
            yg[1], iprint, &sbgnrm, &info, &epsmch);
/*<       if (info .ne. 0) then  >*/
    if (info != 0) {
/*         singular triangular system detected; refresh the lbfgs memory. */
/*<          if(iprint .ge. 1) write (6, 1005) >*/
/*
 1005 format (/, 
     +' Singular triangular system detected;',/,
     +'   refresh the lbfgs memory and restart the iteration.')
*/
        if (*iprint >= 1) {
            printf(" Singular triangular system detected;\n"
                   "   refresh the lbfgs memory and restart the iteration.\n");