lbfgsb.c

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

/* opt/lbfgsb.f -- translated by f2c (version 20050501).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

#ifdef __cplusplus
extern "C" {
#endif
#include "v3p_netlib.h"

#undef abs
#undef min
#undef max
#include <math.h>
#include <stdio.h>
#define abs(x) ((x) >= 0 ? (x) : -(x))
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))

/* Code below must be fixed to pass a real FILE* value to fprintf if
   this macro is defined.  */
#undef LBFGSB_ENABLE_ITERATE_FILE

static void lbfgsb_printf_vec(const char* name,
                              double const* v,
                              integer len)
{
/*<             write (6,1004) 'G =',(g(i), i = 1, n) >*/
/*
  1004 format (/,a4, 1p, 6(1x,d11.4),/,(4x,1p,6(1x,d11.4)))
*/
  integer i;
  printf("%s =", name);
  for(i = 1; i <= len; ++i)
    {
    printf(" %11.4g", v[i]);
    }
  printf("\n");
}

/* Table of constant values */

static doublereal c_b9 = 0.;
static integer c__1 = 1;
static integer c__11 = 11;
static doublereal c_b275 = .001;
static doublereal c_b276 = .9;
static doublereal c_b277 = .1;

/* ================    L-BFGS-B (version 2.1)   ========================== */
/*<    >*/
/* Subroutine */ int setulb_(integer *n, integer *m, doublereal *x, 
        doublereal *l, doublereal *u, integer *nbd, doublereal *f, doublereal 
        *g, doublereal *factr, doublereal *pgtol, doublereal *wa, integer *
        iwa, char *task, integer *iprint, char *csave, logical *lsave, 
        integer *isave, doublereal *dsave)
{
    /* System generated locals */
    integer i__1;

    /* Builtin functions */
    integer s_cmp(char *, char *, ftnlen, ftnlen);

    /* Local variables */
    integer l1, l2, l3, ld, lr, lt, lz, lwa, lsg, lyg, lwn, lss, lws, lwt, 
            lsy, lwy, lyy, lsnd, lsgo, lygo;
    extern /* Subroutine */ int mainlb_(integer *, integer *, doublereal *, 
            doublereal *, doublereal *, integer *, doublereal *, doublereal *,
             doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *, doublereal *, doublereal *, doublereal *, 
            doublereal *, doublereal *, doublereal *, integer *, integer *, 
            integer *, char *, integer *, char *, logical *, integer *, 
            doublereal *, ftnlen, ftnlen);

/*<       character*60     task, csave >*/
/*<       logical          lsave(4) >*/
/*<    >*/
/*<    >*/
/*     ************ */

/*     Subroutine setulb */

/*     This subroutine partitions the working arrays wa and iwa, and */
/*       then uses the limited memory BFGS method to solve the bound */
/*       constrained optimization problem by calling mainlb. */
/*       (The direct method will be used in the subspace minimization.) */

/*     n is an integer variable. */
/*       On entry n is the dimension of the problem. */
/*       On exit n is unchanged. */

/*     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. */

/*     x is a double precision array of dimension n. */
/*       On entry x is an approximation to the solution. */
/*       On exit x is the current approximation. */

/*     l is a double precision array of dimension n. */
/*       On entry l is the lower bound on x. */
/*       On exit l is unchanged. */

/*     u is a double precision array of dimension n. */
/*       On entry u is the upper bound on x. */
/*       On exit u is unchanged. */

/*     nbd is an integer array of dimension n. */
/*       On entry nbd represents the type of bounds imposed on the */
/*         variables, and must be specified as follows: */
/*         nbd(i)=0 if x(i) is unbounded, */
/*                1 if x(i) has only a lower bound, */
/*                2 if x(i) has both lower and upper bounds, and */
/*                3 if x(i) has only an upper bound. */
/*       On exit nbd is unchanged. */

/*     f is a double precision variable. */
/*       On first entry f is unspecified. */
/*       On final exit f is the value of the function at x. */

/*     g is a double precision array of dimension n. */
/*       On first entry g is unspecified. */
/*       On final exit g is the value of the gradient at x. */

/*     factr is a double precision variable. */
/*       On entry factr >= 0 is specified by the user.  The iteration */
/*         will stop when */

/*         (f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr*epsmch */

/*         where epsmch is the machine precision, which is automatically */
/*         generated by the code. Typical values for factr: 1.d+12 for */
/*         low accuracy; 1.d+7 for moderate accuracy; 1.d+1 for extremely */
/*         high accuracy. */
/*       On exit factr is unchanged. */

/*     pgtol is a double precision variable. */
/*       On entry pgtol >= 0 is specified by the user.  The iteration */
/*         will stop when */

/*                 max{|proj g_i | i = 1, ..., n} <= pgtol */

/*         where pg_i is the ith component of the projected gradient. */
/*       On exit pgtol is unchanged. */

/*     wa is a double precision working array of length */
/*       (2mmax + 4)nmax + 12mmax^2 + 12mmax. */

/*     iwa is an integer working array of length 3nmax. */

/*     task is a working string of characters of length 60 indicating */
/*       the current job when entering and quitting this subroutine. */

/*     iprint is an integer variable that must be set by the user. */
/*       It controls the frequency and type of output generated: */
/*        iprint<0    no output is generated; */
/*        iprint=0    print only one line at the last iteration; */
/*        0<iprint<99 print also f and |proj g| every iprint iterations; */
/*        iprint=99   print details of every iteration except n-vectors; */
/*        iprint=100  print also the changes of active set and final x; */
/*        iprint>100  print details of every iteration including x and g; */
/*       When iprint > 0, the file iterate.dat will be created to */
/*                        summarize the iteration. */

/*     csave is a working string of characters of length 60. */

/*     lsave is a logical working array of dimension 4. */
/*       On exit with 'task' = NEW_X, the following information is */
/*                                                             available: */
/*         If lsave(1) = .true.  then  the initial X has been replaced by */
/*                                     its projection in the feasible set; */
/*         If lsave(2) = .true.  then  the problem is constrained; */
/*         If lsave(3) = .true.  then  each variable has upper and lower */
/*                                     bounds; */

/*     isave is an integer working array of dimension 44. */
/*       On exit with 'task' = NEW_X, the following information is */
/*                                                             available: */
/*         isave(22) = the total number of intervals explored in the */
/*                         search of Cauchy points; */
/*         isave(26) = the total number of skipped BFGS updates before */
/*                         the current iteration; */
/*         isave(30) = the number of current iteration; */
/*         isave(31) = the total number of BFGS updates prior the current */
/*                         iteration; */
/*         isave(33) = the number of intervals explored in the search of */
/*                         Cauchy point in the current iteration; */
/*         isave(34) = the total number of function and gradient */
/*                         evaluations; */
/*         isave(36) = the number of function value or gradient */
/*                                  evaluations in the current iteration; */
/*         if isave(37) = 0  then the subspace argmin is within the box; */
/*         if isave(37) = 1  then the subspace argmin is beyond the box; */
/*         isave(38) = the number of free variables in the current */
/*                         iteration; */
/*         isave(39) = the number of active constraints in the current */
/*                         iteration; */
/*         n + 1 - isave(40) = the number of variables leaving the set of */
/*                           active constraints in the current iteration; */
/*         isave(41) = the number of variables entering the set of active */
/*                         constraints in the current iteration. */

/*     dsave is a double precision working array of dimension 29. */
/*       On exit with 'task' = NEW_X, the following information is */
/*                                                             available: */
/*         dsave(1) = current 'theta' in the BFGS matrix; */
/*         dsave(2) = f(x) in the previous iteration; */
/*         dsave(3) = factr*epsmch; */
/*         dsave(4) = 2-norm of the line search direction vector; */
/*         dsave(5) = the machine precision epsmch generated by the code; */
/*         dsave(7) = the accumulated time spent on searching for */
/*                                                         Cauchy points; */
/*         dsave(8) = the accumulated time spent on */
/*                                                 subspace minimization; */
/*         dsave(9) = the accumulated time spent on line search; */
/*         dsave(11) = the slope of the line search function at */
/*                                  the current point of line search; */
/*         dsave(12) = the maximum relative step length imposed in */
/*                                                           line search; */
/*         dsave(13) = the infinity norm of the projected gradient; */
/*         dsave(14) = the relative step length in the line search; */
/*         dsave(15) = the slope of the line search function at */
/*                                 the starting point of the line search; */
/*         dsave(16) = the square of the 2-norm of the line search */
/*                                                      direction vector. */

/*     Subprograms called: */

/*       L-BFGS-B Library ... mainlb. */


/*     References: */

/*       [1] R. H. Byrd, P. Lu, J. Nocedal and C. Zhu, ``A limited */
/*       memory algorithm for bound constrained optimization'', */
/*       SIAM J. Scientific Computing 16 (1995), no. 5, pp. 1190--1208. */

/*       [2] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, ``L-BFGS-B: a */
/*       limited memory FORTRAN code for solving bound constrained */
/*       optimization problems'', Tech. Report, NAM-11, EECS Department, */
/*       Northwestern University, 1994. */

/*       (Postscript files of these papers are available via anonymous */
/*        ftp to eecs.nwu.edu in the directory pub/lbfgs/lbfgs_bcm.) */

/*                           *  *  * */

/*     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. */


/*     ************ */
/*<    >*/
/*<       if (task .eq. 'START') then >*/
    /* Parameter adjustments */
    --iwa;
    --g;
    --nbd;
    --u;
    --l;
    --x;
    --wa;
    --lsave;
    --isave;
    --dsave;

    /* Function Body */
    if (s_cmp(task, "START", (ftnlen)60, (ftnlen)5) == 0) {
/*<          isave(1)  = m*n >*/
        isave[1] = *m * *n;
/*<          isave(2)  = m**2 >*/
/* Computing 2nd power */
        i__1 = *m;
        isave[2] = i__1 * i__1;
/*<          isave(3)  = 4*m**2 >*/
/* Computing 2nd power */
        i__1 = *m;
        isave[3] = i__1 * i__1 << 2;
/*<          isave(4)  = 1 >*/
        isave[4] = 1;
/*<          isave(5)  = isave(4)  + isave(1) >*/
        isave[5] = isave[4] + isave[1];
/*<          isave(6)  = isave(5)  + isave(1) >*/
        isave[6] = isave[5] + isave[1];
/*<          isave(7)  = isave(6)  + isave(2) >*/
        isave[7] = isave[6] + isave[2];
/*<          isave(8)  = isave(7)  + isave(2) >*/
        isave[8] = isave[7] + isave[2];
/*<          isave(9)  = isave(8)  + isave(2) >*/
        isave[9] = isave[8] + isave[2];
/*<          isave(10) = isave(9)  + isave(2) >*/
        isave[10] = isave[9] + isave[2];
/*<          isave(11) = isave(10) + isave(3) >*/
        isave[11] = isave[10] + isave[3];
/*<          isave(12) = isave(11) + isave(3) >*/
        isave[12] = isave[11] + isave[3];
/*<          isave(13) = isave(12) + n >*/
        isave[13] = isave[12] + *n;
/*<          isave(14) = isave(13) + n >*/
        isave[14] = isave[13] + *n;
/*<          isave(15) = isave(14) + n >*/
        isave[15] = isave[14] + *n;
/*<          isave(16) = isave(15) + n >*/
        isave[16] = isave[15] + *n;
/*<          isave(17) = isave(16) + 8*m >*/
        isave[17] = isave[16] + (*m << 3);
/*<          isave(18) = isave(17) + m >*/
        isave[18] = isave[17] + *m;
/*<          isave(19) = isave(18) + m >*/
        isave[19] = isave[18] + *m;
/*<          isave(20) = isave(19) + m    >*/
        isave[20] = isave[19] + *m;
/*<       endif >*/
    }
/*<       l1   = isave(1) >*/
    l1 = isave[1];
/*<       l2   = isave(2) >*/
    l2 = isave[2];
/*<       l3   = isave(3) >*/
    l3 = isave[3];
/*<       lws  = isave(4) >*/
    lws = isave[4];
/*<       lwy  = isave(5) >*/
    lwy = isave[5];
/*<       lsy  = isave(6) >*/
    lsy = isave[6];
/*<       lss  = isave(7) >*/
    lss = isave[7];
/*<       lyy  = isave(8) >*/
    lyy = isave[8];
/*<       lwt  = isave(9) >*/
    lwt = isave[9];
/*<       lwn  = isave(10) >*/
    lwn = isave[10];