GALAHAD QPB package#
purpose#
The qpb package uses a
primal-dual interior-point method to solve a given
non-convex quadratic program.
The aim is to minimize the quadratic objective function
If the matrix \(H\) is positive semi-definite, a global solution is found. However, if \(H\) is indefinite, the procedure may find a (weak second-order) critical point that is not the global solution to the given problem.
See Section 4 of $GALAHAD/doc/qpb.pdf for additional details.
terminolgy#
Any required solution \(x\) necessarily satisfies the primal optimality conditions
method#
Primal-dual interior point methods iterate towards a point that satisfies these conditions by ultimately aiming to satisfy (1), (3) and (5), while ensuring that (2) and (4) are satisfied as strict inequalities at each stage. Appropriate norms of the amounts by which (1), (3) and (5) fail to be satisfied are known as the primal and dual infeasibility, and the violation of complementary slackness, respectively. The fact that (2) and (4) are satisfied as strict inequalities gives such methods their other title, namely interior-point methods.
The problem is solved in two phases. The goal of the first
“initial feasible point” phase is
to find a strictly interior point which is primal feasible, that is that
(1) is satisfied. The package LSQP is used for this
purpose, and offers the options of either accepting the first
strictly feasible point found, or preferably of aiming for the
so-called “analytic center” of the feasible region.
Having found such a suitable initial feasible point, the second “optimality”
phase ensures that (1) remains satisfied while iterating to
satisfy dual feasibility (3) and complementary slackness (5).
The optimality phase proceeds by approximately minimizing a
sequence of barrier functions
Each of the barrier subproblems is solved using a trust-region method. Such a method generates a trial correction step \(\Delta (x, c)\) to the current iterate \((x, c)\) by replacing the nonlinear barrier function locally by a suitable quadratic model, and approximately minimizing this model in the intersection of (1) and a trust region \(\|\Delta (x, c)\| \leq \Delta\) for some appropriate strictly positive trust-region radius \(\Delta\) and norm \(\| \cdot \|\). The step is accepted/rejected and the radius adjusted on the basis of how accurately the model reproduces the value of barrier function at the trial step. If the step proves to be unacceptable, a linesearch is performed along the step to obtain an acceptable new iterate. In practice, the natural primal “Newton” model of the barrier function is frequently less successful than an alternative primal-dual model, and consequently the primal-dual model is usually to be preferred.
Once a barrier subproblem has been solved, extrapolation based on values and derivatives encountered on the central path is optionally used to determine a good starting point for the next subproblem. Traditional Taylor-series extrapolation has been superceded by more accurate Puiseux-series methods as these are particularly suited to deal with degeneracy.
The trust-region subproblem is approximately solved using the combined
conjugate-gradient/Lanczos method implemented in the package GLTR.
Such a method requires a suitable preconditioner, and in our case, the
only flexibility we have is in approximating the model of the
Hessian. Although using a fixed form of preconditioning is sometimes
effective, we have provided the option of an automatic choice, that aims
to balance the cost of applying the preconditioner against the needs for
an accurate solution of the trust-region subproblem. The preconditioner
is applied using the matrix factorization package SBLS, but options
at this stage are to factorize the preconditioner as a whole (the
so-called “augmented system” approach), or to perform a block
elimination first (the “Schur-complement” approach). The latter is
usually to be prefered when a (non-singular) diagonal preconditioner is
used, but may be inefficient if any of the columns of \(A\) is too dense.
In order to make the solution as efficient as possible, the
variables and constraints are reordered internally
by the package QPP prior to solution.
In particular, fixed variables, and
free (unbounded on both sides) constraints are temporarily removed.
reference#
The method is described in detail in
A. R. Conn, N. I. M. Gould, D. Orban and Ph. L. Toint, ``A primal-dual trust-region algorithm for minimizing a non-convex function subject to general inequality and linear equality constraints’’. Mathematical Programming 87 (1999) 215-249.
matrix storage#
unsymmetric storage#
The unsymmetric \(m\) by \(n\) matrix \(A\) may be presented and stored in a variety of convenient input formats.
Dense storage format: The matrix \(A\) is stored as a compact dense matrix by rows, that is, the values of the entries of each row in turn are stored in order within an appropriate real one-dimensional array. In this case, component \(n \ast i + j\) of the storage array A_val will hold the value \(A_{ij}\) for \(1 \leq i \leq m\), \(1 \leq j \leq n\). The string A_type = ‘dense’ should be specified.
Dense by columns storage format: The matrix \(A\) is stored as a compact dense matrix by columns, that is, the values of the entries of each column in turn are stored in order within an appropriate real one-dimensional array. In this case, component \(m \ast j + i\) of the storage array A_val will hold the value \(A_{ij}\) for \(1 \leq i \leq m\), \(1 \leq j \leq n\). The string A_type = ‘dense_by_columns’ should be specified.
Sparse co-ordinate storage format: Only the nonzero entries of the matrices are stored. For the \(l\)-th entry, \(1 \leq l \leq ne\), of \(A\), its row index i, column index j and value \(A_{ij}\), \(1 \leq i \leq m\), \(1 \leq j \leq n\), are stored as the \(l\)-th components of the integer arrays A_row and A_col and real array A_val, respectively, while the number of nonzeros is recorded as A_ne = \(ne\). The string A_type = ‘coordinate’should be specified.
Sparse row-wise storage format: Again only the nonzero entries are stored, but this time they are ordered so that those in row i appear directly before those in row i+1. For the i-th row of \(A\) the i-th component of the integer array A_ptr holds the position of the first entry in this row, while A_ptr(m+1) holds the total number of entries plus one. The column indices j, \(1 \leq j \leq n\), and values \(A_{ij}\) of the nonzero entries in the i-th row are stored in components l = A_ptr(i), \(\ldots\), A_ptr(i+1)-1, \(1 \leq i \leq m\), of the integer array A_col, and real array A_val, respectively. For sparse matrices, this scheme almost always requires less storage than its predecessor. The string A_type = ‘sparse_by_rows’ should be specified.
Sparse column-wise storage format: Once again only the nonzero entries are stored, but this time they are ordered so that those in column j appear directly before those in column j+1. For the j-th column of \(A\) the j-th component of the integer array A_ptr holds the position of the first entry in this column, while A_ptr(n+1) holds the total number of entries plus one. The row indices i, \(1 \leq i \leq m\), and values \(A_{ij}\) of the nonzero entries in the j-th columnsare stored in components l = A_ptr(j), \(\ldots\), A_ptr(j+1)-1, \(1 \leq j \leq n\), of the integer array A_row, and real array A_val, respectively. As before, for sparse matrices, this scheme almost always requires less storage than the co-ordinate format. The string A_type = ‘sparse_by_columns’ should be specified.
symmetric storage#
The symmetric \(n\) by \(n\) matrix \(H\) may also be presented and stored in a variety of formats. But crucially symmetry is exploited by only storing values from the lower triangular part (i.e, those entries that lie on or below the leading diagonal).
Dense storage format: The matrix \(H\) is stored as a compact dense matrix by rows, that is, the values of the entries of each row in turn are stored in order within an appropriate real one-dimensional array. Since \(H\) is symmetric, only the lower triangular part (that is the part \(H_{ij}\) for \(1 \leq j \leq i \leq n\)) need be held. In this case the lower triangle should be stored by rows, that is component \((i-1) * i / 2 + j\) of the storage array H_val will hold the value \(H_{ij}\) (and, by symmetry, \(H_{ji}\)) for \(1 \leq j \leq i \leq n\). The string H_type = ‘dense’ should be specified.
Sparse co-ordinate storage format: Only the nonzero entries of the matrices are stored. For the \(l\)-th entry, \(1 \leq l \leq ne\), of \(H\), its row index i, column index j and value \(H_{ij}\), \(1 \leq j \leq i \leq n\), are stored as the \(l\)-th components of the integer arrays H_row and H_col and real array H_val, respectively, while the number of nonzeros is recorded as H_ne = \(ne\). Note that only the entries in the lower triangle should be stored. The string H_type = ‘coordinate’ should be specified.
Sparse row-wise storage format: Again only the nonzero entries are stored, but this time they are ordered so that those in row i appear directly before those in row i+1. For the i-th row of \(H\) the i-th component of the integer array H_ptr holds the position of the first entry in this row, while H_ptr(n+1) holds the total number of entries plus one. The column indices j, \(1 \leq j \leq i\), and values \(H_{ij}\) of the entries in the i-th row are stored in components l = H_ptr(i), …, H_ptr(i+1)-1 of the integer array H_col, and real array H_val, respectively. Note that as before only the entries in the lower triangle should be stored. For sparse matrices, this scheme almost always requires less storage than its predecessor. The string H_type = ‘sparse_by_rows’ should be specified.
Diagonal storage format: If \(H\) is diagonal (i.e., \(H_{ij} = 0\) for all \(1 \leq i \neq j \leq n\)) only the diagonals entries \(H_{ii}\), \(1 \leq i \leq n\) need be stored, and the first n components of the array H_val may be used for the purpose. The string H_type = ‘diagonal’ should be specified.
Multiples of the identity storage format: If \(H\) is a multiple of the identity matrix, (i.e., \(H = \alpha I\) where \(I\) is the n by n identity matrix and \(\alpha\) is a scalar), it suffices to store \(\alpha\) as the first component of H_val. The string H_type = ‘scaled_identity’ should be specified.
The identity matrix format: If \(H\) is the identity matrix, no values need be stored. The string H_type = ‘identity’ should be specified.
The zero matrix format: The same is true if \(H\) is the zero matrix, but now the string H_type = ‘zero’ or ‘none’ should be specified.
introduction to function calls#
To solve a given problem, functions from the qpb package must be called in the following order:
qpb_initialize - provide default control parameters and set up initial data structures
qpb_read_specfile (optional) - override control values by reading replacement values from a file
qpb_import - set up problem data structures and fixed values
qpb_reset_control (optional) - possibly change control parameters if a sequence of problems are being solved
qpb_solve_qp - solve the quadratic program
qpb_information (optional) - recover information about the solution and solution process
qpb_terminate - deallocate data structures
See the examples section for illustrations of use.
parametric real type T#
Below, the symbol T refers to a parametric real type that may be Float32 (single precision), Float64 (double precision) or, if supported, Float128 (quadruple precision).
callable functions#
function qpb_initialize(T, data, control, status)
Set default control values and initialize private data
Parameters:
data |
holds private internal data |
control |
is a structure containing control information (see qpb_control_type) |
status |
is a scalar variable of type Int32 that gives the exit status from the package. Possible values are (currently):
|
function qpb_read_specfile(T, control, specfile)
Read the content of a specification file, and assign values associated with given keywords to the corresponding control parameters. An in-depth discussion of specification files is available, and a detailed list of keywords with associated default values is provided in $GALAHAD/src/qpb/QPB.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/qpb.pdf for a list of how these keywords relate to the components of the control structure.
Parameters:
control |
is a structure containing control information (see qpb_control_type) |
specfile |
is a one-dimensional array of type Vararg{Cchar} that must give the name of the specification file |
function qpb_import(T, control, data, status, n, m, H_type, H_ne, H_row, H_col, H_ptr, A_type, A_ne, A_row, A_col, A_ptr)
Import problem data into internal storage prior to solution.
Parameters:
control |
is a structure whose members provide control parameters for the remaining procedures (see qpb_control_type) |
data |
holds private internal data |
status |
is a scalar variable of type Int32 that gives the exit status from the package. Possible values are:
|
n |
is a scalar variable of type Int32 that holds the number of variables. |
m |
is a scalar variable of type Int32 that holds the number of general linear constraints. |
H_type |
is a one-dimensional array of type Vararg{Cchar} that specifies the symmetric storage scheme used for the Hessian, \(H\). It should be one of ‘coordinate’, ‘sparse_by_rows’, ‘dense’, ‘diagonal’, ‘scaled_identity’, ‘identity’, ‘zero’ or ‘none’, the latter pair if \(H=0\); lower or upper case variants are allowed. |
H_ne |
is a scalar variable of type Int32 that holds the number of entries in the lower triangular part of \(H\) in the sparse co-ordinate storage scheme. It need not be set for any of the other schemes. |
H_row |
is a one-dimensional array of size H_ne and type Int32 that holds the row indices of the lower triangular part of \(H\) in the sparse co-ordinate storage scheme. It need not be set for any of the other three schemes, and in this case can be C_NULL. |
H_col |
is a one-dimensional array of size H_ne and type Int32 that holds the column indices of the lower triangular part of \(H\) in either the sparse co-ordinate, or the sparse row-wise storage scheme. It need not be set when the dense, diagonal or (scaled) identity storage schemes are used, and in this case can be C_NULL. |
H_ptr |
is a one-dimensional array of size n+1 and type Int32 that holds the starting position of each row of the lower triangular part of \(H\), as well as the total number of entries, in the sparse row-wise storage scheme. It need not be set when the other schemes are used, and in this case can be C_NULL. |
A_type |
is a one-dimensional array of type Vararg{Cchar} that specifies the unsymmetric storage scheme used for the constraint Jacobian, \(A\). It should be one of ‘coordinate’, ‘sparse_by_rows’ or ‘dense; lower or upper case variants are allowed. |
A_ne |
is a scalar variable of type Int32 that holds the number of entries in \(A\) in the sparse co-ordinate storage scheme. It need not be set for any of the other schemes. |
A_row |
is a one-dimensional array of size A_ne and type Int32 that holds the row indices of \(A\) in the sparse co-ordinate storage scheme. It need not be set for any of the other schemes, and in this case can be C_NULL. |
A_col |
is a one-dimensional array of size A_ne and type Int32 that holds the column indices of \(A\) in either the sparse co-ordinate, or the sparse row-wise storage scheme. It need not be set when the dense or diagonal storage schemes are used, and in this case can be C_NULL. |
A_ptr |
is a one-dimensional array of size n+1 and type Int32 that holds the starting position of each row of \(A\), as well as the total number of entries, in the sparse row-wise storage scheme. It need not be set when the other schemes are used, and in this case can be C_NULL. |
function qpb_reset_control(T, control, data, status)
Reset control parameters after import if required.
Parameters:
control |
is a structure whose members provide control parameters for the remaining procedures (see qpb_control_type) |
data |
holds private internal data |
status |
is a scalar variable of type Int32 that gives the exit status from the package. Possible values are:
|
function qpb_solve_qp(T, data, status, n, m, h_ne, H_val, g, f, a_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z, x_stat, c_stat)
Solve the quadratic program when the Hessian \(H\) is available.
Parameters:
data |
holds private internal data |
status |
is a scalar variable of type Int32 that gives the entry and exit status from the package. Possible exit values are:
|
n |
is a scalar variable of type Int32 that holds the number of variables |
m |
is a scalar variable of type Int32 that holds the number of general linear constraints. |
h_ne |
is a scalar variable of type Int32 that holds the number of entries in the lower triangular part of the Hessian matrix \(H\). |
H_val |
is a one-dimensional array of size h_ne and type T that holds the values of the entries of the lower triangular part of the Hessian matrix \(H\) in any of the available storage schemes. |
g |
is a one-dimensional array of size n and type T that holds the linear term \(g\) of the objective function. The j-th component of |
f |
is a scalar of type T that holds the constant term \(f\) of the objective function. |
a_ne |
is a scalar variable of type Int32 that holds the number of entries in the constraint Jacobian matrix \(A\). |
A_val |
is a one-dimensional array of size a_ne and type T that holds the values of the entries of the constraint Jacobian matrix \(A\) in any of the available storage schemes. |
c_l |
is a one-dimensional array of size m and type T that holds the lower bounds \(c^l\) on the constraints \(A x\). The i-th component of |
c_u |
is a one-dimensional array of size m and type T that holds the upper bounds \(c^l\) on the constraints \(A x\). The i-th component of |
x_l |
is a one-dimensional array of size n and type T that holds the lower bounds \(x^l\) on the variables \(x\). The j-th component of |
x_u |
is a one-dimensional array of size n and type T that holds the upper bounds \(x^l\) on the variables \(x\). The j-th component of |
x |
is a one-dimensional array of size n and type T that holds the values \(x\) of the optimization variables. The j-th component of |
c |
is a one-dimensional array of size m and type T that holds the residual \(c(x)\). The i-th component of |
y |
is a one-dimensional array of size n and type T that holds the values \(y\) of the Lagrange multipliers for the general linear constraints. The j-th component of |
z |
is a one-dimensional array of size n and type T that holds the values \(z\) of the dual variables. The j-th component of |
x_stat |
is a one-dimensional array of size n and type Int32 that gives the optimal status of the problem variables. If x_stat(j) is negative, the variable \(x_j\) most likely lies on its lower bound, if it is positive, it lies on its upper bound, and if it is zero, it lies between its bounds. |
c_stat |
is a one-dimensional array of size m and type Int32 that gives the optimal status of the general linear constraints. If c_stat(i) is negative, the constraint value \(a_i^Tx\) most likely lies on its lower bound, if it is positive, it lies on its upper bound, and if it is zero, it lies between its bounds. |
function qpb_information(T, data, inform, status)
Provides output information
Parameters:
data |
holds private internal data |
inform |
is a structure containing output information (see qpb_inform_type) |
status |
is a scalar variable of type Int32 that gives the exit status from the package. Possible values are (currently):
|
function qpb_terminate(T, data, control, inform)
Deallocate all internal private storage
Parameters:
data |
holds private internal data |
control |
is a structure containing control information (see qpb_control_type) |
inform |
is a structure containing output information (see qpb_inform_type) |
available structures#
qpb_control_type structure#
struct qpb_control_type{T} f_indexing::Bool error::Int32 out::Int32 print_level::Int32 start_print::Int32 stop_print::Int32 maxit::Int32 itref_max::Int32 cg_maxit::Int32 indicator_type::Int32 restore_problem::Int32 extrapolate::Int32 path_history::Int32 factor::Int32 max_col::Int32 indmin::Int32 valmin::Int32 infeas_max::Int32 precon::Int32 nsemib::Int32 path_derivatives::Int32 fit_order::Int32 sif_file_device::Int32 infinity::T stop_p::T stop_d::T stop_c::T theta_d::T theta_c::T beta::T prfeas::T dufeas::T muzero::T reduce_infeas::T obj_unbounded::T pivot_tol::T pivot_tol_for_dependencies::T zero_pivot::T identical_bounds_tol::T inner_stop_relative::T inner_stop_absolute::T initial_radius::T mu_min::T inner_fraction_opt::T indicator_tol_p::T indicator_tol_pd::T indicator_tol_tapia::T cpu_time_limit::T clock_time_limit::T remove_dependencies::Bool treat_zero_bounds_as_general::Bool center::Bool primal::Bool puiseux::Bool feasol::Bool array_syntax_worse_than_do_loop::Bool space_critical::Bool deallocate_error_fatal::Bool generate_sif_file::Bool sif_file_name::NTuple{31,Cchar} prefix::NTuple{31,Cchar} lsqp_control::lsqp_control_type{T} fdc_control::fdc_control_type{T} sbls_control::sbls_control_type{T} gltr_control::gltr_control_type{T} fit_control::fit_control_type
detailed documentation#
control derived type as a Julia structure
components#
Bool f_indexing
use C or Fortran sparse matrix indexing
Int32 error
error and warning diagnostics occur on stream error
Int32 out
general output occurs on stream out
Int32 print_level
the level of output required is specified by print_level
Int32 start_print
any printing will start on this iteration
Int32 stop_print
any printing will stop on this iteration
Int32 maxit
at most maxit inner iterations are allowed
Int32 itref_max
the maximum number of iterative refinements allowed
Int32 cg_maxit
the maximum number of CG iterations allowed. If cg_maxit < 0, this number will be reset to the dimension of the system + 1
Int32 indicator_type
specifies the type of indicator function used. Pssible values are
1 primal indicator: constraint active <=> distance to nearest bound <= .indicator_p_tol
2 primal-dual indicator: constraint active <=> distance to nearest bound <= .indicator_tol_pd * size of corresponding multiplier
3 primal-dual indicator: constraint active <=> distance to nearest bound <= .indicator_tol_tapia * distance to same bound at previous iteration
Int32 restore_problem
indicate whether and how much of the input problem should be restored on output. Possible values are
0 nothing restored
1 scalar and vector parameters
2 all parameters
Int32 extrapolate
should extrapolation be used to track the central path? Possible values
0 never
1 after the final major iteration
2 at each major iteration
Int32 path_history
the maximum number of previous path points to use when fitting the data
Int32 factor
the factorization to be used. Possible values are
0 automatic
1 Schur-complement factorization
2 augmented-system factorization
Int32 max_col
the maximum number of nonzeros in a column of A which is permitted with the Schur-complement factorization
Int32 indmin
an initial guess as to the integer workspace required by SBLS
Int32 valmin
an initial guess as to the real workspace required by SBLS
Int32 infeas_max
the number of iterations for which the overall infeasibility of the problem is not reduced by at least a factor .reduce_infeas before the problem is flagged as infeasible (see reduce_infeas)
Int32 precon
the preconditioner to be used for the CG is defined by precon. Possible values are
0 automatic
1 no preconditioner, i.e, the identity within full factorization
2 full factorization
3 band within full factorization
4 diagonal using the barrier terms within full factorization
Int32 nsemib
the semi-bandwidth of a band preconditioner, if appropriate
Int32 path_derivatives
the maximum order of path derivative to use
Int32 fit_order
the order of (Puiseux) series to fit to the path data: <=0 to fit all data
Int32 sif_file_device
specifies the unit number to write generated SIF file describing the current problem
T infinity
any bound larger than infinity in modulus will be regarded as infinite
T stop_p
the required accuracy for the primal infeasibility
T stop_d
the required accuracy for the dual infeasibility
T stop_c
the required accuracy for the complementarity
T theta_d
tolerances used to terminate the inner iteration (for given mu): dual feasibility <= MAX( theta_d * mu ** beta, 0.99 * stop_d ) complementarity <= MAX( theta_c * mu ** beta, 0.99 * stop_d )
T theta_c
see theta_d
T beta
see theta_d
T prfeas
initial primal variables will not be closer than prfeas from their bound
T dufeas
initial dual variables will not be closer than dufeas from their bounds
T muzero
the initial value of the barrier parameter. If muzero is not positive, it will be reset to an appropriate value
T reduce_infeas
if the overall infeasibility of the problem is not reduced by at least a factor reduce_infeas over .infeas_max iterations, the problem is flagged as infeasible (see infeas_max)
T obj_unbounded
if the objective function value is smaller than obj_unbounded, it will be flagged as unbounded from below.
T pivot_tol
the threshold pivot used by the matrix factorization. See the documentation for SBLS for details
T pivot_tol_for_dependencies
the threshold pivot used by the matrix factorization when attempting to detect linearly dependent constraints. See the documentation for FDC for details
T zero_pivot
any pivots smaller than zero_pivot in absolute value will be regarded to zero when attempting to detect linearly dependent constraints
T identical_bounds_tol
any pair of constraint bounds (c_l,c_u) or (x_l,x_u) that are closer than identical_bounds_tol will be reset to the average of their values
T inner_stop_relative
the search direction is considered as an acceptable approximation to the minimizer of the model if the gradient of the model in the preconditioning(inverse) norm is less than max( inner_stop_relative * initial preconditioning(inverse) gradient norm, inner_stop_absolute )
T inner_stop_absolute
see inner_stop_relative
T initial_radius
the initial trust-region radius
T mu_min
start terminal extrapolation when mu reaches mu_min
T inner_fraction_opt
a search direction which gives at least inner_fraction_opt times the optimal model decrease will be found
T indicator_tol_p
if .indicator_type = 1, a constraint/bound will be deemed to be active <=> distance to nearest bound <= .indicator_p_tol
T indicator_tol_pd
if .indicator_type = 2, a constraint/bound will be deemed to be active <=> distance to nearest bound <= .indicator_tol_pd * size of corresponding multiplier
T indicator_tol_tapia
if .indicator_type = 3, a constraint/bound will be deemed to be active <=> distance to nearest bound <= .indicator_tol_tapia * distance to same bound at previous iteration
T cpu_time_limit
the maximum CPU time allowed (-ve means infinite)
T clock_time_limit
the maximum elapsed clock time allowed (-ve means infinite)
Bool remove_dependencies
the equality constraints will be preprocessed to remove any linear dependencies if true
Bool treat_zero_bounds_as_general
any problem bound with the value zero will be treated as if it were a general value if true
Bool center
if .center is true, the algorithm will use the analytic center of the feasible set as its initial feasible point. Otherwise, a feasible point as close as possible to the initial point will be used. We recommend using the analytic center
Bool primal
if .primal, is true, a primal barrier method will be used in place of t primal-dual method
Bool puiseux
If extrapolation is to be used, decide between Puiseux and Taylor series.
Bool feasol
if .feasol is true, the final solution obtained will be perturbed so that variables close to their bounds are moved onto these bounds
Bool array_syntax_worse_than_do_loop
if .array_syntax_worse_than_do_loop is true, f77-style do loops will be used rather than f90-style array syntax for vector operations
Bool space_critical
if .space_critical true, every effort will be made to use as little space as possible. This may result in longer computation time
Bool deallocate_error_fatal
if .deallocate_error_fatal is true, any array/pointer deallocation error will terminate execution. Otherwise, computation will continue
Bool generate_sif_file
if .generate_sif_file is .true. if a SIF file describing the current problem is to be generated
NTuple{31,Cchar} sif_file_name
name of generated SIF file containing input problem
NTuple{31,Cchar} prefix
all output lines will be prefixed by .prefix(2:LEN(TRIM(.prefix))-1) where .prefix contains the required string enclosed in quotes, e.g. “string” or ‘string’
struct lsqp_control_type lsqp_control
control parameters for LSQP
struct fdc_control_type fdc_control
control parameters for FDC
struct sbls_control_type sbls_control
control parameters for SBLS
struct gltr_control_type gltr_control
control parameters for GLTR
struct fit_control_type fit_control
control parameters for FIT
qpb_time_type structure#
struct qpb_time_type{T} total::T preprocess::T find_dependent::T analyse::T factorize::T solve::T phase1_total::T phase1_analyse::T phase1_factorize::T phase1_solve::T clock_total::T clock_preprocess::T clock_find_dependent::T clock_analyse::T clock_factorize::T clock_solve::T clock_phase1_total::T clock_phase1_analyse::T clock_phase1_factorize::T clock_phase1_solve::T
detailed documentation#
time derived type as a Julia structure
components#
T total
the total CPU time spent in the package
T preprocess
the CPU time spent preprocessing the problem
T find_dependent
the CPU time spent detecting linear dependencies
T analyse
the CPU time spent analysing the required matrices prior to factorizatio
T factorize
the CPU time spent factorizing the required matrices
T solve
the CPU time spent computing the search direction
T phase1_total
the total CPU time spent in the initial-point phase of the package
T phase1_analyse
the CPU time spent analysing the required matrices prior to factorizatio in the inital-point phase
T phase1_factorize
the CPU time spent factorizing the required matrices in the inital-point phase
T phase1_solve
the CPU time spent computing the search direction in the inital-point ph
T clock_total
the total clock time spent in the package
T clock_preprocess
the clock time spent preprocessing the problem
T clock_find_dependent
the clock time spent detecting linear dependencies
T clock_analyse
the clock time spent analysing the required matrices prior to factorizat
T clock_factorize
the clock time spent factorizing the required matrices
T clock_solve
the clock time spent computing the search direction
T clock_phase1_total
the total clock time spent in the initial-point phase of the package
T clock_phase1_analyse
the clock time spent analysing the required matrices prior to factorizat in the inital-point phase
T clock_phase1_factorize
the clock time spent factorizing the required matrices in the inital-poi phase
T clock_phase1_solve
the clock time spent computing the search direction in the inital-point
qpb_inform_type structure#
struct qpb_inform_type{T} status::Int32 alloc_status::Int32 bad_alloc::NTuple{81,Cchar} iter::Int32 cg_iter::Int32 factorization_status::Int32 factorization_integer::Int64 factorization_real::Int64 nfacts::Int32 nbacts::Int32 nmods::Int32 obj::T non_negligible_pivot::T feasible::Bool time::qpb_time_type{T} lsqp_inform::lsqp_inform_type{T} fdc_inform::fdc_inform_type{T} sbls_inform::sbls_inform_type{T} gltr_inform::gltr_inform_type{T} fit_inform::fit_inform_type
detailed documentation#
inform derived type as a Julia structure
components#
Int32 status
return status. See QPB_solve for details
Int32 alloc_status
the status of the last attempted allocation/deallocation
NTuple{81,Cchar} bad_alloc
the name of the array for which an allocation/deallocation error occurred
Int32 iter
the total number of iterations required
Int32 cg_iter
the total number of conjugate gradient iterations required
Int32 factorization_status
the return status from the factorization
Int64 factorization_integer
the total integer workspace required for the factorization
Int64 factorization_real
the total real workspace required for the factorization
Int32 nfacts
the total number of factorizations performed
Int32 nbacts
the total number of “wasted” function evaluations during the linesearch
Int32 nmods
the total number of factorizations which were modified to ensure that th matrix was an appropriate preconditioner
T obj
the value of the objective function at the best estimate of the solution determined by QPB_solve
T non_negligible_pivot
the smallest pivot which was not judged to be zero when detecting linear dependent constraints
Bool feasible
is the returned “solution” feasible?
struct qpb_time_type time
timings (see above)
struct lsqp_inform_type lsqp_inform
inform parameters for LSQP
struct fdc_inform_type fdc_inform
inform parameters for FDC
struct sbls_inform_type sbls_inform
inform parameters for SBLS
struct gltr_inform_type gltr_inform
return information from GLTR
struct fit_inform_type fit_inform
return information from FIT
example calls#
This is an example of how to use the package to solve a quadratic program; the code is available in $GALAHAD/src/qpb/Julia/test_qpb.jl . A variety of supported Hessian and constraint matrix storage formats are shown.
# test_qpb.jl
# Simple code to test the Julia interface to QPB
using GALAHAD
using Test
using Printf
using Accessors
using Quadmath
function test_qpb(::Type{T}) where T
# Derived types
data = Ref{Ptr{Cvoid}}()
control = Ref{qpb_control_type{T}}()
inform = Ref{qpb_inform_type{T}}()
# Set problem data
n = 3 # dimension
m = 2 # number of general constraints
H_ne = 3 # Hesssian elements
H_row = Cint[1, 2, 3] # row indices, NB lower triangle
H_col = Cint[1, 2, 3] # column indices, NB lower triangle
H_ptr = Cint[1, 2, 3, 4] # row pointers
H_val = T[1.0, 1.0, 1.0] # values
g = T[0.0, 2.0, 0.0] # linear term in the objective
f = one(T) # constant term in the objective
A_ne = 4 # Jacobian elements
A_row = Cint[1, 1, 2, 2] # row indices
A_col = Cint[1, 2, 2, 3] # column indices
A_ptr = Cint[1, 3, 5] # row pointers
A_val = T[2.0, 1.0, 1.0, 1.0] # values
c_l = T[1.0, 2.0] # constraint lower bound
c_u = T[2.0, 2.0] # constraint upper bound
x_l = T[-1.0, -Inf, -Inf] # variable lower bound
x_u = T[1.0, Inf, 2.0] # variable upper bound
# Set output storage
c = zeros(T, m) # constraint values
x_stat = zeros(Cint, n) # variable status
c_stat = zeros(Cint, m) # constraint status
st = ' '
status = Ref{Cint}()
@printf(" Fortran sparse matrix indexing\n\n")
@printf(" basic tests of qp storage formats\n\n")
for d in 1:7
# Initialize QPB
qpb_initialize(T, data, control, status)
# Set user-defined control options
@reset control[].f_indexing = true # Fortran sparse matrix indexing
# Start from 0
x = T[0.0, 0.0, 0.0]
y = T[0.0, 0.0]
z = T[0.0, 0.0, 0.0]
# sparse co-ordinate storage
if d == 1
st = 'C'
qpb_import(T, control, data, status, n, m,
"coordinate", H_ne, H_row, H_col, C_NULL,
"coordinate", A_ne, A_row, A_col, C_NULL)
qpb_solve_qp(T, data, status, n, m, H_ne, H_val, g, f,
A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
x_stat, c_stat)
end
# sparse by rows
if d == 2
st = 'R'
qpb_import(T, control, data, status, n, m,
"sparse_by_rows", H_ne, C_NULL, H_col, H_ptr,
"sparse_by_rows", A_ne, C_NULL, A_col, A_ptr)
qpb_solve_qp(T, data, status, n, m, H_ne, H_val, g, f,
A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
x_stat, c_stat)
end
# dense
if d == 3
st = 'D'
H_dense_ne = 6 # number of elements of H
A_dense_ne = 6 # number of elements of A
H_dense = T[1.0, 0.0, 1.0, 0.0, 0.0, 1.0]
A_dense = T[2.0, 1.0, 0.0, 0.0, 1.0, 1.0]
qpb_import(T, control, data, status, n, m,
"dense", H_ne, C_NULL, C_NULL, C_NULL,
"dense", A_ne, C_NULL, C_NULL, C_NULL)
qpb_solve_qp(T, data, status, n, m, H_dense_ne, H_dense, g, f,
A_dense_ne, A_dense, c_l, c_u, x_l, x_u,
x, c, y, z, x_stat, c_stat)
end
# diagonal
if d == 4
st = 'L'
qpb_import(T, control, data, status, n, m,
"diagonal", H_ne, C_NULL, C_NULL, C_NULL,
"sparse_by_rows", A_ne, C_NULL, A_col, A_ptr)
qpb_solve_qp(T, data, status, n, m, H_ne, H_val, g, f,
A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
x_stat, c_stat)
end
# scaled identity
if d == 5
st = 'S'
qpb_import(T, control, data, status, n, m,
"scaled_identity", H_ne, C_NULL, C_NULL, C_NULL,
"sparse_by_rows", A_ne, C_NULL, A_col, A_ptr)
qpb_solve_qp(T, data, status, n, m, H_ne, H_val, g, f,
A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
x_stat, c_stat)
end
# identity
if d == 6
st = 'I'
qpb_import(T, control, data, status, n, m,
"identity", H_ne, C_NULL, C_NULL, C_NULL,
"sparse_by_rows", A_ne, C_NULL, A_col, A_ptr)
qpb_solve_qp(T, data, status, n, m, H_ne, H_val, g, f,
A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
x_stat, c_stat)
end
# zero
if d == 7
st = 'Z'
qpb_import(T, control, data, status, n, m,
"zero", H_ne, C_NULL, C_NULL, C_NULL,
"sparse_by_rows", A_ne, C_NULL, A_col, A_ptr)
qpb_solve_qp(T, data, status, n, m, H_ne, H_val, g, f,
A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
x_stat, c_stat)
end
qpb_information(T, data, inform, status)
if inform[].status == 0
@printf("%c:%6i iterations. Optimal objective value = %5.2f status = %1i\n",
st, inform[].iter, inform[].obj, inform[].status)
else
@printf("%c: QPB_solve exit status = %1i\n", st, inform[].status)
end
# @printf("x: ")
# for i = 1:n
# @printf("%f ", x[i])
# end
# @printf("\n")
# @printf("gradient: ")
# for i = 1:n
# @printf("%f ", g[i])
# end
# @printf("\n")
# Delete internal workspace
qpb_terminate(T, data, control, inform)
end
return 0
end
@testset "QPB" begin
@test test_qpb(Float32) == 0
@test test_qpb(Float64) == 0
@test test_qpb(Float128) == 0
end