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 options['transf_buffer_size'], the transformation buffer size, then they are saved in a “history” file, whose name may be chosen by specifying the options['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}}\)

or

\(\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 options['y_sign'] and options['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 the ``options` 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.

introduction to function calls#

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

See the examples section for illustrations of use.

callable functions#

overview of functions provided#

// namespaces

namespace conf;

// typedefs

typedef float spc_;
typedef double rpc_;
typedef int ipc_;

// structs

struct presolve_control_type;
struct presolve_inform_type;

// global functions

void presolve_initialize(
    void **data,
    struct presolve_control_type* control,
    ipc_ *status
);

void presolve_read_specfile(
    struct presolve_control_type* control,
    const char specfile[]
);

void presolve_import_problem(
    struct presolve_control_type* control,
    void **data,
    ipc_ *status,
    ipc_ n,
    ipc_ m,
    const char H_type[],
    ipc_ H_ne,
    const ipc_ H_row[],
    const ipc_ H_col[],
    const ipc_ H_ptr[],
    const rpc_ H_val[],
    const rpc_ g[],
    const rpc_ f,
    const char A_type[],
    ipc_ A_ne,
    const ipc_ A_row[],
    const ipc_ A_col[],
    const ipc_ A_ptr[],
    const rpc_ A_val[],
    const rpc_ c_l[],
    const rpc_ c_u[],
    const rpc_ x_l[],
    const rpc_ x_u[],
    ipc_ *n_out,
    ipc_ *m_out,
    ipc_ *H_ne_out,
    ipc_ *A_ne_out
);

void presolve_transform_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    ipc_ m,
    ipc_ H_ne,
    ipc_ H_col[],
    ipc_ H_ptr[],
    rpc_ H_val[],
    rpc_ g[],
    rpc_* f,
    ipc_ A_ne,
    ipc_ A_col[],
    ipc_ A_ptr[],
    rpc_ A_val[],
    rpc_ c_l[],
    rpc_ c_u[],
    rpc_ x_l[],
    rpc_ x_u[],
    rpc_ y_l[],
    rpc_ y_u[],
    rpc_ z_l[],
    rpc_ z_u[]
);

void presolve_restore_solution(
    void **data,
    ipc_ *status,
    ipc_ n_in,
    ipc_ m_in,
    const rpc_ x_in[],
    const rpc_ c_in[],
    const rpc_ y_in[],
    const rpc_ z_in[],
    ipc_ n,
    ipc_ m,
    rpc_ x[],
    rpc_ c[],
    rpc_ y[],
    rpc_ z[]
);

void presolve_information(
    void **data,
    struct presolve_inform_type* inform,
    ipc_ *status
);

void presolve_terminate(
    void **data,
    struct presolve_control_type* control,
    struct presolve_inform_type* inform
);

typedefs#

typedef float spc_

spc_ is real single precision

typedef double rpc_

rpc_ is the real working precision used, but may be changed to float by defining the preprocessor variable REAL_32 or (if supported) to __real128 using the variable REAL_128.

typedef int ipc_

ipc_ is the default integer word length used, but may be changed to int64_t by defining the preprocessor variable INTEGER_64.

function calls#

void presolve_initialize(
    void **data,
    struct presolve_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 presolve_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 presolve_read_specfile(
    struct presolve_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/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 struct containing control information (see presolve_control_type)

specfile

is a character string containing the name of the specification file

void presolve_import_problem(
    struct presolve_control_type* control,
    void **data,
    ipc_ *status,
    ipc_ n,
    ipc_ m,
    const char H_type[],
    ipc_ H_ne,
    const ipc_ H_row[],
    const ipc_ H_col[],
    const ipc_ H_ptr[],
    const rpc_ H_val[],
    const rpc_ g[],
    const rpc_ f,
    const char A_type[],
    ipc_ A_ne,
    const ipc_ A_row[],
    const ipc_ A_col[],
    const ipc_ A_ptr[],
    const rpc_ A_val[],
    const rpc_ c_l[],
    const rpc_ c_u[],
    const rpc_ x_l[],
    const rpc_ x_u[],
    ipc_ *n_out,
    ipc_ *m_out,
    ipc_ *H_ne_out,
    ipc_ *A_ne_out
)

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

Parameters:

control

is a struct whose members provide control paramters for the remaining prcedures (see presolve_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:

  • 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 ipc_, that holds the number of variables.

m

is a scalar variable of type ipc_, that holds the number of general linear constraints.

H_type

is a one-dimensional array of type char that specifies the symmetric storage scheme used for the Hessian, \(H\). It should be one of ‘coordinate’, ‘sparse_by_rows’, ‘dense’, ‘diagonal’, ‘scaled_identity’, ‘identity’, ‘zero’ or ‘none’, the latter pair if \(H=0\); lower or upper case variants are allowed.

H_ne

is a scalar variable of type ipc_, that holds the number of entries in the lower triangular part of \(H\) in the sparse co-ordinate storage scheme. It need not be set for any of the other schemes.

H_row

is a one-dimensional array of size H_ne and type ipc_, that holds the row indices of the lower triangular part of \(H\) in the sparse co-ordinate storage scheme. It need not be set for any of the other three schemes, and in this case can be NULL.

H_col

is a one-dimensional array of size H_ne and type ipc_, that holds the column indices of the lower triangular part of \(H\) in either the sparse co-ordinate, or the sparse row-wise storage scheme. It need not be set when the dense, diagonal or (scaled) identity storage schemes are used, and in this case can be NULL.

H_ptr

is a one-dimensional array of size n+1 and type ipc_, that holds the starting position of each row of the lower triangular part of \(H\), as well as the total number of entries, in the sparse row-wise storage scheme. It need not be set when the other schemes are used, and in this case can be NULL.

H_val

is a one-dimensional array of size h_ne and type rpc_, that holds the values of the entries of the lower triangular part of the Hessian matrix \(H\) in any of the available storage schemes.

g

is a one-dimensional array of size n and type rpc_, that holds the linear term \(g\) of the objective function. The j-th component of g, j = 0, … , n-1, contains \(g_j\).

f

is a scalar of type rpc_, that holds the constant term \(f\) of the objective function.

A_type

is a one-dimensional array of type char that specifies the unsymmetric storage scheme used for the constraint Jacobian, \(A\). It should be one of ‘coordinate’, ‘sparse_by_rows’ or ‘dense; lower or upper case variants are allowed.

A_ne

is a scalar variable of type ipc_, that holds the number of entries in \(A\) in the sparse co-ordinate storage scheme. It need not be set for any of the other schemes.

A_row

is a one-dimensional array of size A_ne and type ipc_, that holds the row indices of \(A\) in the sparse co-ordinate storage scheme. It need not be set for any of the other schemes, and in this case can be NULL.

A_col

is a one-dimensional array of size A_ne and type ipc_, that holds the column indices of \(A\) in either the sparse co-ordinate, or the sparse row-wise storage scheme. It need not be set when the dense or diagonal storage schemes are used, and in this case can be NULL.

A_ptr

is a one-dimensional array of size n+1 and type ipc_, that holds the starting position of each row of \(A\), as well as the total number of entries, in the sparse row-wise storage scheme. It need not be set when the other schemes are used, and in this case can be NULL.

A_val

is a one-dimensional array of size a_ne and type rpc_, that holds the values of the entries of the constraint Jacobian matrix \(A\) in any of the available storage schemes.

c_l

is a one-dimensional array of size m and type rpc_, that holds the lower bounds \(c^l\) on the constraints \(A x\). The i-th component of c_l, i = 0, … , m-1, contains \(c^l_i\).

c_u

is a one-dimensional array of size m and type rpc_, that holds the upper bounds \(c^l\) on the constraints \(A x\). The i-th component of c_u, i = 0, … , m-1, contains \(c^u_i\).

x_l

is a one-dimensional array of size n and type rpc_, that holds the lower bounds \(x^l\) on the variables \(x\). The j-th component of x_l, j = 0, … , n-1, contains \(x^l_j\).

x_u

is a one-dimensional array of size n and type rpc_, that holds the upper bounds \(x^l\) on the variables \(x\). The j-th component of x_u, j = 0, … , n-1, contains \(x^l_j\).

n_out

is a scalar variable of type ipc_, that holds the number of variables in the transformed problem.

m_out

is a scalar variable of type ipc_, that holds the number of general linear constraints in the transformed problem.

H_ne_out

is a scalar variable of type ipc_, 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 ipc_, that holds the number of entries in \(A\) in the transformed problem.

void presolve_transform_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    ipc_ m,
    ipc_ H_ne,
    ipc_ H_col[],
    ipc_ H_ptr[],
    rpc_ H_val[],
    rpc_ g[],
    rpc_* f,
    ipc_ A_ne,
    ipc_ A_col[],
    ipc_ A_ptr[],
    rpc_ A_val[],
    rpc_ c_l[],
    rpc_ c_u[],
    rpc_ x_l[],
    rpc_ x_u[],
    rpc_ y_l[],
    rpc_ y_u[],
    rpc_ z_l[],
    rpc_ 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 ipc_, 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 ipc_, 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 ipc_, 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 ipc_, 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 ipc_, 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 ipc_, 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 rpc_, 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 rpc_, that holds the the transformed 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 transformed constant term \(f\) of the objective function.

A_ne

is a scalar variable of type ipc_, 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 ipc_, 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 ipc_, 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 rpc_, 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 rpc_, that holds the transformed lower bounds \(c^l\) on the constraints \(A x\). The i-th component of c_l, i = 0, … , m-1, contains \(c^l_i\).

c_u

is a one-dimensional array of size m and type rpc_, that holds the transformed upper bounds \(c^l\) on the constraints \(A x\). The i-th component of c_u, i = 0, … , m-1, contains \(c^u_i\).

x_l

is a one-dimensional array of size n and type rpc_, that holds the transformed 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 transformed upper bounds \(x^l\) on the variables \(x\). The j-th component of x_u, j = 0, … , n-1, contains \(x^l_j\).

y_l

is a one-dimensional array of size m and type rpc_, that holds the implied lower bounds \(y^l\) on the transformed Lagrange multipliers \(y\). The i-th component of y_l, i = 0, … , m-1, contains \(y^l_i\).

y_u

is a one-dimensional array of size m and type rpc_, that holds the implied upper bounds \(y^u\) on the transformed Lagrange multipliers \(y\). The i-th component of y_u, i = 0, … , m-1, contains \(y^u_i\).

z_l

is a one-dimensional array of size m and type rpc_, that holds the implied lower bounds \(y^l\) on the transformed dual variables \(z\). The j-th component of z_l, j = 0, … , n-1, contains \(z^l_i\).

z_u

is a one-dimensional array of size m and type rpc_, that holds the implied upper bounds \(y^u\) on the transformed dual variables \(z\). The j-th component of z_u, j = 0, … , n-1, contains \(z^u_i\).

void presolve_restore_solution(
    void **data,
    ipc_ *status,
    ipc_ n_in,
    ipc_ m_in,
    const rpc_ x_in[],
    const rpc_ c_in[],
    const rpc_ y_in[],
    const rpc_ z_in[],
    ipc_ n,
    ipc_ m,
    rpc_ x[],
    rpc_ c[],
    rpc_ y[],
    rpc_ 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 ipc_, 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 ipc_, 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 ipc_, 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 rpc_, that holds the transformed values \(x\) of the optimization variables. The j-th component of x, j = 0, … , n-1, contains \(x_j\).

c_in

is a one-dimensional array of size m and type rpc_, that holds the transformed residual \(c(x)\). The i-th component of c, j = 0, … , n-1, contains \(c_j(x)\).

y_in

is a one-dimensional array of size n_in and type rpc_, that holds the values \(y\) of the transformed Lagrange multipliers for the general linear constraints. The j-th component of y, j = 0, … , n-1, contains \(y_j\).

z_in

is a one-dimensional array of size n_in and type rpc_, that holds the values \(z\) of the transformed dual variables. The j-th component of z, j = 0, … , n-1, contains \(z_j\).

n

is a scalar variable of type ipc_, 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 ipc_, 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 rpc_, that holds the transformed values \(x\) of the optimization variables. The j-th component of x, j = 0, … , n-1, contains \(x_j\).

c

is a one-dimensional array of size m and type rpc_, that holds the transformed residual \(c(x)\). The i-th component of c, j = 0, … , n-1, contains \(c_j(x)\).

y

is a one-dimensional array of size n and type rpc_, that holds the values \(y\) of the transformed Lagrange multipliers for the general linear constraints. The j-th component of y, j = 0, … , n-1, contains \(y_j\).

z

is a one-dimensional array of size n and type rpc_, that holds the values \(z\) of the transformed dual variables. The j-th component of z, j = 0, … , n-1, contains \(z_j\).

void presolve_information(
    void **data,
    struct presolve_inform_type* inform,
    ipc_ *status
)

Provides output information

Parameters:

data

holds private internal data

inform

is a struct containing output information (see presolve_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 presolve_terminate(
    void **data,
    struct presolve_control_type* control,
    struct presolve_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data

holds private internal data

control

is a struct containing control information (see presolve_control_type)

inform

is a struct containing output information (see presolve_inform_type)

available structures#

presolve_control_type structure#

control derived type as a C struct More…

#include <galahad_presolve.h>

struct presolve_control_type {
    // fields

    bool f_indexing;
    ipc_ termination;
    ipc_ max_nbr_transforms;
    ipc_ max_nbr_passes;
    rpc_ c_accuracy;
    rpc_ z_accuracy;
    rpc_ infinity;
    ipc_ out;
    ipc_ errout;
    ipc_ print_level;
    bool dual_transformations;
    bool redundant_xc;
    ipc_ primal_constraints_freq;
    ipc_ dual_constraints_freq;
    ipc_ singleton_columns_freq;
    ipc_ doubleton_columns_freq;
    ipc_ unc_variables_freq;
    ipc_ dependent_variables_freq;
    ipc_ sparsify_rows_freq;
    ipc_ max_fill;
    ipc_ transf_file_nbr;
    ipc_ transf_buffer_size;
    ipc_ transf_file_status;
    char transf_file_name[31];
    ipc_ y_sign;
    ipc_ inactive_y;
    ipc_ z_sign;
    ipc_ inactive_z;
    ipc_ final_x_bounds;
    ipc_ final_z_bounds;
    ipc_ final_c_bounds;
    ipc_ final_y_bounds;
    ipc_ check_primal_feasibility;
    ipc_ check_dual_feasibility;
    rpc_ pivot_tol;
    rpc_ min_rel_improve;
    rpc_ max_growth_factor;
};

detailed documentation#

control derived type as a C struct

components#

bool f_indexing

use C or Fortran sparse matrix indexing

ipc_ 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.

ipc_ max_nbr_transforms

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

ipc_ max_nbr_passes

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

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

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

rpc_ infinity

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

ipc_ out

The unit number associated with the device used for printout.

ipc_ errout

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

ipc_ 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.

ipc_ 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.

ipc_ 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.

ipc_ 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.

ipc_ 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.

ipc_ 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.

ipc_ 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.

ipc_ 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.

ipc_ 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.

ipc_ transf_file_nbr

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

ipc_ transf_buffer_size

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

ipc_ 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.

ipc_ 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.

ipc_ 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.

ipc_ 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.

ipc_ 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.

ipc_ 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).

ipc_ 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).

ipc_ 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.

ipc_ 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).

ipc_ 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.

ipc_ 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.

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

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

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

#include <galahad_presolve.h>

struct presolve_inform_type {
    // fields

    ipc_ status;
    ipc_ status_continue;
    ipc_ status_continued;
    ipc_ nbr_transforms;
    char message[3][81];
};

detailed documentation#

inform derived type as a C struct

components#

ipc_ 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

ipc_ 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.

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/C/presolvet.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.

/* presolvet.c */
/* Full test for the PRESOLVE C interface using C sparse matrix indexing */

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

int main(void) {

    // Derived types
    void *data;
    struct presolve_control_type control;
    struct presolve_inform_type inform;

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

    // Set output storage
    char st = ' ';
    ipc_ status;

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

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

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

      ipc_ n_trans, m_trans, H_ne_trans, A_ne_trans;

        // Initialize PRESOLVE
        presolve_initialize( &data, &control, &status );

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

        switch(d){
            case 1: // sparse co-ordinate storage
                st = 'C';
                presolve_import_problem( &control, &data, &status, n, m,
                           "coordinate", H_ne, H_row, H_col, NULL, H_val, g, f,
                           "coordinate", A_ne, A_row, A_col, NULL, A_val,
                           c_l, c_u, x_l, x_u,
                           &n_trans, &m_trans, &H_ne_trans, &A_ne_trans );
                break;
            printf(" case %1" i_ipc_ " break\n",d);
            case 2: // sparse by rows
                st = 'R';
                presolve_import_problem( &control, &data, &status, n, m,
                       "sparse_by_rows", H_ne, NULL, H_col, H_ptr, H_val, g, f,
                       "sparse_by_rows", A_ne, 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 );
                break;
            case 3: // dense
                st = 'D';
                ipc_ H_dense_ne = n*(n+1)/2; // number of elements of H
                ipc_ A_dense_ne = m*n; // number of elements of A
                rpc_ H_dense[] = {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};
                rpc_ A_dense[] = {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( &control, &data, &status, n, m,
                            "dense", H_dense_ne, NULL, NULL, NULL, H_dense, g,
                            f, "dense", A_dense_ne, NULL, NULL, NULL, A_dense,
                            c_l, c_u, x_l, x_u,
                            &n_trans, &m_trans, &H_ne_trans, &A_ne_trans );
                break;
            case 4: // diagonal
                st = 'L';
                presolve_import_problem( &control, &data, &status, n, m,
                            "diagonal", n, NULL, NULL, NULL, H_val, g, f,
                            "sparse_by_rows", A_ne, 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 );
                break;

            case 5: // scaled identity
                st = 'S';
                presolve_import_problem( &control, &data, &status, n, m,
                        "scaled_identity", 1, NULL, NULL, NULL, H_val, g, f,
                        "sparse_by_rows", A_ne, 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 );
                break;
            case 6: // identity
                st = 'I';
                presolve_import_problem( &control, &data, &status, n, m,
                            "identity", 0, NULL, NULL, NULL, NULL, g, f,
                            "sparse_by_rows", A_ne, 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 );
                break;
            case 7: // zero
                st = 'Z';
                presolve_import_problem( &control, &data, &status, n, m,
                            "zero", 0, NULL, NULL, NULL, NULL, g, f,
                            "sparse_by_rows", A_ne, 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 );
                break;
            }

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

        presolve_transform_problem( &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 );

        rpc_ x_trans[n_trans]; // transformed variables
        for( ipc_ i = 0; i < n_trans; i++) x_trans[i] = 0.0;
        rpc_ c_trans[m_trans]; // transformed constraints
        for( ipc_ i = 0; i < m_trans; i++) c_trans[i] = 0.0;
        rpc_ y_trans[m_trans]; // transformed Lagrange multipliers
        for( ipc_ i = 0; i < m_trans; i++) y_trans[i] = 0.0;
        rpc_ z_trans[n_trans]; // transformed dual variables
        for( ipc_ i = 0; i < n_trans; i++) z_trans[i] = 0.0;

        rpc_ x[n]; // primal variables
        rpc_ c[m]; // constraint values
        rpc_ y[m]; // Lagrange multipliers
        rpc_ z[n]; // dual variables

        //printf("%c: n_trans, m_trans, n, m = %2" i_ipc_ ", %2" i_ipc_ ", %2" i_ipc_ ", %2" i_ipc_ "\n",
        //           st, n_trans, m_trans, n, m );
        presolve_restore_solution( &data, &status, n_trans, m_trans,
                  x_trans, c_trans, y_trans, z_trans, n, m, x, c, y, z );

        presolve_information( &data, &inform, &status );

        if(inform.status == 0){
            printf("%c:%6" i_ipc_ " transformations, n, m = %2" i_ipc_ ", %2" i_ipc_ ", status = %1" i_ipc_ "\n",
                   st, inform.nbr_transforms, n_trans, m_trans, inform.status);
        }else{
            printf("%c: PRESOLVE_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
        presolve_terminate( &data, &control, &inform );
    }
}

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

/* presolvetf.c */
/* Full test for the PRESOLVE C interface using Fortran sparse matrix indexing */

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

int main(void) {

    // Derived types
    void *data;
    struct presolve_control_type control;
    struct presolve_inform_type inform;

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

    // Set output storage
    char st = ' ';
    ipc_ status;

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

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

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

      ipc_ n_trans, m_trans, H_ne_trans, A_ne_trans;

        // Initialize PRESOLVE
        presolve_initialize( &data, &control, &status );

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

        switch(d){
            case 1: // sparse co-ordinate storage
                st = 'C';
                presolve_import_problem( &control, &data, &status, n, m,
                           "coordinate", H_ne, H_row, H_col, NULL, H_val, g, f,
                           "coordinate", A_ne, A_row, A_col, NULL, A_val,
                           c_l, c_u, x_l, x_u,
                           &n_trans, &m_trans, &H_ne_trans, &A_ne_trans );
                break;
            printf(" case %1" i_ipc_ " break\n",d);
            case 2: // sparse by rows
                st = 'R';
                presolve_import_problem( &control, &data, &status, n, m,
                       "sparse_by_rows", H_ne, NULL, H_col, H_ptr, H_val, g, f,
                       "sparse_by_rows", A_ne, 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 );
                break;
            case 3: // dense
                st = 'D';
                ipc_ H_dense_ne = n*(n+1)/2; // number of elements of H
                ipc_ A_dense_ne = m*n; // number of elements of A
                rpc_ H_dense[] = {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};
                rpc_ A_dense[] = {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( &control, &data, &status, n, m,
                            "dense", H_dense_ne, NULL, NULL, NULL, H_dense, g,
                            f, "dense", A_dense_ne, NULL, NULL, NULL, A_dense,
                            c_l, c_u, x_l, x_u,
                            &n_trans, &m_trans, &H_ne_trans, &A_ne_trans );
                break;
            case 4: // diagonal
                st = 'L';
                presolve_import_problem( &control, &data, &status, n, m,
                            "diagonal", n, NULL, NULL, NULL, H_val, g, f,
                            "sparse_by_rows", A_ne, 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 );
                break;

            case 5: // scaled identity
                st = 'S';
                presolve_import_problem( &control, &data, &status, n, m,
                        "scaled_identity", 1, NULL, NULL, NULL, H_val, g, f,
                        "sparse_by_rows", A_ne, 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 );
                break;
            case 6: // identity
                st = 'I';
                presolve_import_problem( &control, &data, &status, n, m,
                            "identity", 0, NULL, NULL, NULL, NULL, g, f,
                            "sparse_by_rows", A_ne, 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 );
                break;
            case 7: // zero
                st = 'Z';
                presolve_import_problem( &control, &data, &status, n, m,
                            "zero", 0, NULL, NULL, NULL, NULL, g, f,
                            "sparse_by_rows", A_ne, 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 );
                break;
            }

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

        presolve_transform_problem( &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 );

        rpc_ x_trans[n_trans]; // transformed variables
        for( ipc_ i = 0; i < n_trans; i++) x_trans[i] = 0.0;
        rpc_ c_trans[m_trans]; // transformed constraints
        for( ipc_ i = 0; i < m_trans; i++) c_trans[i] = 0.0;
        rpc_ y_trans[m_trans]; // transformed Lagrange multipliers
        for( ipc_ i = 0; i < m_trans; i++) y_trans[i] = 0.0;
        rpc_ z_trans[n_trans]; // transformed dual variables
        for( ipc_ i = 0; i < n_trans; i++) z_trans[i] = 0.0;

        rpc_ x[n]; // primal variables
        rpc_ c[m]; // constraint values
        rpc_ y[m]; // Lagrange multipliers
        rpc_ z[n]; // dual variables

        //printf("%c: n_trans, m_trans, n, m = %2" i_ipc_ ", %2" i_ipc_ ", %2" i_ipc_ ", %2" i_ipc_ "\n",
        //           st, n_trans, m_trans, n, m );
        presolve_restore_solution( &data, &status, n_trans, m_trans,
                  x_trans, c_trans, y_trans, z_trans, n, m, x, c, y, z );

        presolve_information( &data, &inform, &status );

        if(inform.status == 0){
            printf("%c:%6" i_ipc_ " transformations, n, m = %2" i_ipc_ ", %2" i_ipc_ ", status = %1" i_ipc_ "\n",
                   st, inform.nbr_transforms, n_trans, m_trans, inform.status);
        }else{
            printf("%c: PRESOLVE_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
        presolve_terminate( &data, &control, &inform );
    }
}