GALAHAD SCU package#

purpose#

The scu package computes the solution to an extended system of $n + m$ sparse real linear equations in $n + m$ unknowns,

\[\begin{split}\begin{pmatrix}A & B \\ C & D\end{pmatrix} \begin{pmatrix}x_1 \\ x_2\end{pmatrix} = \begin{pmatrix}b_1 \\ b_2\end{pmatrix}\end{split}\]

in the case where the $n$ by $n$ matrix $A$ is nonsingular and solutions to the systems

\[A x = b \;\;\mbox{and}\;\; A^T y = c\]

may be obtained from an external source, such as an existing factorization. The subroutine uses reverse communication to obtain the solution to such smaller systems. The method makes use of the Schur complement matrix

\[S = D - C A^{-1} B.\]

The Schur complement is stored and factorized as a dense matrix and the subroutine is thus appropriate only if there is sufficient storage for this matrix. Special advantage is taken of symmetry and definiteness in the coefficient matrices. Provision is made for introducing additional rows and columns to, and removing existing rows and columns from, the extended.

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/scu.pdf for additional details.

method#

The function scu_factorize forms the Schur complement $S$ of $A$ in the extended matrix by repeated reverse communication to obtain the columns of $A^{-1} B$. The Schur complement or its negative is then factorized into its QR or, if possible, Cholesky factors.

The function scu_solve solves the extended system using the following well-known scheme:

(i) Compute the solution to $A u = b_1$;
(ii) Compute $x_2$ from $S x_2 = b_2 - C u$;
(iii) Compute the solution to $A v = B x_2$; and
(iv) Compute $x_1 = u - v$.

The functions scu__append and scu_delete compute the factorization of the Schur complement after a row and column have been appended to, and removed from, the extended matrix, respectively. The existing factorization is updated to obtain the new one; this is normally more efficient than forming the factorization from scratch.

introduction to function calls#

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

scu_initialize - provide default control parameters and set up initial data structures
scu_read_specfile (optional) - override control values by reading replacement values from a file
scu_form_and_factorize - form and factorize the Schur-complement matrix $S$
scu_solve_system - solve the block system (1)
scu_add_rows_and_cols (optional) - update the factors of the Schur-complement matrix when rows and columns are added to (1).
scu_delete_rows_and_cols (optional) - update the factors of the Schur-complement matrix when rows and columns are removed from (1).
scu_information (optional) - recover information about the solution and solution process
scu_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 scu_control_type;
struct scu_inform_type;

// global functions

void scu_information(void **data, struct scu_inform_type* inform, ipc_ *status);

void scu_terminate(
    void **data,
    struct scu_control_type* control,
    struct scu_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 and structure names#

The function and structure names described below are appropriate for the default real working precision (double) and integer word length (int32_t). To use the functions and structures with different precisions and integer word lengths, an additional suffix must be added to their names (and the arguments set accordingly). The appropriate suffices are:

_s for single precision (float) reals and standard 32-bit (int32_t) integers;

_q for quadruple precision (__real128) reals (if supported) and standard 32-bit (int32_t) integers;

_64 for standard precision (double) reals and 64-bit (int64_t) integers;

_s_64 for single precision (float) reals and 64-bit (int64_t) integers; and

_q_64 for quadruple precision (__real128) reals (if supported) and 64-bit (int64_t) integers.

Thus a call to scu_information below will instead be

void scu_information_s_64(void **data, struct scu_inform_type_s_64* inform,
                          int64_t *status)

if single precision (float) reals and 64-bit (int64_t) integers are required. Thus it is possible to call functions for this package with more that one precision and/or integer word length at same time. An example is provided for the package expo, and the obvious modifications apply equally here.

function calls#

void scu_information(void **data, struct scu_inform_type* inform,
                     ipc_ *status)

Provides output information

Parameters:

data

holds private internal data

inform

is a struct containing output information (see scu_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 scu_terminate(void **data,
                   struct scu_control_type* control,
                   struct scu_inform_type* inform)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a struct containing control information (see scu_control_type)
inform	is a struct containing output information (see scu_inform_type)

available structures#

scu_control_type structure#

#include <galahad_scu.h>

struct scu_control_type {
    // fields

    bool f_indexing;
};

detailed documentation#

control derived type as a C struct

components#

bool f_indexing

use C or Fortran sparse matrix indexing

scu_inform_type structure#

#include <galahad_scu.h>

struct scu_inform_type {
    // fields

    ipc_ status;
    ipc_ alloc_status;
    ipc_ inertia[3];
};

detailed documentation#

inform derived type as a C struct

components#

ipc_ status

return status. A non-zero value indicates an error or a request for further information. See SCU_solve for details.

ipc_ alloc_status

the return status from the last attempted internal workspace array allocation or deallocation. A non-zero value indicates that the allocation or deallocation was unsuccessful, and corresponds to the fortran STAT= value on the user’s system.

ipc_ inertia[3]

the inertia of $S$ when the extended matrix is symmetric. Specifically, inertia(i), i=0,1,2 give the number of positive, negative and zero eigenvalues of $S$ respectively.

Index