GALAHAD CRO package#

purpose#

The cro package provides a crossover from a primal-dual interior-point solution to given convex quadratic program to a basic one in which the matrix of defining active constraints/variables is of full rank. This applies to the problem of minimizing the quadratic objective function

\[q(x) = f + g^T x + \frac{1}{2} x^T H x\]
subject to the general linear constraints and simple bounds
\[c_l \leq A x \leq c_u \;\;\mbox{and} \;\; x_l \leq x \leq x_u,\]
where \(H\) and \(A\) are, respectively, given \(n\) by \(n\) symmetric postive-semi-definite and \(m\) by \(n\) matrices, \(g\) is a vector, \(f\) is a scalar, and any of the components of the vectors \(c_l\), \(c_u\), \(x_l\) or \(x_u\) may be infinite. The method is most suitable for problems involving a large number of unknowns \(x\).

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

terminology#

Any required solution \(x\) 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,\;\;\mbox{(1)}\]
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 component-wise.

method#

Denote the active constraints by \(A_A x = c_A\) and the active bounds by \(I_A x = x_A\). Then any optimal solution satisfies the linear system

\[\begin{split}\left(\begin{array}{ccc} H & - A_A^T & - I^T_A \\ A_A & 0 & 0 \\ I_A & 0 & 0 \end{array}\right) \left(\begin{array}{c} x \\ y_A \\ z_A \end{array}\right) = \left(\begin{array}{c} - g \\ c_A \\ x_A \end{array}\right), \end{split}\]
where \(y_A\) and \(z_A\) are the corresponding active Lagrange multipliers and dual variables respectively. Consequently the difference between any two solutions \((\Delta x, \Delta y, \Delta z)\) must satisfy
\[\begin{split}\left(\begin{array}{ccc} H & - A_A^T & - I^T_A \\ A_A & 0 & 0 \\ I_A & 0 & 0 \end{array}\right) \left(\begin{array}{c} \Delta x \\ \Delta y_A \\ \Delta z_A \end{array}\right) = 0.\;\;\mbox{(2)}\end{split}\]
Thus there can only be multiple solution if the coefficient matrix \(K\) of (2) is singular. The algorithm used in \fullpackagename\ exploits this. The matrix \(K\) is checked for singularity using the package ULS. If \(K\) is non singular, the solution is unique and the solution input by the user provides a linearly independent active set. Otherwise \(K\) is singular, and partitions \(A_A^T = ( A_{AB}^T \;\; A_{AN}^T)\) and \(I_A^T = ( I_{AB}^T \;\; I_{AN}^T)\) are found so that
\[\begin{split}\left(\begin{array}{ccc} H & - A_{AB}^T & - I^T_{AB} \\ A_{AB} & 0 & 0 \\ I_{AB} & 0 & 0 \end{array}\right)\end{split}\]
is non-singular and the non-basic constraints \(A_{AN}^T\) and \(I_{AN}^T\) are linearly dependent on the basic ones \(( A_{AB}^T \;\; I_{AB}^T)\). In this case (2) is equivalent to
\[\begin{split}\left(\begin{array}{ccc} H & - A_{AB}^T & - I^T_{AB} \\ A_{AB} & 0 & 0 \\ I_{AB} & 0 & 0 \end{array}\right) \left(\begin{array}{c} \Delta x \\ \Delta y_{AB} \\ \Delta z_{AB} \end{array}\right) = \left(\begin{array}{c} A_{AN}^T \\ 0 \\ 0 \end{array}\right) \Delta y_{AN} + \left(\begin{array}{c} I^T_{AN} \\ 0 \\ 0 \end{array}\right) \Delta z_{AN} .\;\;\mbox{(3)}\end{split}\]
Thus, starting from the user’s \((x, y, z)\) and with a factorization of the coefficient matrix of (3) found by the package SLS, the alternative solution \((x + \alpha x, y + \alpha y, z + \alpha z)\), featuring \((\Delta x, \Delta y_{AB}, \Delta z_{AB})\) from (3) in which successively one of the components of \(\Delta y_{AN}\) and \(\Delta z_{AN}\) in turn is non zero, is taken. The scalar \(\alpha\) at each stage is chosen to be the largest possible that guarantees (1); this may happen when a non-basic multiplier/dual variable reaches zero, in which case the corresponding constraint is disregarded, or when this happens for a basic multiplier/dual variable, in which case this constraint is exchanged with the non-basic one under consideration and disregarded. The latter corresponds to changing the basic-non-basic partition in (3), and subsequent solutions may be found by updating the factorization of the coefficient matrix in (3) following the basic-non-basic swap using the package SCU.

matrix storage#

The unsymmetric \(m\) by \(n\) matrix \(A\) must be presented and stored in sparse row-wise storage format. For this, only the nonzero entries are stored, and they are ordered so that those in row i appear directly before those in row i+1. For the i-th row of \(A\) the i-th component of the integer array A_ptr holds the position of the first entry in this row, while A_ptr(m) holds the total number of entries. The column indices j, \(0 \leq j \leq n-1\), and values \(A_{ij}\) of the nonzero entries in the i-th row are stored in components l = A_ptr(i), \(\ldots\), A_ptr(i+1)-1, \(0 \leq i \leq m-1\), of the integer array A_col, and real array A_val, respectively.

The symmetric \(n\) by \(n\) matrix \(H\) must also be presented and stored in sparse row-wise storage format. But, crucially, now symmetry is exploited by only storing values from the lower triangular part (i.e, those entries that lie on or below the leading diagonal). As before, only the nonzero entries of the matrices are stored. Only the nonzero entries from the lower triangle are stored, and these are ordered so that those in row i appear directly before those in row i+1. For the i-th row of \(H\) the i-th component of the integer array H_ptr holds the position of the first entry in this row, while H_ptr(n) holds the total number of entries. The column indices j, \(0 \leq j \leq i\), and values \(H_{ij}\) of the entries in the i-th row are stored in components l = H_ptr(i), …, H_ptr(i+1)-1 of the integer array H_col, and real array H_val, respectively.

introduction to function calls#

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

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

See the examples section for illustrations of use.

callable functions#

overview of functions provided#

// typedefs

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

// structs

struct cro_control_type;
struct cro_inform_type;
struct cro_time_type;

// function calls

void cro_initialize(void **data, struct cro_control_type* control, ipc_ *status);
void cro_read_specfile(struct cro_control_type* control, const char specfile[]);

void cro_crossover_solution(
    void **data,
    struct cro_control_type* control,
    struct cro_inform_type* inform,
    ipc_ n,
    ipc_ m,
    ipc_ m_equal,
    ipc_ h_ne,
    const rpc_ H_val[],
    const ipc_ H_col[],
    const ipc_ H_ptr[],
    ipc_ a_ne,
    const rpc_ A_val[],
    const ipc_ A_col[],
    const ipc_ A_ptr[],
    const rpc_ g[],
    const rpc_ c_l[],
    const rpc_ c_u[],
    const rpc_ x_l[],
    const rpc_ x_u[],
    rpc_ x[],
    rpc_ c[],
    rpc_ y[],
    rpc_ z[],
    ipc_ x_stat[],
    ipc_ c_stat[]
);

void cro_terminate(
    void **data,
    struct cro_control_type* control,
    struct cro_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 cro_initialize(void **data, struct cro_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 cro_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 cro_read_specfile(struct cro_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/cro/CRO.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/cro.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 cro_control_type)

specfile

is a character string containing the name of the specification file

 
void cro_crossover_solution(
    void **data,
    struct cro_control_type* control,
    struct cro_inform_type* inform,
    ipc_ n,
    ipc_ m,
    ipc_ m_equal,
    ipc_ h_ne,
    const rpc_ H_val[],
    const ipc_ H_col[],
    const ipc_ H_ptr[],
    ipc_ a_ne,
    const rpc_ A_val[],
    const ipc_ A_col[],
    const ipc_ A_ptr[],
    const rpc_ g[],
    const rpc_ c_l[],
    const rpc_ c_u[],
    const rpc_ x_l[],
    const rpc_ x_u[],
    rpc_ x[],
    rpc_ c[],
    rpc_ y[],
    rpc_ z[],
    ipc_ x_stat[],
    ipc_ c_stat[]
)

Crosover the solution from a primal-dual to a basic one.

Parameters:

control

is a struct whose members provide control paramters for the remaining prcedures (see cro_control_type). The parameter .status is as follows:

data

holds private internal data.

inform

is a struct containing output information (see cro_inform_type). The component .status gives the exit status from the package. Possible values are:

  • 0

    The crossover 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 >= m_equal >= 0 has been violated.

  • -4

    the bound constraints are inconsistent.

  • -5

    the general constraints are likely inconsistent.

  • -9

    an error has occured in SLS_analyse.

  • -10

    an error has occured in SLS_factorize.

  • -11

    an error has occured in SLS_solve.

  • -12

    an error has occured in ULS_factorize.

  • -14

    an error has occured in ULS_solve.

  • -16

    the residuals are large; the factorization may be unsatisfactory.

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.

m_equal

is a scalar variable of type ipc_, that holds the number of general linear equality constraints. Such constraints must occur first in \(A\).

h_ne

is a scalar variable of type ipc_, that holds the number of entries in the lower triangular part of the Hessian matrix \(H\).

H_val

is a one-dimensional array of type rpc_, that holds the values of the entries of the lower triangular part of the Hessian matrix \(H\). The entries are stored by consecutive rows, the order within each row is unimportant.

H_col

is a one-dimensional array of type ipc_, that holds the column indices of the lower triangular part of \(H\), in the same order as those in H_val.

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\). The n+1-st component holds the total number of entries (plus one if fortran indexing is used).

a_ne

is a scalar variable of type ipc_, that holds the number of entries in the constraint Jacobian matrix \(A\).

A_val

is a one-dimensional array of type rpc_, that holds the values of the entries of the constraint Jacobian matrix \(A\). The entries are stored by consecutive rows, the order within each row is unimportant. Equality constraints must be ordered first.

A_col

is a one-dimensional array of size A_ne and type ipc_, that holds the column indices of \(A\) in the same order as those in A_val.

A_ptr

is a one-dimensional array of size m+1 and type ipc_, that holds the starting position of each row of \(A\). The m+1-st component holds the total number of entries (plus one if fortran indexing is used).

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

c_l

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

c_u

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

x_l

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

x_u

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

x

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

c

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

y

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

z

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

x_stat

is a one-dimensional array of size n and type ipc_, that must be set on entry to give the status of the problem variables. If x_stat(j) is negative, the variable \(x_j\) is active on its lower bound, if it is positive, it is active and lies on its upper bound, and if it is zero, it is inactiive and lies between its bounds. On exit, the \(j\) -th component of x_stat is -1 if the variable is basic and active on its lower bound, -2 it is non-basic but active on its lower bound, 1 if it is basic and active on its upper bound, 2 it is non-basic but active on its upper bound, and 0 if it is inactive.

c_stat

is a one-dimensional array of size m and type ipc_, that must be set on entry to give the status of the general linear constraints. If c_stat(i) is negative, the constraint value \(a_i^Tx\) is active on its lower bound, if it is positive, it is active and lies on its upper bound, and if it is zero, it is inactiive and lies between its bounds. On exit, the \(i\) -th component of x_stat is -1 if the constraint is basic and active on its lower bound, -2 it is non-basic but active on its lower bound, 1 if it is basic and active on its upper bound, 2 it is non-basic but active on its upper bound, and 0 if it is inactive.

 
void cro_terminate(
    void **data,
    struct cro_control_type* control,
    struct cro_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data

holds private internal data

control

is a struct containing control information (see cro_control_type)

inform

is a struct containing output information (see cro_inform_type)

available structures#

cro_control_type structure#

#include <galahad_cro.h>

struct cro_control_type {
    // components

    bool f_indexing;
    ipc_ error;
    ipc_ out;
    ipc_ print_level;
    ipc_ max_schur_complement;
    rpc_ infinity;
    rpc_ feasibility_tolerance;
    bool check_io;
    bool refine_solution;
    bool space_critical;
    bool deallocate_error_fatal;
    char symmetric_linear_solver[31];
    char unsymmetric_linear_solver[31];
    char prefix[31];
    struct sls_control_type sls_control;
    struct sbls_control_type sbls_control;
    struct uls_control_type uls_control;
    struct ir_control_type ir_control;
};

detailed documentation#

control derived type as a C struct

components#

bool f_indexing

use C or Fortran sparse matrix indexing

ipc_ error

error and warning diagnostics occur on stream error

ipc_ out

general output occurs on stream out

ipc_ print_level

the level of output required is specified by print_level

ipc_ max_schur_complement

the maximum permitted size of the Schur complement before a refactorization is performed

rpc_ infinity

any bound larger than infinity in modulus will be regarded as infinite

rpc_ feasibility_tolerance

feasibility tolerance for KKT violation

bool check_io

if .check_io is true, the input (x,y,z) will be fully tested for consistency

bool refine_solution

if .refine solution is true, attempt to satisfy the KKT conditions as accurately as possible

bool space_critical

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

bool deallocate_error_fatal

if .deallocate_error_fatal is true, any array/pointer deallocation error will terminate execution. Otherwise, computation will continue

char symmetric_linear_solver[31]

the name of the symmetric-indefinite linear equation solver used. Possible choices are currently: ‘sils’, ‘ma27’, ‘ma57’, ‘ma77’, ‘ma86’, ‘ma97’, ‘ssids’, ‘mumps’, ‘pardiso’, ‘mkl_pardiso’, ‘pastix’, ‘wsmp’, and ‘sytr’, although only ‘sytr’ and, for OMP 4.0-compliant compilers, ‘ssids’ are installed by default; others are easily installed (see README.external). More details of the capabilities of each solver are provided in the documentation for galahad_sls.

char unsymmetric_linear_solver[31]

the name of the unsymmetric linear equation solver used. Possible choices are currently: ‘gls’, ‘ma48’ and ‘getr’, although only ‘getr’ is installed by default; others are easily installed (see README.external). More details of the capabilities of each solver are provided in the documentation for galahad_uls.

char prefix[31]

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

struct sls_control_type sls_control

control parameters for SLS

struct sbls_control_type sbls_control

control parameters for SBLS

struct uls_control_type uls_control

control parameters for ULS

struct ir_control_type ir_control

control parameters for iterative refinement

cro_time_type structure#

#include <galahad_cro.h>

struct cro_time_type {
    // components

    spc_ total;
    spc_ analyse;
    spc_ factorize;
    spc_ solve;
    rpc_ clock_total;
    rpc_ clock_analyse;
    rpc_ clock_factorize;
    rpc_ clock_solve;
};

detailed documentation#

time derived type as a C struct

components#

spc_ total

the total CPU time spent in the package

spc_ analyse

the CPU time spent reordering the matrix prior to factorization

spc_ factorize

the CPU time spent factorizing the required matrices

spc_ solve

the CPU time spent computing corrections

rpc_ clock_total

the total clock time spent in the package

rpc_ clock_analyse

the clock time spent analysing the required matrices prior to factorizat

rpc_ clock_factorize

the clock time spent factorizing the required matrices

rpc_ clock_solve

the clock time spent computing corrections

cro_inform_type structure#

#include <galahad_cro.h>

struct cro_inform_type {
    // components

    ipc_ status;
    ipc_ alloc_status;
    char bad_alloc[81];
    ipc_ dependent;
    struct cro_time_type time;
    struct sls_inform_type sls_inform;
    struct sbls_inform_type sbls_inform;
    struct uls_inform_type uls_inform;
    ipc_ scu_status;
    struct scu_inform_type scu_inform;
    struct ir_inform_type ir_inform;
};

detailed documentation#

inform derived type as a C struct

components#

ipc_ status

return status. See CRO_solve for details

ipc_ alloc_status

the status of the last attempted allocation/deallocation

char bad_alloc[81]

the name of the array for which an allocation/deallocation error occurred

ipc_ dependent

the number of dependent active constraints

struct cro_time_type time

timings (see above)

struct sls_inform_type sls_inform

information from SLS

struct sbls_inform_type sbls_inform

information from SBLS

struct uls_inform_type uls_inform

information from ULS

ipc_ scu_status

information from SCU

struct scu_inform_type scu_inform

see scu_status

struct ir_inform_type ir_inform

information from IR

Index

example calls#

This is an example of how to use the package to crossover from a primal-dual QP solution to a basic one; the code is available in $GALAHAD/src/cro/C/crot.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.

/* crot.c */
/* Full test for the CRO C interface using C sparse matrix indexing */

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

int main(void) {

    // Derived types
    void *data;
    struct cro_control_type control;
    struct cro_inform_type inform;


    // Set problem dimensions
    ipc_ n = 11; // dimension
    ipc_ m = 3; // number of general constraints
    ipc_ m_equal = 1; // number of equality constraints

    //  describe the objective function

    ipc_ H_ne = 21;
    rpc_ H_val[] = {1.0,0.5,1.0,0.5,1.0,0.5,1.0,0.5,1.0,0.5,
                        1.0,0.5,1.0,0.5,1.0,0.5,1.0,0.5,1.0,0.5,1.0};
    ipc_ H_col[] = {0,0,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10};
    ipc_ H_ptr[] = {0,1,3,5,7,9,11,13,15,17,19,21};
    rpc_ g[] = {0.5,-0.5,-1.0,-1.0,-1.0, -1.0,-1.0,-1.0,-1.0,-1.0,-0.5};

    //  describe constraints

    ipc_ A_ne = 30;
    rpc_ A_val[] = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
                        1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
                        1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
    ipc_ A_col[] = {0,1,2,3,4,5,6,7,8,9,10,2,3,4,5,6,7,8,9,10,
                   1,2,3,4,5,6,7,8,9,10};
    ipc_ A_ptr[] = {0,11,20,30};
    rpc_ c_l[] = {10.0,9.0,-INFINITY};
    rpc_ c_u[] = {10.0,INFINITY,10.0};
    rpc_ x_l[] = {0.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
    rpc_ x_u[] = {INFINITY,INFINITY,INFINITY,INFINITY,INFINITY,INFINITY,
                      INFINITY,INFINITY,INFINITY,INFINITY,INFINITY};

    // provide optimal variables, Lagrange multipliers and dual variables
    rpc_ x[] = {0.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.0,1.0};
    rpc_ c[] = {10.0,9.0,10.0};
    rpc_ y[] = { -1.0,1.5,-2.0};
    rpc_ z[] = {2.0,4.0,2.5,2.5,2.5,2.5,2.5,2.5,2.5,2.5,2.5};
    // provide interior-point constraint and variable status
    ipc_ c_stat[] = {-1,-1,1};
    ipc_ x_stat[] = {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};

    // Set output storage

    ipc_ status;

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

    // Initialize CRO
    cro_initialize( &data, &control, &status );

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

    // crossover the solution
    cro_crossover_solution( &data, &control, &inform,
                            n, m, m_equal,
                            H_ne, H_val, H_col, H_ptr,
                            A_ne, A_val, A_col, A_ptr,
                            g, c_l, c_u, x_l, x_u, x, c, y, z,
                            x_stat, c_stat );

    printf(" CRO_crossover exit status = %1" i_ipc_ "\n", inform.status);

    // Delete internal workspace
    cro_terminate( &data, &control, &inform );

}

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

/* crotf.c */
/* Full test for the CRO C interface using C sparse matrix indexing */

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

int main(void) {

    // Derived types
    void *data;
    struct cro_control_type control;
    struct cro_inform_type inform;

    // Set problem dimensions
    ipc_ n = 11; // dimension
    ipc_ m = 3; // number of general constraints
    ipc_ m_equal = 1; // number of equality constraints

    //  describe the objective function

    ipc_ H_ne = 21;
    rpc_ H_val[] = {1.0,0.5,1.0,0.5,1.0,0.5,1.0,0.5,1.0,0.5,
                        1.0,0.5,1.0,0.5,1.0,0.5,1.0,0.5,1.0,0.5,1.0};
    ipc_ H_col[] = {1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,11};
    ipc_ H_ptr[] = {1,2,4,6,8,10,12,14,16,18,20,22};
    rpc_ g[] = {0.5,-0.5,-1.0,-1.0,-1.0, -1.0,-1.0,-1.0,-1.0,-1.0,-0.5};

    //  describe constraints

    ipc_ A_ne = 30;
    rpc_ A_val[] = {1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
                        1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
                        1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
    ipc_ A_col[] = {1,2,3,4,5,6,7,8,9,10,11,3,4,5,6,7,8,9,10,11,
                   2,3,4,5,6,7,8,9,10,11 };
    ipc_ A_ptr[] = {1,12,21,31};
    rpc_ c_l[] = {10.0,9.0,-INFINITY};
    rpc_ c_u[] = {10.0,INFINITY,10.0};
    rpc_ x_l[] = {0.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0};
    rpc_ x_u[] = {INFINITY,INFINITY,INFINITY,INFINITY,INFINITY,INFINITY,
                      INFINITY,INFINITY,INFINITY,INFINITY,INFINITY};

    // provide optimal variables, Lagrange multipliers and dual variables
    rpc_ x[] = {0.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.0,1.0};
    rpc_ c[] = {10.0,9.0,10.0};
    rpc_ y[] = { -1.0,1.5,-2.0};
    rpc_ z[] = {2.0,4.0,2.5,2.5,2.5,2.5,2.5,2.5,2.5,2.5,2.5};
    // provide interior-point constraint and variable status
    ipc_ c_stat[] = {-1,-1,1};
    ipc_ x_stat[] = {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};

    // Set output storage

    ipc_ status;

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

    // Initialize CRO
    cro_initialize( &data, &control, &status );

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

    // crossover the solution
    cro_crossover_solution( &data, &control, &inform,
                            n, m, m_equal,
                            H_ne, H_val, H_col, H_ptr,
                            A_ne, A_val, A_col, A_ptr,
                            g, c_l, c_u, x_l, x_u, x, c, y, z,
                            x_stat, c_stat );

    printf(" CRO_crossover exit status = %1" i_ipc_ "\n", inform.status);

    // Delete internal workspace
    cro_terminate( &data, &control, &inform );

}