"""
    setulb(n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave)
This function wraps L-BFGS-B's fortran subroutine `setulb` which 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`.

Note, `nmax` is the dimension of the largest problem to be solved, `mmax` is the maximum number of limited memory corrections.

# Arguments
- `n::Ref{Cint}`: the dimension of the problem.
- `m::Ref{Cint}`: the maximum number of variable metric corrections used to define the limited memory matrix.
- `x::Vector{Cdouble}`: an approximation to the solution, of length n.
- `l::Vector{Cdouble}`: the lower bound on x, of length n.
- `u::Vector{Cdouble}`: the upper bound on x, of length n.
- `nbd::Vector{Cint}`: represents the type of bounds imposed on the variables, and must be specified as follows:
    1. nbd[i]=0 if x[i] is unbounded;
    2. nbd[i]=1 if x(i) has only a lower bound;
    3. nbd[i]=2 if x[i] has both lower and upper bounds;
    4. nbd[i]=3 if x[i] has only an upper bound.
- `f::Ref{Cdouble}`: the value of the function at x(on first entry f is unspecified).
- `g::Vector{Cdouble}`: the value of the gradient at x(on first entry g is unspecified).
- `factr::Ref{Cdouble}`: 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.e12 for low accuracy;
    * 1.e7 for moderate accuracy;
    * 1.e1 for extremely high accuracy.
- `pgtol::Ref{Cdouble}`: the iteration will stop when max{|proj g_i | i = 1, ..., n} <= pgtol where pg_i is the ith component of the projected gradient.
- `wa::Vector{Cdouble}`: the working array of length (2mmax + 5)nmax + 12mmax^2 + 12mmax.
- `iwa::Vector{Cint}`: the integer working array of length 3nmax.
- `task::Vector{Cuchar}`: a working string of characters of length 60 indicating the current job when entering and quitting this function.
- `iprint::Ref{Cint}`: it controls the frequency and type of output generated:
    1. iprint<0 no output is generated;
    2. iprint=0 print only one line at the last iteration;
    3. 0<iprint<99 print also f and |proj g| every iprint iterations;
    4. iprint=99 print details of every iteration except n-vectors;
    5. iprint=100 print also the changes of active set and final x;
    6. iprint>100 print details of every iteration including x and g.
    7. when iprint > 0, the file iterate.dat will be created to summarize the iteration.
- `csave::Vector{Cuchar}`: a working string of characters of length 60.
- `lsave::Vector{Bool}`: a logical working array of dimension 4. if `task` = NEW_X, the following information is available:
    1. if lsave[1] == true then the initial X has been replaced by its projection in the feasible set;
    2. if lsave[2] == true then the problem is constrained;
    3. if lsave[3] == true then each variable has upper and lower bounds.
- `isave::Vector{Cint}`: an integer working array of dimension 44. if `task` = NEW_X, the following information is available:
    1. isave[22] = the total number of intervals explored in the search of Cauchy points;
    2. isave[26] = the total number of skipped BFGS updates before the current iteration;
    3. isave[30] = the number of current iteration;
    4. isave[31] = the total number of BFGS updates prior the current iteration;
    5. isave[33] = the number of intervals explored in the search of Cauchy point in the current iteration;
    6. isave[34] = the total number of function and gradient evaluations;
    7. isave[36] = the number of function value or gradient evaluations in the current iteration;
    8. if isave[37] = 0  then the subspace argmin is within the box;
    9. if isave[37] = 1  then the subspace argmin is beyond the box;
    10. isave[38] = the number of free variables in the current iteration;
    11. isave[39] = the number of active constraints in the current iteration;
    12. n + 1 - isave[40] = the number of variables leaving the set of active constraints in the current iteration;
    13. isave[41] = the number of variables entering the set of active constraints in the current iteration.
- `dsave::Vector{Cdouble}`: a double precision working array of dimension 29. if `task` = NEW_X, the following information is available:
    1. dsave[1] = current 'theta' in the BFGS matrix;
    2. dsave[2] = f(x) in the previous iteration;
    3. dsave[3] = factr*epsmch;
    4. dsave[4] = 2-norm of the line search direction vector;
    5. dsave[5] = the machine precision epsmch generated by the code;
    6. dsave[7] = the accumulated time spent on searching for Cauchy points;
    7. dsave[8] = the accumulated time spent on subspace minimization;
    8. dsave[9] = the accumulated time spent on line search;
    9. dsave[11] = the slope of the line search function at the current point of line search;
    10. dsave[12] = the maximum relative step length imposed in line search;
    11. dsave[13] = the infinity norm of the projected gradient;
    12. dsave[14] = the relative step length in the line search;
    13. dsave[15] = the slope of the line search function at the starting point of the line search;
    14. dsave[16] = the square of the 2-norm of the line search direction vector.
"""
function setulb(n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave)
    global liblbfgsb
    hdl = Libdl.dlopen_e(liblbfgsb)
    @assert hdl != C_NULL "Could not open $liblbfgsb"
    setulb_ = Libdl.dlsym_e(hdl, :setulb_)
    @assert setulb_ != C_NULL "Could not find `setulb` within $liblbfgsb"
    ccall(setulb_, Cvoid, (Ref{Cint}, Ref{Cint}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
          Ref{Cint}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
          Ref{Cint}, Ref{Cuchar}, Ref{Cint}, Ref{Cuchar}, Ref{Cint}, Ref{Cint}, Ref{Cdouble}, Csize_t, Csize_t),
          n, m, x, l, u, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave, 60, 60)
end

"""
    timer(x)
The double precision cpu timing subroutine in the L-BFGS-B code.
"""
function timer(x)
    global liblbfgsb
    hdl = Libdl.dlopen_e(liblbfgsb)
    @assert hdl != C_NULL "Could not open $liblbfgsb"
    timer_ = Libdl.dlsym_e(hdl, :timer_)
    @assert timer_ != C_NULL "Could not find `timer` within $liblbfgsb"
    ccall(timer_, Cvoid, (Ref{Cdouble},), x)
end
