GALAHAD QPA package#
purpose#
The qpa
package uses a working-set method to solve
non-convex quadratic programs in various guises.
The first is the l\({}_1\) quadratic programming problem
that aims to minimize
The package may also be used to solve the standard quadratic programming problem whose aim is to minimize \(q(x)\) subject to the general linear constraints and simple bounds
Similarly, the package is capable of solving the bound-constrained l\({}_1\) quadratic programming problem whose intention is to minimize \(q(x) + \rho_g v_g(x),\) subject to the above simple bound constraints by automatically adjusting \(\rho_b\) in \(f(x;\rho_g,\rho_b)\).
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/qpa.pdf for additional details.
N.B. In many cases, the alternative quadratic programming package qpb
is faster, and thus to be preferred.
terminolgy#
Any required solution \(x\) for the standard quadratic programming problem necessarily satisfies the primal optimality conditions
method#
At the \(k\)-th iteration of the method, an improvement to the value of the merit function \(m(x, \rho_g, \rho_b ) = q(x) + \rho_g v_g(x) + \rho_b v_b(x)\) at \(x = x^{(k)}\) is sought. This is achieved by first computing a search direction \(s^{(k)}\), and then setting \(x^{(k+1)} = x^{(k)} + \alpha^{(k)} s^{(k)}\), where the stepsize \(\alpha^{(k)}\) is chosen as the first local minimizer of \(\phi ( \alpha ) = m( x^{(k)} + \alpha s^{(k)} , \rho_g, \rho_b )\) as \(\alpha\) incesases from zero. The stepsize calculation is straightforward, and exploits the fact that \(\phi ( \alpha )\) is a piecewise quadratic function of \(\alpha\).
The search direction is defined by a subset of the “active” terms in \(v(x)\), i.e., those for which \(a_i^T x = c_i^l\) or \(c_i^u\) (for \(i=1,\ldots,m\)) or \(x_j = x_j^l\) or \(x_j^u\) (for (\(j=1,\ldots,n\)). The “working” set \(W^{(k)}\) is chosen from the active terms, and is such that its members have linearly independent gradients. The search direction \(s^{(k)}\) is chosen as an approximate solution of the equality-constrained quadratic program
Preconditioning of the conjugate gradient iteration requires the solution of one or more linear systems of the form
SLS
. This
reference matrix may be factorized as a whole (the so-called “augmented
system” approach), or by performing a block elimination first (the
“Schur-complement” approach). The latter is usually to be preferred when
a (non-singular) diagonal preconditioner is used, but may be inefficient
if any of the columns of \(A^{(l)}\) is too dense. Subsequent iterations
within the current major iteration obtain solutions to (3) via the
factors of \(K^{(l)}\) and an appropriate (dense) Schur complement,
obtained from the package SCU
. The major iteration terminates once
the space required to hold the factors of the (growing) Schur complement
exceeds a given threshold.
The working set changes by (a) adding an active term encountered during the determination of the stepsize, or (b) the removal of a term if \(s = 0\) solves (1)–(2). The decision on which to remove in the latter case is based upon the expected decrease upon the removal of an individual term, and this information is available from the magnitude and sign of the components of the auxiliary vector \(u\) computed in (3). At optimality, the components of \(u\) for \(a_i\) terms will all lie between \(0\) and \(\rho_g\) — and those for the other terms between \(0\) and \(\rho_b\) — and any violation of this rule indicates further progress is possible. The components of \(u\) corresonding to the terms involving \(a_i^T x\) are sometimes known as Lagrange multipliers (or generalized gradients) and denoted by \(y\), while those for the remaining \(x_j\) terms are dual variables and denoted by \(z\).
To solve the standard quadratic programming problem, a sequence of \(\ell_1\)-quadratic programs are solved, each with a larger value of \(\rho_g\) and/or \(\rho_b\) than its predecessor. The required solution has been found once the infeasibilities \(v_g(x)\) and \(v_b(x)\) have been reduced to zero at the solution of the \(\ell_1\)-problem for the given \(\rho_g\) and \(\rho_b\).
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
N. I. M. Gould and Ph. L. Toint ``An iterative working-set method for large-scale non-convex quadratic programming’’. Applied Numerical Mathematics 43(1–2) (2002) 109–128.
matrix 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 \(0 \leq i \leq m-1\), \(0 \leq j \leq n-1\). 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 \(0 \leq i \leq m-1\), \(0 \leq j \leq n-1\). 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, \(0 \leq l \leq ne-1\), of \(A\), its row index i, column index j and value \(A_{ij}\), \(0 \leq i \leq m-1\), \(0 \leq j \leq n-1\), 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) holds the total number of entries. The column indices j, \(0 \leq j \leq n-1\), 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, \(0 \leq i \leq m-1\), 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) holds the total number of entries. The row indices i, \(0 \leq i \leq m-1\), 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, \(0 \leq j \leq n-1\), 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.
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 \(0 \leq j \leq i \leq n-1\)) need be held. In this case the lower triangle should be stored by rows, that is component \(i * i / 2 + j\) of the storage array H_val will hold the value \(H_{ij}\) (and, by symmetry, \(H_{ji}\)) for \(0 \leq j \leq i \leq n-1\). 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, \(0 \leq l \leq ne-1\), of \(H\), its row index i, column index j and value \(H_{ij}\), \(0 \leq j \leq i \leq n-1\), 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) holds the total number of entries. The column indices j, \(0 \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 \(0 \leq i \neq j \leq n-1\)) only the diagonals entries \(H_{ii}\), \(0 \leq i \leq n-1\) 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 qpa package must be called in the following order:
qpa_initialize - provide default control parameters and set up initial data structures
qpa_read_specfile (optional) - override control values by reading replacement values from a file
qpa_import - set up problem data structures and fixed values
qpa_reset_control (optional) - possibly change control parameters if a sequence of problems are being solved
solve the problem by calling one of
qpa_solve_qp - solve the quadratic program (2)-(4)
qpa_solve_l1qp - solve the l1 quadratic program (1)
qpa_solve_bcl1qp - solve the bound constrained l1 quadratic program (4)-(5)
qpa_information (optional) - recover information about the solution and solution process
qpa_terminate - deallocate data structures
See the examples section for illustrations of use.
callable functions#
overview of functions provided#
// typedefs typedef float spc_; typedef double rpc_; typedef int ipc_; // structs struct qpa_control_type; struct qpa_inform_type; struct qpa_time_type; // function calls void qpa_initialize(void **data, struct qpa_control_type* control, ipc_ *status); void qpa_read_specfile(struct qpa_control_type* control, const char specfile[]); void qpa_import( struct qpa_control_type* control, void **data, ipc_ *status, ipc_ n, ipc_ m, const char H_type[], ipc_ H_ne, const ipc_ H_row[], const ipc_ H_col[], const ipc_ H_ptr[], const char A_type[], ipc_ A_ne, const ipc_ A_row[], const ipc_ A_col[], const ipc_ A_ptr[] ); void qpa_reset_control( struct qpa_control_type* control, void **data, ipc_ *status ); void qpa_solve_qp( void **data, ipc_ *status, ipc_ n, ipc_ m, ipc_ h_ne, const rpc_ H_val[], const rpc_ g[], const rpc_ f, ipc_ a_ne, const rpc_ A_val[], const rpc_ c_l[], const rpc_ c_u[], const rpc_ x_l[], const rpc_ x_u[], rpc_ x[], rpc_ c[], rpc_ y[], rpc_ z[], ipc_ x_stat[], ipc_ c_stat[] ); void qpa_solve_l1qp( void **data, ipc_ *status, ipc_ n, ipc_ m, ipc_ h_ne, const rpc_ H_val[], const rpc_ g[], const rpc_ f, const rpc_ rho_g, const rpc_ rho_b, ipc_ a_ne, const rpc_ A_val[], const rpc_ c_l[], const rpc_ c_u[], const rpc_ x_l[], const rpc_ x_u[], rpc_ x[], rpc_ c[], rpc_ y[], rpc_ z[], ipc_ x_stat[], ipc_ c_stat[] ); void qpa_solve_bcl1qp( void **data, ipc_ *status, ipc_ n, ipc_ m, ipc_ h_ne, const rpc_ H_val[], const rpc_ g[], const rpc_ f, const rpc_ rho_g, ipc_ a_ne, const rpc_ A_val[], const rpc_ c_l[], const rpc_ c_u[], const rpc_ x_l[], const rpc_ x_u[], rpc_ x[], rpc_ c[], rpc_ y[], rpc_ z[], ipc_ x_stat[], ipc_ c_stat[] ); void qpa_information(void **data, struct qpa_inform_type* inform, ipc_ *status); void qpa_terminate( void **data, struct qpa_control_type* control, struct qpa_inform_type* inform );
typedefs#
typedef float spc_
spc_
is real single precision
typedef double rpc_
rpc_
is the real working precision used, but may be changed to float
by
defining the preprocessor variable REAL_32
or (if supported) to
__real128
using the variable REAL_128
.
typedef int ipc_
ipc_
is the default integer word length used, but may be changed to
int64_t
by defining the preprocessor variable INTEGER_64
.
function calls#
void qpa_initialize(void **data, struct qpa_control_type* control, ipc_ *status)
Set default control values and initialize private data
Parameters:
data |
holds private internal data |
control |
is a struct containing control information (see qpa_control_type) |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are (currently):
|
void qpa_read_specfile(struct qpa_control_type* control, const char 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/qpa/QPA.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/qpa.pdf for a list of how these keywords relate to the components of the control structure.
Parameters:
control |
is a struct containing control information (see qpa_control_type) |
specfile |
is a character string containing the name of the specification file |
void qpa_import( struct qpa_control_type* control, void **data, ipc_ *status, ipc_ n, ipc_ m, const char H_type[], ipc_ H_ne, const ipc_ H_row[], const ipc_ H_col[], const ipc_ H_ptr[], const char A_type[], ipc_ A_ne, const ipc_ A_row[], const ipc_ A_col[], const ipc_ A_ptr[] )
Import problem data into internal storage prior to solution.
Parameters:
control |
is a struct whose members provide control paramters for the remaining prcedures (see qpa_control_type) |
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are:
|
n |
is a scalar variable of type ipc_, that holds the number of variables. |
m |
is a scalar variable of type ipc_, that holds the number of general linear constraints. |
H_type |
is a one-dimensional array of type char 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 ipc_, 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 ipc_, 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 NULL. |
H_col |
is a one-dimensional array of size H_ne and type ipc_, 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 NULL. |
H_ptr |
is a one-dimensional array of size n+1 and type ipc_, 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 NULL. |
A_type |
is a one-dimensional array of type char 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 ipc_, 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 ipc_, 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 NULL. |
A_col |
is a one-dimensional array of size A_ne and type ipc_, 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 NULL. |
A_ptr |
is a one-dimensional array of size n+1 and type ipc_, 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 NULL. |
void qpa_reset_control( struct qpa_control_type* control, void **data, ipc_ *status )
Reset control parameters after import if required.
Parameters:
control |
is a struct whose members provide control paramters for the remaining prcedures (see qpa_control_type) |
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are:
|
void qpa_solve_qp( void **data, ipc_ *status, ipc_ n, ipc_ m, ipc_ h_ne, const rpc_ H_val[], const rpc_ g[], const rpc_ f, ipc_ a_ne, const rpc_ A_val[], const rpc_ c_l[], const rpc_ c_u[], const rpc_ x_l[], const rpc_ x_u[], rpc_ x[], rpc_ c[], rpc_ y[], rpc_ z[], ipc_ x_stat[], ipc_ c_stat[] )
Solve the quadratic program (2)-(4).
Parameters:
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the entry and exit status from the package. Possible exit values are:
|
n |
is a scalar variable of type ipc_, that holds the number of variables |
m |
is a scalar variable of type ipc_, that holds the number of general linear constraints. |
h_ne |
is a scalar variable of type ipc_, 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 rpc_, 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 rpc_, that holds the linear term \(g\) of the objective function. The j-th component of g, j = 0, … , n-1, contains \(g_j\). |
f |
is a scalar of type rpc_, that holds the constant term \(f\) of the objective function. |
a_ne |
is a scalar variable of type ipc_, 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 rpc_, 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 rpc_, that holds the lower bounds \(c^l\) on the constraints \(A x\). The i-th component of c_l, i = 0, … , m-1, contains \(c^l_i\). |
c_u |
is a one-dimensional array of size m and type rpc_, that holds the upper bounds \(c^l\) on the constraints \(A x\). The i-th component of c_u, i = 0, … , m-1, contains \(c^u_i\). |
x_l |
is a one-dimensional array of size n and type rpc_, that holds the lower bounds \(x^l\) on the variables \(x\). The j-th component of x_l, j = 0, … , n-1, contains \(x^l_j\). |
x_u |
is a one-dimensional array of size n and type rpc_, that holds the upper bounds \(x^l\) on the variables \(x\). The j-th component of x_u, j = 0, … , n-1, contains \(x^l_j\). |
x |
is a one-dimensional array of size n and type rpc_, that holds the values \(x\) of the optimization variables. The j-th component of x, j = 0, … , n-1, contains \(x_j\). |
c |
is a one-dimensional array of size m and type rpc_, that holds the residual \(c(x)\). The i-th component of c, j = 0, … , n-1, contains \(c_j(x)\). |
y |
is a one-dimensional array of size n and type rpc_, that holds the values \(y\) of the Lagrange multipliers for the general linear constraints. The j-th component of y, j = 0, … , n-1, contains \(y_j\). |
z |
is a one-dimensional array of size n and type rpc_, that holds the values \(z\) of the dual variables. The j-th component of z, j = 0, … , n-1, contains \(z_j\). |
x_stat |
is a one-dimensional array of size n and type ipc_, that gives the current 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. On entry, if control.cold_start = 0, x_stat should be set as above to provide a guide to the initial working set. |
c_stat |
is a one-dimensional array of size m and type ipc_, that gives the current 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. On entry, if control.cold_start = 0, c_stat should be set as above to provide a guide to the initial working set. |
void qpa_solve_l1qp( void **data, ipc_ *status, ipc_ n, ipc_ m, ipc_ h_ne, const rpc_ H_val[], const rpc_ g[], const rpc_ f, const rpc_ rho_g, const rpc_ rho_b, ipc_ a_ne, const rpc_ A_val[], const rpc_ c_l[], const rpc_ c_u[], const rpc_ x_l[], const rpc_ x_u[], rpc_ x[], rpc_ c[], rpc_ y[], rpc_ z[], ipc_ x_stat[], ipc_ c_stat[] )
Solve the l_1 quadratic program (1).
Parameters:
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the entry and exit status from the package. Possible exit values are:
|
n |
is a scalar variable of type ipc_, that holds the number of variables |
m |
is a scalar variable of type ipc_, that holds the number of general linear constraints. |
h_ne |
is a scalar variable of type ipc_, 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 rpc_, 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 rpc_, that holds the linear term \(g\) of the objective function. The j-th component of g, j = 0, … , n-1, contains \(g_j\). |
f |
is a scalar of type rpc_, that holds the constant term \(f\) of the objective function. |
rho_g |
is a scalar of type rpc_, that holds the parameter \(\rho_g\) associated with the linear constraints. |
rho_b |
is a scalar of type rpc_, that holds the parameter \(\rho_b\) associated with the simple bound constraints. |
a_ne |
is a scalar variable of type ipc_, 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 rpc_, 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 rpc_, that holds the lower bounds \(c^l\) on the constraints \(A x\). The i-th component of c_l, i = 0, … , m-1, contains \(c^l_i\). |
c_u |
is a one-dimensional array of size m and type rpc_, that holds the upper bounds \(c^l\) on the constraints \(A x\). The i-th component of c_u, i = 0, … , m-1, contains \(c^u_i\). |
x_l |
is a one-dimensional array of size n and type rpc_, that holds the lower bounds \(x^l\) on the variables \(x\). The j-th component of x_l, j = 0, … , n-1, contains \(x^l_j\). |
x_u |
is a one-dimensional array of size n and type rpc_, that holds the upper bounds \(x^l\) on the variables \(x\). The j-th component of x_u, j = 0, … , n-1, contains \(x^l_j\). |
x |
is a one-dimensional array of size n and type rpc_, that holds the values \(x\) of the optimization variables. The j-th component of x, j = 0, … , n-1, contains \(x_j\). |
c |
is a one-dimensional array of size m and type rpc_, that holds the residual \(c(x)\). The i-th component of c, j = 0, … , n-1, contains \(c_j(x)\). |
y |
is a one-dimensional array of size n and type rpc_, that holds the values \(y\) of the Lagrange multipliers for the general linear constraints. The j-th component of y, j = 0, … , n-1, contains \(y_j\). |
z |
is a one-dimensional array of size n and type rpc_, that holds the values \(z\) of the dual variables. The j-th component of z, j = 0, … , n-1, contains \(z_j\). |
x_stat |
is a one-dimensional array of size n and type ipc_, that gives the current 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. On entry, if control.cold_start = 0, x_stat should be set as above to provide a guide to the initial working set. |
c_stat |
is a one-dimensional array of size m and type ipc_, that gives the current 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. On entry, if control.cold_start = 0, c_stat should be set as above to provide a guide to the initial working set. |
void qpa_solve_bcl1qp( void **data, ipc_ *status, ipc_ n, ipc_ m, ipc_ h_ne, const rpc_ H_val[], const rpc_ g[], const rpc_ f, const rpc_ rho_g, ipc_ a_ne, const rpc_ A_val[], const rpc_ c_l[], const rpc_ c_u[], const rpc_ x_l[], const rpc_ x_u[], rpc_ x[], rpc_ c[], rpc_ y[], rpc_ z[], ipc_ x_stat[], ipc_ c_stat[] )
Solve the bound-constrained l_1 quadratic program (4)-(5)
Parameters:
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the entry and exit status from the package. Possible exit values are:
|
n |
is a scalar variable of type ipc_, that holds the number of variables |
m |
is a scalar variable of type ipc_, that holds the number of general linear constraints. |
h_ne |
is a scalar variable of type ipc_, 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 rpc_, 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 rpc_, that holds the linear term \(g\) of the objective function. The j-th component of g, j = 0, … , n-1, contains \(g_j\). |
f |
is a scalar of type rpc_, that holds the constant term \(f\) of the objective function. |
rho_g |
is a scalar of type rpc_, that holds the parameter \(\rho_g\) associated with the linear constraints. |
a_ne |
is a scalar variable of type ipc_, 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 rpc_, 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 rpc_, that holds the lower bounds \(c^l\) on the constraints \(A x\). The i-th component of c_l, i = 0, … , m-1, contains \(c^l_i\). |
c_u |
is a one-dimensional array of size m and type rpc_, that holds the upper bounds \(c^l\) on the constraints \(A x\). The i-th component of c_u, i = 0, … , m-1, contains \(c^u_i\). |
x_l |
is a one-dimensional array of size n and type rpc_, that holds the lower bounds \(x^l\) on the variables \(x\). The j-th component of x_l, j = 0, … , n-1, contains \(x^l_j\). |
x_u |
is a one-dimensional array of size n and type rpc_, that holds the upper bounds \(x^l\) on the variables \(x\). The j-th component of x_u, j = 0, … , n-1, contains \(x^l_j\). |
x |
is a one-dimensional array of size n and type rpc_, that holds the values \(x\) of the optimization variables. The j-th component of x, j = 0, … , n-1, contains \(x_j\). |
c |
is a one-dimensional array of size m and type rpc_, that holds the residual \(c(x)\). The i-th component of c, j = 0, … , n-1, contains \(c_j(x)\). |
y |
is a one-dimensional array of size n and type rpc_, that holds the values \(y\) of the Lagrange multipliers for the general linear constraints. The j-th component of y, j = 0, … , n-1, contains \(y_j\). |
z |
is a one-dimensional array of size n and type rpc_, that holds the values \(z\) of the dual variables. The j-th component of z, j = 0, … , n-1, contains \(z_j\). |
x_stat |
is a one-dimensional array of size n and type ipc_, that gives the current 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. On entry, if control.cold_start = 0, x_stat should be set as above to provide a guide to the initial working set. |
c_stat |
is a one-dimensional array of size m and type ipc_, that gives the current 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. On entry, if control.cold_start = 0, c_stat should be set as above to provide a guide to the initial working set. |
void qpa_information(void **data, struct qpa_inform_type* inform, ipc_ *status)
Provides output information
Parameters:
data |
holds private internal data |
inform |
is a struct containing output information (see qpa_inform_type) |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are (currently):
|
void qpa_terminate( void **data, struct qpa_control_type* control, struct qpa_inform_type* inform )
Deallocate all internal private storage
Parameters:
data |
holds private internal data |
control |
is a struct containing control information (see qpa_control_type) |
inform |
is a struct containing output information (see qpa_inform_type) |
available structures#
qpa_control_type structure#
#include <galahad_qpa.h> struct qpa_control_type { // components bool f_indexing; ipc_ error; ipc_ out; ipc_ print_level; ipc_ start_print; ipc_ stop_print; ipc_ maxit; ipc_ factor; ipc_ max_col; ipc_ max_sc; ipc_ indmin; ipc_ valmin; ipc_ itref_max; ipc_ infeas_check_interval; ipc_ cg_maxit; ipc_ precon; ipc_ nsemib; ipc_ full_max_fill; ipc_ deletion_strategy; ipc_ restore_problem; ipc_ monitor_residuals; ipc_ cold_start; ipc_ sif_file_device; rpc_ infinity; rpc_ feas_tol; rpc_ obj_unbounded; rpc_ increase_rho_g_factor; rpc_ infeas_g_improved_by_factor; rpc_ increase_rho_b_factor; rpc_ infeas_b_improved_by_factor; rpc_ pivot_tol; rpc_ pivot_tol_for_dependencies; rpc_ zero_pivot; rpc_ inner_stop_relative; rpc_ inner_stop_absolute; rpc_ multiplier_tol; rpc_ cpu_time_limit; rpc_ clock_time_limit; bool treat_zero_bounds_as_general; bool solve_qp; bool solve_within_bounds; bool randomize; bool array_syntax_worse_than_do_loop; bool space_critical; bool deallocate_error_fatal; bool generate_sif_file; char symmetric_linear_solver[31]; char sif_file_name[31]; char prefix[31]; bool each_interval; struct sls_control_type sls_control; };
detailed documentation#
control derived type as a C struct
components#
bool f_indexing
use C or Fortran sparse matrix indexing
ipc_ error
error and warning diagnostics occur on stream error
ipc_ out
general output occurs on stream out
ipc_ print_level
the level of output required is specified by print_level
ipc_ start_print
any printing will start on this iteration
ipc_ stop_print
any printing will stop on this iteration
ipc_ maxit
at most maxit inner iterations are allowed
ipc_ factor
the factorization to be used. Possible values are 0 automatic 1 Schur-complement factorization 2 augmented-system factorization
ipc_ max_col
the maximum number of nonzeros in a column of A which is permitted with the Schur-complement factorization
ipc_ max_sc
the maximum permitted size of the Schur complement before a refactorization is performed
ipc_ indmin
an initial guess as to the integer workspace required by SLS (OBSOLETE)
ipc_ valmin
an initial guess as to the real workspace required by SLS (OBSOLETE)
ipc_ itref_max
the maximum number of iterative refinements allowed (OBSOLETE)
ipc_ infeas_check_interval
the infeasibility will be checked for improvement every infeas_check_interval iterations (see infeas_g_improved_by_factor and infeas_b_improved_by_factor below)
ipc_ 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
ipc_ 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
ipc_ nsemib
the semi-bandwidth of a band preconditioner, if appropriate
ipc_ full_max_fill
if the ratio of the number of nonzeros in the factors of the reference matrix to the number of nonzeros in the matrix itself exceeds full_max_fill, and the preconditioner is being selected automatically (precon = 0), a banded approximation will be used instead
ipc_ deletion_strategy
the constraint deletion strategy to be used. Possible values are:
0 most violated of all 1 LIFO (last in, first out) k LIFO(k) most violated of the last k in LIFO
ipc_ 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
ipc_ monitor_residuals
the frequency at which residuals will be monitored
ipc_ cold_start
indicates whether a cold or warm start should be made. Possible values are
0 warm start - the values set in C_stat and B_stat indicate which constraints will be included in the initial working set. 1 cold start from the value set in X; constraints active at X will determine the initial working set. 2 cold start with no active constraints 3 cold start with only equality constraints active 4 cold start with as many active constraints as possible
ipc_ sif_file_device
specifies the unit number to write generated SIF file describing the current problem
rpc_ infinity
any bound larger than infinity in modulus will be regarded as infinite
rpc_ feas_tol
any constraint violated by less than feas_tol will be considered to be satisfied
rpc_ obj_unbounded
if the objective function value is smaller than obj_unbounded, it will be flagged as unbounded from below.
rpc_ increase_rho_g_factor
if the problem is currently infeasible and solve_qp (see below) is .TRUE. the current penalty parameter for the general constraints will be increased by increase_rho_g_factor when needed
rpc_ infeas_g_improved_by_factor
if the infeasibility of the general constraints has not dropped by a fac of infeas_g_improved_by_factor over the previous infeas_check_interval iterations, the current corresponding penalty parameter will be increase
rpc_ increase_rho_b_factor
if the problem is currently infeasible and solve_qp or solve_within_boun (see below) is .TRUE., the current penalty parameter for the simple boun constraints will be increased by increase_rho_b_factor when needed
rpc_ infeas_b_improved_by_factor
if the infeasibility of the simple bounds has not dropped by a factor of infeas_b_improved_by_factor over the previous infeas_check_interval iterations, the current corresponding penalty parameter will be increase
rpc_ pivot_tol
the threshold pivot used by the matrix factorization. See the documentation for SLS for details (OBSOLE
rpc_ pivot_tol_for_dependencies
the threshold pivot used by the matrix factorization when attempting to detect linearly dependent constraints.
rpc_ zero_pivot
any pivots smaller than zero_pivot in absolute value will be regarded to zero when attempting to detect linearly dependent constraints (OBSOLE
rpc_ 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 )
rpc_ inner_stop_absolute
see inner_stop_relative
rpc_ multiplier_tol
any dual variable or Lagrange multiplier which is less than multiplier_t outside its optimal interval will be regarded as being acceptable when checking for optimality
rpc_ cpu_time_limit
the maximum CPU time allowed (-ve means infinite)
rpc_ clock_time_limit
the maximum elapsed clock time allowed (-ve means infinite)
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 solve_qp
if solve_qp is .TRUE., the value of prob.rho_g and prob.rho_b will be increased as many times as are needed to ensure that the output solution is feasible, and thus aims to solve the quadratic program (2)-(4)
bool solve_within_bounds
if solve_within_bounds is .TRUE., the value of prob.rho_b will be increased as many times as are needed to ensure that the output solution is feasible with respect to the simple bounds, and thus aims to solve the bound-constrained quadratic program (4)-(5)
bool randomize
if randomize is .TRUE., the constraint bounds will be perturbed by small random quantities during the first stage of the solution process. Any randomization will ultimately be removed. Randomization helps when solving degenerate problems
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 (OBSOLETE)
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
char symmetric_linear_solver[31]
the name of the symmetric-indefinite linear equation solver used. Possible choices are currently: ‘sils’, ‘ma27’, ‘ma57’, ‘ma77’, ‘ma86’, ‘ma97’, ‘ssids’, ‘mumps’, ‘pardiso’, ‘mkl_pardiso’, ‘pastix’, ‘wsmp’, and ‘sytr’, although only ‘sytr’ and, for OMP 4.0-compliant compilers, ‘ssids’ are installed by default; others are easily installed (see README.external). More details of the capabilities of each solver are provided in the documentation for galahad_sls.
char sif_file_name[31]
definite linear equation solver
name of generated SIF file containing input problem
char prefix[31]
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’
bool each_interval
component specifically for parametric problems (not used at present)
struct sls_control_type sls_control
control parameters for SLS
qpa_time_type structure#
#include <galahad_qpa.h> struct qpa_time_type { // components rpc_ total; rpc_ preprocess; rpc_ analyse; rpc_ factorize; rpc_ solve; rpc_ clock_total; rpc_ clock_preprocess; rpc_ clock_analyse; rpc_ clock_factorize; rpc_ clock_solve; };
detailed documentation#
time derived type as a C struct
components#
rpc_ total
the total CPU time spent in the package
rpc_ preprocess
the CPU time spent preprocessing the problem
rpc_ analyse
the CPU time spent analysing the required matrices prior to factorizatio
rpc_ factorize
the CPU time spent factorizing the required matrices
rpc_ solve
the CPU time spent computing the search direction
rpc_ clock_total
the total clock time spent in the package
rpc_ clock_preprocess
the clock time spent preprocessing the problem
rpc_ clock_analyse
the clock time spent analysing the required matrices prior to factorizat
rpc_ clock_factorize
the clock time spent factorizing the required matrices
rpc_ clock_solve
the clock time spent computing the search direction
qpa_inform_type structure#
#include <galahad_qpa.h> struct qpa_inform_type { // components ipc_ status; ipc_ alloc_status; char bad_alloc[81]; ipc_ major_iter; ipc_ iter; ipc_ cg_iter; ipc_ factorization_status; int64_t factorization_integer; int64_t factorization_real; ipc_ nfacts; ipc_ nmods; ipc_ num_g_infeas; ipc_ num_b_infeas; rpc_ obj; rpc_ infeas_g; rpc_ infeas_b; rpc_ merit; struct qpa_time_type time; struct sls_inform_type sls_inform; };
detailed documentation#
inform derived type as a C struct
components#
ipc_ status
return status. See QPA_solve for details
ipc_ alloc_status
the status of the last attempted allocation/deallocation
char bad_alloc[81]
the name of the array for which an allocation/deallocation error occurred
ipc_ major_iter
the total number of major iterations required
ipc_ iter
the total number of iterations required
ipc_ cg_iter
the total number of conjugate gradient iterations required
ipc_ factorization_status
the return status from the factorization
int64_t factorization_integer
the total integer workspace required for the factorization
int64_t factorization_real
the total real workspace required for the factorization
ipc_ nfacts
the total number of factorizations performed
ipc_ nmods
the total number of factorizations which were modified to ensure that th matrix was an appropriate preconditioner
ipc_ num_g_infeas
the number of infeasible general constraints
ipc_ num_b_infeas
the number of infeasible simple-bound constraints
rpc_ obj
the value of the objective function at the best estimate of the solution determined by QPA_solve
rpc_ infeas_g
the 1-norm of the infeasibility of the general constraints
rpc_ infeas_b
the 1-norm of the infeasibility of the simple-bound constraints
rpc_ merit
the merit function value = obj + rho_g * infeas_g + rho_b * infeas_b
struct qpa_time_type time
timings (see above)
struct sls_inform_type sls_inform
inform parameters for SLS
example calls#
This is an example of how to use the package to solve a quadratic program; the code is available in $GALAHAD/src/qpa/C/qpat.c . A variety of supported Hessian and constraint matrix storage formats are shown.
Notice that C-style indexing is used, and that this is flagged by setting
control.f_indexing
to false
. The floating-point type rpc_
is set in galahad_precision.h
to double
by default, but to float
if the preprocessor variable SINGLE
is defined. Similarly, the integer
type ipc_
from galahad_precision.h
is set to int
by default,
but to int64_t
if the preprocessor variable INTEGER_64
is defined.
/* qpat.c */
/* Full test for the QPA C interface using C sparse matrix indexing */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_qpa.h"
#ifdef REAL_128
#include <quadmath.h>
#endif
int main(void) {
// Derived types
void *data;
struct qpa_control_type control;
struct qpa_inform_type inform;
// Set problem data
ipc_ n = 3; // dimension
ipc_ m = 2; // number of general constraints
ipc_ H_ne = 3; // Hesssian elements
ipc_ H_row[] = {0, 1, 2 }; // row indices, NB lower triangle
ipc_ H_col[] = {0, 1, 2}; // column indices, NB lower triangle
ipc_ H_ptr[] = {0, 1, 2, 3}; // row pointers
rpc_ H_val[] = {1.0, 1.0, 1.0 }; // values
rpc_ g[] = {0.0, 2.0, 0.0}; // linear term in the objective
rpc_ f = 1.0; // constant term in the objective
rpc_ rho_g = 0.1; // penalty paramter for general constraints
rpc_ rho_b = 0.1; // penalty paramter for simple bound constraints
ipc_ A_ne = 4; // Jacobian elements
ipc_ A_row[] = {0, 0, 1, 1}; // row indices
ipc_ A_col[] = {0, 1, 1, 2}; // column indices
ipc_ A_ptr[] = {0, 2, 4}; // row pointers
rpc_ A_val[] = {2.0, 1.0, 1.0, 1.0 }; // values
rpc_ c_l[] = {1.0, 2.0}; // constraint lower bound
rpc_ c_u[] = {2.0, 2.0}; // constraint upper bound
rpc_ x_l[] = {-1.0, - INFINITY, - INFINITY}; // variable lower bound
rpc_ x_u[] = {1.0, INFINITY, 2.0}; // variable upper bound
// Set output storage
rpc_ c[m]; // constraint values
ipc_ x_stat[n]; // variable status
ipc_ c_stat[m]; // constraint status
char st = ' ';
ipc_ status;
printf(" C sparse matrix indexing\n\n");
printf(" basic tests of qp storage formats\n\n");
for( ipc_ d=1; d <= 7; d++){
// Initialize QPA
qpa_initialize( &data, &control, &status );
strcpy(control.symmetric_linear_solver, "sytr ") ;
// Set user-defined control options
control.f_indexing = false; // C sparse matrix indexing
// Start from 0
rpc_ x[] = {0.0,0.0,0.0};
rpc_ y[] = {0.0,0.0};
rpc_ z[] = {0.0,0.0,0.0};
switch(d){
case 1: // sparse co-ordinate storage
st = 'C';
qpa_import( &control, &data, &status, n, m,
"coordinate", H_ne, H_row, H_col, NULL,
"coordinate", A_ne, A_row, A_col, NULL );
qpa_solve_qp( &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 );
break;
printf(" case %1" i_ipc_ " break\n",d);
case 2: // sparse by rows
st = 'R';
qpa_import( &control, &data, &status, n, m,
"sparse_by_rows", H_ne, NULL, H_col, H_ptr,
"sparse_by_rows", A_ne, NULL, A_col, A_ptr );
qpa_solve_qp( &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 );
break;
case 3: // dense
st = 'D';
ipc_ H_dense_ne = 6; // number of elements of H
ipc_ A_dense_ne = 6; // number of elements of A
rpc_ H_dense[] = {1.0, 0.0, 1.0, 0.0, 0.0, 1.0};
rpc_ A_dense[] = {2.0, 1.0, 0.0, 0.0, 1.0, 1.0};
qpa_import( &control, &data, &status, n, m,
"dense", H_ne, NULL, NULL, NULL,
"dense", A_ne, NULL, NULL, NULL );
qpa_solve_qp( &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 );
break;
case 4: // diagonal
st = 'L';
qpa_import( &control, &data, &status, n, m,
"diagonal", H_ne, NULL, NULL, NULL,
"sparse_by_rows", A_ne, NULL, A_col, A_ptr );
qpa_solve_qp( &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 );
break;
case 5: // scaled identity
st = 'S';
qpa_import( &control, &data, &status, n, m,
"scaled_identity", H_ne, NULL, NULL, NULL,
"sparse_by_rows", A_ne, NULL, A_col, A_ptr );
qpa_solve_qp( &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 );
break;
case 6: // identity
st = 'I';
qpa_import( &control, &data, &status, n, m,
"identity", H_ne, NULL, NULL, NULL,
"sparse_by_rows", A_ne, NULL, A_col, A_ptr );
qpa_solve_qp( &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 );
break;
case 7: // zero
st = 'Z';
qpa_import( &control, &data, &status, n, m,
"zero", H_ne, NULL, NULL, NULL,
"sparse_by_rows", A_ne, NULL, A_col, A_ptr );
qpa_solve_qp( &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 );
break;
}
qpa_information( &data, &inform, &status );
if(inform.status == 0){
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_f.h
#include "galahad_pquad_f.h"
#else
printf("%c:%6" i_ipc_ " iterations. Optimal objective "
"value = %.2f status = %1" i_ipc_ "\n",
st, inform.iter, inform.obj, inform.status);
#endif
}else{
printf("%c: QPA_solve exit status = %1" i_ipc_ "\n",
st, inform.status);
}
//printf("x: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", x[i]);
//printf("\n");
//printf("gradient: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", g[i]);
//printf("\n");
// Delete internal workspace
qpa_terminate( &data, &control, &inform );
}
printf("\n basic tests of l_1 qp storage formats\n\n");
qpa_initialize( &data, &control, &status );
strcpy(control.symmetric_linear_solver, "sytr ") ;
// Set user-defined control options
control.f_indexing = false; // C sparse matrix indexing
// Start from 0
rpc_ x[] = {0.0,0.0,0.0};
rpc_ y[] = {0.0,0.0};
rpc_ z[] = {0.0,0.0,0.0};
// solve the l_1qp problem
qpa_import( &control, &data, &status, n, m,
"coordinate", H_ne, H_row, H_col, NULL,
"coordinate", A_ne, A_row, A_col, NULL );
qpa_solve_l1qp( &data, &status, n, m, H_ne, H_val, g, f, rho_g, rho_b,
A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
x_stat, c_stat );
qpa_information( &data, &inform, &status );
if(inform.status == 0){
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_f.h
#include "galahad_pquad_f.h"
#else
printf("%c:%6" i_ipc_ " iterations. Optimal objective "
"value = %.2f status = %1" i_ipc_ "\n",
st, inform.iter, inform.obj, inform.status);
#endif
}else{
printf("%c: QPA_solve exit status = %1" i_ipc_ "\n", st, inform.status);
}
// Start from 0
for( ipc_ i=0; i <= n-1; i++) x[i] = 0.0;
for( ipc_ i=0; i <= m-1; i++) y[i] = 0.0;
for( ipc_ i=0; i <= n-1; i++) z[i] = 0.0;
// solve the bound constrained l_1qp problem
qpa_import( &control, &data, &status, n, m,
"coordinate", H_ne, H_row, H_col, NULL,
"coordinate", A_ne, A_row, A_col, NULL );
qpa_solve_bcl1qp( &data, &status, n, m, H_ne, H_val, g, f, rho_g,
A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
x_stat, c_stat );
qpa_information( &data, &inform, &status );
if(inform.status == 0){
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_f.h
#include "galahad_pquad_f.h"
#else
printf("%c:%6" i_ipc_ " iterations. Optimal objective "
"value = %.2f status = %1" i_ipc_ "\n",
st, inform.iter, inform.obj, inform.status);
#endif
}else{
printf("%c: QPA_solve exit status = %1" i_ipc_ "\n", st, inform.status);
}
// Delete internal workspace
qpa_terminate( &data, &control, &inform );
}
This is the same example, but now fortran-style indexing is used; the code is available in $GALAHAD/src/qpa/C/qpatf.c .
/* qpatf.c */
/* Full test for the QPA C interface using Fortran sparse matrix indexing */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_qpa.h"
#ifdef REAL_128
#include <quadmath.h>
#endif
int main(void) {
// Derived types
void *data;
struct qpa_control_type control;
struct qpa_inform_type inform;
// Set problem data
ipc_ n = 3; // dimension
ipc_ m = 2; // number of general constraints
ipc_ H_ne = 3; // Hesssian elements
ipc_ H_row[] = {1, 2, 3 }; // row indices, NB lower triangle
ipc_ H_col[] = {1, 2, 3}; // column indices, NB lower triangle
ipc_ H_ptr[] = {1, 2, 3, 4}; // row pointers
rpc_ H_val[] = {1.0, 1.0, 1.0 }; // values
rpc_ g[] = {0.0, 2.0, 0.0}; // linear term in the objective
rpc_ f = 1.0; // constant term in the objective
rpc_ rho_g = 0.1; // penalty paramter for general constraints
rpc_ rho_b = 0.1; // penalty paramter for simple bound constraints
ipc_ A_ne = 4; // Jacobian elements
ipc_ A_row[] = {1, 1, 2, 2}; // row indices
ipc_ A_col[] = {1, 2, 2, 3}; // column indices
ipc_ A_ptr[] = {1, 3, 5}; // row pointers
rpc_ A_val[] = {2.0, 1.0, 1.0, 1.0 }; // values
rpc_ c_l[] = {1.0, 2.0}; // constraint lower bound
rpc_ c_u[] = {2.0, 2.0}; // constraint upper bound
rpc_ x_l[] = {-1.0, - INFINITY, - INFINITY}; // variable lower bound
rpc_ x_u[] = {1.0, INFINITY, 2.0}; // variable upper bound
// Set output storage
rpc_ c[m]; // constraint values
ipc_ x_stat[n]; // variable status
ipc_ c_stat[m]; // constraint status
char st = ' ';
ipc_ status;
printf(" Fortran sparse matrix indexing\n\n");
printf(" basic tests of qp storage formats\n\n");
for( ipc_ d=1; d <= 7; d++){
// Initialize QPA
qpa_initialize( &data, &control, &status );
strcpy(control.symmetric_linear_solver, "sytr ") ;
// Set user-defined control options
control.f_indexing = true; // Fortran sparse matrix indexing
// Start from 0
rpc_ x[] = {0.0,0.0,0.0};
rpc_ y[] = {0.0,0.0};
rpc_ z[] = {0.0,0.0,0.0};
switch(d){
case 1: // sparse co-ordinate storage
st = 'C';
qpa_import( &control, &data, &status, n, m,
"coordinate", H_ne, H_row, H_col, NULL,
"coordinate", A_ne, A_row, A_col, NULL );
qpa_solve_qp( &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 );
break;
printf(" case %1" i_ipc_ " break\n",d);
case 2: // sparse by rows
st = 'R';
qpa_import( &control, &data, &status, n, m,
"sparse_by_rows", H_ne, NULL, H_col, H_ptr,
"sparse_by_rows", A_ne, NULL, A_col, A_ptr );
qpa_solve_qp( &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 );
break;
case 3: // dense
st = 'D';
ipc_ H_dense_ne = 6; // number of elements of H
ipc_ A_dense_ne = 6; // number of elements of A
rpc_ H_dense[] = {1.0, 0.0, 1.0, 0.0, 0.0, 1.0};
rpc_ A_dense[] = {2.0, 1.0, 0.0, 0.0, 1.0, 1.0};
qpa_import( &control, &data, &status, n, m,
"dense", H_ne, NULL, NULL, NULL,
"dense", A_ne, NULL, NULL, NULL );
qpa_solve_qp( &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 );
break;
case 4: // diagonal
st = 'L';
qpa_import( &control, &data, &status, n, m,
"diagonal", H_ne, NULL, NULL, NULL,
"sparse_by_rows", A_ne, NULL, A_col, A_ptr );
qpa_solve_qp( &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 );
break;
case 5: // scaled identity
st = 'S';
qpa_import( &control, &data, &status, n, m,
"scaled_identity", H_ne, NULL, NULL, NULL,
"sparse_by_rows", A_ne, NULL, A_col, A_ptr );
qpa_solve_qp( &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 );
break;
case 6: // identity
st = 'I';
qpa_import( &control, &data, &status, n, m,
"identity", H_ne, NULL, NULL, NULL,
"sparse_by_rows", A_ne, NULL, A_col, A_ptr );
qpa_solve_qp( &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 );
break;
case 7: // zero
st = 'Z';
qpa_import( &control, &data, &status, n, m,
"zero", H_ne, NULL, NULL, NULL,
"sparse_by_rows", A_ne, NULL, A_col, A_ptr );
qpa_solve_qp( &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 );
break;
}
qpa_information( &data, &inform, &status );
if(inform.status == 0){
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_f.h
#include "galahad_pquad_f.h"
#else
printf("%c:%6" i_ipc_ " iterations. Optimal objective "
"value = %.2f status = %1" i_ipc_ "\n",
st, inform.iter, inform.obj, inform.status);
#endif
}else{
printf("%c: QPA_solve exit status = %1" i_ipc_ "\n",
st, inform.status);
}
//printf("x: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", x[i]);
//printf("\n");
//printf("gradient: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", g[i]);
//printf("\n");
// Delete internal workspace
qpa_terminate( &data, &control, &inform );
}
printf("\n basic tests of l_1 qp storage formats\n\n");
qpa_initialize( &data, &control, &status );
strcpy(control.symmetric_linear_solver, "sytr ") ;
// Set user-defined control options
control.f_indexing = true; // Fortran sparse matrix indexing
// Start from 0
rpc_ x[] = {0.0,0.0,0.0};
rpc_ y[] = {0.0,0.0};
rpc_ z[] = {0.0,0.0,0.0};
// solve the l_1qp problem
qpa_import( &control, &data, &status, n, m,
"coordinate", H_ne, H_row, H_col, NULL,
"coordinate", A_ne, A_row, A_col, NULL );
qpa_solve_l1qp( &data, &status, n, m, H_ne, H_val, g, f, rho_g, rho_b,
A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
x_stat, c_stat );
qpa_information( &data, &inform, &status );
if(inform.status == 0){
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_f.h
#include "galahad_pquad_f.h"
#else
printf("%c:%6" i_ipc_ " iterations. Optimal objective "
"value = %.2f status = %1" i_ipc_ "\n",
st, inform.iter, inform.obj, inform.status);
#endif
}else{
printf("%c: QPA_solve exit status = %1" i_ipc_ "\n", st, inform.status);
}
// Start from 0
for( ipc_ i=0; i <= n-1; i++) x[i] = 0.0;
for( ipc_ i=0; i <= m-1; i++) y[i] = 0.0;
for( ipc_ i=0; i <= n-1; i++) z[i] = 0.0;
// solve the bound constrained l_1qp problem
qpa_import( &control, &data, &status, n, m,
"coordinate", H_ne, H_row, H_col, NULL,
"coordinate", A_ne, A_row, A_col, NULL );
qpa_solve_bcl1qp( &data, &status, n, m, H_ne, H_val, g, f, rho_g,
A_ne, A_val, c_l, c_u, x_l, x_u, x, c, y, z,
x_stat, c_stat );
qpa_information( &data, &inform, &status );
if(inform.status == 0){
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_f.h
#include "galahad_pquad_f.h"
#else
printf("%c:%6" i_ipc_ " iterations. Optimal objective "
"value = %.2f status = %1" i_ipc_ "\n",
st, inform.iter, inform.obj, inform.status);
#endif
}else{
printf("%c: QPA_solve exit status = %1" i_ipc_ "\n", st, inform.status);
}
// Delete internal workspace
qpa_terminate( &data, &control, &inform );
}