GALAHAD FDC package#

purpose#

Given an under-determined set of linear equations/constraints $a_i^T x = b_i^{}$, $i = 1, \ldots, m$ involving $n \geq m$ unknowns $x$, the fdc package determines whether the constraints are consistent, and if so how many of the constraints are dependent; a list of dependent constraints, that is, those which may be removed without changing the solution set, will be found and the remaining $a_i$ will be linearly independent. Full advantage is taken of any zero coefficients in the matrix $A$ whose columns are the vectors $a_i^T$.

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

method#

A choice of two methods is available. In the first, the matrix

\[\begin{split}K = \begin{pmatrix}\alpha I & A^T \\ A & 0 \end{pmatrix}\end{split}\]

is formed and factorized for some small $\alpha > 0$ using the SLS package — the factors $K = P L D L^T P^T$ are used to determine whether $A$ has dependent rows. In particular, in exact arithmetic dependencies in $A$ will correspond to zero pivots in the block diagonal matrix $D$.

The second choice of method finds factors $A = P L U Q$ of the rectangular matrix $A$ using the ULS package. In this case, dependencies in $A$ will be reflected in zero diagonal entries in $U$ in exact arithmetic.

The factorization in either case may also be used to determine whether the system is consistent.

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.

introduction to function calls#

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

fdc_initialize - provide default control parameters and set up initial data structures
fdc_read_specfile (optional) - override control values by reading replacement values from a file
fdc_find_dependent_rows - find the number of dependent rows and, if there are any, whether the constraints are independent
fdc_terminate - deallocate data structures

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 fdc_control_type;
struct fdc_inform_type;
struct fdc_time_type;

// global functions

void fdc_initialize(void **data, struct fdc_control_type* control, ipc_ *status);
void fdc_read_specfile(struct fdc_control_type* control, const char specfile[]);

void fdc_find_dependent_rows(
    struct fdc_control_type* control,
    void **data,
    struct fdc_inform_type* inform,
    ipc_ *status,
    ipc_ m,
    ipc_ n,
    ipc_ A_ne,
    const ipc_ A_col[],
    const ipc_ A_ptr[],
    const rpc_ A_val[],
    const rpc_ b[],
    ipc_ *n_depen,
    ipc_ depen[]
);

void fdc_terminate(
    void **data,
    struct fdc_control_type* control,
    struct fdc_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 fdc_initialize(void **data, struct fdc_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 fdc_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 fdc_read_specfile(struct fdc_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/fdc/FDC.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/fdc.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 fdc_control_type)
specfile	is a character string containing the name of the specification file

void fdc_find_dependent_rows(
    struct fdc_control_type* control,
    void **data,
    struct fdc_inform_type* inform,
    ipc_ *status,
    ipc_ m,
    ipc_ n,
    ipc_ A_ne,
    const ipc_ A_col[],
    const ipc_ A_ptr[],
    const rpc_ A_val[],
    const rpc_ b[],
    ipc_ *n_depen,
    ipc_ depen[]
)

Find dependent rows and, if any, check if $A x = b$ is consistent

Parameters:

control	is a struct containing control information (see fdc_control_type)
data	holds private internal data
inform	is a struct containing output information (see fdc_inform_type)
status	is a scalar variable of type ipc_, that gives the entry and exit status from the package. Possible exit values are: 0 The run was successful. -1 An allocation error occurred. A message indicating the offending array is written on unit control.error, and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -2 A deallocation error occurred. A message indicating the offending array is written on unit control.error and the returned allocation status and a string containing the name of the offending array are held in inform.alloc_status and inform.bad_alloc respectively. -3 The restrictions n > 0 and m > 0 or requirement that a type contains its relevant string ‘dense’, ‘coordinate’, ‘sparse_by_rows’, ‘diagonal’, ‘scaled_identity’, ‘identity’, ‘zero’ or ‘none’ has been violated. -5 The constraints appear to be inconsistent. -9 The analysis phase of the factorization failed; the return status from the factorization package is given in the component inform.factor_status -10 The factorization failed; the return status from the factorization package is given in the component inform.factor_status.
m	is a scalar variable of type ipc_, that holds the number of rows of $A$.
n	is a scalar variable of type ipc_, that holds the number of columns of $A$.
A_ne	is a scalar variable of type ipc_, that holds the number of nonzero entries in $A$.
A_col	is a one-dimensional array of size A_ne and type ipc_, that holds the column indices of $A$ in a row-wise storage scheme. The nonzeros must be ordered so that those in row i appear directly before those in row i+1, the order within each row is unimportant.
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.
A_val	is a one-dimensional array of size a_ne and type rpc_, that holds the values of the entries of the $A$ ordered as in A_col and A_ptr.
b	is a one-dimensional array of size m and type rpc_, that holds the linear term $b$ in the constraints. The i-th component of b, i = 0, … , m-1, contains $b_i$.
n_depen	is a scalar variable of type ipc_, that holds the number of dependent constraints, if any.
depen	is a one-dimensional array of size m and type ipc_, whose first n_depen components contain the indices of dependent constraints.

void fdc_terminate(
    void **data,
    struct fdc_control_type* control,
    struct fdc_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a struct containing control information (see fdc_control_type)
inform	is a struct containing output information (see fdc_inform_type)

available structures#

fdc_control_type structure#

#include <galahad_fdc.h>

struct fdc_control_type {
    // fields

    bool f_indexing;
    ipc_ error;
    ipc_ out;
    ipc_ print_level;
    ipc_ indmin;
    ipc_ valmin;
    rpc_ pivot_tol;
    rpc_ zero_pivot;
    rpc_ max_infeas;
    bool use_sls;
    bool scale;
    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 uls_control_type uls_control;
};

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

unit for monitor output

ipc_ print_level

controls level of diagnostic output

ipc_ indmin

initial estimate of integer workspace for sls (obsolete)

ipc_ valmin

initial estimate of real workspace for sls (obsolete)

rpc_ pivot_tol

the relative pivot tolerance (obsolete)

rpc_ zero_pivot

the absolute pivot tolerance used (obsolete)

rpc_ max_infeas

the largest permitted residual

bool use_sls

choose whether SLS or ULS is used to determine dependencies

bool scale

should the rows of A be scaled to have unit infinity norm or should no scaling be applied

bool space_critical

if space is critical, ensure allocated arrays are no bigger than needed

bool deallocate_error_fatal

exit if any deallocation fails

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 uls_control_type uls_control

control parameters for ULS

fdc_time_type structure#

#include <galahad_fdc.h>

struct fdc_time_type {
    // fields

    rpc_ total;
    rpc_ analyse;
    rpc_ factorize;
    rpc_ clock_total;
    rpc_ clock_analyse;
    rpc_ clock_factorize;
};

detailed documentation#

time derived type as a C struct

components#

rpc_ total

the total CPU time spent in the package

rpc_ analyse

the CPU time spent analysing the required matrices prior to factorization

rpc_ factorize

the CPU time spent factorizing the required matrices

rpc_ clock_total

the total clock time spent in the package

rpc_ clock_analyse

the clock time spent analysing the required matrices prior to factorization

rpc_ clock_factorize

the clock time spent factorizing the required matrices

fdc_inform_type structure#

#include <galahad_fdc.h>

struct fdc_inform_type {
    // fields

    ipc_ status;
    ipc_ alloc_status;
    char bad_alloc[81];
    ipc_ factorization_status;
    int64_t factorization_integer;
    int64_t factorization_real;
    rpc_ non_negligible_pivot;
    struct fdc_time_type time;
    struct sls_inform_type sls_inform;
    struct uls_inform_type uls_inform;
};

detailed documentation#

inform derived type as a C struct

components#

ipc_ status

return status. See FDC_find_dependent 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_ factorization_status

the return status from the factorization

int64_t factorization_integer

the total integer workspace required for the factorization

int64_t factorization_real

the total real workspace required for the factorization

rpc_ non_negligible_pivot

the smallest pivot which was not judged to be zero when detecting linear dependent constraints

struct fdc_time_type time

timings (see above)

struct sls_inform_type sls_inform

SLS inform type.

struct uls_inform_type uls_inform

ULS inform type.

Index

example calls#

This is an example of how to use the package to find a subset of independent linear constraints; the code is available in $GALAHAD/src/fdc/C/fdct.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.

/* fdct.c */
/* Full test for the FDC C interface using C sparse matrix indexing */

#include <stdio.h>
#include <math.h>
#include <string.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_fdc.h"

int main(void) {

    // Derived types
    void *data;
    struct fdc_control_type control;
    struct fdc_inform_type inform;

    // Set problem data
    ipc_ m = 3; // number of rows
    ipc_ n = 4; // number of columns
    ipc_ A_ne = 10; // number of nonzeros
    ipc_ A_col[] = {0, 1, 2, 3, 0, 1, 2, 3, 1, 3}; // column indices
    ipc_ A_ptr[] = {0, 4, 8, 10}; // row pointers
    rpc_ A_val[] = {1.0, 2.0, 3.0, 4.0, 2.0, -4.0, 6.0, -8.0, 5.0, 10.0};
    rpc_ b[] = {5.0, 10.0, 0.0};

    // Set output storage
    ipc_ depen[m]; // dependencies, if any
    ipc_ n_depen;
    ipc_ status;

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

    // Initialize FDC
    fdc_initialize( &data, &control, &status );

    // Set user-defined control options
    control.f_indexing = false; // C sparse matrix indexing
    control.use_sls = true;
    strcpy(control.symmetric_linear_solver, "sytr ");

    // Start from 0

    fdc_find_dependent_rows( &control, &data, &inform, &status, m, n, A_ne,
                             A_col, A_ptr, A_val, b, &n_depen, depen );

    if(status == 0){
      if(n_depen == 0){
        printf("FDC_find_dependent - no dependent rows, status = %1" i_ipc_ "\n",
               status);
      }else{
        printf("FDC_find_dependent - dependent rows(s):" );
        for( ipc_ i = 0; i < n_depen; i++) printf(" %" i_ipc_ "", depen[i]);
        printf(", status = %" i_ipc_ "\n", status);
      }
    }else{
        printf("FDC_find_dependent - exit status = %1" i_ipc_ "\n", status);
    }

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

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

/* fdctf.c */
/* Full test for the FDC C interface using Fortran sparse matrix indexing */

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

int main(void) {

    // Derived types
    void *data;
    struct fdc_control_type control;
    struct fdc_inform_type inform;

    // Set problem data
    ipc_ m = 3; // number of rows
    ipc_ n = 4; // number of columns
    ipc_ A_ne = 10; // number of nonzeros
    ipc_ A_col[] = {1, 2, 3, 4, 1, 2, 3, 4, 2, 4}; // column indices
    ipc_ A_ptr[] = {1, 5, 9, 11}; // row pointers
    rpc_ A_val[] = {1.0, 2.0, 3.0, 4.0, 2.0, -4.0, 6.0, -8.0, 5.0, 10.0};
    rpc_ b[] = {5.0, 10.0, 0.0};

    // Set output storage
    ipc_ depen[m]; // dependencies, if any
    ipc_ n_depen;
    ipc_ status;

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

    // Initialize FDC
    fdc_initialize( &data, &control, &status );

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

    // Start from 0

    fdc_find_dependent_rows( &control, &data, &inform, &status, m, n, A_ne,
                             A_col, A_ptr, A_val, b, &n_depen, depen );

    if(status == 0){
      if(n_depen == 0){
        printf("FDC_find_dependent - no dependent rows, status = %" i_ipc_ "\n",
               status);
      }else{
        printf("FDC_find_dependent - dependent rows(s):" );
        for( ipc_ i = 0; i < n_depen; i++) printf(" %" i_ipc_ "", depen[i]);
        printf(", status = %" i_ipc_ "\n", status);
      }
    }else{
        printf("FDC_find_dependent - exit status = %1" i_ipc_ "\n", status);
    }

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