GALAHAD ULS package#

purpose#

The uls package solves dense or sparse unsymmetric systems of linear equations using variants of Gaussian elimination. Given a sparse matrix $A = \{ a_{ij} \}_{m \times n}$, and an $n$-vector $b$, this function solves the systems $A x = b$ or $A^T x = b$. Both square ($m=n$) and rectangular ($m\neq n$) matrices are handled; one of an infinite class of solutions for consistent systems will be returned whenever $A$ is not of full rank.

The method provides a common interface to a variety of well-known solvers from HSL and elsewhere. Currently supported solvers include MA28/GLS and HSL_MA48 from {HSL}, as well as GETR from LAPACK. Note that, with the exception of he Netlib reference LAPACK code, the solvers themselves do not form part of this package and must be obtained/linked to separately. Dummy instances are provided for solvers that are unavailable. Also note that additional flexibility may be obtained by calling the solvers directly rather that via this package.

terminology#

The solvers used each produce an $L U$ factorization of $A$, where $L$ and U are permuted lower and upper triangular matrices (respectively). It is convenient to write this factorization in the form

\[A = P_R L U P_C,\]

where $P_R$ and $P_C$ are row and column permutation matrices, respectively.

supported solvers#

The key features of the external solvers supported by uls are given in the following table:

External solver characteristics#
solver	factorization	out-of-core	parallelised
`GLS/MA28`	sparse	no	no
`HSL_MA48`	sparse	no	no
`GETR`	dense	no	with parallel LAPACK

method#

Variants of sparse Gaussian elimination are used. See Section 4 of $GALAHAD/doc/uls.pdf for a brief description of the method employed and other details.

The solver GLS is available as part of GALAHAD and relies on the HSL Archive package MA33. To obtain HSL Archive packages, see

http://hsl.rl.ac.uk/archive/ .

The solver HSL_MA48 is part of HSL 2011. To obtain HSL 2011 packages, see

http://hsl.rl.ac.uk .

The solver GETR is available as S/DGETRF/S as part of LAPACK. Reference versions are provided by GALAHAD, but for good performance machined-tuned versions should be used.

The methods used are described in the user-documentation for

HSL 2011, A collection of Fortran codes for large-scale scientific computation (2011).

http://www.hsl.rl.ac.uk .

matrix storage#

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

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

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

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

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

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

introduction to function calls#

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

uls_initialize - provide default control parameters and set up initial data structures
uls_read_specfile (optional) - override control values by reading replacement values from a file
uls_factorize_matrix - set up matrix data structures, analyse the structure to choose a suitable order for factorization, and then factorize the matrix $A$
uls_reset_control (optional) - possibly change control parameters if a sequence of problems are being solved
uls_solve_system - solve the linear system of equations $Ax=b$ or $A^Tx=b$
uls_information (optional) - recover information about the solution and solution process
uls_terminate - deallocate data structures

See the examples section for illustrations of use.

callable functions#

overview of functions provided#

// typedefs

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

// structs

struct uls_control_type;
struct uls_inform_type;

// global functions

void uls_initialize(
    const char solver[],
    void **data,
    struct uls_control_type* control,
    ipc_ *status
);

void uls_read_specfile(struct uls_control_type* control, const char specfile[]);

void uls_factorize_matrix(
    struct uls_control_type* control,
    void **data,
    ipc_ *status,
    ipc_ m,
    ipc_ n,
    const char type[],
    ipc_ ne,
    const rpc_ val[],
    const ipc_ row[],
    const ipc_ col[],
    const ipc_ ptr[]
);

void uls_reset_control(
    struct uls_control_type* control,
    void **data,
    ipc_ *status
);

void uls_solve_system(
    void **data,
    ipc_ *status,
    ipc_ m,
    ipc_ n,
    rpc_ sol[],
    bool trans
);

void uls_information(void **data, struct uls_inform_type* inform, ipc_ *status);

void uls_terminate(
    void **data,
    struct uls_control_type* control,
    struct uls_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 uls_initialize(
    const char solver[],
    void **data,
    struct uls_control_type* control,
    ipc_ *status
)

Set default control values and initialize private data

Select solver, set default control values and initialize private data

Parameters:

solver	is a one-dimensional array of type char that specifies the solver package that should be used to factorize the matrix $A$. It should be one of ‘gls’, ‘ma28’, ‘ma48 or ‘getr’; lower or upper case variants are allowed. Only ‘getr’ is available by default, but others are easily installed (see README.external).
data	holds private internal data
control	is a struct containing control information (see uls_control_type)
status	is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are: 0 The initialization was successful. -26 The requested solver is not available**

void uls_read_specfile(struct uls_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/uls/ULS.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/uls.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 uls_control_type)
specfile	is a character string containing the name of the specification file

void uls_factorize_matrix(
    struct uls_control_type* control,
    void **data,
    ipc_ *status,
    ipc_ m,
    ipc_ n,
    const char type[],
    ipc_ ne,
    const rpc_ val[],
    const ipc_ row[],
    const ipc_ col[],
    const ipc_ ptr[]
)

Import matrix data into internal storage prior to solution, analyse the sparsity patern, and subsequently factorize the matrix

Parameters:

control	is a struct whose members provide control paramters for the remaining prcedures (see uls_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, analysis and factorization were conducted successfully. -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 and m> 0 or requirement that the matrix type must contain the relevant string ‘dense’, ‘coordinate’ or ‘sparse_by_rows has been violated. -26 The requested solver is not available. -29 This option is not available with this solver. -32 More than control.max integer factor size words of internal integer storage are required for in-core factorization. -50 A solver-specific error occurred; check the solver-specific information component of inform along with the solver’s documentation for more details.
m	is a scalar variable of type ipc_, that holds the number of rows in the unsymmetric matrix $A$.
n	is a scalar variable of type ipc_, that holds the number of columns in the unsymmetric matrix $A$.
type	is a one-dimensional array of type char that specifies the unsymmetric storage scheme used for the matrix $A$. It should be one of ‘coordinate’, ‘sparse_by_rows’ or ‘dense’; lower or upper case variants are allowed.
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.
val	is a one-dimensional array of size ne and type rpc_, that holds the values of the entries of the matrix $A$ in any of the supported storage schemes.
row	is a one-dimensional array of size ne and type ipc_, that holds the row indices of the matrix $A$ 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.
col	is a one-dimensional array of size ne and type ipc_, that holds the column indices of the matrix $A$ in either the sparse co-ordinate, or the sparse row-wise storage scheme. It need not be set when the dense storage schemes is used, and in this case can be NULL.
ptr	is a one-dimensional array of size m+1 and type ipc_, that holds the starting position of each row of the matrix $A$, as well as the total number of entries, in the sparse row-wise storage scheme. It need not be set when the other schemes are used, and in this case can be NULL.

void uls_reset_control(
    struct uls_control_type* control,
    void **data,
    ipc_ *status
)

Reset control parameters after import if required.

Parameters:

control

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

void uls_solve_system(
    void **data,
    ipc_ *status,
    ipc_ m,
    ipc_ n,
    rpc_ sol[],
    bool trans
)

Solve the linear system $Ax=b$ or $A^Tx=b$.

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 required solution was obtained. -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. -34 The package PARDISO failed; check the solver-specific information components inform.pardiso iparm and inform.pardiso_dparm along with PARDISO’s documentation for more details. -35 The package WSMP failed; check the solver-specific information components inform.wsmp_iparm and inform.wsmp dparm along with WSMP’s documentation for more details.
m	is a scalar variable of type ipc_, that holds the number of rows in the unsymmetric matrix $A$.
n	is a scalar variable of type ipc_, that holds the number of columns in the unsymmetric matrix $A$.
sol	is a one-dimensional array of size n and type double. On entry, it must hold the vector $b$. On a successful exit, its contains the solution $x$.
trans	is a scalar variable of type bool, that specifies whether to solve the equation $A^Tx=b$ (trans=true) or $Ax=b$ (trans=false).

void uls_information(void **data, struct uls_inform_type* inform, ipc_ *status)

Provides output information

Parameters:

data

holds private internal data

inform

is a struct containing output information (see uls_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 uls_terminate(
    void **data,
    struct uls_control_type* control,
    struct uls_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a struct containing control information (see uls_control_type)
inform	is a struct containing output information (see uls_inform_type)

available structures#

uls_control_type structure#

#include <galahad_uls.h>

struct uls_control_type {
    // fields

    bool f_indexing;
    ipc_ error;
    ipc_ warning;
    ipc_ out;
    ipc_ print_level;
    ipc_ print_level_solver;
    ipc_ initial_fill_in_factor;
    ipc_ min_real_factor_size;
    ipc_ min_integer_factor_size;
    int64_t max_factor_size;
    ipc_ blas_block_size_factorize;
    ipc_ blas_block_size_solve;
    ipc_ pivot_control;
    ipc_ pivot_search_limit;
    ipc_ minimum_size_for_btf;
    ipc_ max_iterative_refinements;
    bool stop_if_singular;
    rpc_ array_increase_factor;
    rpc_ switch_to_full_code_density;
    rpc_ array_decrease_factor;
    rpc_ relative_pivot_tolerance;
    rpc_ absolute_pivot_tolerance;
    rpc_ zero_tolerance;
    rpc_ acceptable_residual_relative;
    rpc_ acceptable_residual_absolute;
    char prefix[31];
};

detailed documentation#

control derived type as a C struct

components#

bool f_indexing

use C or Fortran sparse matrix indexing

ipc_ error

unit for error messages

ipc_ warning

unit for warning messages

ipc_ out

unit for monitor output

ipc_ print_level

controls level of diagnostic output

ipc_ print_level_solver

controls level of diagnostic output from external solver

ipc_ initial_fill_in_factor

prediction of factor by which the fill-in will exceed the initial number of nonzeros in $A$

ipc_ min_real_factor_size

initial size for real array for the factors and other data

ipc_ min_integer_factor_size

initial size for integer array for the factors and other data

int64_t max_factor_size

maximum size for real array for the factors and other data

ipc_ blas_block_size_factorize

level 3 blocking in factorize

ipc_ blas_block_size_solve

level 2 and 3 blocking in solve

ipc_ pivot_control

pivot control:

1 Threshold Partial Pivoting is desired
2 Threshold Rook Pivoting is desired
3 Threshold Complete Pivoting is desired
4 Threshold Symmetric Pivoting is desired
5 Threshold Diagonal Pivoting is desired

ipc_ pivot_search_limit

number of rows/columns pivot selection restricted to (0 = no restriction)

ipc_ minimum_size_for_btf

the minimum permitted size of blocks within the block-triangular form

ipc_ max_iterative_refinements

maximum number of iterative refinements allowed

bool stop_if_singular

stop if the matrix is found to be structurally singular

rpc_ array_increase_factor

factor by which arrays sizes are to be increased if they are too small

rpc_ switch_to_full_code_density

switch to full code when the density exceeds this factor

rpc_ array_decrease_factor

if previously allocated internal workspace arrays are greater than array_decrease_factor times the currently required sizes, they are reset to current requirements

rpc_ relative_pivot_tolerance

pivot threshold

rpc_ absolute_pivot_tolerance

any pivot small than this is considered zero

rpc_ zero_tolerance

any entry smaller than this in modulus is reset to zero

rpc_ acceptable_residual_relative

refinement will cease as soon as the residual $\|Ax-b\|$ falls below max( acceptable_residual_relative \* $\|b\|$, acceptable_residual_absolute )

rpc_ acceptable_residual_absolute

see acceptable_residual_relative

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’

uls_inform_type structure#

#include <galahad_uls.h>

struct uls_inform_type {
    // fields

    ipc_ status;
    ipc_ alloc_status;
    char bad_alloc[81];
    ipc_ more_info;
    int64_t out_of_range;
    int64_t duplicates;
    int64_t entries_dropped;
    int64_t workspace_factors;
    ipc_ compresses;
    int64_t entries_in_factors;
    ipc_ rank;
    ipc_ structural_rank;
    ipc_ pivot_control;
    ipc_ iterative_refinements;
    bool alternative;
    char solver[21];
    struct gls_ainfo gls_ainfo;
    struct gls_finfo gls_finfo;
    struct gls_sinfo gls_sinfo;
    struct ma48_ainfo ma48_ainfo;
    struct ma48_finfo ma48_finfo;
    struct ma48_sinfo ma48_sinfo;
    ipc_ lapack_error;
};

detailed documentation#

inform derived type as a C struct

components#

ipc_ status

reported return status:

0

success
-1

allocation error
-2

deallocation error
-3

matrix data faulty (m < 1, n < 1, ne < 0)
-26

unknown solver
-29

unavailable option
-31

input order is not a permutation or is faulty in some other way
-32

error with integer workspace
-33

error with real workspace
-50

solver-specific error; see the solver’s info parameter

ipc_ alloc_status

STAT value after allocate failure.

char bad_alloc[81]

name of array which provoked an allocate failure

ipc_ more_info

further information on failure

int64_t out_of_range

number of indices out-of-range

int64_t duplicates

number of duplicates

int64_t entries_dropped

number of entries dropped during the factorization

int64_t workspace_factors

predicted or actual number of reals and integers to hold factors

ipc_ compresses

number of compresses of data required

int64_t entries_in_factors

number of entries in factors

ipc_ rank

estimated rank of the matrix

ipc_ structural_rank

structural rank of the matrix

ipc_ pivot_control

pivot control:

1

Threshold Partial Pivoting has been used
2

Threshold Rook Pivoting has been used
3

Threshold Complete Pivoting has been desired
4

Threshold Symmetric Pivoting has been desired
5

Threshold Diagonal Pivoting has been desired

ipc_ iterative_refinements

number of iterative refinements performed

bool alternative

has an “alternative” y: A^T y = 0 and yT b > 0 been found when trying to solve A x = b ?

char solver[21]

name of external solver used to factorize and solve

struct gls_ainfo gls_ainfo

the output arrays from GLS

struct gls_finfo gls_finfo

see gls_ainfo

struct gls_sinfo gls_sinfo

see gls_ainfo

struct ma48_ainfo ma48_ainfo

the output arrays from HSL’s MA48

struct ma48_finfo ma48_finfo

see ma48_ainfo

struct ma48_sinfo ma48_sinfo

see ma48_ainfo

ipc_ lapack_error

the LAPACK error return code

Index

example calls#

This is an example of how to use the package to solve a linear system; the code is available in $GALAHAD/src/uls/C/ulst.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.

/* ulst.c */
/* Full test for the ULS C interface using C sparse matrix indexing */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <float.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_uls.h"

ipc_ maxabsarray(rpc_ a[], ipc_ n, rpc_ *maxabs);

int main(void) {

    // Derived types
    void *data;
    struct uls_control_type control;
    struct uls_inform_type inform;

    // Set problem data
    ipc_ m = 5; // column dimension of A
    ipc_ n = 5; // column dimension of A
    ipc_ ne = 7; // number of entries of A
    ipc_ dense_ne = 25; // number of elements of A as a dense matrix
    ipc_ row[] = {0, 1, 1, 2, 2, 3, 4}; // row indices
    ipc_ col[] = {0, 0, 4, 1, 2, 2, 3}; // column indices
    ipc_ ptr[] = {0, 1, 3, 5, 6, 7}; // pointers to indices
    rpc_ val[] = {2.0, 3.0, 6.0, 4.0, 1.0, 5.0, 1.0}; // values
    rpc_ dense[] = {2.0, 0.0, 0.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 6.0,
                      0.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0,
                      0.0, 0.0, 0.0, 1.0, 0.0};
    rpc_ rhs[] = {2.0, 33.0, 11.0, 15.0, 4.0};
    rpc_ rhst[] = {8.0, 12.0, 23.0, 5.0, 12.0};
    rpc_ sol[] = {1.0, 2.0, 3.0, 4.0, 5.0};
    ipc_ i, status;
    rpc_ x[n];
    rpc_ error[n];
    _Bool trans;

    rpc_ norm_residual;
    rpc_ good_x = pow( DBL_EPSILON, 0.3333 );

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

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

    printf(" storage          RHS   refine   RHST  refine\n");

    for( ipc_ d=1; d <= 3; d++){
        // Initialize ULS - use the gls solver
        uls_initialize( "getr", &data, &control, &status );

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

        switch(d){ // import matrix data and factorize
            case 1: // sparse co-ordinate storage
                printf(" coordinate     ");
                uls_factorize_matrix( &control, &data, &status, m, n,
                                    "coordinate", ne, val, row, col, NULL );
                break;
            case 2: // sparse by rows
                printf(" sparse by rows ");
                uls_factorize_matrix( &control, &data, &status, m, n,
                                    "sparse_by_rows", ne, val, NULL, col, ptr );
                break;
            case 3: // dense
                printf(" dense          ");
                uls_factorize_matrix( &control, &data, &status, m, n, "dense",
                                      dense_ne, dense, NULL, NULL, NULL );
                break;
            }

        // Set right-hand side and solve the system A x = b
        for(i=0; i<n; i++) x[i] = rhs[i];
        trans = false;
        uls_solve_system( &data, &status, m, n, x, trans );
        uls_information( &data, &inform, &status );

        if(inform.status == 0){
          for(i=0; i<n; i++) error[i] = x[i]-sol[i];
          status = maxabsarray( error, n, &norm_residual );
          if(norm_residual < good_x){
            printf("   ok  ");
          }else{
            printf("  fail ");
          }
        }else{
            printf(" ULS_solve exit status = %1" i_ipc_ "\n", inform.status);
        }
        // printf("sol: ");
        // for( ipc_ i = 0; i < n; i++) printf("%f ", x[i]);

        // resolve, this time using iterative refinement
        control.max_iterative_refinements = 1;
        uls_reset_control( &control, &data, &status );
        for(i=0; i<n; i++) x[i] = rhs[i];
        uls_solve_system( &data, &status, m, n, x, trans );
        uls_information( &data, &inform, &status );

        if(inform.status == 0){
          for(i=0; i<n; i++) error[i] = x[i]-sol[i];
          status = maxabsarray( error, n, &norm_residual );
          if(norm_residual < good_x){
            printf("    ok  ");
          }else{
            printf("   fail ");
          }
        }else{
            printf(" ULS_solve exit status = %1" i_ipc_ "\n", inform.status);
        }

        // Set right-hand side and solve the system A^T x = b
        for(i=0; i<n; i++) x[i] = rhst[i];
        trans = true;
        uls_solve_system( &data, &status, m, n, x, trans );
        uls_information( &data, &inform, &status );

        if(inform.status == 0){
          for(i=0; i<n; i++) error[i] = x[i]-sol[i];
          status = maxabsarray( error, n, &norm_residual );
          if(norm_residual < good_x){
            printf("   ok  ");
          }else{
            printf("  fail ");
          }
        }else{
            printf(" ULS_solve exit status = %1" i_ipc_ "\n", inform.status);
        }
        // printf("sol: ");
        // for( ipc_ i = 0; i < n; i++) printf("%f ", x[i]);

        // resolve, this time using iterative refinement
        control.max_iterative_refinements = 1;
        uls_reset_control( &control, &data, &status );
        for(i=0; i<n; i++) x[i] = rhst[i];
        uls_solve_system( &data, &status, m, n, x, trans );
        uls_information( &data, &inform, &status );

        if(inform.status == 0){
          for(i=0; i<n; i++) error[i] = x[i]-sol[i];
          status = maxabsarray( error, n, &norm_residual );
          if(norm_residual < good_x){
            printf("    ok  ");
          }else{
            printf("   fail ");
          }
        }else{
            printf(" ULS_solve exit status = %1" i_ipc_ "\n", inform.status);
        }

        // Delete internal workspace
        uls_terminate( &data, &control, &inform );
        printf("\n");
    }
}

ipc_ maxabsarray(rpc_ a[], ipc_ n, rpc_ *maxabs)
 {
    ipc_ i;
    rpc_ b, max;
    max=abs(a[0]);
    for(i=1; i<n; i++)
    {
        b = fabs(a[i]);
    if(max<b)
          max=b;
    }
    *maxabs=max;
    return 0;
 }

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

/* ulstf.c */
/* Full test for the ULS C interface using Fortran sparse matrix indexing */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <float.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_uls.h"

ipc_ maxabsarray(rpc_ a[], ipc_ n, rpc_ *maxabs);

int main(void) {

    // Derived types
    void *data;
    struct uls_control_type control;
    struct uls_inform_type inform;

    // Set problem data
    ipc_ m = 5; // column dimension of A
    ipc_ n = 5; // column dimension of A
    ipc_ ne = 7; // number of entries of A
    ipc_ dense_ne = 25; // number of elements of A as a dense matrix
    ipc_ row[] = {1, 2, 2, 3, 3, 4, 5}; // row indices
    ipc_ col[] = {1, 1, 5, 2, 3, 3, 4}; // column indices
    ipc_ ptr[] = {1, 2, 4, 6, 7, 8}; // pointers to indices
    rpc_ val[] = {2.0, 3.0, 6.0, 4.0, 1.0, 5.0, 1.0}; // values
    rpc_ dense[] = {2.0, 0.0, 0.0, 0.0, 0.0, 3.0, 0.0, 0.0, 0.0, 6.0,
                      0.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0, 0.0,
                      0.0, 0.0, 0.0, 1.0, 0.0};
    rpc_ rhs[] = {2.0, 33.0, 11.0, 15.0, 4.0};
    rpc_ rhst[] = {8.0, 12.0, 23.0, 5.0, 12.0};
    rpc_ sol[] = {1.0, 2.0, 3.0, 4.0, 5.0};
    ipc_ i, status;
    rpc_ x[n];
    rpc_ error[n];
    _Bool trans;

    rpc_ norm_residual;
    rpc_ good_x = pow( DBL_EPSILON, 0.3333 );

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

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

    printf(" storage          RHS   refine   RHST  refine\n");

    for( ipc_ d=1; d <= 3; d++){
        // Initialize ULS - use the gls solver
        uls_initialize( "getr", &data, &control, &status );

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

        switch(d){ // import matrix data and factorize
            case 1: // sparse co-ordinate storage
                printf(" coordinate     ");
                uls_factorize_matrix( &control, &data, &status, m, n,
                                    "coordinate", ne, val, row, col, NULL );
                break;
            case 2: // sparse by rows
                printf(" sparse by rows ");
                uls_factorize_matrix( &control, &data, &status, m, n,
                                    "sparse_by_rows", ne, val, NULL, col, ptr );
                break;
            case 3: // dense
                printf(" dense          ");
                uls_factorize_matrix( &control, &data, &status, m, n,
                                     "dense", dense_ne, dense, NULL, NULL, NULL );
                break;
            }

        // Set right-hand side and solve the system A x = b
        for(i=0; i<n; i++) x[i] = rhs[i];
        trans = false;
        uls_solve_system( &data, &status, m, n, x, trans );
        uls_information( &data, &inform, &status );

        if(inform.status == 0){
          for(i=0; i<n; i++) error[i] = x[i]-sol[i];
          status = maxabsarray( error, n, &norm_residual );
          if(norm_residual < good_x){
            printf("   ok  ");
          }else{
            printf("  fail ");
          }
        }else{
            printf(" ULS_solve exit status = %1" i_ipc_ "\n", inform.status);
        }
        // printf("sol: ");
        // for( ipc_ i = 0; i < n; i++) printf("%f ", x[i]);

        // resolve, this time using iterative refinement
        control.max_iterative_refinements = 1;
        uls_reset_control( &control, &data, &status );
        for(i=0; i<n; i++) x[i] = rhs[i];
        uls_solve_system( &data, &status, m, n, x, trans );
        uls_information( &data, &inform, &status );

        if(inform.status == 0){
          for(i=0; i<n; i++) error[i] = x[i]-sol[i];
          status = maxabsarray( error, n, &norm_residual );
          if(norm_residual < good_x){
            printf("    ok  ");
          }else{
            printf("   fail ");
          }
        }else{
            printf(" ULS_solve exit status = %1" i_ipc_ "\n", inform.status);
        }

        // Set right-hand side and solve the system A^T x = b
        for(i=0; i<n; i++) x[i] = rhst[i];
        trans = true;
        uls_solve_system( &data, &status, m, n, x, trans );
        uls_information( &data, &inform, &status );

        if(inform.status == 0){
          for(i=0; i<n; i++) error[i] = x[i]-sol[i];
          status = maxabsarray( error, n, &norm_residual );
          if(norm_residual < good_x){
            printf("   ok  ");
          }else{
            printf("  fail ");
          }
        }else{
            printf(" ULS_solve exit status = %1" i_ipc_ "\n", inform.status);
        }
        // printf("sol: ");
        // for( ipc_ i = 0; i < n; i++) printf("%f ", x[i]);

        // resolve, this time using iterative refinement
        control.max_iterative_refinements = 1;
        uls_reset_control( &control, &data, &status );
        for(i=0; i<n; i++) x[i] = rhst[i];
        uls_solve_system( &data, &status, m, n, x, trans );
        uls_information( &data, &inform, &status );

        if(inform.status == 0){
          for(i=0; i<n; i++) error[i] = x[i]-sol[i];
          status = maxabsarray( error, n, &norm_residual );
          if(norm_residual < good_x){
            printf("    ok  ");
          }else{
            printf("   fail ");
          }
        }else{
            printf(" ULS_solve exit status = %1" i_ipc_ "\n", inform.status);
        }

        // Delete internal workspace
        uls_terminate( &data, &control, &inform );
        printf("\n");
    }
}

ipc_ maxabsarray(rpc_ a[], ipc_ n, rpc_ *maxabs)
 {
    ipc_ i;
    rpc_ b, max;
    max=abs(a[0]);
    for(i=1; i<n; i++)
    {
        b = fabs(a[i]);
    if(max<b)
          max=b;
    }
    *maxabs=max;
    return 0;
 }