BSC#

purpose#

The bsc package takes given matrices $A$ and (diagonal) $D$ , and builds the Schur complement $S = A D A^{T}$ in sparse co-ordinate (and optionally sparse column) format(s). Full advantage is taken of any zero coefficients in the matrix $A$ .

See Section 4 of $GALAHAD/doc/bsc.pdf for a brief description of the method employed and other details.

matrix storage#

The unsymmetric $m$ by $n$ matrix $A$ may be presented and stored in a variety of convenient input formats.

Dense storage format: The matrix $A$ is stored as a compact dense matrix by rows, that is, the values of the entries of each row in turn are stored in order within an appropriate real one-dimensional array. In this case, component $n * i + j$ of the storage array A_val will hold the value $A_{i j}$ for $0 \leq i \leq m - 1$ , $0 \leq j \leq n - 1$ . The string A_type = ‘dense’ should be specified.

Sparse co-ordinate storage format: Only the nonzero entries of the matrices are stored. For the $l$ -th entry, $0 \leq l \leq n e - 1$ , of $A$ , its row index i, column index j and value $A_{i j}$ , $0 \leq i \leq m - 1$ , $0 \leq j \leq n - 1$ , are stored as the $l$ -th components of the integer arrays A_row and A_col and real array A_val, respectively, while the number of nonzeros is recorded as A_ne = $n e$ . The string A_type = ‘coordinate’should be specified.

Sparse row-wise storage format: Again only the nonzero entries are stored, but this time 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_{i j}$ of the nonzero entries in the i-th row are stored in components l = A_ptr(i), $\dots$ , A_ptr(i+1)-1, $0 \leq i \leq m - 1,$ of the integer array A_col, and real array A_val, respectively. For sparse matrices, this scheme almost always requires less storage than its predecessor. The string A_type = ‘sparse_by_rows’ should be specified.

Sparse column-wise storage format: Once again only the nonzero entries are stored, but this time they are ordered so that those in column j appear directly before those in column j+1. For the j-th column of $A$ the j-th component of the integer array A_ptr holds the position of the first entry in this column, while A_ptr(n) holds the total number of entries. The row indices i, $0 \leq i \leq m - 1$ , and values $A_{i j}$ of the nonzero entries in the j-th columns are stored in components l = A_ptr(j), $\dots$ , A_ptr(j+1)-1, $0 \leq j \leq n - 1$ , of the integer array A_row, and real array A_val, respectively. As before, for sparse matrices, this scheme almost always requires less storage than the co-ordinate format. The string A_type = ‘sparse_by_columns’ should be specified.

The symmetric $n$ by $n$ Schur complement matrix $S$ may be returned in a couple of formats. But now crucially symmetry is exploited by only storing values from the lower triangular part (i.e, those entries that lie on or below the leading diagonal).

Sparse co-ordinate storage format: Only the nonzero entries of the matrices are stored. For the $l$ -th entry, $0 \leq l \leq n e - 1$ , of $S$ , its row index i, column index j and value $S_{i j}$ , $0 \leq j \leq i \leq n - 1$ , are stored as the $l$ -th components of the integer arrays S_row and S_col and real array S_val, respectively, while the number of nonzeros is recorded as S_ne = $n e$ . Note that only the entries in the lower triangle will be returned.

Sparse row-wise storage format: Again only the nonzero entries are stored, but this time they are ordered so that those in row i appear directly before those in row i+1. For the i-th row of $S$ the i-th component of the integer array S_ptr holds the position of the first entry in this row, while S_ptr(n) holds the total number of entries. The column indices j, $0 \leq j \leq i$ , and values $S_{i j}$ of the entries in the i-th row are stored in components l = S_ptr(i), …, S_ptr(i+1)-1, $0 \leq i \leq n - 1,$ of the integer array S_col, and real array S_val, respectively. Note that as before only the entries in the lower triangle will be stored. For sparse matrices, this scheme almost always requires less storage than its predecessor.

functions#

bsc.initialize()#

Set default option values and initialize private data.

Returns:

optionsdict

dictionary containing default control options:

errorint
error and warning diagnostics occur on stream error.

outint
general output occurs on stream out.

print_levelint
the level of output required is specified by print_level.

max_colint
maximum permitted number of nonzeros in a column of \f $A \f$ ; -ve means unlimit.

new_aint
how much has $A$ changed since it was last accessed:

0

unchanged.

1

values changed.

2

structure changed.

3

structure changed but values not required.

extra_space_sint
how much extra space is to be allocated in $S$ above that needed to hold the Schur complement.

s_also_by_columnbool
should s.ptr also be set to indicate the first entry in each column of $S$ .

space_criticalbool
if space_critical is True, every effort will be made to use as little space as possible. This may result in longer computation time.

deallocate_error_fatalbool
if deallocate_error_fatal is True, any array/pointer deallocation error will terminate execution. Otherwise, computation will continue.

prefixstr
all output lines will be prefixed by the string contained in quotes within prefix, e.g. ‘word’ (note the qutoes) will result in the prefix word.

bsc.load(m, n, A_type, A_ne, A_row, A_col, A_ptr, options=None)#

Import the structure of $A$ to build that of $S$ .

Parameters:

mint
holds the number of rows of $A$ .

nint
holds the number of columns of $A$ .

A_typestring
specifies the unsymmetric storage scheme used for the matrix $A$ . It should be one of ‘coordinate’, ‘sparse_by_rows’ or ‘dense’; lower or upper case variants are allowed.

A_neint
holds the number of entries in $A$ in the sparse co-ordinate storage scheme. It need not be set for any of the other two schemes.

A_rowndarray(A_ne)
holds the row indices of $A$ in the sparse co-ordinate storage scheme. It need not be set for any of the other two schemes, and in this case can be None.

A_colndarray(A_ne)
holds the column indices of $A$ in either the sparse co-ordinate, or the sparse row-wise storage scheme. It need not be set when the dense storage scheme is used, and in this case can be None.

A_ptrndarray(m+1)
holds the starting position of each row of $A$ , as well as the total number of entries, in the sparse row-wise storage scheme. It need not be set when the other schemes are used, and in this case can be None.

optionsdict, optional
dictionary of control options (see bsc.initialize).

Returns:

S_neint
holds the number of entries in $S$ .

bsc.form(m, n, A_ne, A_val, S_ne, D)#

Form the Schur complement matrix $S = A D A^{T}$ .

Parameters:

mint
holds the number of rows of $A$ .

nint
holds the number of columns of $A$ .

A_neint
holds the number of entries in the matrix $A$ .

A_valndarray(A_ne)
holds the values of the nonzeros in the matrix $A$ in the same order as specified in the sparsity pattern in bsc.load.

S_neint
holds the number of entries in the matrix $S$ , as returned by bsc.load.

Dndarray(n)
holds the values of diagonal matrix $D$ . If $D$ is the identity matrix, it can take the value None to save storage.

Returns:

S_rowndarray(S_ne)
holds the row indices of $S$ in the sparse co-ordinate storage scheme.

S_colndarray(S_ne)
holds the column indices of $S$ in either the sparse co-ordinate, or the sparse row-wise storage scheme.

S_ptrndarray(n+1)
holds the starting position of each row of $S$ , as well as the total number of entries, in the sparse row-wise storage scheme.

S_valndarray(S_ne)
holds the values of the nonzeros in the matrix $S$ .

[optional] bsc.information()

Provide optional output information.

Returns:

informdict

dictionary containing output information:

statusint
the return status. Possible values are:

0

The call was successful.

-1

An allocation error occurred. A message indicating the offending array is written on unit options[‘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 options[‘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 restriction n > 0 or m > 0 or requirement that type contains its relevant string ‘dense’, ‘coordinate’ or ‘sparse_by_rows’ has been violated.

alloc_statusint
the status of the last attempted allocation/deallocation.

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

max_col_aint
the maximum number of entries in a column of $A$ .

exceeds_max_colint
the number of columns of $A$ that have more than control.max_col entries.

timefloat
the total CPU time spent in the package.

clock_timefloat
the total clock time spent in the package.

bsc.finalize()#

Deallocate all internal private storage.

example code#

from galahad import bsc
import numpy as np
np.set_printoptions(precision=2,suppress=True,floatmode='fixed')
print("\n** python test: bsc")

# set parameters
n = 4
m = 3

# describe matrix
A_type = 'coordinate'
A_ne = 6
A_row = np.array([0,0,1,1,2,2])
A_col = np.array([0,1,2,3,0,3])
A_ptr = None
A_val = np.array([1.0,1.0,1.0,1.0,1.0,1.0])
D = np.array([1.0,2.0,3.0,4.0])

# allocate internal data and set default options
options = bsc.initialize()

# set some non-default options
options['print_level'] = 0
#print("options:", options)

# load data (and optionally non-default options)
S_ne = bsc.load(m, n, A_type, A_ne, A_row, A_col, A_ptr, options)
print(" S_ne:",S_ne)

# form S = A D A'
S_row, S_col, S_ptr, S_val = bsc.form(m, n, A_ne, A_val, S_ne, D)
print(" S_row:",S_row)
print(" S_col:",S_col)
print(" S_ptr:",S_ptr)
print(" S_val:",S_val)

# get information
inform = bsc.information()
print('** bsc exit status:', inform['status'])

# deallocate internal data

bsc.terminate()

This example code is available in $GALAHAD/src/bsc/Python/test_bsc.py .