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+1) holds the total number of entries plus one. The column indices j, $1 \leq j \leq n$, and values $A_{ij}$ of the nonzero entries in the i-th row are stored in components l = A_ptr(i), $\ldots$, A_ptr(i+1)-1, $1 \leq i \leq m$, of the integer array A_col, and real array A_val, respectively.

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.

parametric real type T#

Below, the symbol T refers to a parametric real type that may be Float32 (single precision) or Float64 (double precision).

callable functions#

    function fdc_initialize(T, data, control, status)

Set default control values and initialize private data

Parameters:

data

holds private internal data

control

is a structure containing control information (see fdc_control_type)

status

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

0
The import was successful.

    function fdc_read_specfile(T, control, specfile)

Read the content of a specification file, and assign values associated with given keywords to the corresponding control parameters. An in-depth discussion of specification files is available, and a detailed list of keywords with associated default values is provided in $GALAHAD/src/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 structure containing control information (see fdc_control_type)
specfile	is a one-dimensional array of type Vararg{Cchar} that must give the name of the specification file

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

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

Parameters:

control	is a structure containing control information (see fdc_control_type)
data	holds private internal data
inform	is a structure containing output information (see fdc_inform_type)
status	is a scalar variable of type Int32 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 Int32 that holds the number of rows of $A$.
n	is a scalar variable of type Int32 that holds the number of columns of $A$.
A_ne	is a scalar variable of type Int32 that holds the number of nonzero entries in $A$.
A_col	is a one-dimensional array of size A_ne and type Int32 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 Int32 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 T 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 T that holds the linear term $b$ in the constraints. The i-th component of `b`, i = 1, … , m, contains $b_i$.
n_depen	is a scalar variable of type Int32 that holds the number of dependent constraints, if any.
depen	is a one-dimensional array of size m and type Int32 whose first n_depen components contain the indices of dependent constraints.

    function fdc_terminate(T, data, control, inform)

Deallocate all internal private storage

Parameters:

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

available structures#

fdc_control_type structure#

    struct fdc_control_type{T}
      f_indexing::Bool
      error::Int32
      out::Int32
      print_level::Int32
      indmin::Int32
      valmin::Int32
      pivot_tol::T
      zero_pivot::T
      max_infeas::T
      use_sls::Bool
      scale::Bool
      space_critical::Bool
      deallocate_error_fatal::Bool
      symmetric_linear_solver::NTuple{31,Cchar}
      unsymmetric_linear_solver::NTuple{31,Cchar}
      prefix::NTuple{31,Cchar}
      sls_control::sls_control_type{T}
      uls_control::uls_control_type{T}

detailed documentation#

control derived type as a Julia structure

components#

Bool f_indexing

use C or Fortran sparse matrix indexing

Int32 error

unit for error messages

Int32 out

unit for monitor output

Int32 print_level

controls level of diagnostic output

Int32 indmin

initial estimate of integer workspace for sls (obsolete)

Int32 valmin

initial estimate of real workspace for sls (obsolete)

T pivot_tol

the relative pivot tolerance (obsolete)

T zero_pivot

the absolute pivot tolerance used (obsolete)

T 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]

symmetric (indefinite) linear equation solver

char unsymmetric_linear_solver[31]

unsymmetric linear equation solver

NTuple{31,Cchar} prefix

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#

    struct fdc_time_type{T}
      total::T
      analyse::T
      factorize::T
      clock_total::T
      clock_analyse::T
      clock_factorize::T

detailed documentation#

time derived type as a Julia structure

components#

T total

the total CPU time spent in the package

T analyse

the CPU time spent analysing the required matrices prior to factorization

T factorize

the CPU time spent factorizing the required matrices

T clock_total

the total clock time spent in the package

T clock_analyse

the clock time spent analysing the required matrices prior to factorization

T clock_factorize

the clock time spent factorizing the required matrices

fdc_inform_type structure#

    struct fdc_inform_type{T}
      status::Int32
      alloc_status::Int32
      bad_alloc::NTuple{81,Cchar}
      factorization_status::Int32
      factorization_integer::Int64
      factorization_real::Int64
      non_negligible_pivot::T
      time::fdc_time_type{T}
      sls_inform::sls_inform_type{T}
      uls_inform::uls_inform_type{T}

detailed documentation#

inform derived type as a Julia structure

components#

Int32 status

return status. See FDC_find_dependent for details

Int32 alloc_status

the status of the last attempted allocation/deallocation

NTuple{81,Cchar} bad_alloc

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

Int32 factorization_status

the return status from the factorization

Int64 factorization_integer

the total integer workspace required for the factorization

Int64 factorization_real

the total real workspace required for the factorization

T 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/Julia/test_fdc.jl . A variety of supported Hessian and constraint matrix storage formats are shown.

# test_fdc.jl
# Simple code to test the Julia interface to FDC

using GALAHAD
using Test
using Printf
using Accessors

function test_fdc(::Type{T}) where T
  # Derived types
  data = Ref{Ptr{Cvoid}}()
  control = Ref{fdc_control_type{T}}()
  inform = Ref{fdc_inform_type{T}}()

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

  # Set output storage
  depen = zeros(Cint, m) # dependencies, if any
  n_depen = Ref{Cint}()
  status = Ref{Cint}()

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

  # Initialize FDC
  fdc_initialize(T, data, control, status)

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

  # Start from 0
  fdc_find_dependent_rows(T, 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\n", status[])
    else
      @printf("FDC_find_dependent - dependent rows(s):")
      for i in 1:n_depen[]
        @printf(" %i", depen[i])
      end
      @printf(", status = %i\n", status[])
    end
  else
    @printf("FDC_find_dependent - exit status = %1i\n", status[])
  end

  # Delete internal workspace
  fdc_terminate(T, data, control, inform)
  return 0
end

@testset "FDC" begin
  @test test_fdc(Float32) == 0
  @test test_fdc(Float64) == 0
end