GALAHAD DQP package#
purpose#
The dqp
package uses a
dual gradient-projection method to solve a given
stricly-convex quadratic program.
The aim is to minimize the quadratic objective function
See Section 4 of $GALAHAD/doc/dqp.pdf for additional details.
terminology#
Any required solution \(x\) necessarily satisfies the primal optimality conditions
In the shifted-least-distance case, \(g\) is shifted by \(-W^2 x^0\), and \(H = W^2\), where \(W\) is the diagonal matrix whose entries are the \(w_j\).
method#
Dual gradient-projection methods solve the quadratic programmimg problem by instead solving the dual quadratic program
Both phases require the solution of sparse systems of symmetric linear
equations, and these are handled by the matrix factorization package
SBLS
or the conjugate-gradient package GLTR
. The systems are
commonly singular, and this leads to a requirement to find the Fredholm
Alternative for the given matrix and its right-hand side. In the
non-singular case, there is an option to update existing factorizations
using the “Schur-complement” approach given by the package SCU
.
Optionally, the problem may be pre-processed temporarily to eliminate dependent
constraints using the package FDC
. This may improve the
performance of the subsequent iteration.
reference#
The basic algorithm is described in
N. I. M. Gould and D. P. Robinson, ``A dual gradient-projection method for large-scale strictly-convex quadratic problems’’, Computational Optimization and Applications 67(1) (2017) 1-38.
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 dqp package must be called in the following order:
dqp_initialize - provide default control parameters and set up initial data structures
dqp_read_specfile (optional) - override control values by reading replacement values from a file
dqp_import - set up problem data structures and fixed values
dqp_reset_control (optional) - possibly change control parameters if a sequence of problems are being solved
solve the problem by calling one of
dqp_solve_qp - solve the quadratic program
dqp_solve_sldqp - solve the shifted least-distance problem
dqp_information (optional) - recover information about the solution and solution process
dqp_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 dqp_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 dqp_control_type) |
status |
is a scalar variable of type Int32 that gives the exit status from the package. Possible values are (currently):
|
function dqp_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/dqp/DQP.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/dqp.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 dqp_control_type) |
specfile |
is a one-dimensional array of type Vararg{Cchar} that must give the name of the specification file |
function dqp_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 dqp_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’, or ‘identity’; 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 dqp_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 dqp_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 dqp_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 dqp_solve_sldqp(T, data, status, n, m, w, x0, g, f, a_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z, x_stat, c_stat)
Solve the shifted least-distance quadratic program
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. |
w |
is a one-dimensional array of size n and type T that holds the values of the weights \(w\). |
x0 |
is a one-dimensional array of size n and type T that holds the values of the shifts \(x^0\). |
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 dqp_information(T, data, inform, status)
Provides output information
Parameters:
data |
holds private internal data |
inform |
is a structure containing output information (see dqp_inform_type) |
status |
is a scalar variable of type Int32 that gives the exit status from the package. Possible values are (currently):
|
function dqp_terminate(T, data, control, inform)
Deallocate all internal private storage
Parameters:
data |
holds private internal data |
control |
is a structure containing control information (see dqp_control_type) |
inform |
is a structure containing output information (see dqp_inform_type) |
available structures#
dqp_control_type structure#
struct dqp_control_type{T} f_indexing::Bool error::Int32 out::Int32 print_level::Int32 start_print::Int32 stop_print::Int32 print_gap::Int32 dual_starting_point::Int32 maxit::Int32 max_sc::Int32 cauchy_only::Int32 arc_search_maxit::Int32 cg_maxit::Int32 explore_optimal_subspace::Int32 restore_problem::Int32 sif_file_device::Int32 qplib_file_device::Int32 rho::T infinity::T stop_abs_p::T stop_rel_p::T stop_abs_d::T stop_rel_d::T stop_abs_c::T stop_rel_c::T stop_cg_relative::T stop_cg_absolute::T cg_zero_curvature::T max_growth::T identical_bounds_tol::T cpu_time_limit::T clock_time_limit::T initial_perturbation::T perturbation_reduction::T final_perturbation::T factor_optimal_matrix::Bool remove_dependencies::Bool treat_zero_bounds_as_general::Bool exact_arc_search::Bool subspace_direct::Bool subspace_alternate::Bool subspace_arc_search::Bool space_critical::Bool deallocate_error_fatal::Bool generate_sif_file::Bool generate_qplib_file::Bool symmetric_linear_solver::NTuple{31,Cchar} definite_linear_solver::NTuple{31,Cchar} unsymmetric_linear_solver::NTuple{31,Cchar} sif_file_name::NTuple{31,Cchar} qplib_file_name::NTuple{31,Cchar} prefix::NTuple{31,Cchar} fdc_control::fdc_control_type{T} sls_control::sls_control_type{T} sbls_control::sbls_control_type{T} gltr_control::gltr_control_type{T}
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 print_gap
printing will only occur every print_gap iterations
Int32 dual_starting_point
which starting point should be used for the dual problem
-1 user supplied comparing primal vs dual variables
0 user supplied
1 minimize linearized dual
2 minimize simplified quadratic dual
3 all free (= all active primal costraints)
4 all fixed on bounds (= no active primal costraints)
Int32 maxit
at most maxit inner iterations are allowed
Int32 max_sc
the maximum permitted size of the Schur complement before a refactorization is performed (used in the case where there is no Fredholm Alternative, 0 = refactor every iteration)
Int32 cauchy_only
a subspace step will only be taken when the current Cauchy step has changed no more than than cauchy_only active constraints; the subspace step will always be taken if cauchy_only < 0
Int32 arc_search_maxit
how many iterations are allowed per arc search (-ve = as many as require
Int32 cg_maxit
how many CG iterations to perform per DQP iteration (-ve reverts to n+1)
Int32 explore_optimal_subspace
once a potentially optimal subspace has been found, investigate it
0 as per an ordinary subspace
1 by increasing the maximum number of allowed CG iterations
2 by switching to a direct method
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 sif_file_device
specifies the unit number to write generated SIF file describing the current problem
Int32 qplib_file_device
specifies the unit number to write generated QPLIB file describing the current problem
T rho
the penalty weight, rho. The general constraints are not enforced explicitly, but instead included in the objective as a penalty term weighted by rho when rho > 0. If rho <= 0, the general constraints are explicit (that is, there is no penalty term in the objective function)
T infinity
any bound larger than infinity in modulus will be regarded as infinite
T stop_abs_p
the required absolute and relative accuracies for the primal infeasibilies
T stop_rel_p
see stop_abs_p
T stop_abs_d
the required absolute and relative accuracies for the dual infeasibility
T stop_rel_d
see stop_abs_d
T stop_abs_c
the required absolute and relative accuracies for the complementarity
T stop_rel_c
see stop_abs_c
T stop_cg_relative
the CG iteration will be stopped as soon as the current norm of the preconditioned gradient is smaller than max( stop_cg_relative * initial preconditioned gradient, stop_cg_absolute )
T stop_cg_absolute
see stop_cg_relative
T cg_zero_curvature
threshold below which curvature is regarded as zero if CG is used
T max_growth
maximum growth factor allowed without a refactorization
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 cpu_time_limit
the maximum CPU time allowed (-ve means infinite)
T clock_time_limit
the maximum elapsed clock time allowed (-ve means infinite)
T initial_perturbation
the initial penalty weight (for DLP only)
T perturbation_reduction
the penalty weight reduction factor (for DLP only)
T final_perturbation
the final penalty weight (for DLP only)
Bool factor_optimal_matrix
are the factors of the optimal augmented matrix required? (for DLP only)
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 exact_arc_search
if .exact_arc_search is true, an exact piecewise arc search will be performed. Otherwise an ineaxt search using a backtracing Armijo strategy will be employed
Bool subspace_direct
if .subspace_direct is true, the subspace step will be calculated using a direct (factorization) method, while if it is false, an iterative (conjugate-gradient) method will be used.
Bool subspace_alternate
if .subspace_alternate is true, the subspace step will alternate between a direct (factorization) method and an iterative (GLTR conjugate-gradient) method. This will override .subspace_direct
Bool subspace_arc_search
if .subspace_arc_search is true, a piecewise arc search will be performed along the subspace step. Otherwise the search will stop at the firstconstraint encountered
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
Bool generate_qplib_file
if .generate_qplib_file is .true. if a QPLIB file describing the current problem is to be generated
char symmetric_linear_solver[31]
indefinite linear equation solver set in symmetric_linear_solver
char definite_linear_solver[31]
definite linear equation solver
char unsymmetric_linear_solver[31]
unsymmetric linear equation solver
NTuple{31,Cchar} sif_file_name
name of generated SIF file containing input problem
NTuple{31,Cchar} qplib_file_name
name of generated QPLIB 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 fdc_control_type fdc_control
control parameters for FDC
struct sls_control_type sls_control
control parameters for SLS
struct sbls_control_type sbls_control
control parameters for SBLS
struct gltr_control_type gltr_control
control parameters for GLTR
dqp_time_type structure#
struct dqp_time_type{T} total::T preprocess::T find_dependent::T analyse::T factorize::T solve::T search::T clock_total::T clock_preprocess::T clock_find_dependent::T clock_analyse::T clock_factorize::T clock_solve::T clock_search::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 factorization
T factorize
the CPU time spent factorizing the required matrices
T solve
the CPU time spent computing the search direction
T search
the CPU time spent in the linesearch
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 factorization
T clock_factorize
the clock time spent factorizing the required matrices
T clock_solve
the clock time spent computing the search direction
T clock_search
the clock time spent in the linesearch
dqp_inform_type structure#
struct dqp_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 threads::Int32 obj::T primal_infeasibility::T dual_infeasibility::T complementary_slackness::T non_negligible_pivot::T feasible::Bool checkpointsIter::NTuple{16,Cint} checkpointsTime::NTuple{16,T} time::dqp_time_type{T} fdc_inform::fdc_inform_type{T} sls_inform::sls_inform_type{T} sbls_inform::sbls_inform_type{T} gltr_inform::gltr_inform_type{T} scu_status::Int32 scu_inform::scu_inform_type rpd_inform::rpd_inform_type
detailed documentation#
inform derived type as a Julia structure
components#
Int32 status
return status. See DQP_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 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 threads
the number of threads used
T obj
the value of the objective function at the best estimate of the solution determined by DQP_solve
T primal_infeasibility
the value of the primal infeasibility
T dual_infeasibility
the value of the dual infeasibility
T complementary_slackness
the value of the complementary slackness
T non_negligible_pivot
the smallest pivot that was not judged to be zero when detecting linearly dependent constraints
Bool feasible
is the returned “solution” feasible?
Int32 checkpointsIter[16]
checkpoints(i) records the iteration at which the criticality measures first fall below \(10^{-i-1}\), i = 0, …, 15 (-1 means not achieved)
T checkpointsTime[16]
see checkpointsIter
struct dqp_time_type time
timings (see above)
struct fdc_inform_type fdc_inform
inform parameters for FDC
struct sls_inform_type sls_inform
inform parameters for SLS
struct sbls_inform_type sbls_inform
inform parameters for SBLS
struct gltr_inform_type gltr_inform
return information from GLTR
Int32 scu_status
inform parameters for SCU
struct scu_inform_type scu_inform
see scu_status
struct rpd_inform_type rpd_inform
inform parameters for RPD
example calls#
This is an example of how to use the package to solve a given convex quadratic program; the code is available in $GALAHAD/src/dqp/Julia/test_dqp.jl . A variety of supported Hessian and constraint matrix storage formats are shown.
# test_dqp.jl
# Simple code to test the Julia interface to DQP
using GALAHAD
using Test
using Printf
using Accessors
using Quadmath
function test_dqp(::Type{T}) where T
# Derived types
data = Ref{Ptr{Cvoid}}()
control = Ref{dqp_control_type{T}}()
inform = Ref{dqp_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:6
# Initialize DQP
dqp_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'
dqp_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)
dqp_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'
dqp_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)
dqp_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]
dqp_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)
dqp_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'
dqp_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)
dqp_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'
dqp_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)
dqp_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'
dqp_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)
dqp_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
dqp_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: DQP_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
dqp_terminate(T, data, control, inform)
end
# test shifted least-distance interface
for d in 1:1
# Initialize DQP
dqp_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]
# Set shifted least-distance data
w = T[1.0, 1.0, 1.0]
x_0 = T[0.0, 0.0, 0.0]
# sparse co-ordinate storage
if d == 1
st = 'W'
dqp_import(T, control, data, status, n, m,
"shifted_least_distance", H_ne, C_NULL, C_NULL, C_NULL,
"coordinate", A_ne, A_row, A_col, C_NULL)
dqp_solve_sldqp(T, data, status, n, m, w, x_0, g, f,
A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
x_stat, c_stat)
end
dqp_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: DQP_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
dqp_terminate(T, data, control, inform)
end
return 0
end
@testset "DQP" begin
@test test_dqp(Float32) == 0
@test test_dqp(Float64) == 0
@test test_dqp(Float128) == 0
end