GALAHAD PRESOLVE package#

purpose#

The presolve package transforms linear and quadratic programming data so that the resulting problem is easier to solve. This reduced problem may then be passed to an appropriate solver. Once the reduced problem has been solved, it is then a restored to recover the solution for the original formulation.

The package applies presolving techniques to linear programs, whose aim is to minimize the linear objective function

\[\ell(x) = f + g^T x,\]

or quadratic programs, which target the quadratic objective function

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

subject to the general linear constraints

\[c_i^l \leq a_i^{T}x \leq c_i^u, \;\;\; i = 1, \ldots , m,\]

and simple bounds

\[x_j^l \leq x_j^{ } \leq x_j^u , \;\;\; j = 1, \ldots , n,\]

where the scalar $f$, the $n$-dimensional vectors $g$, $x^l$ and $x^u$, the $m$-dimensional vectors $c^l$ and $c^u$, the $n \times n$ symmetric matrix $H$ and the $m \times n$ matrix $A$ (whose rows are the vectors $a_i^T$) are given. Furthermore, bounds on the Lagrange multipliers $y$ associated with the general linear constraints and on the dual variables $z$ associated with the simple bound constraints

\[y_i^l \leq y_i \leq y_i^u, \;\;\; i = 1, \ldots , m,\]

and

\[z_i^l \leq z_i \leq z_i^u, \;\;\; i = 1, \ldots , n,\]

are also provided, where the $m$-dimensional vectors $y^l$ and $y^u$, as well as the $n$-dimensional vectors $x^l$ and $x^u$ are given. Any component of $c^l$, $c^u$, $x^l$, $x^u$, $y^l$, $y^u$, $z^l$ or $z^u$ may be infinite.

Currently only the options and inform dictionaries are exposed; these are provided and used by other GALAHAD packages with Python interfaces. Please contact us if you would like full functionality!

See Section 4 of $GALAHAD/doc/presolve.pdf for additional details.

method#

The required solution $x$ of the problem necessarily satisfies the primal optimality conditions

\[A x = c\]

and

\[c^l \leq c \leq c^u, \;\; x^l \leq x \leq x^u,\]

the dual optimality conditions

\[H x + g = A^{T} y + z, \;\; y = y^l + y^u \;\; \mbox{and} \;\; z = z^l + z^u,\]

and

\[y^l \geq 0 , \;\; y^u \leq 0 , \;\; z^l \geq 0 \;\; \mbox{and} \;\; z^u \leq 0,\]

and the complementary slackness conditions

\[ ( A x - c^l )^{T} y^l = 0, \;\; ( A x - c^u )^{T} y^u = 0, \;\; (x -x^l )^{T} z^l = 0 \;\;\mbox{and}\;\; (x -x^u )^{T} z^u = 0, \]

where the vectors $y$ and $z$ are known as the Lagrange multipliers for the general linear constraints, and the dual variables for the bounds, respectively, and where the vector inequalities hold componentwise.

The purpose of presolving is to exploit these equations in order to reduce the problem to the standard form defined as follows:

$*$ The variables are ordered so that their bounds appear in the order

variable bound order#
free			$x$
non-negativity	0	$\leq$	$x$
lower	$x^l$	$\leq$	$x$
range	$x^l$	$\leq$	$x$	$\leq$	$x^u$
upper			$x$	$\leq$	$x^u$
non-positivity			$x$	$\leq$	0

Fixed variables are removed. Within each category, the variables are further ordered so that those with non-zero diagonal Hessian entries occur before the remainder.

$*$ The constraints are ordered so that their bounds appear in the order

constraint bound order#
non-negativity	0	$\leq$	$A x$
equality	$c^l$	$=$	$A x$
lower	$c^l$	$\leq$	$A x$
range	$c^l$	$\leq$	$A x$	$\leq$	$c^u$
upper			$A x$	$\leq$	$c^u$
non-positivity			$A x$	$\leq$	0

Free constraints are removed.

$*$ In addition, constraints may be removed or bounds tightened, to reduce the size of the feasible region or simplify the problem if this is possible, and bounds may be tightened on the dual variables and the multipliers associated with the problem.

The presolving algorithm proceeds by applying a (potentially long) series of simple transformations to the problem, each transformation introducing a further simplification of the problem. These involve the removal of empty and singleton rows, the removal of redundant and forcing primal constraints, the tightening of primal and dual bounds, the exploitation of linear singleton, linear doubleton and linearly unconstrained columns, the merging dependent variables, row sparsification and split equalities. Transformations are applied in successive passes, each pass involving the following actions:

$*$ remove empty and singletons rows,
$*$ try to eliminate variables that are linearly unconstrained,
$*$ attempt to exploit the presence of linear singleton columns,
$*$ attempt to exploit the presence of linear doubleton columns,
$*$ complete the analysis of the dual constraints,
$*$ remove empty and singletons rows,
$*$ possibly remove dependent variables,
$*$ analyze the primal constraints,
$*$ try to make $A$ sparser by combining its rows,
$*$ check the current status of the variables, dual variables and multipliers.

All these transformations are applied to the structure of the original problem, which is only permuted to standard form after all transformations are completed. Note that the Hessian and Jacobian of the resulting reduced problem are always stored in sparse row-wise format. The reduced problem is then solved by a quadratic or linear programming solver, thus ensuring sufficiently small primal-dual feasibility and complementarity. Finally, the solution of the simplified problem is re-translated in the variables/constraints/format of the original problem formulation by a restoration phase.

If the number of problem transformations exceeds control.transf_buffer_size, the transformation buffer size, then they are saved in a “history” file, whose name may be chosen by specifying the control.transf_file_name parameter. When this is the case, this file is subsequently reread by presolve_restore. It must not be altered by the user.

At the overall level, the presolving process follows one of the two sequences:

$\boxed{\mbox{initialize}}$ $\rightarrow$ $\bigg[$ $\boxed{\mbox{apply transformations}}$ $\rightarrow$ (solve problem) $\rightarrow$ $\boxed{\mbox{restore}}$ $\bigg]$ $\rightarrow$ $\boxed{\mbox{terminate}}$

$\boxed{\mbox{initialize}}$ $\rightarrow$ $\bigg[$ $\boxed{\mbox{read specfile}}$ $\rightarrow$ $\boxed{\mbox{apply transformations}}$ $\rightarrow$ (solve problem) $\rightarrow$ $\boxed{\mbox{restore}}$ $\bigg]$ $\rightarrow$ $\boxed{\mbox{terminate}}$

where the procedure’s control parameter may be modified by reading the specfile, and where (solve problem) indicates that the reduced problem is solved. Each of the “boxed” steps in these sequences corresponds to calling a specific function in the package. In the above diagrams, brackated subsequence of steps means that they can be repeated with problem having the same structure. The value of the new_problem_structure argument must be true on entry to presolve_apply on the first time it is used in this repeated subsequence. Such a subsequence must be terminated by a call to presolve_terminate before presolving is applied to a problem with a different structure.

Note that the values of the multipliers and dual variables (and thus of their respective bounds) depend on the functional form assumed for the Lagrangian function associated with the problem. This form is given by

\[L( x, y, z ) = q( x ) - {\tt y_{sign}} * y^T ( A x - c ) - {\tt z_{sign}} * z,\]

(considering only active constraints $A x = c$), where the parameters ${\tt y_{sign}}$ and ${\tt z_{sign}}$ are +1 or -1 and can be chosen by the user using control.y_sign and control.z_sign. Thus, if ${\tt y_{sign} = +1}$, the multipliers associated to active constraints originally posed as inequalities are non-negative if the inequality is a lower bound and non-positive if it is an upper bound. Obvioulsy they are not constrained in sign for constraints originally posed as equalities. These sign conventions are reversed if ${\tt y_{sign} = -1}$. Similarly, if ${\tt z_{sign} = +1}$, the dual variables associated to active bounds are non-negative if the original bound is an lower bound, non-positive if it is an upper bound, or unconstrained in sign if the variables is fixed; and this convention is reversed in ${\tt z_{sign} = -1}$. The values of ${\tt z_{sign}}$ and ${\tt y_{sign}}$ may be chosen by setting the corresponding components of control to 1 or -1.

references#

The algorithm is described in more detail in

N. I. M. Gould and Ph. L. Toint, ``Presolving for quadratic programming’’. Mathematical Programming 100(1) (2004) 95–132.

matrix storage#

unsymmetric storage#

The unsymmetric $m$ by $n$ matrix $A$ may be presented and stored in a variety of convenient input formats.

Dense storage format: The matrix $A$ is stored as a compact dense matrix by rows, that is, the values of the entries of each row in turn are stored in order within an appropriate real one-dimensional array. In this case, component $n \ast i + j$ of the storage array A_val will hold the value $A_{ij}$ for $1 \leq i \leq m$, $1 \leq j \leq n$. The string A_type = ‘dense’ should be specified.

Dense by columns storage format: The matrix $A$ is stored as a compact dense matrix by columns, that is, the values of the entries of each column in turn are stored in order within an appropriate real one-dimensional array. In this case, component $m \ast j + i$ of the storage array A_val will hold the value $A_{ij}$ for $1 \leq i \leq m$, $1 \leq j \leq n$. The string A_type = ‘dense_by_columns’ should be specified.

Sparse co-ordinate storage format: Only the nonzero entries of the matrices are stored. For the $l$-th entry, $1 \leq l \leq ne$, of $A$, its row index i, column index j and value $A_{ij}$, $1 \leq i \leq m$, $1 \leq j \leq n$, are stored as the $l$-th components of the integer arrays A_row and A_col and real array A_val, respectively, while the number of nonzeros is recorded as A_ne = $ne$. The string A_type = ‘coordinate’should be specified.

Sparse row-wise storage format: Again only the nonzero entries are stored, but this time they are ordered so that those in row i appear directly before those in row i+1. For the i-th row of $A$ the i-th component of the integer array A_ptr holds the position of the first entry in this row, while A_ptr(m+1) holds the total number of entries plus one. The column indices j, $1 \leq j \leq n$, and values $A_{ij}$ of the nonzero entries in the i-th row are stored in components l = A_ptr(i), $\ldots$, A_ptr(i+1)-1, $1 \leq i \leq m$, of the integer array A_col, and real array A_val, respectively. For sparse matrices, this scheme almost always requires less storage than its predecessor. The string A_type = ‘sparse_by_rows’ should be specified.

Sparse column-wise storage format: Once again only the nonzero entries are stored, but this time they are ordered so that those in column j appear directly before those in column j+1. For the j-th column of $A$ the j-th component of the integer array A_ptr holds the position of the first entry in this column, while A_ptr(n+1) holds the total number of entries plus one. The row indices i, $1 \leq i \leq m$, and values $A_{ij}$ of the nonzero entries in the j-th columnsare stored in components l = A_ptr(j), $\ldots$, A_ptr(j+1)-1, $1 \leq j \leq n$, of the integer array A_row, and real array A_val, respectively. As before, for sparse matrices, this scheme almost always requires less storage than the co-ordinate format. The string A_type = ‘sparse_by_columns’ should be specified.

symmetric storage#

The symmetric $n$ by $n$ matrix $H$ may also be presented and stored in a variety of formats. But crucially symmetry is exploited by only storing values from the lower triangular part (i.e, those entries that lie on or below the leading diagonal).

Dense storage format: The matrix $H$ is stored as a compact dense matrix by rows, that is, the values of the entries of each row in turn are stored in order within an appropriate real one-dimensional array. Since $H$ is symmetric, only the lower triangular part (that is the part $H_{ij}$ for $1 \leq j \leq i \leq n$) need be held. In this case the lower triangle should be stored by rows, that is component $(i-1) * i / 2 + j$ of the storage array H_val will hold the value $H_{ij}$ (and, by symmetry, $H_{ji}$) for $1 \leq j \leq i \leq n$. The string H_type = ‘dense’ should be specified.

Sparse co-ordinate storage format: Only the nonzero entries of the matrices are stored. For the $l$-th entry, $1 \leq l \leq ne$, of $H$, its row index i, column index j and value $H_{ij}$, $1 \leq j \leq i \leq n$, are stored as the $l$-th components of the integer arrays H_row and H_col and real array H_val, respectively, while the number of nonzeros is recorded as H_ne = $ne$. Note that only the entries in the lower triangle should be stored. The string H_type = ‘coordinate’ should be specified.

Sparse row-wise storage format: Again only the nonzero entries are stored, but this time they are ordered so that those in row i appear directly before those in row i+1. For the i-th row of $H$ the i-th component of the integer array H_ptr holds the position of the first entry in this row, while H_ptr(n+1) holds the total number of entries plus one. The column indices j, $1 \leq j \leq i$, and values $H_{ij}$ of the entries in the i-th row are stored in components l = H_ptr(i), …, H_ptr(i+1)-1 of the integer array H_col, and real array H_val, respectively. Note that as before only the entries in the lower triangle should be stored. For sparse matrices, this scheme almost always requires less storage than its predecessor. The string H_type = ‘sparse_by_rows’ should be specified.

Diagonal storage format: If $H$ is diagonal (i.e., $H_{ij} = 0$ for all $1 \leq i \neq j \leq n$) only the diagonals entries $H_{ii}$, $1 \leq i \leq n$ need be stored, and the first n components of the array H_val may be used for the purpose. The string H_type = ‘diagonal’ should be specified.

Multiples of the identity storage format: If $H$ is a multiple of the identity matrix, (i.e., $H = \alpha I$ where $I$ is the n by n identity matrix and $\alpha$ is a scalar), it suffices to store $\alpha$ as the first component of H_val. The string H_type = ‘scaled_identity’ should be specified.

The identity matrix format: If $H$ is the identity matrix, no values need be stored. The string H_type = ‘identity’ should be specified.

The zero matrix format: The same is true if $H$ is the zero matrix, but now the string H_type = ‘zero’ or ‘none’ should be specified.

introduction to function calls#

To solve a given problem, functions from the presolve package must be called in the following order:

presolve_initialize - provide default control parameters and set up initial data structures
presolve_read_specfile (optional) - override control values by reading replacement values from a file
presolve_import_problem - import the problem data and report the dimensions of the transformed problem
presolve_transform_problem - apply the presolve algorithm to transform the data
presolve_restore_solution - restore the solution from that of the transformed problem
presolve_information (optional) - recover information about the solution and solution process
presolve_terminate - deallocate data structures

See the examples section for illustrations of use.

parametric real type T and integer type INT#

Below, the symbol T refers to a parametric real type that may be Float32 (single precision), Float64 (double precision) or, if supported, Float128 (quadruple precision). The symbol INT refers to a parametric integer type that may be Int32 (32-bit integer) or Int64 (64-bit integer).

callable functions#

    function presolve_initialize(T, INT, data, control, status)

Set default control values and initialize private data

Parameters:

data

holds private internal data

control

is a structure containing control information (see presolve_control_type)

status

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

0
The initialization was successful.

    function presolve_read_specfile(T, INT, control, specfile)

Read the content of a specification file, and assign values associated with given keywords to the corresponding control parameters. An in-depth discussion of specification files is available, and a detailed list of keywords with associated default values is provided in $GALAHAD/src/presolve/PRESOLVE.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/presolve.pdf for a list of how these keywords relate to the components of the control structure.

Parameters:

control	is a structure containing control information (see presolve_control_type)
specfile	is a one-dimensional array of type Vararg{Cchar} that must give the name of the specification file

    function presolve_import_problem(T, INT, control, data, status, n, m,
                                     H_type, H_ne, H_row, H_col, H_ptr,
                                     H_val,  g, f, A_type, A_ne, A_row,
                                     A_col, A_ptr, A_val, c_l, c_u,
                                     x_l, x_u, n_out, m_out, H_ne_out,
                                     A_ne_out)

Import the initial data, and apply the presolve algorithm to report crucial characteristics of the transformed variant

Parameters:

control	is a structure whose members provide control parameters for the remaining procedures (see presolve_control_type)
data	holds private internal data
status	is a scalar variable of type INT that gives the exit status from the package. Possible values are: 0 The import 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 restrictions n > 0 or m > 0 or requirement that a type contains its relevant string ‘dense’, ‘coordinate’, ‘sparse_by_rows’ or ‘diagonal’ has been violated. -23 An entry from the strict upper triangle of $H$ has been specified.
n	is a scalar variable of type INT that holds the number of variables.
m	is a scalar variable of type INT that holds the number of general linear constraints.
H_type	is a one-dimensional array of type Vararg{Cchar} that specifies the symmetric storage scheme used for the Hessian, $H$. It should be one of ‘coordinate’, ‘sparse_by_rows’, ‘dense’, ‘diagonal’, ‘scaled_identity’, ‘identity’, ‘zero’ or ‘none’, the latter pair if $H=0$; lower or upper case variants are allowed.
H_ne	is a scalar variable of type INT 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 INT that holds the row indices of the lower triangular part of $H$ in the sparse co-ordinate storage scheme. It need not be set for any of the other three schemes, and in this case can be C_NULL.
H_col	is a one-dimensional array of size H_ne and type INT that holds the column indices of the lower triangular part of $H$ in either the sparse co-ordinate, or the sparse row-wise storage scheme. It need not be set when the dense, diagonal or (scaled) identity storage schemes are used, and in this case can be C_NULL.
H_ptr	is a one-dimensional array of size n+1 and type INT that holds the starting position of each row of the lower triangular part of $H$, as well as the total number of entries, in the sparse row-wise storage scheme. It need not be set when the other schemes are used, and in this case can be C_NULL.
H_val	is a one-dimensional array of size h_ne and type T that holds the values of the entries of the lower triangular part of the Hessian matrix $H$ in any of the available storage schemes.
g	is a one-dimensional array of size n and type T that holds the linear term $g$ of the objective function. The j-th component of `g`, j = 1, … , n, contains $g_j$.
f	is a scalar of type T that holds the constant term $f$ of the objective function.
A_type	is a one-dimensional array of type Vararg{Cchar} that specifies the unsymmetric storage scheme used for the constraint Jacobian, $A$. It should be one of ‘coordinate’, ‘sparse_by_rows’ or ‘dense; lower or upper case variants are allowed.
A_ne	is a scalar variable of type INT 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 INT that holds the row indices of $A$ in the sparse co-ordinate storage scheme. It need not be set for any of the other schemes, and in this case can be C_NULL.
A_col	is a one-dimensional array of size A_ne and type INT that holds the column indices of $A$ in either the sparse co-ordinate, or the sparse row-wise storage scheme. It need not be set when the dense or diagonal storage schemes are used, and in this case can be C_NULL.
A_ptr	is a one-dimensional array of size n+1 and type INT that holds the starting position of each row of $A$, as well as the total number of entries, in the sparse row-wise storage scheme. It need not be set when the other schemes are used, and in this case can be C_NULL.
A_val	is a one-dimensional array of size a_ne and type T that holds the values of the entries of the constraint Jacobian matrix $A$ in any of the available storage schemes.
c_l	is a one-dimensional array of size m and type T that holds the lower bounds $c^l$ on the constraints $A x$. The i-th component of `c_l`, i = 1, … , m, contains $c^l_i$.
c_u	is a one-dimensional array of size m and type T that holds the upper bounds $c^l$ on the constraints $A x$. The i-th component of `c_u`, i = 1, … , m, contains $c^u_i$.
x_l	is a one-dimensional array of size n and type T that holds the lower bounds $x^l$ on the variables $x$. The j-th component of `x_l`, j = 1, … , n, contains $x^l_j$.
x_u	is a one-dimensional array of size n and type T that holds the upper bounds $x^l$ on the variables $x$. The j-th component of `x_u`, j = 1, … , n, contains $x^l_j$.
n_out	is a scalar variable of type INT that holds the number of variables in the transformed problem.
m_out	is a scalar variable of type INT that holds the number of general linear constraints in the transformed problem.
H_ne_out	is a scalar variable of type INT that holds the number of entries in the lower triangular part of $H$ in the transformed problem.
A_ne_out	is a scalar variable of type INT that holds the number of entries in $A$ in the transformed problem.

    function presolve_transform_problem(T, INT, data, status, n, m, H_ne, H_col,
                                        H_ptr, H_val, g, f, A_ne, A_col,
                                        A_ptr, A_val, c_l, c_u, x_l, x_u,
                                        y_l, y_u, z_l, z_u)

Apply the presolve algorithm to simplify the input problem, and output the transformed variant

Parameters:

data	holds private internal data
status	is a scalar variable of type INT that gives the exit status from the package. Possible values are: 0 The import 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 input values n, m, A_ne or H_ne do not agree with those output as necessary from presolve_import_problem.
n	is a scalar variable of type INT that holds the number of variables in the transformed problem. This must match the value n_out from the last call to presolve_import_problem.
m	is a scalar variable of type INT that holds the number of general linear constraints. This must match the value m_out from the last call to presolve_import_problem.
H_ne	is a scalar variable of type INT that holds the number of entries in the lower triangular part of the transformed $H$. This must match the value H_ne_out from the last call to presolve_import_problem.
H_col	is a one-dimensional array of size H_ne and type INT that holds the column indices of the lower triangular part of the transformed $H$ in the sparse row-wise storage scheme.
H_ptr	is a one-dimensional array of size n+1 and type INT that holds the starting position of each row of the lower triangular part of the transformed $H$ in the sparse row-wise storage scheme.
H_val	is a one-dimensional array of size h_ne and type T that holds the values of the entries of the lower triangular part of the the transformed Hessian matrix $H$ in the sparse row-wise storage scheme.
g	is a one-dimensional array of size n and type T that holds the the transformed linear term $g$ of the objective function. The j-th component of `g`, j = 1, … , n, contains $g_j$.
f	is a scalar of type T that holds the transformed constant term $f$ of the objective function.
A_ne	is a scalar variable of type INT that holds the number of entries in the transformed $A$. This must match the value A_ne_out from the last call to presolve_import_problem.
A_col	is a one-dimensional array of size A_ne and type INT that holds the column indices of the transformed $A$ in the sparse row-wise storage scheme.
A_ptr	is a one-dimensional array of size n+1 and type INT that holds the starting position of each row of the transformed $A$, as well as the total number of entries, in the sparse row-wise storage scheme.
A_val	is a one-dimensional array of size a_ne and type T that holds the values of the entries of the transformed constraint Jacobian matrix $A$ in the sparse row-wise storage scheme.
c_l	is a one-dimensional array of size m and type T that holds the transformed lower bounds $c^l$ on the constraints $A x$. The i-th component of `c_l`, i = 1, … , m, contains $c^l_i$.
c_u	is a one-dimensional array of size m and type T that holds the transformed upper bounds $c^l$ on the constraints $A x$. The i-th component of `c_u`, i = 1, … , m, contains $c^u_i$.
x_l	is a one-dimensional array of size n and type T that holds the transformed lower bounds $x^l$ on the variables $x$. The j-th component of `x_l`, j = 1, … , n, contains $x^l_j$.
x_u	is a one-dimensional array of size n and type T that holds the transformed upper bounds $x^l$ on the variables $x$. The j-th component of `x_u`, j = 1, … , n, contains $x^l_j$.
y_l	is a one-dimensional array of size m and type T that holds the implied lower bounds $y^l$ on the transformed Lagrange multipliers $y$. The i-th component of `y_l`, i = 1, … , m, contains $y^l_i$.
y_u	is a one-dimensional array of size m and type T that holds the implied upper bounds $y^u$ on the transformed Lagrange multipliers $y$. The i-th component of `y_u`, i = 1, … , m, contains $y^u_i$.
z_l	is a one-dimensional array of size m and type T that holds the implied lower bounds $y^l$ on the transformed dual variables $z$. The j-th component of `z_l`, j = 1, … , n, contains $z^l_i$.
z_u	is a one-dimensional array of size m and type T that holds the implied upper bounds $y^u$ on the transformed dual variables $z$. The j-th component of `z_u`, j = 1, … , n, contains $z^u_i$.

    function presolve_restore_solution(T, INT, data, status, n_in, m_in, x_in,
                                        c_in, y_in, z_in, n, m, x, c, y, z)

Given the solution (x_in,c_in,y_in,z_in) to the transformed problem, restore to recover the solution (x,c,y,z) to the original

Parameters:

data	holds private internal data
status	is a scalar variable of type INT that gives the exit status from the package. Possible values are: 0 The import 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 input values n, m, n_in and m_in do not agree with those input to and output as necessary from presolve_import_problem.
n_in	is a scalar variable of type INT that holds the number of variables in the transformed problem. This must match the value n_out from the last call to presolve_import_problem.
m_in	is a scalar variable of type INT that holds the number of general linear constraints. This must match the value m_out from the last call to presolve_import_problem.
x_in	is a one-dimensional array of size n_in and type T that holds the transformed values $x$ of the optimization variables. The j-th component of `x`, j = 1, … , n, contains $x_j$.
c_in	is a one-dimensional array of size m and type T that holds the transformed residual $c(x)$. The i-th component of `c`, j = 1, … , m, contains $c_j(x)$.
y_in	is a one-dimensional array of size n_in and type T that holds the values $y$ of the transformed Lagrange multipliers for the general linear constraints. The j-th component of `y`, j = 1, … , m, contains $y_j$.
z_in	is a one-dimensional array of size n_in and type T that holds the values $z$ of the transformed dual variables. The j-th component of `z`, j = 1, … , n, contains $z_j$.
n	is a scalar variable of type INT that holds the number of variables in the transformed problem. This must match the value n as input to presolve_import_problem.
m	is a scalar variable of type INT that holds the number of general linear constraints. This must match the value m as input to presolve_import_problem.
x	is a one-dimensional array of size n and type T that holds the transformed values $x$ of the optimization variables. The j-th component of `x`, j = 1, … , n, contains $x_j$.
c	is a one-dimensional array of size m and type T that holds the transformed residual $c(x)$. The i-th component of `c`, j = 1, … , m, contains $c_j(x)$.
y	is a one-dimensional array of size n and type T that holds the values $y$ of the transformed Lagrange multipliers for the general linear constraints. The j-th component of `y`, j = 1, … , m, contains $y_j$.
z	is a one-dimensional array of size n and type T that holds the values $z$ of the transformed dual variables. The j-th component of `z`, j = 1, … , n, contains $z_j$.

    function presolve_information(T, INT, data, inform, status)

Provides output information

Parameters:

data

holds private internal data

inform

is a structure containing output information (see presolve_inform_type)

status

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

0
The values were recorded successfully

    function presolve_terminate(T, INT, data, control, inform)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a structure containing control information (see presolve_control_type)
inform	is a structure containing output information (see presolve_inform_type)

available structures#

presolve_control_type structure#

control derived type as a Julia structureure More…

    struct presolve_control_type{T,INT}
      f_indexing::Bool
      termination::INT
      max_nbr_transforms::INT
      max_nbr_passes::INT
      c_accuracy::T
      z_accuracy::T
      infinity::T
      out::INT
      errout::INT
      print_level::INT
      dual_transformations::Bool
      redundant_xc::Bool
      primal_constraints_freq::INT
      dual_constraints_freq::INT
      singleton_columns_freq::INT
      doubleton_columns_freq::INT
      unc_variables_freq::INT
      dependent_variables_freq::INT
      sparsify_rows_freq::INT
      max_fill::INT
      transf_file_nbr::INT
      transf_buffer_size::INT
      transf_file_status::INT
      transf_file_name::NTuple{31,Cchar}
      y_sign::INT
      inactive_y::INT
      z_sign::INT
      inactive_z::INT
      final_x_bounds::INT
      final_z_bounds::INT
      final_c_bounds::INT
      final_y_bounds::INT
      check_primal_feasibility::INT
      check_dual_feasibility::INT
      pivot_tol::T
      min_rel_improve::T
      max_growth_factor::T

detailed documentation#

control derived type as a Julia structure

components#

Bool f_indexing

use C or Fortran sparse matrix indexing

INT termination

Determines the strategy for terminating the presolve analysis. Possible values are:

1 presolving is continued as long as one of the sizes of the problem (n, m, a_ne, or h_ne) is being reduced;
2 presolving is continued as long as problem transformations remain possible. NOTE: the maximum number of analysis passes (control.max_nbr_passes) and the maximum number of problem transformations (control.max_nbr_transforms) set an upper limit on the presolving effort irrespective of the choice of control.termination. The only effect of this latter parameter is to allow for early termination.

INT max_nbr_transforms

The maximum number of problem transformations, cumulated over all calls to presolve.

INT max_nbr_passes

The maximum number of analysis passes for problem analysis during a single call of presolve_transform_problem.

T c_accuracy

The relative accuracy at which the general linear constraints are satisfied at the exit of the solver. Note that this value is not used before the restoration of the problem.

T z_accuracy

The relative accuracy at which the dual feasibility constraints are satisfied at the exit of the solver. Note that this value is not used before the restoration of the problem.

T infinity

The value beyond which a number is deemed equal to plus infinity (minus infinity being defined as its opposite)

INT out

The unit number associated with the device used for printout.

INT errout

The unit number associated with the device used for error ouput.

INT print_level

The level of printout requested by the user. Can take the values:

0 no printout is produced
1 only reports the major steps in the analysis
2 reports the identity of each problem transformation
3 reports more details
4 reports lots of information.
5 reports a completely silly amount of information

Bool dual_transformations

true if dual transformations of the problem are allowed. Note that this implies that the reduced problem is solved accurately (for the dual feasibility condition to hold) as to be able to restore the problem to the original constraints and variables. false prevents dual transformations to be applied, thus allowing for inexact solution of the reduced problem. The setting of this control parameter overides that of get_z, get_z_bounds, get_y, get_y_bounds, dual_constraints_freq, singleton_columns_freq, doubleton_columns_freq, z_accuracy, check_dual_feasibility.

Bool redundant_xc

true if the redundant variables and constraints (i.e. variables that do not appear in the objective function and appear with a consistent sign in the constraints) are to be removed with their associated constraints before other transformations are attempted.

INT primal_constraints_freq

The frequency of primal constraints analysis in terms of presolving passes. A value of j = 2 indicates that primal constraints are analyzed every 2 presolving passes. A zero value indicates that they are never analyzed.

INT dual_constraints_freq

The frequency of dual constraints analysis in terms of presolving passes. A value of j = 2 indicates that dual constraints are analyzed every 2 presolving passes. A zero value indicates that they are never analyzed.

INT singleton_columns_freq

The frequency of singleton column analysis in terms of presolving passes. A value of j = 2 indicates that singleton columns are analyzed every 2 presolving passes. A zero value indicates that they are never analyzed.

INT doubleton_columns_freq

The frequency of doubleton column analysis in terms of presolving passes. A value of j indicates that doubleton columns are analyzed every 2 presolving passes. A zero value indicates that they are never analyzed.

INT unc_variables_freq

The frequency of the attempts to fix linearly unconstrained variables, expressed in terms of presolving passes. A value of j = 2 indicates that attempts are made every 2 presolving passes. A zero value indicates that no attempt is ever made.

INT dependent_variables_freq

The frequency of search for dependent variables in terms of presolving passes. A value of j = 2 indicates that dependent variables are searched for every 2 presolving passes. A zero value indicates that they are never searched for.

INT sparsify_rows_freq

The frequency of the attempts to make A sparser in terms of presolving passes. A value of j = 2 indicates that attempts are made every 2 presolving passes. A zero value indicates that no attempt is ever made.

INT max_fill

The maximum percentage of fill in each row of A. Note that this is a row-wise measure: globally fill never exceeds the storage initially used for A, no matter how large control.max_fill is chosen. If max_fill is negative, no limit is put on row fill.

INT transf_file_nbr

The unit number to be associated with the file(s) used for saving problem transformations on a disk file.

INT transf_buffer_size

The number of transformations that can be kept in memory at once (that is without being saved on a disk file).

INT transf_file_status

The exit status of the file where problem transformations are saved:

0 the file is not deleted after program termination
1 the file is not deleted after program termination

char transf_file_name[31]

The name of the file (to be) used for storing problem transformation on disk. NOTE: this parameter must be identical for all calls to presolve following presolve_read_specfile. It can then only be changed after calling presolve_terminate.

INT y_sign

Determines the convention of sign used for the multipliers associated with the general linear constraints.

1 All multipliers corresponding to active inequality constraints are non-negative for lower bound constraints and non-positive for upper bounds constraints.
-1 All multipliers corresponding to active inequality constraints are non-positive for lower bound constraints and non-negative for upper bounds constraints.

INT inactive_y

Determines whether or not the multipliers corresponding to constraints that are inactive at the unreduced point corresponding to the reduced point on input to presolve_restore_solution must be set to zero. Possible values are: associated with the general linear constraints.

0 All multipliers corresponding to inactive inequality constraints are forced to zero, possibly at the expense of deteriorating the dual feasibility condition.
1 Multipliers corresponding to inactive inequality constraints are left unaltered.

INT z_sign

Determines the convention of sign used for the dual variables associated with the bound constraints.

1 All dual variables corresponding to active lower bounds are non-negative, and non-positive for active upper bounds.
-1 All dual variables corresponding to active lower bounds are non-positive, and non-negative for active upper bounds.

INT inactive_z

Determines whether or not the dual variables corresponding to bounds that are inactive at the unreduced point corresponding to the reduced point on input to presolve_restore_solution must be set to zero. Possible values are: associated with the general linear constraints.

0: All dual variables corresponding to inactive bounds are forced to zero, possibly at the expense of deteriorating the dual feasibility condition.
1 Dual variables corresponding to inactive bounds are left unaltered.

INT final_x_bounds

The type of final bounds on the variables returned by the package. This parameter can take the values:

0 the final bounds are the tightest bounds known on the variables (at the risk of being redundant with other constraints, which may cause degeneracy);
1 the best known bounds that are known to be non-degenerate. This option implies that an additional real workspace of size 2 * n must be allocated.
2 the loosest bounds that are known to keep the problem equivalent to the original problem. This option also implies that an additional real workspace of size 2 * n must be allocated.

NOTE: this parameter must be identical for all calls to presolve (except presolve_initialize).

INT final_z_bounds

The type of final bounds on the dual variables returned by the package. This parameter can take the values:

0 the final bounds are the tightest bounds known on the dual variables (at the risk of being redundant with other constraints, which may cause degeneracy);
1 the best known bounds that are known to be non-degenerate. This option implies that an additional real workspace of size 2 * n must be allocated.
2 the loosest bounds that are known to keep the problem equivalent to the original problem. This option also implies that an additional real workspace of size 2 * n must be allocated.

NOTE: this parameter must be identical for all calls to presolve (except presolve_initialize).

INT final_c_bounds

The type of final bounds on the constraints returned by the package. This parameter can take the values:

0 the final bounds are the tightest bounds known on the constraints (at the risk of being redundant with other constraints, which may cause degeneracy);
1 the best known bounds that are known to be non-degenerate. This option implies that an additional real workspace of size 2 * m must be allocated.
2 the loosest bounds that are known to keep the problem equivalent to the original problem. This option also implies that an additional real workspace of size 2 * n must be allocated.

NOTES: 1) This parameter must be identical for all calls to presolve (except presolve_initialize). 2) If different from 0, its value must be identical to that of control.final_x_bounds.

INT final_y_bounds

The type of final bounds on the multipliers returned by the package. This parameter can take the values:

0 the final bounds are the tightest bounds known on the multipliers (at the risk of being redundant with other constraints, which may cause degeneracy);
1 the best known bounds that are known to be non-degenerate. This option implies that an additional real workspace of size 2 * m must be allocated.
2 the loosest bounds that are known to keep the problem equivalent to the original problem. This option also implies that an additional real workspace of size 2 * n must be allocated.

NOTE: this parameter must be identical for all calls to presolve (except presolve_initialize).

INT check_primal_feasibility

The level of feasibility check (on the values of x) at the start of the restoration phase. This parameter can take the values:

0 no check at all;
1 the primal constraints are recomputed at x and a message issued if the computed value does not match the input value, or if it is out of bounds (if control.print_level >= 2);
2 the same as for 1, but presolve is terminated if an incompatibilty is detected.

INT check_dual_feasibility

The level of dual feasibility check (on the values of x, y and z) at the start of the restoration phase. This parameter can take the values:

0 no check at all;
1 the dual feasibility condition is recomputed at ( x, y, z ) and a message issued if the computed value does not match the input value (if control.print_level >= 2);
2 the same as for 1, but presolve is terminated if an incompatibilty is detected. The last two values imply the allocation of an additional real workspace vector of size equal to the number of variables in the reduced problem.

T pivot_tol

The relative pivot tolerance above which pivoting is considered as numerically stable in transforming the coefficient matrix A. A zero value corresponds to a totally unsafeguarded pivoting strategy (potentially unstable).

T min_rel_improve

The minimum relative improvement in the bounds on x, y and z for a tighter bound on these quantities to be accepted in the course of the analysis. More formally, if lower is the current value of the lower bound on one of the x, y or z, and if new_lower is a tentative tighter lower bound on the same quantity, it is only accepted if.

new_lower >= lower + tol * MAX( 1, ABS( lower ) ),

where

tol = control.min_rel_improve.

Similarly, a tentative tighter upper bound new_upper only replaces the current upper bound upper if

new_upper <= upper - tol * MAX( 1, ABS( upper ) ).

Note that this parameter must exceed the machine precision significantly.

T max_growth_factor

The maximum growth factor (in absolute value) that is accepted between the maximum data item in the original problem and any data item in the reduced problem. If a transformation results in this bound being exceeded, the transformation is skipped.

presolve_inform_type structure#

    struct presolve_inform_type{INT}
      status::INT
      status_continue::INT
      status_continued::INT
      nbr_transforms::INT
      message::NTuple{3,NTuple{81,Cchar}}

detailed documentation#

inform derived type as a Julia structure

components#

INT status

The presolve exit condition. It can take the following values (symbol in parentheses is the related Fortran code):

0

(OK) successful exit;
1

(MAX_NBR_TRANSF) the maximum number of problem transformation has been reached NOTE: this exit is not really an error, since the problem can nevertheless be permuted and solved. It merely signals that further problem reduction could possibly be obtained with a larger value of the parameter control.max_nbr_transforms
-1

(MEMORY_FULL) memory allocation failed
-2 (FILE_NOT_OPENED) a file intended for saving problem
transformations could not be opened;
-3 (COULD_NOT_WRITE) an IO error occurred while saving
transformations on the relevant disk file;
-4

(TOO_FEW_BITS_PER_BYTE) an integer contains less than NBRH + 1 bits.
-21

(PRIMAL_INFEASIBLE) the problem is primal infeasible;
-22

(DUAL_INFEASIBLE) the problem is dual infeasible;
-23

(WRONG_G_DIMENSION) the dimension of the gradient is incompatible with the problem dimension;
-24

(WRONG_HVAL_DIMENSION) the dimension of the vector containing the entries of the Hessian is erroneously specified;
-25

(WRONG_HPTR_DIMENSION) the dimension of the vector containing the addresses of the first entry of each Hessian row is erroneously specified;
-26

(WRONG_HCOL_DIMENSION) the dimension of the vector containing the column indices of the nonzero Hessian entries is erroneously specified;
-27

(WRONG_HROW_DIMENSION) the dimension of the vector containing the row indices of the nonzero Hessian entries is erroneously specified;
-28

(WRONG_AVAL_DIMENSION) the dimension of the vector containing the entries of the Jacobian is erroneously specified;
-29

(WRONG_APTR_DIMENSION) the dimension of the vector containing the addresses of the first entry of each Jacobian row is erroneously specified;
-30 (WRONG_ACOL_DIMENSION) the dimension of the vector containing
the column indices of the nonzero Jacobian entries is erroneously specified;
-31

(WRONG_AROW_DIMENSION) the dimension of the vector containing the row indices of the nonzero Jacobian entries is erroneously specified;
-32

(WRONG_X_DIMENSION) the dimension of the vector of variables is incompatible with the problem dimension;
-33

(WRONG_XL_DIMENSION) the dimension of the vector of lower bounds on the variables is incompatible with the problem dimension;
-34

(WRONG_XU_DIMENSION) the dimension of the vector of upper bounds on the variables is incompatible with the problem dimension;
-35

(WRONG_Z_DIMENSION) the dimension of the vector of dual variables is incompatible with the problem dimension;
-36

(WRONG_ZL_DIMENSION) the dimension of the vector of lower bounds on the dual variables is incompatible with the problem dimension;
-37

(WRONG_ZU_DIMENSION) the dimension of the vector of upper bounds on the dual variables is incompatible with the problem dimension;
-38

(WRONG_C_DIMENSION) the dimension of the vector of constraints values is incompatible with the problem dimension;
-39

(WRONG_CL_DIMENSION) the dimension of the vector of lower bounds on the constraints is incompatible with the problem dimension;
-40 (WRONG_CU_DIMENSION) the dimension of the vector of upper
bounds on the constraints is incompatible with the problem dimension;
-41

(WRONG_Y_DIMENSION) the dimension of the vector of multipliers values is incompatible with the problem dimension;
-42

(WRONG_YL_DIMENSION) the dimension of the vector of lower bounds on the multipliers is incompatible with the problem dimension;
-43

(WRONG_YU_DIMENSION) the dimension of the vector of upper bounds on the multipliers is incompatible with the problem dimension;
-44

(STRUCTURE_NOT_SET) the problem structure has not been set or has been cleaned up before an attempt to analyze;
-45

(PROBLEM_NOT_ANALYZED) the problem has not been analyzed before an attempt to permute it;
-46

(PROBLEM_NOT_PERMUTED) the problem has not been permuted or fully reduced before an attempt to restore it
-47

(H_MISSPECIFIED) the column indices of a row of the sparse Hessian are not in increasing order, in that they specify an entry above the diagonal;
-48

(CORRUPTED_SAVE_FILE) one of the files containing saved problem transformations has been corrupted between writing and reading;
-49

(WRONG_XS_DIMENSION) the dimension of the vector of variables’ status is incompatible with the problem dimension;
-50 (WRONG_CS_DIMENSION) the dimension of the vector of
constraints’ status is incompatible with the problem dimension;
-52

(WRONG_N) the problem does not contain any (active) variable;
-53

(WRONG_M) the problem contains a negative number of constraints;
-54

(SORT_TOO_LONG) the vectors are too long for the sorting routine;
-55

(X_OUT_OF_BOUNDS) the value of a variable that is obtained by substitution from a constraint is incoherent with the variable’s bounds. This may be due to a relatively loose accuracy on the linear constraints. Try to increase control.c_accuracy.
-56

(X_NOT_FEASIBLE) the value of a constraint that is obtained by recomputing its value on input of presolve_restore_solution from the current x is incompatible with its declared value or its bounds. This may caused the restored problem to be infeasible.
-57

(Z_NOT_FEASIBLE) the value of a dual variable that is obtained by recomputing its value on input to presolve_restore_solution (assuming dual feasibility) from the current values of $(x, y, z)$ is incompatible with its declared value. This may caused the restored problem to be infeasible or suboptimal.
-58

(Z_CANNOT_BE_ZEROED) a dual variable whose value is nonzero because the corresponding primal is at an artificial bound cannot be zeroed while maintaining dual feasibility (on restoration). This can happen when $( x, y, z)$ on input of RESTORE are not (sufficiently) optimal.
-60

(UNRECOGNIZED_KEYWORD) a keyword was not recognized in the analysis of the specification file
-61

(UNRECOGNIZED_VALUE) a value was not recognized in the analysis of the specification file
-63

(G_NOT_ALLOCATED) the vector G has not been allocated although it has general values
-64

(C_NOT_ALLOCATED) the vector C has not been allocated although m > 0
-65

(AVAL_NOT_ALLOCATED) the vector A.val has not been allocated although m > 0
-66

(APTR_NOT_ALLOCATED) the vector A.ptr has not been allocated although m > 0 and A is stored in row-wise sparse format
-67

(ACOL_NOT_ALLOCATED) the vector A.col has not been allocated although m > 0 and A is stored in row-wise sparse format or sparse coordinate format
-68

(AROW_NOT_ALLOCATED) the vector A.row has not been allocated although m > 0 and A is stored in sparse coordinate format
-69 (HVAL_NOT_ALLOCATED) the vector H.val has not been allocated
although H.ne > 0
-70

(HPTR_NOT_ALLOCATED) the vector H.ptr has not been allocated although H.ne > 0 and H is stored in row-wise sparse format
-71

(HCOL_NOT_ALLOCATED) the vector H.col has not been allocated although H.ne > 0 and H is stored in row-wise sparse format or sparse coordinate format
-72

(HROW_NOT_ALLOCATED) the vector H.row has not been allocated although H.ne > 0 and A is stored in sparse coordinate format
-73

(WRONG_ANE) incompatible value of A_ne
-74

(WRONG_HNE) incompatible value of H_ne

INT nbr_transforms

The final number of problem transformations, as reported to the user at exit.

char message[3][81]

A few lines containing a description of the exit condition on exit of PRESOLVE, typically including more information than indicated in the description of control.status above. It is printed out on device errout at the end of execution if control.print_level >= 1.

Index

example calls#

This is an example of how to use the package to presolve a quadratic program; the code is available in $GALAHAD/src/presolve/Julia/test_presolve.jl . A variety of supported Hessian and constraint matrix storage formats are shown.

# test_presolve.jl
# Simple code to test the Julia interface to PRESOLVE

using GALAHAD
using Test
using Printf
using Accessors
using Quadmath

function test_presolve(::Type{T}, ::Type{INT}) where {T,INT}
  # Derived types
  data = Ref{Ptr{Cvoid}}()
  control = Ref{presolve_control_type{T,INT}}()
  inform = Ref{presolve_inform_type{INT}}()

  # Set problem data
  n = INT(6)  # dimension
  m = INT(5)  # number of general constraints
  H_ne = INT(1)  # Hesssian elements
  H_row = INT[1]  # row indices, NB lower triangle
  H_col = INT[1]  # column indices, NB lower triangle
  H_ptr = INT[1, 2, 2, 2, 2, 2, 2]  # row pointers
  H_val = T[1.0]  # values
  g = T[1.0, 1.0, 1.0, 1.0, 1.0, 1.0]  # linear term in the objective
  f = one(T)  # constant term in the objective
  A_ne = INT(8)  # Jacobian elements
  A_row = INT[3, 3, 3, 4, 4, 5, 5, 5]  # row indices
  A_col = INT[3, 4, 5, 3, 6, 4, 5, 6]  # column indices
  A_ptr = INT[1, 1, 1, 4, 6, 9]  # row pointers
  A_val = T[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]  # values
  c_l = T[0.0, 0.0, 2.0, 1.0, 3.0]  # constraint lower bound
  c_u = T[1.0, 1.0, 3.0, 3.0, 3.0]  # constraint upper bound
  x_l = T[-3.0, 0.0, 0.0, 0.0, 0.0, 0.0]  # variable lower bound
  x_u = T[3.0, 1.0, 1.0, 1.0, 1.0, 1.0]  # variable upper bound

  # Set output storage
  st = ' '
  status = Ref{INT}()

  @printf(" Fortran sparse matrix indexing\n\n")
  @printf(" basic tests of qp storage formats\n\n")
  for d in 1:7
    n_trans = Ref{INT}()
    m_trans = Ref{INT}()
    H_ne_trans = Ref{INT}()
    A_ne_trans = Ref{INT}()

    # Initialize PRESOLVE
    presolve_initialize(T, INT, data, control, status)

    # Set user-defined control options
    @reset control[].f_indexing = true # Fortran sparse matrix indexing

    # sparse co-ordinate storage
    if d == 1
      st = 'C'
      presolve_import_problem(T, INT, control, data, status, n, m,
                              "coordinate", H_ne, H_row, H_col, C_NULL, H_val, g, f,
                              "coordinate", A_ne, A_row, A_col, C_NULL, A_val,
                              c_l, c_u, x_l, x_u,
                              n_trans, m_trans, H_ne_trans, A_ne_trans)
    end

    # sparse by rows
    if d == 2
      st = 'R'
      presolve_import_problem(T, INT, control, data, status, n, m,
                              "sparse_by_rows", H_ne, C_NULL, H_col, H_ptr, H_val, g, f,
                              "sparse_by_rows", A_ne, C_NULL, A_col, A_ptr, A_val,
                              c_l, c_u, x_l, x_u,
                              n_trans, m_trans, H_ne_trans, A_ne_trans)
    end

    # dense
    if d == 3
      st = 'D'
      H_dense_ne = div(n * (n + 1), 2) # number of elements of H
      A_dense_ne = m * n # number of elements of A
      H_dense = T[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
      A_dense = T[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                        0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0,
                        0.0, 1.0, 1.0, 1.0]
      presolve_import_problem(T, INT, control, data, status, n, m,
                              "dense", H_dense_ne, C_NULL, C_NULL, C_NULL, H_dense, g,
                              f, "dense", A_dense_ne, C_NULL, C_NULL, C_NULL, A_dense,
                              c_l, c_u, x_l, x_u,
                              n_trans, m_trans, H_ne_trans, A_ne_trans)
    end

    # diagonal
    if d == 4
      st = 'L'
      presolve_import_problem(T, INT, control, data, status, n, m,
                              "diagonal", n, C_NULL, C_NULL, C_NULL, H_val, g, f,
                              "sparse_by_rows", A_ne, C_NULL, A_col, A_ptr, A_val,
                              c_l, c_u, x_l, x_u,
                              n_trans, m_trans, H_ne_trans, A_ne_trans)
    end

    # scaled identity
    if d == 5
      st = 'S'
      presolve_import_problem(T, INT, control, data, status, n, m,
                              "scaled_identity", 1, C_NULL, C_NULL, C_NULL, H_val, g, f,
                              "sparse_by_rows", A_ne, C_NULL, A_col, A_ptr, A_val,
                              c_l, c_u, x_l, x_u,
                              n_trans, m_trans, H_ne_trans, A_ne_trans)
    end

    # identity
    if d == 6
      st = 'I'
      presolve_import_problem(T, INT, control, data, status, n, m,
                              "identity", 0, C_NULL, C_NULL, C_NULL, C_NULL, g, f,
                              "sparse_by_rows", A_ne, C_NULL, A_col, A_ptr, A_val,
                              c_l, c_u, x_l, x_u,
                              n_trans, m_trans, H_ne_trans, A_ne_trans)
    end

    # zero
    if d == 7
      st = 'Z'
      presolve_import_problem(T, INT, control, data, status, n, m,
                              "zero", 0, C_NULL, C_NULL, C_NULL, C_NULL, g, f,
                              "sparse_by_rows", A_ne, C_NULL, A_col, A_ptr, A_val,
                              c_l, c_u, x_l, x_u,
                              n_trans, m_trans, H_ne_trans, A_ne_trans)
    end

    # @printf("%c: n, m, h_ne, a_ne = %2i, %2i, %2i, %2i\n", st, n_trans, m_trans, H_ne_trans, A_ne_trans)
    f_trans = Ref{T}() # transformed constant term in the objective
    H_ptr_trans = zeros(INT, n_trans[] + 1) # transformed Hessian row pointers
    H_col_trans = zeros(INT, H_ne_trans[]) # transformed Hessian column indices
    H_val_trans = zeros(T, H_ne_trans[]) # transformed Hessian values
    g_trans = zeros(T, n_trans[]) # transformed gradient
    A_ptr_trans = zeros(INT, m_trans[] + 1) # transformed Jacobian row pointers
    A_col_trans = zeros(INT, A_ne_trans[]) # transformed Jacobian column indices
    A_val_trans = zeros(T, A_ne_trans[]) # transformed Jacobian values
    x_l_trans = zeros(T, n_trans[]) # transformed lower variable bounds
    x_u_trans = zeros(T, n_trans[]) # transformed upper variable bounds
    c_l_trans = zeros(T, m_trans[]) # transformed lower constraint bounds
    c_u_trans = zeros(T, m_trans[]) # transformed upper constraint bounds
    y_l_trans = zeros(T, m_trans[]) # transformed lower multiplier bounds
    y_u_trans = zeros(T, m_trans[]) # transformed upper multiplier bounds
    z_l_trans = zeros(T, n_trans[]) # transformed lower dual variable bounds
    z_u_trans = zeros(T, n_trans[]) # transformed upper dual variable bounds

    presolve_transform_problem(T, INT, data, status, n_trans[], m_trans[],
                               H_ne_trans[], H_col_trans, H_ptr_trans,
                               H_val_trans, g_trans, f_trans, A_ne_trans[],
                               A_col_trans, A_ptr_trans, A_val_trans,
                               c_l_trans, c_u_trans, x_l_trans, x_u_trans,
                               y_l_trans, y_u_trans, z_l_trans, z_u_trans)

    x_trans = zeros(T, n_trans[]) # transformed variables
    for i in 1:n_trans[]
      x_trans[i] = 0.0
    end
    c_trans = zeros(T, m_trans[]) # transformed constraints
    for i in 1:n_trans[]
      c_trans[i] = 0.0
    end
    y_trans = zeros(T, m_trans[]) # transformed Lagrange multipliers
    for i in 1:n_trans[]
      y_trans[i] = 0.0
    end
    z_trans = zeros(T, n_trans[]) # transformed dual variables
    for i in 1:n_trans[]
      z_trans[i] = 0.0
    end

    x = zeros(T, n) # primal variables
    c = zeros(T, m) # constraint values
    y = zeros(T, m) # Lagrange multipliers
    z = zeros(T, n) # dual variables

    # @printf("%c: n_trans, m_trans, n, m = %2i, %2i, %2i, %2i\n", st, n_trans, m_trans, n, m)
    presolve_restore_solution(T, INT, data, status, n_trans[], m_trans[],
                              x_trans, c_trans, y_trans, z_trans, n, m, x, c, y, z)

    presolve_information(T, INT, data, inform, status)

    if inform[].status == 0
      @printf("%c:%6i transformations, n, m = %2i, %2i, status = %1i\n", st,
              inform[].nbr_transforms, n_trans[], m_trans[], inform[].status)
    else
      @printf("%c: PRESOLVE_solve exit status = %1i\n", st, inform[].status)
    end

    # @printf("x: ")
    # for i = 1:n
    #   @printf("%f ", x[i])
    # end
    # @printf("\n")
    # @printf("gradient: ")
    # for i = 1:n
    #   @printf("%f ", g[i])
    # end
    # @printf("\n")

    # Delete internal workspace
    presolve_terminate(T, INT, data, control, inform)
  end

  return 0
end

for (T, INT, libgalahad) in ((Float32 , Int32, GALAHAD.libgalahad_single      ),
                             (Float32 , Int64, GALAHAD.libgalahad_single_64   ),
                             (Float64 , Int32, GALAHAD.libgalahad_double      ),
                             (Float64 , Int64, GALAHAD.libgalahad_double_64   ),
                             (Float128, Int32, GALAHAD.libgalahad_quadruple   ),
                             (Float128, Int64, GALAHAD.libgalahad_quadruple_64))
  if isfile(libgalahad)
    @testset "PRESOLVE -- $T -- $INT" begin
      @test test_presolve(T, INT) == 0
    end
  end
end

free			\(x\)
non-negativity	0	\(\leq\)	\(x\)
lower	\(x^l\)	\(\leq\)	\(x\)
range	\(x^l\)	\(\leq\)	\(x\)	\(\leq\)	\(x^u\)
upper			\(x\)	\(\leq\)	\(x^u\)
non-positivity			\(x\)	\(\leq\)	0

non-negativity	0	\(\leq\)	\(A x\)
equality	\(c^l\)	\(=\)	\(A x\)
lower	\(c^l\)	\(\leq\)	\(A x\)
range	\(c^l\)	\(\leq\)	\(A x\)	\(\leq\)	\(c^u\)
upper			\(A x\)	\(\leq\)	\(c^u\)
non-positivity			\(A x\)	\(\leq\)	0