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.

parametric real type T#

Below, the symbol T refers to a parametric real type that may be Float32 (single precision) or Float64 (double precision). Calable functions as described are with T as Float64, but variants (with the additional suffix _s, e.g., scu_initialize_s) are available with T as Float32.

callable functions#

    function scu_information(data, inform, status)

Provides output information

Parameters:

data

holds private internal data

inform

is a structure containing output information (see scu_inform_type)

status

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

  • 0

    The values were recorded successfully

 
    function scu_terminate(data, control, inform)

Deallocate all internal private storage

Parameters:

data

holds private internal data

control

is a structure containing control information (see scu_control_type)

inform

is a structure containing output information (see scu_inform_type)

available structures#

scu_control_type structure#

    struct scu_control_type
      f_indexing::Bool

detailed documentation#

control derived type as a Julia structure

components#

Bool f_indexing

use C or Fortran sparse matrix indexing

scu_inform_type structure#

    struct scu_inform_type
      status::Int32
      alloc_status::Int32
      inertia::NTuple{3,Cint}

detailed documentation#

inform derived type as a Julia structure

components#

Int32 status

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

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

Int32 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