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 SINGLE.

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