scitbx::lbfgsb::raw Namespace Reference

C++ port of the raw FORTRAN interface of L-BFGS-B Version 2.1. More...


Classes

class  ref1
 Emulation of 1-dimensional FORTRAN arrays with offset 1. More...
class  ref2
 Emulation of 2-dimensional FORTRAN arrays with offset 1. More...

Functions

template<typename T>
max3 (T const &v0, T const &v1, T const &v2)
 Emulation of the intrinsic function max for three arguments.
template<typename FloatType>
void timer (FloatType &ttime)
 Current user time for the process.
template<typename ElementType>
void write_ref1 (std::string const &label, ref1< ElementType > const &a, const char *fmt=" %11.4E")
 Emulation of write statement with implicit loop.
template<typename FloatType>
void dcopy (int const &n, ref1< FloatType > const &dx, int const &incx, ref1< FloatType > const &dy, int const &incy)
 copies a vector, x, to a vector, y.
template<typename FloatType>
void dscal (int const &n, FloatType const &da, ref1< FloatType > const &dx, int const &incx)
 scales a vector by a constant.
template<typename FloatType, typename SizeType>
void daxpy (SizeType const &n, FloatType const &da, const FloatType *dx, FloatType *dy)
 constant times a vector plus a vector.
template<typename FloatType>
void projgr (int const &n, ref1< FloatType > const &l, ref1< FloatType > const &u, ref1< int > const &nbd, ref1< FloatType > const &x, ref1< FloatType > const &g, FloatType &sbgnrm)
 This subroutine computes the infinity norm of the projected gradient.
template<typename FloatType>
void dtrsl (ref2< FloatType > const &t, int const &, int const &n, ref1< FloatType > const &b, int const &job, int &info)
 Solution of t * x = b.
template<typename FloatType>
void bmv (int const &m, ref2< FloatType > const &sy, ref2< FloatType > const &wt, int const &col, ref1< FloatType > const &v, ref1< FloatType > const &p, int &info)
 Product of the 2m x 2m middle matrix in the compact L-BFGS formula of B and a 2m vector v.
template<typename FloatType>
void hpsolb (int const &n, ref1< FloatType > const &t, ref1< int > const &iorder, int const &iheap)
 Sorting of t.
template<typename FloatType>
void errclb (int const &n, int const &m, FloatType const &factr, ref1< FloatType > const &l, ref1< FloatType > const &u, ref1< int > const &nbd, std::string &task, int &info, int &k)
 Check validity of input data.
template<typename FloatType>
void prn1lb (int const &n, int const &m, ref1< FloatType > const &l, ref1< FloatType > const &u, ref1< FloatType > const &x, int const &iprint, int const &, FloatType const &epsmch)
 Printing function 1.
template<typename FloatType>
void prn2lb (int const &, ref1< FloatType > const &x, FloatType const &f, ref1< FloatType > const &g, int const &iprint, int const &, int const &iter, int const &nfgv, int const &nact, FloatType const &sbgnrm, int const &nint, std::string &word, int const &iword, int const &iback, FloatType const &stp, FloatType const &xstep)
 Printing function 2.
template<typename FloatType>
void prn3lb (int const &n, ref1< FloatType > const &x, FloatType const &f, std::string const &task, int const &iprint, int const &info, int const &, int const &iter, int const &nfgv, int const &nintol, int const &nskip, int const &nact, FloatType const &sbgnrm, FloatType const &time, int const &nint, std::string const &word, int const &iback, FloatType const &stp, FloatType const &xstep, int const &k, FloatType const &cachyt, FloatType const &sbtime, FloatType const &lnscht)
 Printing function 3.
template<typename FloatType>
void active (int const &n, ref1< FloatType > const &l, ref1< FloatType > const &u, ref1< int > const &nbd, ref1< FloatType > const &x, ref1< int > const &iwhere, int const &iprint, bool &prjctd, bool &cnstnd, bool &boxed)
 Initializes iwhere and projects the initial x to the feasible set if necessary.
template<typename FloatType>
bool cauchy_loop (int const &n, ref1< FloatType > const &x, ref1< FloatType > const &l, ref1< FloatType > const &u, ref1< int > const &iorder, ref1< int > const &iwhere, ref1< FloatType > const &t, ref1< FloatType > const &d, ref1< FloatType > const &xcp, int const &m, ref2< FloatType > const &wy, ref2< FloatType > const &ws, ref2< FloatType > const &sy, ref2< FloatType > const &wt, FloatType const &theta, int const &col, int const &head, ref1< FloatType > const &p, ref1< FloatType > const &c, ref1< FloatType > const &wbp, ref1< FloatType > const &v, int &nint, int const &iprint, int &info, FloatType const &epsmch, FloatType const &bkmin, bool const &bnded, int const &col2, FloatType &dtm, FloatType &f1, FloatType &f2, FloatType &f2_org, int const &ibkmin, int const &nbreak, FloatType &tsum)
 Fragment from subroutine cauchy.
template<typename FloatType>
void cauchy (int const &n, ref1< FloatType > const &x, ref1< FloatType > const &l, ref1< FloatType > const &u, ref1< int > const &nbd, ref1< FloatType > const &g, ref1< int > const &iorder, ref1< int > const &iwhere, ref1< FloatType > const &t, ref1< FloatType > const &d, ref1< FloatType > const &xcp, int const &m, ref2< FloatType > const &wy, ref2< FloatType > const &ws, ref2< FloatType > const &sy, ref2< FloatType > const &wt, FloatType const &theta, int const &col, int const &head, ref1< FloatType > const &p, ref1< FloatType > const &c, ref1< FloatType > const &wbp, ref1< FloatType > const &v, int &nint, ref1< FloatType > const &, ref1< FloatType > const &, int const &iprint, FloatType const &sbgnrm, int &info, FloatType const &epsmch)
 Computation of the generalized Cauchy point (GCP).
void freev (int const &n, int &nfree, ref1< int > const &index, int &nenter, int &ileave, ref1< int > const &indx2, ref1< int > const &iwhere, bool &wrk, bool const &updatd, bool const &cnstnd, int const &iprint, int const &iter)
 Determination of the index set of free and active variables at the GCP.
template<typename FloatType>
void dpofa (ref2< FloatType > const &a, int const &, int const &n, int &info)
 Factorization of a symmetric positive definite matrix.
template<typename FloatType>
void formk (int const &n, int const &nsub, ref1< int > const &ind, int const &nenter, int const &ileave, ref1< int > const &indx2, int const &iupdat, bool const &updatd, ref2< FloatType > const &wn, ref2< FloatType > const &wn1, int const &m, ref2< FloatType > const &ws, ref2< FloatType > const &wy, ref2< FloatType > const &sy, FloatType const &theta, int const &col, int const &head, int &info)
 LEL^T factorization.
template<typename FloatType>
void cmprlb (int const &n, int const &m, ref1< FloatType > const &x, ref1< FloatType > const &g, ref2< FloatType > const &ws, ref2< FloatType > const &wy, ref2< FloatType > const &sy, ref2< FloatType > const &wt, ref1< FloatType > const &z, ref1< FloatType > const &r, ref1< FloatType > const &wa, ref1< int > const &index, FloatType const &theta, int const &col, int const &head, int const &nfree, bool const &cnstnd, int &info)
 Computation of -Z'B(xcp-xk)-Z'g.
template<typename FloatType>
void subsm (int const &, int const &m, int const &nsub, ref1< int > const &ind, ref1< FloatType > const &l, ref1< FloatType > const &u, ref1< int > const &nbd, ref1< FloatType > const &x, ref1< FloatType > const &d, ref2< FloatType > const &ws, ref2< FloatType > const &wy, FloatType const &theta, int const &col, int const &head, int &iword, ref1< FloatType > const &wv, ref2< FloatType > const &wn, int const &iprint, int &info)
 Approximate solution of the subspace problem.
template<typename FloatType>
void dcstep (FloatType &stx, FloatType &fx, FloatType &dx, FloatType &sty, FloatType &fy, FloatType &dy, FloatType &stp, FloatType const &fp, FloatType const &dp, bool &brackt, FloatType const &stpmin, FloatType const &stpmax)
 Computation of a safeguarded step for a search procedure.
template<typename FloatType>
void dcsrch (FloatType const &f, FloatType const &g, FloatType &stp, FloatType const &ftol, FloatType const &gtol, FloatType const &xtol, FloatType const &stpmin, FloatType const &stpmax, std::string &task, ref1< int > const &isave, ref1< FloatType > const &dsave)
 Determination of a step that satisfies a sufficient decrease condition and a curvature condition.
template<typename FloatType>
void lnsrlb (int const &n, ref1< FloatType > const &l, ref1< FloatType > const &u, ref1< int > const &nbd, ref1< FloatType > const &x, FloatType const &f, FloatType &fold, FloatType &gd, FloatType &gdold, ref1< FloatType > const &g, ref1< FloatType > const &d, ref1< FloatType > const &r, ref1< FloatType > const &t, ref1< FloatType > const &z, FloatType &stp, FloatType &dnorm, FloatType &dtd, FloatType &xstep, FloatType &stpmx, int const &iter, int &ifun, int &iback, int &nfgv, int &info, std::string &task, bool const &boxed, bool const &cnstnd, std::string &csave, ref1< int > const &isave, ref1< FloatType > const &dsave, bool enable_stp_init)
 Line search.
template<typename FloatType>
void matupd (int const &n, int const &m, ref2< FloatType > const &ws, ref2< FloatType > const &wy, ref2< FloatType > const &sy, ref2< FloatType > const &ss, ref1< FloatType > const &d, ref1< FloatType > const &r, int &itail, int const &iupdat, int &col, int &head, FloatType &theta, FloatType const &rr, FloatType const &dr, FloatType const &stp, FloatType const &dtd)
 Updates matrices WS and WY, and forms the middle matrix in B.
template<typename FloatType>
void formt (int const &m, ref2< FloatType > const &wt, ref2< FloatType > const &sy, ref2< FloatType > const &ss, int const &col, FloatType const &theta, int &info)
 Forms the upper half of the pos. def. and symm. T.
template<typename FloatType>
void mainlb (int const &n, int const &m, ref1< FloatType > const &x, ref1< FloatType > const &l, ref1< FloatType > const &u, ref1< int > const &nbd, FloatType &f, ref1< FloatType > const &g, FloatType const &factr, FloatType const &pgtol, ref2< FloatType > const &ws, ref2< FloatType > const &wy, ref2< FloatType > const &sy, ref2< FloatType > const &ss, ref2< FloatType > const &, ref2< FloatType > const &wt, ref2< FloatType > const &wn, ref2< FloatType > const &snd, ref1< FloatType > const &z, ref1< FloatType > const &r, ref1< FloatType > const &d, ref1< FloatType > const &t, ref1< FloatType > const &wa, ref1< FloatType > const &sg, ref1< FloatType > const &, ref1< FloatType > const &yg, ref1< FloatType > const &, ref1< int > const &index, ref1< int > const &iwhere, ref1< int > const &indx2, std::string &task, int const &iprint, std::string &csave, ref1< bool > const &lsave, ref1< int > const &isave, ref1< FloatType > const &dsave, bool enable_stp_init)
 Solution of bound constrained optimization problems.
template<typename FloatType>
void setulb (int const &n, int const &m, ref1< FloatType > const &x, ref1< FloatType > const &l, ref1< FloatType > const &u, ref1< int > const &nbd, FloatType &f, ref1< FloatType > const &g, FloatType const &factr, FloatType const &pgtol, ref1< FloatType > const &wa, ref1< int > const &iwa, std::string &task, int const &iprint, std::string &csave, ref1< bool > const &lsave, ref1< int > const &isave, ref1< FloatType > const &dsave, bool enable_stp_init=false)
 Main L-BFGS-B interface function.


Detailed Description

C++ port of the raw FORTRAN interface of L-BFGS-B Version 2.1.

Original FORTRAN distribution:

http://www.ece.northwestern.edu/~nocedal/lbfgsb.html

Written by Ciyou Zhu in collaboration with R.H. Byrd, P. Lu-Chen and J. Nocedal.

C++ port by Ralf W. Grosse-Kunstleve.


Function Documentation

void scitbx::lbfgsb::raw::active ( int const &  n,
ref1< FloatType > const &  l,
ref1< FloatType > const &  u,
ref1< int > const &  nbd,
ref1< FloatType > const &  x,
ref1< int > const &  iwhere,
int const &  iprint,
bool &  prjctd,
bool &  cnstnd,
bool &  boxed 
) [inline]

Initializes iwhere and projects the initial x to the feasible set if necessary.

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.

Referenced by mainlb().

void scitbx::lbfgsb::raw::bmv ( int const &  m,
ref2< FloatType > const &  sy,
ref2< FloatType > const &  wt,
int const &  col,
ref1< FloatType > const &  v,
ref1< FloatType > const &  p,
int &  info 
) [inline]

Product of the 2m x 2m middle matrix in the compact L-BFGS formula of B and a 2m vector v.

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 upper triangular matrix J' which is the Cholesky factor of (thetaS'S+LD^(-1)L'). On exit wt is unchanged.

col is an integer variable. On entry col specifies the number of s-vectors (or y-vectors) stored in the compact L-BFGS formula. On exit col is unchanged.

v is a double precision array of dimension 2col. On entry v specifies vector v. On exit v is unchanged.

p is a double precision array of dimension 2col. On entry p is unspecified. On exit p is the product Mv.

info is an integer variable. On entry info is unspecified. On exit info = 0 for normal return, = nonzero for abnormal return when the system to be solved by dtrsl is singular.

Subprograms called:

Linpack ... dtrsl.

* *

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.

References dtrsl(), ref1::get1(), ref2::get2(), SCITBX_ASSERT, and const_ref::size().

Referenced by cauchy(), cauchy_loop(), and cmprlb().

void scitbx::lbfgsb::raw::cauchy ( int const &  n,
ref1< FloatType > const &  x,
ref1< FloatType > const &  l,
ref1< FloatType > const &  u,
ref1< int > const &  nbd,
ref1< FloatType > const &  g,
ref1< int > const &  iorder,
ref1< int > const &  iwhere,
ref1< FloatType > const &  t,
ref1< FloatType > const &  d,
ref1< FloatType > const &  xcp,
int const &  m,
ref2< FloatType > const &  wy,
ref2< FloatType > const &  ws,
ref2< FloatType > const &  sy,
ref2< FloatType > const &  wt,
FloatType const &  theta,
int const &  col,
int const &  head,
ref1< FloatType > const &  p,
ref1< FloatType > const &  c,
ref1< FloatType > const &  wbp,
ref1< FloatType > const &  v,
int &  nint,
ref1< FloatType > const &  ,
ref1< FloatType > const &  ,
int const &  iprint,
FloatType const &  sbgnrm,
int &  info,
FloatType const &  epsmch 
) [inline]

Computation of the generalized Cauchy point (GCP).

For given x, l, u, g (with sbgnrm > 0), and a limited memory BFGS matrix B defined in terms of matrices WY, WS, WT, and scalars head, col, and theta, this subroutine computes the generalized Cauchy point (GCP), defined as the first local minimizer of the quadratic

Q(x + s) = g's + 1/2 s'Bs

along the projected gradient direction P(x-tg,l,u). The routine returns the GCP in xcp.

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

x is a double precision array of dimension n. On entry x is the starting point for the GCP computation. On exit x is unchanged.

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

u is a double precision array of dimension n. On entry u is the upper bound of 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.

g is a double precision array of dimension n. On entry g is the gradient of f(x). g must be a nonzero vector. On exit g is unchanged.

iorder is an integer working array of dimension n. iorder will be used to store the breakpoints in the piecewise linear path and free variables encountered. On exit, iorder(1),...,iorder(nleft) are indices of breakpoints which have not been encountered; iorder(nleft+1),...,iorder(nbreak) are indices of encountered breakpoints; and iorder(nfree),...,iorder(n) are indices of variables which have no bound constraits along the search direction.

iwhere is an integer array of dimension n. On entry iwhere indicates only the permanently fixed (iwhere=3) or free (iwhere= -1) components of x. On exit iwhere records the status of the current x variables. iwhere(i)=-3 if x(i) is free and has bounds, but is not moved 0 if x(i) is free and has bounds, and is moved 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) -1 if x(i) is always free, i.e., it has no bounds.

t is a double precision working array of dimension n. t will be used to store the break points.

d is a double precision array of dimension n used to store the Cauchy direction P(x-tg)-x.

xcp is a double precision array of dimension n used to return the GCP on exit.

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.

ws, wy, sy, and wt are double precision arrays. On entry they store information that defines the limited memory BFGS matrix: ws(n,m) stores S, a set of s-vectors; wy(n,m) stores Y, a set of y-vectors; sy(m,m) stores S'Y; wt(m,m) stores the Cholesky factorization of (theta*S'S+LD^(-1)L'). On exit these arrays are unchanged.

theta is a double precision variable. On entry theta is the scaling factor specifying B_0 = theta I. On exit theta is unchanged.

col is an integer variable. On entry col is the actual number of variable metric corrections stored so far. On exit col is unchanged.

head is an integer variable. On entry head is the location of the first s-vector (or y-vector) in S (or Y). On exit col is unchanged.

p is a double precision working array of dimension 2m. p will be used to store the vector p = W^(T)d.

c is a double precision working array of dimension 2m. c will be used to store the vector c = W^(T)(xcp-x).

wbp is a double precision working array of dimension 2m. wbp will be used to store the row of W corresponding to a breakpoint.

v is a double precision working array of dimension 2m.

nint is an integer variable. On exit nint records the number of quadratic segments explored in searching for the GCP.

sg and yg are double precision arrays of dimension m. On entry sg and yg store S'g and Y'g correspondingly. On exit they are unchanged.

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.

sbgnrm is a double precision variable. On entry sbgnrm is the norm of the projected gradient at x. On exit sbgnrm is unchanged.

info is an integer variable. On entry info is 0. On exit info = 0 for normal return, = nonzero for abnormal return when the the system used in routine bmv is singular.

Subprograms called:

L-BFGS-B Library ... hpsolb, bmv.

Linpack ... dscal dcopy, daxpy.

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: FORTRAN Subroutines for Large Scale Bound Constrained Optimization'' 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.

References ref::begin(), bmv(), cauchy_loop(), daxpy(), dcopy(), dscal(), ref1::get1(), SCITBX_ASSERT, const_ref::size(), and write_ref1().

Referenced by mainlb().

void scitbx::lbfgsb::raw::cmprlb ( int const &  n,
int const &  m,
ref1< FloatType > const &  x,
ref1< FloatType > const &  g,
ref2< FloatType > const &  ws,
ref2< FloatType > const &  wy,
ref2< FloatType > const &  sy,
ref2< FloatType > const &  wt,
ref1< FloatType > const &  z,
ref1< FloatType > const &  r,
ref1< FloatType > const &  wa,
ref1< int > const &  index,
FloatType const &  theta,
int const &  col,
int const &  head,
int const &  nfree,
bool const &  cnstnd,
int &  info 
) [inline]

Computation of -Z'B(xcp-xk)-Z'g.

This subroutine computes r=-Z'B(xcp-xk)-Z'g by using wa(2m+1)=W'(xcp-x) from subroutine cauchy.

Subprograms called:

L-BFGS-B Library ... bmv.

* *

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.

References bmv(), ref1::get1(), SCITBX_ASSERT, and const_ref::size().

Referenced by mainlb().

void scitbx::lbfgsb::raw::dcopy ( int const &  n,
ref1< FloatType > const &  dx,
int const &  incx,
ref1< FloatType > const &  dy,
int const &  incy 
) [inline]

copies a vector, x, to a vector, y.

uses unrolled loops for increments equal to one. jack dongarra, linpack, 3/11/78.

Referenced by cauchy(), formk(), lnsrlb(), mainlb(), and matupd().

void scitbx::lbfgsb::raw::dcsrch ( FloatType const &  f,
FloatType const &  g,
FloatType &  stp,
FloatType const &  ftol,
FloatType const &  gtol,
FloatType const &  xtol,
FloatType const &  stpmin,
FloatType const &  stpmax,
std::string &  task,
ref1< int > const &  isave,
ref1< FloatType > const &  dsave 
) [inline]

Determination of a step that satisfies a sufficient decrease condition and a curvature condition.

This subroutine finds a step that satisfies a sufficient decrease condition and a curvature condition.

Each call of the subroutine updates an interval with endpoints stx and sty. The interval is initially chosen so that it contains a minimizer of the modified function

psi(stp) = f(stp) - f(0) - ftol*stp*f'(0).

If psi(stp) <= 0 and f'(stp) >= 0 for some step, then the interval is chosen so that it contains a minimizer of f.

The algorithm is designed to find a step that satisfies the sufficient decrease condition

f(stp) <= f(0) + ftol*stp*f'(0),

and the curvature condition

abs(f'(stp)) <= gtol*abs(f'(0)).

If ftol is less than gtol and if, for example, the function is bounded below, then there is always a step which satisfies both conditions.

If no step can be found that satisfies both conditions, then the algorithm stops with a warning. In this case stp only satisfies the sufficient decrease condition.

A typical invocation of dcsrch has the following outline:

task = 'START' continue call dcsrch( ... ) if (task .eq. 'FG') then Evaluate the function and the gradient at stp goto 10 end if

NOTE: The user must no alter work arrays between calls.

The subroutine statement is

subroutine dcsrch(f,g,stp,ftol,gtol,xtol,stpmin,stpmax, task,isave,dsave) where

f is a double precision variable. On initial entry f is the value of the function at 0. On subsequent entries f is the value of the function at stp. On exit f is the value of the function at stp.

g is a double precision variable. On initial entry g is the derivative of the function at 0. On subsequent entries g is the derivative of the function at stp. On exit g is the derivative of the function at stp.

stp is a double precision variable. On entry stp is the current estimate of a satisfactory step. On initial entry, a positive initial estimate must be provided. On exit stp is the current estimate of a satisfactory step if task = 'FG'. If task = 'CONV' then stp satisfies the sufficient decrease and curvature condition.

ftol is a double precision variable. On entry ftol specifies a nonnegative tolerance for the sufficient decrease condition. On exit ftol is unchanged.

gtol is a double precision variable. On entry gtol specifies a nonnegative tolerance for the curvature condition. On exit gtol is unchanged.

xtol is a double precision variable. On entry xtol specifies a nonnegative relative tolerance for an acceptable step. The subroutine exits with a warning if the relative difference between sty and stx is less than xtol. On exit xtol is unchanged.

stpmin is a double precision variable. On entry stpmin is a nonnegative lower bound for the step. On exit stpmin is unchanged.

stpmax is a double precision variable. On entry stpmax is a nonnegative upper bound for the step. On exit stpmax is unchanged.

task is a character variable of length at least 60. On initial entry task must be set to 'START'. On exit task indicates the required action:

If task(1:2) = 'FG' then evaluate the function and derivative at stp and call dcsrch again.

If task(1:4) = 'CONV' then the search is successful.

If task(1:4) = 'WARN' then the subroutine is not able to satisfy the convergence conditions. The exit value of stp contains the best point found during the search.

If task(1:5) = 'ERROR' then there is an error in the input arguments.

On exit with convergence, a warning or an error, the variable task contains additional information.

isave is an integer work array of dimension 2.

dsave is a double precision work array of dimension 13.

Subprograms called

MINPACK-2 ... dcstep

MINPACK-1 Project. June 1983. Argonne National Laboratory. Jorge J. More' and David J. Thuente.

MINPACK-2 Project. October 1993. Argonne National Laboratory and University of Minnesota. Brett M. Averick, Richard G. Carter, and Jorge J. More'.

References dcstep(), SCITBX_ASSERT, and const_ref::size().

Referenced by lnsrlb().

void scitbx::lbfgsb::raw::dcstep ( FloatType &  stx,
FloatType &  fx,
FloatType &  dx,
FloatType &  sty,
FloatType &  fy,
FloatType &  dy,
FloatType &  stp,
FloatType const &  fp,
FloatType const &  dp,
bool &  brackt,
FloatType const &  stpmin,
FloatType const &  stpmax 
) [inline]

Computation of a safeguarded step for a search procedure.

This subroutine computes a safeguarded step for a search procedure and updates an interval that contains a step that satisfies a sufficient decrease and a curvature condition.

The parameter stx contains the step with the least function value. If brackt is set to .true. then a minimizer has been bracketed in an interval with endpoints stx and sty. The parameter stp contains the current step. The subroutine assumes that if brackt is set to .true. then

min(stx,sty) < stp < max(stx,sty),

and that the derivative at stx is negative in the direction of the step.

The subroutine statement is

subroutine dcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,brackt, stpmin,stpmax)

where

stx is a double precision variable. On entry stx is the best step obtained so far and is an endpoint of the interval that contains the minimizer. On exit stx is the updated best step.

fx is a double precision variable. On entry fx is the function at stx. On exit fx is the function at stx.

dx is a double precision variable. On entry dx is the derivative of the function at stx. The derivative must be negative in the direction of the step, that is, dx and stp - stx must have opposite signs. On exit dx is the derivative of the function at stx.

sty is a double precision variable. On entry sty is the second endpoint of the interval that contains the minimizer. On exit sty is the updated endpoint of the interval that contains the minimizer.

fy is a double precision variable. On entry fy is the function at sty. On exit fy is the function at sty.

dy is a double precision variable. On entry dy is the derivative of the function at sty. On exit dy is the derivative of the function at the exit sty.

stp is a double precision variable. On entry stp is the current step. If brackt is set to .true. then on input stp must be between stx and sty. On exit stp is a new trial step.

fp is a double precision variable. On entry fp is the function at stp On exit fp is unchanged.

dp is a double precision variable. On entry dp is the the derivative of the function at stp. On exit dp is unchanged.

brackt is an logical variable. On entry brackt specifies if a minimizer has been bracketed. Initially brackt must be set to .false. On exit brackt specifies if a minimizer has been bracketed. When a minimizer is bracketed brackt is set to .true.

stpmin is a double precision variable. On entry stpmin is a lower bound for the step. On exit stpmin is unchanged.

stpmax is a double precision variable. On entry stpmax is an upper bound for the step. On exit stpmax is unchanged.

MINPACK-1 Project. June 1983 Argonne National Laboratory. Jorge J. More' and David J. Thuente.

MINPACK-2 Project. October 1993. Argonne National Laboratory and University of Minnesota. Brett M. Averick and Jorge J. More'.

References max3().

Referenced by dcsrch().

void scitbx::lbfgsb::raw::dpofa ( ref2< FloatType > const &  a,
int const &  ,
int const &  n,
int &  info 
) [inline]

Factorization of a symmetric positive definite matrix.

dpofa factors a double precision symmetric positive definite matrix.

dpofa is usually called by dpoco, but it can be called directly with a saving in time if rcond is not needed. (time for dpoco) = (1 + 18/n)*(time for dpofa) .

on entry

a double precision(lda, n) the symmetric matrix to be factored. only the diagonal and upper triangle are used.

lda integer the leading dimension of the array a .

n integer the order of the matrix a .

on return

a an upper triangular matrix r so that a = trans(r)*r where trans(r) is the transpose. the strict lower triangle is unaltered. if info .ne. 0 , the factorization is not complete.

info integer = 0 for normal return. = k signals an error condition. the leading minor of order k is not positive definite.

linpack. this version dated 08/14/78 . cleve moler, university of new mexico, argonne national lab.

subroutines and functions

blas ddot fortran sqrt

References SCITBX_ASSERT, and const_ref::size().

Referenced by formk(), and formt().

void scitbx::lbfgsb::raw::dscal ( int const &  n,
FloatType const &  da,
ref1< FloatType > const &  dx,
int const &  incx 
) [inline]

scales a vector by a constant.

uses unrolled loops for increment equal to one. jack dongarra, linpack, 3/11/78. modified 3/93 to return if incx .le. 0.

Referenced by cauchy(), and mainlb().

void scitbx::lbfgsb::raw::dtrsl ( ref2< FloatType > const &  t,
int const &  ,
int const &  n,
ref1< FloatType > const &  b,
int const &  job,
int &  info 
) [inline]

Solution of t * x = b.

dtrsl solves systems of the form

t * x = b or trans(t) * x = b

where t is a triangular matrix of order n. here trans(t) denotes the transpose of the matrix t.

on entry

t double precision(ldt,n) t contains the matrix of the system. the zero elements of the matrix are not referenced, and the corresponding elements of the array can be used to store other information.

ldt integer ldt is the leading dimension of the array t.

n integer n is the order of the system.

b double precision(n). b contains the right hand side of the system.

job integer job specifies what kind of system is to be solved. if job is

00 solve t*x=b, t lower triangular, 01 solve t*x=b, t upper triangular, 10 solve trans(t)*x=b, t lower triangular, 11 solve trans(t)*x=b, t upper triangular.

on return

b b contains the solution, if info .eq. 0. otherwise b is unaltered.

info integer info contains zero if the system is nonsingular. otherwise info contains the index of the first zero diagonal element of t.

linpack. this version dated 08/14/78 . g. w. stewart, university of maryland, argonne national lab.

subroutines and functions

blas daxpy,ddot fortran mod

References daxpy(), SCITBX_ASSERT, and const_ref::size().

Referenced by bmv(), formk(), and subsm().

void scitbx::lbfgsb::raw::errclb ( int const &  n,
int const &  m,
FloatType const &  factr,
ref1< FloatType > const &  l,
ref1< FloatType > const &  u,
ref1< int > const &  nbd,
std::string &  task,
int &  info,
int &  k 
) [inline]

Check validity of input data.

This subroutine checks the validity of the input data.

* *

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.

Referenced by mainlb().

void scitbx::lbfgsb::raw::formk ( int const &  n,
int const &  nsub,
ref1< int > const &  ind,
int const &  nenter,
int const &  ileave,
ref1< int > const &  indx2,
int const &  iupdat,
bool const &  updatd,
ref2< FloatType > const &  wn,
ref2< FloatType > const &  wn1,
int const &  m,
ref2< FloatType > const &  ws,
ref2< FloatType > const &  wy,
ref2< FloatType > const &  sy,
FloatType const &  theta,
int const &  col,
int const &  head,
int &  info 
) [inline]

LEL^T factorization.

This subroutine forms 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] The matrix K can be shown to be equal to the matrix M^[-1]N occurring in section 5.1 of [1], as well as to the matrix Mbar^[-1] Nbar in section 5.3.

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

nsub is an integer variable On entry nsub is the number of subspace variables in free set. On exit nsub is not changed.

ind is an integer array of dimension nsub. On entry ind specifies the indices of subspace variables. On exit ind is unchanged.

nenter is an integer variable. On entry nenter is the number of variables entering the free set. On exit nenter is unchanged.

ileave is an integer variable. On entry indx2(ileave),...,indx2(n) are the variables leaving the free set. On exit ileave is unchanged.

indx2 is an integer array of dimension n. On entry indx2(1),...,indx2(nenter) are the variables entering the free set, while indx2(ileave),...,indx2(n) are the variables leaving the free set. On exit indx2 is unchanged.

iupdat is an integer variable. On entry iupdat is the total number of BFGS updates made so far. On exit iupdat is unchanged.

updatd is a logical variable. On entry 'updatd' is true if the L-BFGS matrix is updatd. On exit 'updatd' is unchanged.

wn is a double precision array of dimension 2m x 2m. On entry wn is unspecified. On exit the upper triangle of wn stores the LEL^T factorization of the 2*col x 2*col indefinite matrix [-D -Y'ZZ'Y/theta L_a'-R_z' ] [L_a -R_z theta*S'AA'S ]

wn1 is a double precision array of dimension 2m x 2m. On entry wn1 stores the lower triangular part of [Y' ZZ'Y L_a'+R_z'] [L_a+R_z S'AA'S ] in the previous iteration. On exit wn1 stores the corresponding updated matrices. The purpose of wn1 is just to store these inner products so they can be easily updated and inserted into wn.

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.

ws, wy, sy, and wtyy are double precision arrays; theta is a double precision variable; col is an integer variable; head is an integer variable. On entry they store the information defining the limited memory BFGS matrix: ws(n,m) stores S, a set of s-vectors; wy(n,m) stores Y, a set of y-vectors; sy(m,m) stores S'Y; wtyy(m,m) stores the Cholesky factorization of (theta*S'S+LD^(-1)L') theta is the scaling factor specifying B_0 = theta I; col is the number of variable metric corrections stored; head is the location of the 1st s- (or y-) vector in S (or Y). On exit they are unchanged.

info is an integer variable. On entry info is unspecified. On exit info = 0 for normal return; = -1 when the 1st Cholesky factorization failed; = -2 when the 2st Cholesky factorization failed.

Subprograms called:

Linpack ... dcopy, dpofa, dtrsl.

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.

References dcopy(), dpofa(), dtrsl(), ref2::get1(), ref2::get2(), ref2::get2soft(), SCITBX_ASSERT, and const_ref::size().

Referenced by mainlb().

void scitbx::lbfgsb::raw::formt ( int const &  m,
ref2< FloatType > const &  wt,
ref2< FloatType > const &  sy,
ref2< FloatType > const &  ss,
int const &  col,
FloatType const &  theta,
int &  info 
) [inline]

Forms the upper half of the pos. def. and symm. T.

This subroutine forms the upper half of the pos. def. and symm. T = theta*SS + L*D^(-1)*L', stores T in the upper triangle of the array wt, and performs the Cholesky factorization of T to produce J*J', with J' stored in the upper triangle of wt.

Subprograms called:

Linpack ... dpofa.

* *

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.

References dpofa(), ref2::get2(), SCITBX_ASSERT, and const_ref::size().

Referenced by mainlb().

void scitbx::lbfgsb::raw::freev ( int const &  n,
int &  nfree,
ref1< int > const &  index,
int &  nenter,
int &  ileave,
ref1< int > const &  indx2,
ref1< int > const &  iwhere,
bool &  wrk,
bool const &  updatd,
bool const &  cnstnd,
int const &  iprint,
int const &  iter 
) [inline]

Determination of the index set of free and active variables at the GCP.

This subroutine counts the entering and leaving variables when iter > 0, and finds the index set of free and active variables at the GCP.

cnstnd is a logical variable indicating whether bounds are present

index is an integer array of dimension n for i=1,...,nfree, index(i) are the indices of free variables for i=nfree+1,...,n, index(i) are the indices of bound variables On entry after the first iteration, index gives the free variables at the previous iteration. On exit it gives the free variables based on the determination in cauchy using the array iwhere.

indx2 is an integer array of dimension n On entry indx2 is unspecified. On exit with iter>0, indx2 indicates which variables have changed status since the previous iteration. For i= 1,...,nenter, indx2(i) have changed from bound to free. For i= ileave+1,...,n, indx2(i) have changed from free to bound.

* *

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.

References SCITBX_ASSERT, and const_ref::size().

Referenced by mainlb().

void scitbx::lbfgsb::raw::hpsolb ( int const &  n,
ref1< FloatType > const &  t,
ref1< int > const &  iorder,
int const &  iheap 
) [inline]

Sorting of t.

This subroutine sorts out the least element of t, and puts the remaining elements of t in a heap.

n is an integer variable. On entry n is the dimension of the arrays t and iorder. On exit n is unchanged.

t is a double precision array of dimension n. On entry t stores the elements to be sorted, On exit t(n) stores the least elements of t, and t(1) to t(n-1) stores the remaining elements in the form of a heap.

iorder is an integer array of dimension n. On entry iorder(i) is the index of t(i). On exit iorder(i) is still the index of t(i), but iorder may be permuted in accordance with t.

iheap is an integer variable specifying the task. On entry iheap should be set as follows: iheap .eq. 0 if t(1) to t(n) is not in the form of a heap, iheap .ne. 0 if otherwise. On exit iheap is unchanged.

References: Algorithm 232 of CACM (J. W. J. Williams): HEAPSORT.

* *

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.

References SCITBX_ASSERT, and const_ref::size().

Referenced by cauchy_loop().

void scitbx::lbfgsb::raw::lnsrlb ( int const &  n,
ref1< FloatType > const &  l,
ref1< FloatType > const &  u,
ref1< int > const &  nbd,
ref1< FloatType > const &  x,
FloatType const &  f,
FloatType &  fold,
FloatType &  gd,
FloatType &  gdold,
ref1< FloatType > const &  g,
ref1< FloatType > const &  d,
ref1< FloatType > const &  r,
ref1< FloatType > const &  t,
ref1< FloatType > const &  z,
FloatType &  stp,
FloatType &  dnorm,
FloatType &  dtd,
FloatType &  xstep,
FloatType &  stpmx,
int const &  iter,
int &  ifun,
int &  iback,
int &  nfgv,
int &  info,
std::string &  task,
bool const &  boxed,
bool const &  cnstnd,
std::string &  csave,
ref1< int > const &  isave,
ref1< FloatType > const &  dsave,
bool  enable_stp_init 
) [inline]

Line search.

This subroutine calls subroutine dcsrch from the Minpack2 library to perform the line search. Subroutine dscrch is safeguarded so that all trial points lie within the feasible region.

Subprograms called:

Minpack2 Library ... dcsrch.

Linpack ... dtrsl, ddot.

* *

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.

References ref::begin(), dcopy(), dcsrch(), SCITBX_ASSERT, and const_ref::size().

Referenced by mainlb().

void scitbx::lbfgsb::raw::mainlb ( int const &  n,
int const &  m,
ref1< FloatType > const &  x,
ref1< FloatType > const &  l,
ref1< FloatType > const &  u,
ref1< int > const &  nbd,
FloatType &  f,
ref1< FloatType > const &  g,
FloatType const &  factr,
FloatType const &  pgtol,
ref2< FloatType > const &  ws,
ref2< FloatType > const &  wy,
ref2< FloatType > const &  sy,
ref2< FloatType > const &  ss,
ref2< FloatType > const &  ,
ref2< FloatType > const &  wt,
ref2< FloatType > const &  wn,
ref2< FloatType > const &  snd,
ref1< FloatType > const &  z,
ref1< FloatType > const &  r,
ref1< FloatType > const &  d,
ref1< FloatType > const &  t,
ref1< FloatType > const &  wa,
ref1< FloatType > const &  sg,
ref1< FloatType > const &  ,
ref1< FloatType > const &  yg,
ref1< FloatType > const &  ,
ref1< int > const &  index,
ref1< int > const &  iwhere,
ref1< int > const &  indx2,
std::string &  task,
int const &  iprint,
std::string &  csave,
ref1< bool > const &  lsave,
ref1< int > const &  isave,
ref1< FloatType > const &  dsave,
bool  enable_stp_init 
) [inline]

Solution of bound constrained optimization problems.

This subroutine solves bound constrained optimization problems by using the compact formula of the limited memory BFGS updates.

n is an integer variable. On entry n is the number of variables. On exit n is unchanged.

m is an integer variable. On entry m is the maximum number of variable metric corrections allowed in 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 of x. On exit l is unchanged.

u is a double precision array of dimension n. On entry u is the upper bound of 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, 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. 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.

ws, wy, sy, and wt are double precision working arrays used to store the following information defining the limited memory BFGS matrix: ws, of dimension n x m, stores S, the matrix of s-vectors; wy, of dimension n x m, stores Y, the matrix of y-vectors; sy, of dimension m x m, stores S'Y; ss, of dimension m x m, stores S'S; yy, of dimension m x m, stores Y'Y; wt, of dimension m x m, stores the Cholesky factorization of (theta*S'S+LD^(-1)L'); see eq. (2.26) in [3].

wn is a double precision working array of dimension 2m x 2m used to store 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]

snd is a double precision working array of dimension 2m x 2m used to store the lower triangular part of N = [Y' ZZ'Y L_a'+R_z'] [L_a +R_z S'AA'S ]

z(n),r(n),d(n),t(n),wa(8*m) are double precision working arrays. z is used at different times to store the Cauchy point and the Newton point.

sg(m),sgo(m),yg(m),ygo(m) are double precision working arrays.

index is an integer working array of dimension n. In subroutine freev, index is used to store the free and fixed variables at the Generalized Cauchy Point (GCP).

iwhere is an integer working array of dimension n used to record the status of the vector x for GCP computation. iwhere(i)=0 or -3 if x(i) is free and has bounds, 1 if x(i) is fixed at l(i), and l(i) .ne. u(i) 2 if x(i) is fixed at u(i), and u(i) .ne. l(i) 3 if x(i) is always fixed, i.e., u(i)=x(i)=l(i) -1 if x(i) is always free, i.e., no bounds on it.

indx2 is an integer working array of dimension n. Within subroutine cauchy, indx2 corresponds to the array iorder. In subroutine freev, a list of variables entering and leaving the free set is stored in indx2, and it is passed on to subroutine formk with this information.

task is a working string of characters of length 60 indicating the current job when entering and leaving 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.

isave is an integer working array of dimension 23.

dsave is a double precision working array of dimension 29.

Subprograms called

L-BFGS-B Library ... cauchy, subsm, lnsrlb, formk,

errclb, prn1lb, prn2lb, prn3lb, active, projgr,

freev, cmprlb, matupd, formt.

Minpack2 Library ... timer, dpmeps.

Linpack Library ... dcopy, ddot.

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: FORTRAN Subroutines for Large Scale Bound Constrained Optimization'' Tech. Report, NAM-11, EECS Department, Northwestern University, 1994.

[3] R. Byrd, J. Nocedal and R. Schnabel "Representations of Quasi-Newton Matrices and their use in Limited Memory Methods'', Mathematical Programming 63 (1994), no. 4, pp. 129-156.

(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.

References active(), ref::begin(), cauchy(), cmprlb(), dcopy(), dscal(), errclb(), formk(), formt(), freev(), floating_point_epsilon::get(), ref1::get1(), ref2::get2(), lnsrlb(), matupd(), max3(), prn1lb(), prn2lb(), prn3lb(), projgr(), SCITBX_ASSERT, const_ref::size(), subsm(), and timer().

Referenced by setulb().

void scitbx::lbfgsb::raw::matupd ( int const &  n,
int const &  m,
ref2< FloatType > const &  ws,
ref2< FloatType > const &  wy,
ref2< FloatType > const &  sy,
ref2< FloatType > const &  ss,
ref1< FloatType > const &  d,
ref1< FloatType > const &  r,
int &  itail,
int const &  iupdat,
int &  col,
int &  head,
FloatType &  theta,
FloatType const &  rr,
FloatType const &  dr,
FloatType const &  stp,
FloatType const &  dtd 
) [inline]

Updates matrices WS and WY, and forms the middle matrix in B.

This subroutine updates matrices WS and WY, and forms the middle matrix in B.

Subprograms called:

Linpack ... dcopy, ddot.

* *

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.

References ref::begin(), dcopy(), ref2::get1(), SCITBX_ASSERT, and const_ref::size().

Referenced by mainlb().

void scitbx::lbfgsb::raw::prn1lb ( int const &  n,
int const &  m,
ref1< FloatType > const &  l,
ref1< FloatType > const &  u,
ref1< FloatType > const &  x,
int const &  iprint,
int const &  ,
FloatType const &  epsmch 
) [inline]

Printing function 1.

This subroutine prints the input data, initial point, upper and lower bounds of each variable, machine precision, as well as the headings of the output.

* *

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.

References write_ref1().

Referenced by mainlb().

void scitbx::lbfgsb::raw::prn2lb ( int const &  ,
ref1< FloatType > const &  x,
FloatType const &  f,
ref1< FloatType > const &  g,
int const &  iprint,
int const &  ,
int const &  iter,
int const &  nfgv,
int const &  nact,
FloatType const &  sbgnrm,
int const &  nint,
std::string &  word,
int const &  iword,
int const &  iback,
FloatType const &  stp,
FloatType const &  xstep 
) [inline]

Printing function 2.

This subroutine prints out new information after a successful line search.

* *

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.

References write_ref1().

Referenced by mainlb().

void scitbx::lbfgsb::raw::prn3lb ( int const &  n,
ref1< FloatType > const &  x,
FloatType const &  f,
std::string const &  task,
int const &  iprint,
int const &  info,
int const &  ,
int const &  iter,
int const &  nfgv,
int const &  nintol,
int const &  nskip,
int const &  nact,
FloatType const &  sbgnrm,
FloatType const &  time,
int const &  nint,
std::string const &  word,
int const &  iback,
FloatType const &  stp,
FloatType const &  xstep,
int const &  k,
FloatType const &  cachyt,
FloatType const &  sbtime,
FloatType const &  lnscht 
) [inline]

Printing function 3.

This subroutine prints out information when either a built-in convergence test is satisfied or when an error message is generated.

* *

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.

References write_ref1().

Referenced by mainlb().

void scitbx::lbfgsb::raw::projgr ( int const &  n,
ref1< FloatType > const &  l,
ref1< FloatType > const &  u,
ref1< int > const &  nbd,
ref1< FloatType > const &  x,
ref1< FloatType > const &  g,
FloatType &  sbgnrm 
) [inline]

This subroutine computes the infinity norm of the projected gradient.

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.

Referenced by mainlb().

void scitbx::lbfgsb::raw::setulb ( int const &  n,
int const &  m,
ref1< FloatType > const &  x,
ref1< FloatType > const &  l,
ref1< FloatType > const &  u,
ref1< int > const &  nbd,
FloatType &  f,
ref1< FloatType > const &  g,
FloatType const &  factr,
FloatType const &  pgtol,
ref1< FloatType > const &  wa,
ref1< int > const &  iwa,
std::string &  task,
int const &  iprint,
std::string &  csave,
ref1< bool > const &  lsave,
ref1< int > const &  isave,
ref1< FloatType > const &  dsave,
bool  enable_stp_init = false 
) [inline]

Main L-BFGS-B interface function.

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.

References ref1::get1(), ref1::get2(), mainlb(), SCITBX_ASSERT, and const_ref::size().

Referenced by minimizer::process().

void scitbx::lbfgsb::raw::subsm ( int const &  ,
int const &  m,
int const &  nsub,
ref1< int > const &  ind,
ref1< FloatType > const &  l,
ref1< FloatType > const &  u,
ref1< int > const &  nbd,
ref1< FloatType > const &  x,
ref1< FloatType > const &  d,
ref2< FloatType > const &  ws,
ref2< FloatType > const &  wy,
FloatType const &  theta,
int const &  col,
int const &  head,
int &  iword,
ref1< FloatType > const &  wv,
ref2< FloatType > const &  wn,
int const &  iprint,
int &  info 
) [inline]

Approximate solution of the subspace problem.

Given xcp, l, u, r, an index set that specifies the active set at xcp, and an l-BFGS matrix B (in terms of WY, WS, SY, WT, head, col, and theta), this subroutine computes an approximate solution of the subspace problem

(P) min Q(x) = r'(x-xcp) + 1/2 (x-xcp)' B (x-xcp)

subject to l<=x<=u x_i=xcp_i for all i in A(xcp)

along the subspace unconstrained Newton direction

d = -(Z'BZ)^(-1) r.

The formula for the Newton direction, given the L-BFGS matrix and the Sherman-Morrison formula, is

d = (1/theta)r + (1/theta*2) Z'WK^(-1)W'Z r.

where K = [-D -Y'ZZ'Y/theta L_a'-R_z' ] [L_a -R_z theta*S'AA'S ]

Note that this procedure for computing d differs from that described in [1]. One can show that the matrix K is equal to the matrix M^[-1]N in that paper.

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.

nsub is an integer variable. On entry nsub is the number of free variables. On exit nsub is unchanged.

ind is an integer array of dimension nsub. On entry ind specifies the coordinate indices of free variables. On exit ind is unchanged.

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

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

nbd is a 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.

x is a double precision array of dimension n. On entry x specifies the Cauchy point xcp. On exit x(i) is the minimizer of Q over the subspace of free variables.

d is a double precision array of dimension n. On entry d is the reduced gradient of Q at xcp. On exit d is the Newton direction of Q.

ws and wy are double precision arrays; theta is a double precision variable; col is an integer variable; head is an integer variable. On entry they store the information defining the limited memory BFGS matrix: ws(n,m) stores S, a set of s-vectors; wy(n,m) stores Y, a set of y-vectors; theta is the scaling factor specifying B_0 = theta I; col is the number of variable metric corrections stored; head is the location of the 1st s- (or y-) vector in S (or Y). On exit they are unchanged.

iword is an integer variable. On entry iword is unspecified. On exit iword specifies the status of the subspace solution. iword = 0 if the solution is in the box, 1 if some bound is encountered.

wv is a double precision working array of dimension 2m.

wn is a double precision array of dimension 2m x 2m. On entry the upper triangle of wn stores 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] On exit wn is unchanged.

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.

info is an integer variable. On entry info is unspecified. On exit info = 0 for normal return, = nonzero for abnormal return when the matrix K is ill-conditioned.

Subprograms called:

Linpack dtrsl.

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.

* *

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.

References dtrsl(), ref1::get1(), ref2::get2(), SCITBX_ASSERT, const_ref::size(), and write_ref1().

Referenced by mainlb().


Generated on Thu Jun 19 15:35:08 2014 for cctbx by  doxygen 1.5.6