GALAHAD BQP package#

purpose#

The bqp package uses a preconditioned, projected-gradient method to solve a given bound-constrained convex quadratic program. The aim is to minimize the quadratic objective function

\[q(x) = f + g^T x + \frac{1}{2} x^T H x,\]

subject to the simple bounds

\[x_l \leq x \leq x_u,\]

where $H$ is a given $n$ by $n$ symmetric postive-semi-definite matrix, $g$ is a vector, $f$ is a scalar, and any of the components of the vectors $x_l$ or $x_u$ may be infinite. The method offers the choice of direct and iterative solution of the key regularization subproblems, and is most suitable for problems involving a large number of unknowns $x$.

See Section 4 of $GALAHAD/doc/bqp.pdf for a brief description of the method employed and other details.

terminology#

Any required solution $x$ necessarily satisfies the primal optimality conditions

\[x_l \leq x \leq x_u,\]

the dual optimality conditions

\[H x + g = z, \;\; z = z_l + z_u, z_l \geq 0 \;\;\mbox{and}\;\; z_u \leq 0,\]

and the complementary slackness conditions

\[(x -x_l )^{T} z_l = 0 \;\;\mbox{and}\;\;(x -x_u )^{T} z_u = 0,\]

where the vector $z$ is known as the dual variables for the bounds, and where the vector inequalities hold component-wise.

method#

Projected-gradient methods iterate towards a point that satisfies these conditions by ultimately aiming to satisfy $H x + g = z$ and $z = z_l + z_u$, while satifying the remaining optimality conditions at each stage. Appropriate norms of the amounts by which the optimality conditions fail to be satisfied are known as the primal and dual infeasibility, and the violation of complementary slackness, respectively.

The method is iterative. Each iteration proceeds in two stages. Firstly, the so-called generalized Cauchy point for the quadratic objective is found. (The purpose of this point is to ensure that the algorithm converges and that the set of bounds which are satisfied as equations at the solution is rapidly identified.) Thereafter an improvement to the objective is sought using either a direct-matrix or truncated conjugate-gradient algorithm.

reference#

This is a specialised version of the method presented in

A. R. Conn, N. I. M. Gould and Ph. L. Toint, Global convergence of a class of trust region algorithms for optimization with simple bounds. SIAM Journal on Numerical Analysis 25 (1988) 433-460.

matrix storage#

The symmetric $n$ by $n$ matrix $H$ may 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 bqp package must be called in the following order:

bqp_initialize - provide default control parameters and set up initial data structures
bqp_read_specfile (optional) - override control values by reading replacement values from a file
set up problem data structures and fixed values by caling one of
- bqp_import - in the case that $H$ is explicitly available
- bqp_import_without_h - in the case that only the effect of applying $H$ to a vector is possible
bqp_reset_control (optional) - possibly change control parameters if a sequence of problems are being solved
solve the problem by calling one of
- bqp_solve_given_h - solve the problem using values of $H$
- bqp_solve_reverse_h_prod - solve the problem by returning to the caller for products of $H$ with specified vectors
bqp_information (optional) - recover information about the solution and solution process
bqp_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 bqp_control_type;
struct bqp_inform_type;
struct bqp_time_type;

// function calls

void bqp_initialize(void **data, struct bqp_control_type* control, ipc_ *status);
void bqp_read_specfile(struct bqp_control_type* control, const char specfile[]);

void bqp_import(
    struct bqp_control_type* control,
    void **data,
    ipc_ *status,
    ipc_ n,
    const char H_type[],
    ipc_ ne,
    const ipc_ H_row[],
    const ipc_ H_col[],
    const ipc_ H_ptr[]
);

void bqp_import_without_h(
    struct bqp_control_type* control,
    void **data,
    ipc_ *status,
    ipc_ n
);

void bqp_reset_control(
    struct bqp_control_type* control,
    void **data,
    ipc_ *status
);

void bqp_solve_given_h(
    void **data,
    ipc_ *status,
    ipc_ n,
    ipc_ h_ne,
    const rpc_ H_val[],
    const rpc_ g[],
    const rpc_ f,
    const rpc_ x_l[],
    const rpc_ x_u[],
    rpc_ x[],
    rpc_ z[],
    ipc_ x_stat[]
);

void bqp_solve_reverse_h_prod(
    void **data,
    ipc_ *status,
    ipc_ n,
    const rpc_ g[],
    const rpc_ f,
    const rpc_ x_l[],
    const rpc_ x_u[],
    rpc_ x[],
    rpc_ z[],
    ipc_ x_stat[],
    rpc_ v[],
    const rpc_ prod[],
    ipc_ nz_v[],
    ipc_ *nz_v_start,
    ipc_ *nz_v_end,
    const ipc_ nz_prod[],
    ipc_ nz_prod_end
);

void bqp_information(void **data, struct bqp_inform_type* inform, ipc_ *status);

void bqp_terminate(
    void **data,
    struct bqp_control_type* control,
    struct bqp_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 SINGLE.

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 bqp_initialize(void **data, struct bqp_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 bqp_control_type)

status

is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are (currently):

0
The initialization was successful.

void bqp_read_specfile(struct bqp_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/bqp/BQP.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/bqp.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 bqp_control_type)
specfile	is a character string containing the name of the specification file

void bqp_import(
    struct bqp_control_type* control,
    void **data,
    ipc_ *status,
    ipc_ n,
    const char H_type[],
    ipc_ ne,
    const ipc_ H_row[],
    const ipc_ H_col[],
    const ipc_ H_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 bqp_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: 1 The import was successful, and the package is ready for the solve phase -1 An allocation error occurred. A message indicating the offending array is written on unit control.error, and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -2 A deallocation error occurred. A message indicating the offending array is written on unit control.error and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -3 The restriction n > 0 or requirement that type contains its relevant string ‘dense’, ‘coordinate’, ‘sparse_by_rows’ or ‘diagonal’ has been violated.
n	is a scalar variable of type ipc_, that holds the number of variables.
H_type	is a one-dimensional array of type char that specifies the symmetric storage scheme used for the Hessian. It should be one of ‘coordinate’, ‘sparse_by_rows’, ‘dense’, ‘diagonal’ or ‘absent’, the latter if access to the Hessian is via matrix-vector products; lower or upper case variants are allowed.
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 three schemes.
H_row	is a one-dimensional array of size 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 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 or diagonal 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

void bqp_import_without_h(
    struct bqp_control_type* control,
    void **data,
    ipc_ *status,
    ipc_ n
)

Import problem data into internal storage prior to solution.

Parameters:

control	is a struct whose members provide control paramters for the remaining prcedures (see bqp_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: 1 The import was successful, and the package is ready for the solve phase -1 An allocation error occurred. A message indicating the offending array is written on unit control.error, and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -2 A deallocation error occurred. A message indicating the offending array is written on unit control.error and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -3. The restriction n > 0 has been violated.
n	is a scalar variable of type ipc_, that holds the number of variables.

void bqp_reset_control(
    struct bqp_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 bqp_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:

1
The import was successful, and the package is ready for the solve phase

void bqp_solve_given_h(
    void **data,
    ipc_ *status,
    ipc_ n,
    ipc_ h_ne,
    const rpc_ H_val[],
    const rpc_ g[],
    const rpc_ f,
    const rpc_ x_l[],
    const rpc_ x_u[],
    rpc_ x[],
    rpc_ z[],
    ipc_ x_stat[]
)

Solve the bound-constrained quadratic program when the Hessian $H$ is available.

Parameters:

data	holds private internal data
status	is a scalar variable of type ipc_, that gives the entry and exit status from the package. On initial entry, status must be set to 1. Possible exit values are: 0 The run was successful. -1 An allocation error occurred. A message indicating the offending array is written on unit control.error, and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -2 A deallocation error occurred. A message indicating the offending array is written on unit control.error and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -3 The restriction n > 0 or requirement that a type contains its relevant string ‘dense’, ‘coordinate’, ‘sparse_by_rows’ or ‘diagonal’ has been violated. -4 The simple-bound constraints are inconsistent. -9 The analysis phase of the factorization failed; the return status from the factorization package is given in the component inform.factor_status -10 The factorization failed; the return status from the factorization package is given in the component inform.factor_status. -11 The solution of a set of linear equations using factors from the factorization package failed; the return status from the factorization package is given in the component inform.factor_status. -16 The problem is so ill-conditioned that further progress is impossible. -17 The step is too small to make further impact. -18 Too many iterations have been performed. This may happen if control.maxit is too small, but may also be symptomatic of a badly scaled problem. -19 The CPU time limit has been reached. This may happen if control.cpu_time_limit is too small, but may also be symptomatic of a badly scaled problem. -20 The Hessian matrix $H$ appears to be indefinite. specified. -23 An entry from the strict upper triangle of $H$ has been
n	is a scalar variable of type ipc_, that holds the number of variables
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.
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$.
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 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.

void bqp_solve_reverse_h_prod(
    void **data,
    ipc_ *status,
    ipc_ n,
    const rpc_ g[],
    const rpc_ f,
    const rpc_ x_l[],
    const rpc_ x_u[],
    rpc_ x[],
    rpc_ z[],
    ipc_ x_stat[],
    rpc_ v[],
    const rpc_ prod[],
    ipc_ nz_v[],
    ipc_ *nz_v_start,
    ipc_ *nz_v_end,
    const ipc_ nz_prod[],
    ipc_ nz_prod_end
)

Solve the bound-constrained quadratic program when the products of the Hessian $H$ with specified vectors may be computed by the calling program.

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: The run was successful. -1 An allocation error occurred. A message indicating the offending array is written on unit control.error, and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -2 A deallocation error occurred. A message indicating the offending array is written on unit control.error and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -3 The restriction n > 0 or requirement that a type contains its relevant string ‘dense’, ‘coordinate’, ‘sparse_by_rows’ or ‘diagonal’ has been violated. -4 The simple-bound constraints are inconsistent. -9 The analysis phase of the factorization failed; the return status from the factorization package is given in the component inform.factor_status -10 The factorization failed; the return status from the factorization package is given in the component inform.factor_status. -11 The solution of a set of linear equations using factors from the factorization package failed; the return status from the factorization package is given in the component inform.factor_status. -16 The problem is so ill-conditioned that further progress is impossible. -17 The step is too small to make further impact. -18 Too many iterations have been performed. This may happen if control.maxit is too small, but may also be symptomatic of a badly scaled problem. -19 The CPU time limit has been reached. This may happen if control.cpu_time_limit is too small, but may also be symptomatic of a badly scaled problem. -20 The Hessian matrix $H$ appears to be indefinite. specified. -23 An entry from the strict upper triangle of $H$ has been specified. 2 The product $Hv$ of the Hessian $H$ with a given output vector $v$ is required from the user. The vector $v$ will be stored in v and the product $Hv$ must be returned in prod, and bqp_solve_reverse_h_prod re-entered with all other arguments unchanged. 3 The product $Hv$ of the Hessian H with a given output vector $v$ is required from the user. Only components nz_v[nz_v_start-1:nz_v_end-1] of the vector $v$ stored in v are nonzero. The resulting product $Hv$ must be placed in prod, and bqp_solve_reverse_h_prod re-entered with all other arguments unchanged. 4 The product $Hv$ of the Hessian H with a given output vector $v$ is required from the user. Only components nz_v[nz_v_start-1:nz_v_end-1] of the vector $v$ stored in v are nonzero. The resulting nonzeros in the product $Hv$ must be placed in their appropriate comnponents of prod, while a list of indices of the nonzeros placed in nz_prod[0 : nz_prod_end-1]. bqp_solve_reverse_h_prod should then be re-entered with all other arguments unchanged. Typically v will be very sparse (i.e., nz_p_end-nz_p_start will be small).
n	is a scalar variable of type ipc_, that holds the number of variables
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.
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$.
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 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.
v	is a one-dimensional array of size n and type rpc_, that is used for reverse communication (see status=2-4 above for details)
prod	is a one-dimensional array of size n and type rpc_, that is used for reverse communication (see status=2-4 above for details)
nz_v	is a one-dimensional array of size n and type ipc_, that is used for reverse communication (see status=3-4 above for details)
nz_v_start	is a scalar of type ipc_, that is used for reverse communication (see status=3-4 above for details)
nz_v_end	is a scalar of type ipc_, that is used for reverse communication (see status=3-4 above for details)
nz_prod	is a one-dimensional array of size n and type ipc_, that is used for reverse communication (see status=4 above for details)
nz_prod_end	is a scalar of type ipc_, that is used for reverse communication (see status=4 above for details)

void bqp_information(void **data, struct bqp_inform_type* inform, ipc_ *status)

Provides output information

Parameters:

data

holds private internal data

inform

is a struct containing output information (see bqp_inform_type)

status

is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are (currently):

0
The values were recorded successfully

void bqp_terminate(
    void **data,
    struct bqp_control_type* control,
    struct bqp_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a struct containing control information (see bqp_control_type)
inform	is a struct containing output information (see bqp_inform_type)

available structures#

bqp_control_type structure#

#include <galahad_bqp.h>

struct bqp_control_type {
    // components

    bool f_indexing;
    ipc_ error;
    ipc_ out;
    ipc_ print_level;
    ipc_ start_print;
    ipc_ stop_print;
    ipc_ print_gap;
    ipc_ maxit;
    ipc_ cold_start;
    ipc_ ratio_cg_vs_sd;
    ipc_ change_max;
    ipc_ cg_maxit;
    ipc_ sif_file_device;
    rpc_ infinity;
    rpc_ stop_p;
    rpc_ stop_d;
    rpc_ stop_c;
    rpc_ identical_bounds_tol;
    rpc_ stop_cg_relative;
    rpc_ stop_cg_absolute;
    rpc_ zero_curvature;
    rpc_ cpu_time_limit;
    bool exact_arcsearch;
    bool space_critical;
    bool deallocate_error_fatal;
    bool generate_sif_file;
    char sif_file_name[31];
    char prefix[31];
    struct sbls_control_type sbls_control;
};

detailed documentation#

control derived type as a C struct

components#

bool f_indexing

use C or Fortran sparse matrix indexing

ipc_ error

unit number for error and warning diagnostics

ipc_ out

general output unit number

ipc_ print_level

the level of output required

ipc_ start_print

on which iteration to start printing

ipc_ stop_print

on which iteration to stop printing

ipc_ print_gap

how many iterations between printing

ipc_ maxit

how many iterations to perform (-ve reverts to HUGE(1)-1)

ipc_ cold_start

cold_start should be set to 0 if a warm start is required (with variable assigned according to B_stat, see below), and to any other value if the values given in prob.X suffice

ipc_ ratio_cg_vs_sd

the ratio of how many iterations use CG rather steepest descent

ipc_ change_max

the maximum number of per-iteration changes in the working set permitted when allowing CG rather than steepest descent

ipc_ cg_maxit

how many CG iterations to perform per BQP iteration (-ve reverts to n+1)

ipc_ sif_file_device

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_ stop_p

the required accuracy for the primal infeasibility

rpc_ stop_d

the required accuracy for the dual infeasibility

rpc_ stop_c

the required accuracy for the complementary slackness

rpc_ identical_bounds_tol

any pair of constraint bounds (x_l,x_u) that are closer than i dentical_bounds_tol will be reset to the average of their values

rpc_ 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)

rpc_ stop_cg_absolute

see stop_cg_relative

rpc_ zero_curvature

threshold below which curvature is regarded as zero

rpc_ cpu_time_limit

the maximum CPU time allowed (-ve = no limit)

bool exact_arcsearch

exact_arcsearch is true if an exact arcsearch is required, and false if approximation suffices

bool space_critical

if space_critical is true, every effort will be made to use as little space as possible. This may result in longer computation times

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, a SIF file describing the current problem will be generated

char sif_file_name[31]

name (max 30 characters) of generated SIF file containing input problem

char prefix[31]

all output lines will be prefixed by a string (max 30 characters) prefix(2:LEN(TRIM(.prefix))-1) where prefix contains the required string enclosed in quotes, e.g. “string” or ‘string’

struct sbls_control_type sbls_control

control parameters for SBLS

bqp_time_type structure#

#include <galahad_bqp.h>

struct bqp_time_type {
    // components

    spc_ total;
    spc_ analyse;
    spc_ factorize;
    spc_ solve;
};

detailed documentation#

time derived type as a C struct

components#

spc_ total

total time

spc_ analyse

time for the analysis phase

spc_ factorize

time for the factorization phase

spc_ solve

time for the linear solution phase

bqp_inform_type structure#

#include <galahad_bqp.h>

struct bqp_inform_type {
    // components

    ipc_ status;
    ipc_ alloc_status;
    ipc_ factorization_status;
    ipc_ iter;
    ipc_ cg_iter;
    rpc_ obj;
    rpc_ norm_pg;
    char bad_alloc[81];
    struct bqp_time_type time;
    struct sbls_inform_type sbls_inform;
};

detailed documentation#

inform derived type as a C struct

components#

ipc_ status

reported return status:

0

success
-1

allocation error
-2

deallocation error
-3

matrix data faulty (.n < 1, .ne < 0)
-20

alegedly +ve definite matrix is not

ipc_ alloc_status

Fortran STAT value after allocate failure.

ipc_ factorization_status

status return from factorization

ipc_ iter

number of iterations required

ipc_ cg_iter

number of CG iterations required

rpc_ obj

current value of the objective function

rpc_ norm_pg

current value of the projected gradient

char bad_alloc[81]

name of array which provoked an allocate failure

struct bqp_time_type time

times for various stages

struct sbls_inform_type sbls_inform

inform values from SBLS

Index

example calls#

This is an example of how to use the package to solve a bound-constrained QP; the code is available in $GALAHAD/src/bqp/C/bqpt.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.

/* bqpt.c */
/* Full test for the BQP C interface using C sparse matrix indexing */

#include <stdio.h>
#include <math.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_bqp.h"

int main(void) {

    // Derived types
    void *data;
    struct bqp_control_type control;
    struct bqp_inform_type inform;

    // Set problem data
    ipc_ n = 10; // dimension
    ipc_ H_ne = 2 * n - 1; // Hesssian elements, NB lower triangle
    ipc_ H_dense_ne = n * ( n + 1 ) / 2; // dense Hessian elements
    ipc_ H_row[H_ne]; // row indices,
    ipc_ H_col[H_ne]; // column indices
    ipc_ H_ptr[n+1];  // row pointers
    rpc_ H_val[H_ne]; // values
    rpc_ H_dense[H_dense_ne]; // dense values
    rpc_ H_diag[n];   // diagonal values
    rpc_ g[n];  // linear term in the objective
    rpc_ f = 1.0;  // constant term in the objective
    rpc_ x_l[n]; // variable lower bound
    rpc_ x_u[n]; // variable upper bound
    rpc_ x[n]; // variables
    rpc_ z[n]; // dual variables

    // Set output storage
    ipc_ x_stat[n]; // variable status
    char st = ' ';
    ipc_ i, l, status;

    g[0] = 2.0;
    for( ipc_ i = 1; i < n; i++) g[i] = 0.0;
    x_l[0] = -1.0;
    for( ipc_ i = 1; i < n; i++) x_l[i] = - INFINITY;
    x_u[0] = 1.0;
    x_u[1] = INFINITY;
    for( ipc_ i = 2; i < n; i++) x_u[i] = 2.0;

    // H = tridiag(2,1), H_dense = diag(2)

    l = 0 ;
    H_ptr[0] = l;
    H_row[l] = 0; H_col[l] = 0; H_val[l] = 2.0;
    for( ipc_ i = 1; i < n; i++)
    {
      l = l + 1;
      H_ptr[i] = l;
      H_row[l] = i; H_col[l] = i - 1; H_val[l] = 1.0;
      l = l + 1;
      H_row[l] = i; H_col[l] = i; H_val[l] = 2.0;
    }
    H_ptr[n] = l + 1;

    l = - 1 ;
    for( ipc_ i = 0; i < n; i++)
    {
      H_diag[i] = 2.0;
      for( ipc_ j = 0; j <= i; j++)
      {
        l = l + 1;
        if ( j < i - 1 ) {
          H_dense[l] = 0.0;
        }
        else if ( j == i - 1 ) {
          H_dense[l] = 1.0;
        }
        else {
          H_dense[l] = 2.0;
        }
      }
    }

    printf(" C sparse matrix indexing\n\n");

    printf(" basic tests of bqp storage formats\n\n");

    for( ipc_ d=1; d <= 4; d++){

        // Initialize BQP
        bqp_initialize( &data, &control, &status );

        // Set user-defined control options
        control.f_indexing = false; // C sparse matrix indexing

        // Start from 0
        for( ipc_ i = 0; i < n; i++) x[i] = 0.0;
        for( ipc_ i = 0; i < n; i++) z[i] = 0.0;

        switch(d){
            case 1: // sparse co-ordinate storage
                st = 'C';
                bqp_import( &control, &data, &status, n,
                            "coordinate", H_ne, H_row, H_col, NULL );
                bqp_solve_given_h( &data, &status, n, H_ne, H_val, g, f,
                                   x_l, x_u, x, z, x_stat );
                break;
            printf(" case %1" i_ipc_ " break\n",d);
            case 2: // sparse by rows
                st = 'R';
                bqp_import( &control, &data, &status, n,
                             "sparse_by_rows", H_ne, NULL, H_col, H_ptr );
                bqp_solve_given_h( &data, &status, n, H_ne, H_val, g, f,
                                   x_l, x_u, x, z, x_stat );
                break;
            case 3: // dense
                st = 'D';
                bqp_import( &control, &data, &status, n,
                             "dense", H_dense_ne, NULL, NULL, NULL );
                bqp_solve_given_h( &data, &status, n, H_dense_ne, H_dense,
                                   g, f, x_l, x_u, x, z, x_stat );
                break;
            case 4: // diagonal
                st = 'L';
                bqp_import( &control, &data, &status, n,
                             "diagonal", H_ne, NULL, NULL, NULL );
                bqp_solve_given_h( &data, &status, n, n, H_diag, g, f,
                                   x_l, x_u, x, z, x_stat );
                break;
            }
        bqp_information( &data, &inform, &status );

        if(inform.status == 0){
            printf("%c:%6" i_ipc_ " iterations. Optimal objective value = %5.2f status = %1" i_ipc_ "\n",
                   st, inform.iter, inform.obj, inform.status);
        }else{
            printf("%c: BQP_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
        bqp_terminate( &data, &control, &inform );
    }

    printf("\n tests reverse-communication options\n\n");

    // reverse-communication input/output
    ipc_ nz_v_start, nz_v_end, nz_prod_end;
    ipc_ nz_v[n], nz_prod[n], mask[n];
    rpc_ v[n], prod[n];

    nz_prod_end = 0;

    // Initialize BQP
    bqp_initialize( &data, &control, &status );
   // control.print_level = 1;

   // Set user-defined control options
    control.f_indexing = false; // C sparse matrix indexing

    // Start from 0
    for( ipc_ i = 0; i < n; i++) x[i] = 0.0;
    for( ipc_ i = 0; i < n; i++) z[i] = 0.0;

    st = 'I';
    for( ipc_ i = 0; i < n; i++) mask[i] = 0;
    bqp_import_without_h( &control, &data, &status, n ) ;
    while(true){ // reverse-communication loop
        bqp_solve_reverse_h_prod( &data, &status, n, g, f, x_l, x_u,
                                  x, z, x_stat, v, prod,
                                  nz_v, &nz_v_start, &nz_v_end,
                                  nz_prod, nz_prod_end );
        if(status == 0){ // successful termination
            break;
        }else if(status < 0){ // error exit
            break;
        }else if(status == 2){ // evaluate Hv
          prod[0] = 2.0 * v[0] + v[1];
          for( ipc_ i = 1; i < n-1; i++) prod[i] = 2.0 * v[i] + v[i-1] + v[i+1];
          prod[n-1] = 2.0 * v[n-1] + v[n-2];
        }else if(status == 3){ // evaluate Hv for sparse v
          for( ipc_ i = 0; i < n; i++) prod[i] = 0.0;
          for( ipc_ l = nz_v_start - 1; l < nz_v_end; l++){
             i = nz_v[l];
             if (i > 0) prod[i-1] = prod[i-1] + v[i];
             prod[i] = prod[i] + 2.0 * v[i];
             if (i < n-1) prod[i+1] = prod[i+1] + v[i];
           }
        }else if(status == 4){ // evaluate sarse Hv for sparse v
          nz_prod_end = 0;
          for( ipc_ l = nz_v_start - 1; l < nz_v_end; l++){
             i = nz_v[l];
             if (i > 0){
               if (mask[i-1] == 0){
                 mask[i-1] = 1;
                 nz_prod[nz_prod_end] = i - 1;
                 nz_prod_end = nz_prod_end + 1;
                 prod[i-1] = v[i];
               }else{
                 prod[i-1] = prod[i-1] + v[i];
               }
             }
             if (mask[i] == 0){
               mask[i] = 1;
               nz_prod[nz_prod_end] = i;
               nz_prod_end = nz_prod_end + 1;
               prod[i] = 2.0 * v[i];
             }else{
               prod[i] = prod[i] + 2.0 * v[i];
             }
             if (i < n-1){
               if (mask[i+1] == 0){
                 mask[i+1] = 1;
                 nz_prod[nz_prod_end] = i + 1;
                 nz_prod_end = nz_prod_end + 1;
                 prod[i+1] = prod[i+1] + v[i];
               }else{
                 prod[i+1] = prod[i+1] + v[i];
               }
             }
          }
          for( ipc_ l = 0; l < nz_prod_end; l++) mask[nz_prod[l]] = 0;
        }else{
            printf(" the value %1" i_ipc_ " of status should not occur\n", status);
            break;
        }
    }

    // Record solution information
    bqp_information( &data, &inform, &status );

    // Print solution details
    if(inform.status == 0){
        printf("%c:%6" i_ipc_ " iterations. Optimal objective value = %5.2f status = %1" i_ipc_ "\n",
               st, inform.iter, inform.obj, inform.status);
    }else{
        printf("%c: BQP_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
    bqp_terminate( &data, &control, &inform );
}

This is the same example, but now fortran-style indexing is used; the code is available in $GALAHAD/src/bqp/C/bqptf.c .

/* bqptf.c */
/* Full test for the BQP C interface using fortran sparse matrix indexing */

#include <stdio.h>
#include <math.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_bqp.h"

int main(void) {

    // Derived types
    void *data;
    struct bqp_control_type control;
    struct bqp_inform_type inform;

    // Set problem data
    ipc_ n = 10; // dimension
    ipc_ H_ne = 2 * n - 1; // Hesssian elements, NB lower triangle
    ipc_ H_dense_ne = n * ( n + 1 ) / 2; // dense Hessian elements
    ipc_ H_row[H_ne]; // row indices,
    ipc_ H_col[H_ne]; // column indices
    ipc_ H_ptr[n+1];  // row pointers
    rpc_ H_val[H_ne]; // values
    rpc_ H_dense[H_dense_ne]; // dense values
    rpc_ H_diag[n];   // diagonal values
    rpc_ g[n];  // linear term in the objective
    rpc_ f = 1.0;  // constant term in the objective
    rpc_ x_l[n]; // variable lower bound
    rpc_ x_u[n]; // variable upper bound
    rpc_ x[n]; // variables
    rpc_ z[n]; // dual variables

    // Set output storage
    ipc_ x_stat[n]; // variable status
    char st = ' ';
    ipc_ i, l, status;

    g[0] = 2.0;
    for( ipc_ i = 1; i < n; i++) g[i] = 0.0;
    x_l[0] = -1.0;
    for( ipc_ i = 1; i < n; i++) x_l[i] = - INFINITY;
    x_u[0] = 1.0;
    x_u[1] = INFINITY;
    for( ipc_ i = 2; i < n; i++) x_u[i] = 2.0;

    // H = tridiag(2,1), H_dense = diag(2)

    l = 0 ;
    H_ptr[0] = l + 1;
    H_row[l] = 1; H_col[l] = 1; H_val[l] = 2.0;
    for( ipc_ i = 1; i < n; i++)
    {
      l = l + 1;
      H_ptr[i] = l + 1;
      H_row[l] = i + 1; H_col[l] = i; H_val[l] = 1.0;
      l = l + 1;
      H_row[l] = i + 1; H_col[l] = i + 1; H_val[l] = 2.0;
    }
    H_ptr[n] = l + 2;

    l = - 1;
    for( ipc_ i = 0; i < n; i++)
    {
      H_diag[i] = 2.0;
      for( ipc_ j = 0; j <= i; j++)
      {
        l = l + 1;
        if ( j < i - 1 ) {
          H_dense[l] = 0.0;
        }
        else if ( j == i - 1 ) {
          H_dense[l] = 1.0;
        }
        else {
          H_dense[l] = 2.0;
        }
      }
    }

    printf(" fortran sparse matrix indexing\n\n");

    printf(" basic tests of bqp storage formats\n\n");

    for( ipc_ d=1; d <= 4; d++){

        // Initialize BQP
        bqp_initialize( &data, &control, &status );

        // Set user-defined control options
        control.f_indexing = true; // fortran sparse matrix indexing

        // Start from 0
        for( ipc_ i = 0; i < n; i++) x[i] = 0.0;
        for( ipc_ i = 0; i < n; i++) z[i] = 0.0;

        switch(d){
            case 1: // sparse co-ordinate storage
                st = 'C';
                bqp_import( &control, &data, &status, n,
                            "coordinate", H_ne, H_row, H_col, NULL );
                bqp_solve_given_h( &data, &status, n, H_ne, H_val, g, f,
                                   x_l, x_u, x, z, x_stat );
                break;
            printf(" case %1" i_ipc_ " break\n",d);
            case 2: // sparse by rows
                st = 'R';
                bqp_import( &control, &data, &status, n,
                             "sparse_by_rows", H_ne, NULL, H_col, H_ptr );
                bqp_solve_given_h( &data, &status, n, H_ne, H_val, g, f,
                                   x_l, x_u, x, z, x_stat );
                break;
            case 3: // dense
                st = 'D';
                bqp_import( &control, &data, &status, n,
                             "dense", H_dense_ne, NULL, NULL, NULL );
                bqp_solve_given_h( &data, &status, n, H_dense_ne, H_dense,
                                   g, f, x_l, x_u, x, z, x_stat );
                break;
            case 4: // diagonal
                st = 'L';
                bqp_import( &control, &data, &status, n,
                             "diagonal", H_ne, NULL, NULL, NULL );
                bqp_solve_given_h( &data, &status, n, n, H_diag, g, f,
                                   x_l, x_u, x, z, x_stat );
                break;
            }
        bqp_information( &data, &inform, &status );

        if(inform.status == 0){
            printf("%c:%6" i_ipc_ " iterations. Optimal objective value = %5.2f"
                   " status = %1" i_ipc_ "\n",
                   st, inform.iter, inform.obj, inform.status);
        }else{
            printf("%c: BQP_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
        bqp_terminate( &data, &control, &inform );
    }

    printf("\n tests reverse-communication options\n\n");

    // reverse-communication input/output
    ipc_ nz_v_start, nz_v_end, nz_prod_end;
    ipc_ nz_v[n], nz_prod[n], mask[n];
    rpc_ v[n], prod[n];

    nz_prod_end = 0;

    // Initialize BQP
    bqp_initialize( &data, &control, &status );
   // control.print_level = 1;

   // Set user-defined control options
    control.f_indexing = true; // fortran sparse matrix indexing

    // Start from 0
    for( ipc_ i = 0; i < n; i++) x[i] = 0.0;
    for( ipc_ i = 0; i < n; i++) z[i] = 0.0;

    st = 'I';
    for( ipc_ i = 0; i < n; i++) mask[i] = 0;
    bqp_import_without_h( &control, &data, &status, n ) ;
    while(true){ // reverse-communication loop
        bqp_solve_reverse_h_prod( &data, &status, n, g, f, x_l, x_u,
                                  x, z, x_stat, v, prod,
                                  nz_v, &nz_v_start, &nz_v_end,
                                  nz_prod, nz_prod_end );
        if(status == 0){ // successful termination
            break;
        }else if(status < 0){ // error exit
            break;
        }else if(status == 2){ // evaluate Hv
          prod[0] = 2.0 * v[0] + v[1];
          for( ipc_ i = 1; i < n-1; i++) prod[i] = 2.0 * v[i] + v[i-1] + v[i+1];
          prod[n-1] = 2.0 * v[n-1] + v[n-2];
        }else if(status == 3){ // evaluate Hv for sparse v
          for( ipc_ i = 0; i < n; i++) prod[i] = 0.0;
          for( ipc_ l = nz_v_start - 1; l < nz_v_end; l++){
             i = nz_v[l]-1;
             if (i > 0) prod[i-1] = prod[i-1] + v[i];
             prod[i] = prod[i] + 2.0 * v[i];
             if (i < n-1) prod[i+1] = prod[i+1] + v[i];
           }
        }else if(status == 4){ // evaluate sarse Hv for sparse v
          nz_prod_end = 0;
          for( ipc_ l = nz_v_start - 1; l < nz_v_end; l++){
             i = nz_v[l]-1;
             if (i > 0){
               if (mask[i-1] == 0){
                 mask[i-1] = 1;
                 nz_prod[nz_prod_end] = i - 1;
                 nz_prod_end = nz_prod_end + 1;
                 prod[i-1] = v[i];
               }else{
                 prod[i-1] = prod[i-1] + v[i];
               }
             }
             if (mask[i] == 0){
               mask[i] = 1;
               nz_prod[nz_prod_end] = i;
               nz_prod_end = nz_prod_end + 1;
               prod[i] = 2.0 * v[i];
             }else{
               prod[i] = prod[i] + 2.0 * v[i];
             }
             if (i < n-1){
               if (mask[i+1] == 0){
                 mask[i+1] = 1;
                 nz_prod[nz_prod_end] = i + 1;
                 nz_prod_end = nz_prod_end + 1;
               }
               prod[i+1] = prod[i+1] + v[i];
             }
          }
          for( ipc_ l = 0; l < nz_prod_end; l++) mask[nz_prod[l]] = 0;
        }else{
            printf(" the value %1" i_ipc_ " of status should not occur\n", status);
            break;
        }
    }

    // Record solution information
    bqp_information( &data, &inform, &status );

    // Print solution details
    if(inform.status == 0){
        printf("%c:%6" i_ipc_ " iterations. Optimal objective value = %5.2f"
               " status = %1" i_ipc_ "\n",
               st, inform.iter, inform.obj, inform.status);
    }else{
        printf("%c: BQP_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
    bqp_terminate( &data, &control, &inform );
}

data	holds private internal data
status	is a scalar variable of type ipc_, that gives the entry and exit status from the package. On initial entry, status must be set to 1. Possible exit values are: 0 The run was successful. -1 An allocation error occurred. A message indicating the offending array is written on unit control.error, and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -2 A deallocation error occurred. A message indicating the offending array is written on unit control.error and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -3 The restriction n > 0 or requirement that a type contains its relevant string ‘dense’, ‘coordinate’, ‘sparse_by_rows’ or ‘diagonal’ has been violated. -4 The simple-bound constraints are inconsistent. -9 The analysis phase of the factorization failed; the return status from the factorization package is given in the component inform.factor_status -10 The factorization failed; the return status from the factorization package is given in the component inform.factor_status. -11 The solution of a set of linear equations using factors from the factorization package failed; the return status from the factorization package is given in the component inform.factor_status. -16 The problem is so ill-conditioned that further progress is impossible. -17 The step is too small to make further impact. -18 Too many iterations have been performed. This may happen if control.maxit is too small, but may also be symptomatic of a badly scaled problem. -19 The CPU time limit has been reached. This may happen if control.cpu_time_limit is too small, but may also be symptomatic of a badly scaled problem. -20 The Hessian matrix \(H\) appears to be indefinite. specified. -23 An entry from the strict upper triangle of \(H\) has been
n	is a scalar variable of type ipc_, that holds the number of variables
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.
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\).
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 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.