GALAHAD PSLS package#
purpose#
Given a sparse symmetric matrix \(A = \{ a_{ij} \}_{n \times n}\),
the psls package builds a suitable symmetric, positive definite—or
diagonally dominant—preconditioner \(P\) of \(A\) or a symmetric
sub-matrix thereof. The matrix \(A\) need not be definite. Facilities
are provided to apply the preconditioner to a given vector, and to
remove rows and columns (symmetrically) from the initial preconditioner
without a full re-factorization.
method#
The basic preconditioners are described in detail in
A. R. Conn, N. I. M. Gould and Ph. L. Toint, LANCELOT. A fortran package for large-scale nonlinear optimization (release A). Springer Verlag Series in Computational Mathematics 17, Berlin (1992), Section 3.3.10,
along with the more modern versions implements in ICFS due to
C.-J. Lin and J. J. More’, ``Incomplete Cholesky factorizations with limited memory’’. SIAM Journal on Scientific Computing 21 (1999) 21-45,
and in HSL_MI28 described by
J. A. Scott and M. Tuma, ``HSL_MI28: an efficient and robust limited-memory incomplete Cholesky factorization code’’. ACM Transactions on Mathematical Software 40(4) (2014), Article 24.
The factorization methods used by the package sls in conjunction
with some preconditioners are described in the documentation to that
package. The key features of the external solvers supported by sls
are given in the following table.
solver |
factorization |
indefinite \(A\) |
out-of-core |
parallelised |
|---|---|---|---|---|
|
multifrontal |
yes |
no |
no |
|
multifrontal |
yes |
no |
no |
|
multifrontal |
yes |
yes |
OpenMP core |
|
left-looking |
yes |
no |
OpenMP fully |
|
left-looking |
no |
no |
OpenMP fully |
|
multifrontal |
yes |
no |
OpenMP core |
|
multifrontal |
yes |
no |
CUDA core |
|
multifrontal |
yes |
optionally |
MPI |
|
left-right-looking |
yes |
no |
OpenMP fully |
|
left-right-looking |
yes |
optionally |
OpenMP fully |
|
left-right-looking |
yes |
no |
OpenMP fully |
|
left-right-looking |
yes |
no |
OpenMP fully |
|
dense |
no |
no |
with parallel LAPACK |
|
dense |
yes |
no |
with parallel LAPACK |
|
dense band |
no |
no |
with parallel LAPACK |
Note that, with the exception of SSIDS and the Netlib
reference LAPACK codes,
the solvers themselves do not form part of this package and
must be obtained/linked to separately. See the documentation for sls
for more details and references.
Dummy instances are provided for solvers that are unavailable.
Also note that additional flexibility may be obtained by calling the
solvers directly rather that via this package.
Orderings to reduce the bandwidth, as implemented in HSL’s MC61, are due to
J. K. Reid and J. A. Scott, ``Ordering symmetric sparse matrices for small profile and wavefront’’. International Journal for Numerical Methods in Engineering 45 (1999) 1737-1755.
If a subset of the rows and columns are specified, the remaining rows/columns
are removed before processing. Any subsequent removal of rows and columns
is achieved using the Schur-complement updating package scu
unless a complete re-factorization is likely more efficient.
matrix storage#
The symmetric \(n\) by \(n\) matrix \(A\) may be presented and stored in a variety of formats. But 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).
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. Since \(A\) is symmetric, only the lower triangular part (that is the part \(A_{ij}\) for \(0 \leq j \leq i \leq n-1\)) need be held. In this case the lower triangle should be stored by rows, that is component \(i * i / 2 + j\) of the storage array A_val will hold the value \(A_{ij}\) (and, by symmetry, \(A_{ji}\)) for \(0 \leq j \leq i \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 ne-1\), of \(A\), its row index i, column index j and value \(A_{ij}\), \(0 \leq j \leq i \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 = \(ne\). Note that only the entries in the lower triangle should be stored. 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(n) holds the total number of entries. The column indices j, \(0 \leq j \leq i\), and values \(A_{ij}\) of the entries in the i-th row are stored in components l = A_ptr(i), …, A_ptr(i+1)-1 of the integer array A_col, and real array A_val, respectively. Note that as before only the entries in the lower triangle should be stored. For sparse matrices, this scheme almost always requires less storage than its predecessor. The string A_type = ‘sparse_by_rows’ should be specified.
Diagonal storage format: If \(A\) is diagonal (i.e., \(A_{ij} = 0\) for all \(0 \leq i \neq j \leq n-1\)) only the diagonals entries \(A_{ii}\), \(0 \leq i \leq n-1\) need be stored, and the first n components of the array A_val may be used for the purpose. The string A_type = ‘diagonal’ should be specified.
Multiples of the identity storage format: If \(A\) is a multiple of the identity matrix, (i.e., \(H = \alpha I\) where \(I\) is the n by n identity matrix and \(\alpha\) is a scalar), it suffices to store \(\alpha\) as the first component of A_val. The string A_type = ‘scaled_identity’ should be specified.
The identity matrix format: If \(A\) is the identity matrix, no values need be stored. The string A_type = ‘identity’ should be specified.
The zero matrix format: The same is true if \(A\) is the zero matrix, but now the string A_type = ‘zero’ or ‘none’ should be specified.
introduction to function calls#
To solve a given problem, functions from the psls package must be called in the following order:
psls_initialize - provide default control parameters and set up initial data structures
psls_read_specfile (optional) - override control values by reading replacement values from a file
psls_import - set up matrix data structures for \(A\) prior to solution
psls_reset_control (optional) - possibly change control parameters if a sequence of problems are being solved
one of
psls_form_preconditioner - form and factorize a preconditioner \(P\) of the matrix \(A\)
psls_form_subset_preconditioner - form and factorize a preconditioner \(P\) of a symmetric submatrix of the matrix \(A\)
psls_update_preconditioner (optional) - update the preconditioner \(P\) when rows (amd columns) are removed
psls_apply_preconditioner - solve the linear system of equations \(Px=b\)
psls_information (optional) - recover information about the preconditioner and solution process
psls_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 psls_control_type; struct psls_inform_type; struct psls_time_type; // global functions void psls_initialize( void **data, struct psls_control_type* control, ipc_ *status ); void psls_read_specfile( struct psls_control_type* control, const char specfile[] ); void psls_import( struct psls_control_type* control, void **data, ipc_ *status, ipc_ n, const char type[], ipc_ ne, const ipc_ row[], const ipc_ col[], const ipc_ ptr[] ); void psls_reset_control( struct psls_control_type* control, void **data, ipc_ *status ); void psls_form_preconditioner( void **data, ipc_ *status, ipc_ ne, const rpc_ val[] ); void psls_form_subset_preconditioner( void **data, ipc_ *status, ipc_ ne, const rpc_ val[], ipc_ n_sub, const ipc_ sub[] ); void psls_update_preconditioner( void **data, ipc_ *status, ipc_ ne, const rpc_ val[], ipc_ n_del, const ipc_ del[] ); void psls_apply_preconditioner(void **data, ipc_ *status, ipc_ n, rpc_ sol[]); void psls_information(void **data, struct psls_inform_type* inform, ipc_ *status); void psls_terminate( void **data, struct psls_control_type* control, struct psls_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 calls#
void psls_initialize( void **data, struct psls_control_type* control, ipc_ *status )
Set default control values and initialize private data
Parameters:
data |
holds private internal data |
control |
is a struct containing control information (see psls_control_type) |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are (currently):
|
void psls_read_specfile( struct psls_control_type* control, const char 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/psls/PSLS.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/psls.pdf for a list of how these keywords relate to the components of the control structure.
Parameters:
control |
is a struct containing control information (see psls_control_type) |
specfile |
is a character string containing the name of the specification file |
void psls_import( struct psls_control_type* control, void **data, ipc_ *status, ipc_ n, const char type[], ipc_ ne, const ipc_ row[], const ipc_ col[], const ipc_ ptr[] )
Import structural matrix data into internal storage prior to solution.
Parameters:
control |
is a struct whose members provide control paramters for the remaining prcedures (see psls_control_type) |
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are:
|
n |
is a scalar variable of type ipc_, that holds the number of rows in the symmetric matrix \(A\). |
type |
is a one-dimensional array of type char that specifies the symmetric 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. |
ne |
is a scalar variable of type ipc_, that holds the number of entries in the lower triangular part of \(A\) in the sparse co-ordinate storage scheme. It need not be set for any of the other schemes. |
row |
is a one-dimensional array of size ne and type ipc_, that holds the row indices of the lower triangular part of \(A\) in the sparse co-ordinate storage scheme. It need not be set for any of the other three schemes, and in this case can be NULL. |
col |
is a one-dimensional array of size ne and type ipc_, that holds the column indices of the lower triangular part 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 NULL. |
ptr |
is a one-dimensional array of size n+1 and type ipc_, that holds the starting position of each row of the lower triangular part 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 NULL. |
void psls_reset_control( struct psls_control_type* control, void **data, ipc_ *status )
Reset control parameters after import if required.
Parameters:
control |
is a struct whose members provide control paramters for the remaining prcedures (see psls_control_type) |
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are:
|
void psls_form_preconditioner( void **data, ipc_ *status, ipc_ ne, const rpc_ val[] )
Form and factorize a preconditioner \(P\) of the matrix \(A\).
Parameters:
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are:
|
ne |
is a scalar variable of type ipc_, that holds the number of entries in the lower triangular part of the symmetric matrix \(A\). |
val |
is a one-dimensional array of size ne and type rpc_, that holds the values of the entries of the lower triangular part of the symmetric matrix \(A\) in any of the supported storage schemes. |
void psls_form_subset_preconditioner( void **data, ipc_ *status, ipc_ ne, const rpc_ val[], ipc_ n_sub, const ipc_ sub[] )
Form and factorize a \(P\) preconditioner of a symmetric submatrix of the matrix \(A\).
Parameters:
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are:
|
ne |
is a scalar variable of type ipc_, that holds the number of entries in the lower triangular part of the symmetric matrix \(A\). |
val |
is a one-dimensional array of size ne and type rpc_, that holds the values of the entries of the lower triangular part of the symmetric matrix \(A\) in any of the supported storage schemes. |
n_sub |
is a scalar variable of type ipc_, that holds the number of rows (and columns) of the required submatrix of \(A\). |
sub |
is a one-dimensional array of size n_sub and type ipc_, that holds the indices of the rows of required submatrix. |
void psls_update_preconditioner( void **data, ipc_ *status, ipc_ ne, const rpc_ val[], ipc_ n_del, const ipc_ del[] )
Update the preconditioner \(P\) when rows (amd columns) are removed.
Parameters:
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are:
|
ne |
is a scalar variable of type ipc_, that holds the number of entries in the lower triangular part of the symmetric matrix \(A\). |
val |
is a one-dimensional array of size ne and type rpc_, that holds the values of the entries of the lower triangular part of the symmetric matrix \(A\) in any of the supported storage schemes. |
n_del |
is a scalar variable of type ipc_, that holds the number of rows (and columns) of (sub) matrix that are to be deleted. |
del |
is a one-dimensional array of size n_fix and type ipc_, that holds the indices of the rows that are to be deleted. |
void psls_apply_preconditioner(void **data, ipc_ *status, ipc_ n, rpc_ sol[])
Solve the linear system \(Px=b\).
Parameters:
data |
holds private internal data |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are:
|
n |
is a scalar variable of type ipc_, that holds the number of entries in the vectors \(b\) and \(x\). |
sol |
is a one-dimensional array of size n and type double. On entry, it must hold the vector \(b\). On a successful exit, its contains the solution \(x\). Any component corresponding to rows/columns not in the initial subset recorded by psls_form_subset_preconditioner, or in those subsequently deleted by psls_update_preconditioner, will not be altered. |
void psls_information(void **data, struct psls_inform_type* inform, ipc_ *status)
Provide output information
Parameters:
data |
holds private internal data |
inform |
is a struct containing output information (see psls_inform_type) |
status |
is a scalar variable of type ipc_, that gives the exit status from the package. Possible values are (currently):
|
void psls_terminate( void **data, struct psls_control_type* control, struct psls_inform_type* inform )
Deallocate all internal private storage
Parameters:
data |
holds private internal data |
control |
is a struct containing control information (see psls_control_type) |
inform |
is a struct containing output information (see psls_inform_type) |
available structures#
psls_control_type structure#
#include <galahad_psls.h> struct psls_control_type { // fields bool f_indexing; ipc_ error; ipc_ out; ipc_ print_level; ipc_ preconditioner; ipc_ semi_bandwidth; ipc_ scaling; ipc_ ordering; ipc_ max_col; ipc_ icfs_vectors; ipc_ mi28_lsize; ipc_ mi28_rsize; rpc_ min_diagonal; bool new_structure; bool get_semi_bandwidth; bool get_norm_residual; bool space_critical; bool deallocate_error_fatal; char definite_linear_solver[31]; char prefix[31]; struct sls_control_type sls_control; struct mi28_control mi28_control; };
detailed documentation#
control derived type as a C struct
components#
bool f_indexing
use C or Fortran sparse matrix indexing
ipc_ error
unit for error messages
ipc_ out
unit for monitor output
ipc_ print_level
controls level of diagnostic output
ipc_ preconditioner
which preconditioner to use:
<0 no preconditioning occurs, \(P = I\)
0 the preconditioner is chosen automatically (forthcoming, and currently defaults to 1).
1 \(A\) is replaced by the diagonal, \(P\) = diag( max(\(A\), .min_diagonal ) ).
2 \(A\) is replaced by the band \(P\) = band(\(A\)) with semi-bandwidth .semi_bandwidth.
3 \(A\) is replaced by the reordered band \(P\) = band( order(\(A\)) ) with semi-bandwidth .semi_bandwidth, where order is chosen by the HSL package MC61 to move entries closer to the diagonal.
4 \(P\) is a full factorization of \(A\) using Schnabel-Eskow modifications, in which small or negative diagonals are made sensibly positive during the factorization.
5 \(P\) is a full factorization of \(A\) due to Gill, Murray, Ponceleon and Saunders, in which an indefinite factorization is altered to give a positive definite one.
6 \(P\) is an incomplete Cholesky factorization of \(A\) using the package ICFS due to Lin and More’.
7 \(P\) is an incomplete factorization of \(A\) implemented as HSL_MI28 from HSL.
8 \(P\) is an incomplete factorization of \(A\) due to Munskgaard (forthcoming).
>8 treated as 0.
N.B. Options 3-8 may require additional external software that is not part of the package, and that must be obtained separately.
ipc_ semi_bandwidth
the semi-bandwidth for band(H) when .preconditioner = 2,3
ipc_ scaling
not used at present
ipc_ ordering
see scaling
ipc_ max_col
maximum number of nonzeros in a column of \(A\) for Schur-complement factorization to accommodate newly deleted rpws and columns
ipc_ icfs_vectors
number of extra vectors of length n required by the Lin-More’ incomplete Cholesky preconditioner when .preconditioner = 6
ipc_ mi28_lsize
the maximum number of fill entries within each column of the incomplete factor L computed by HSL_MI28 when .preconditioner = 7. In general, increasing mi28_lsize improve the quality of the preconditioner but increases the time to compute and then apply the preconditioner. Values less than 0 are treated as 0
ipc_ mi28_rsize
the maximum number of entries within each column of the strictly lower triangular matrix \(R\) used in the computation of the preconditioner by HSL_MI28 when .preconditioner = 7. Rank-1 arrays of size .mi28_rsize \* n are allocated internally to hold \(R\). Thus the amount of memory used, as well as the amount of work involved in computing the preconditioner, depends on mi28_rsize. Setting mi28_rsize > 0 generally leads to a higher quality preconditioner than using mi28_rsize = 0, and choosing mi28_rsize >= mi28_lsize is generally recommended
rpc_ min_diagonal
the minimum permitted diagonal in diag(max(H,.min_diagonal))
bool new_structure
set new_structure true if the storage structure for the input matrix has changed, and false if only the values have changed
bool get_semi_bandwidth
set get_semi_bandwidth true if the semi-bandwidth of the submatrix is to be calculated
bool get_norm_residual
set get_norm_residual true if the residual when applying the preconditioner are to be calculated
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 definite_linear_solver[31]
the name of the definite linear equation solver used. Possible choices are currently: ‘sils’, ‘ma27’, ‘ma57’, ‘ma77’, ‘ma86’, ‘ma87’, ‘ma97’, ‘ssids’, ‘mumps’, ‘pardiso’, ‘mkl_pardiso’, ‘pastix’, ‘wsmp’, ‘potr’, ‘sytr’ and ‘pbtr’, although only ‘potr’, ‘sytr’, ‘pbtr’ and, for OMP 4.0-compliant compilers, ‘ssids’ are installed by default; others are easily installed (see README.external). More details of the capabilities of each solver are provided in the documentation for galahad_sls.
char prefix[31]
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 mi28_control mi28_control
control parameters for HSL_MI28
psls_time_type structure#
#include <galahad_psls.h> struct psls_time_type { // fields spc_ total; spc_ assemble; spc_ analyse; spc_ factorize; spc_ solve; spc_ update; rpc_ clock_total; rpc_ clock_assemble; rpc_ clock_analyse; rpc_ clock_factorize; rpc_ clock_solve; rpc_ clock_update; };
detailed documentation#
time derived type as a C struct
components#
spc_ total
total time
spc_ assemble
time to assemble the preconditioner prior to factorization
spc_ analyse
time for the analysis phase
spc_ factorize
time for the factorization phase
spc_ solve
time for the linear solution phase
spc_ update
time to update the factorization
rpc_ clock_total
total clock time spent in the package
rpc_ clock_assemble
clock time to assemble the preconditioner prior to factorization
rpc_ clock_analyse
clock time for the analysis phase
rpc_ clock_factorize
clock time for the factorization phase
rpc_ clock_solve
clock time for the linear solution phase
rpc_ clock_update
clock time to update the factorization
psls_inform_type structure#
#include <galahad_psls.h> struct psls_inform_type { // fields ipc_ status; ipc_ alloc_status; ipc_ analyse_status; ipc_ factorize_status; ipc_ solve_status; int64_t factorization_integer; int64_t factorization_real; ipc_ preconditioner; ipc_ semi_bandwidth; ipc_ reordered_semi_bandwidth; ipc_ out_of_range; ipc_ duplicates; ipc_ upper; ipc_ missing_diagonals; ipc_ semi_bandwidth_used; ipc_ neg1; ipc_ neg2; bool perturbed; rpc_ fill_in_ratio; rpc_ norm_residual; char bad_alloc[81]; ipc_ mc61_info[10]; rpc_ mc61_rinfo[15]; struct psls_time_type time; struct sls_inform_type sls_inform; struct mi28_info mi28_info; };
detailed documentation#
inform derived type as a C struct
components#
ipc_ status
reported return status:
0
success
-1
allocation error
-2
deallocation error
-3
matrix data faulty (.n < 1, .ne < 0)
ipc_ alloc_status
STAT value after allocate failure.
ipc_ analyse_status
status return from factorization
ipc_ factorize_status
status return from factorization
ipc_ solve_status
status return from solution phase
int64_t factorization_integer
number of integer words to hold factors
int64_t factorization_real
number of real words to hold factors
ipc_ preconditioner
code for the actual preconditioner used (see control.preconditioner)
ipc_ semi_bandwidth
the actual semi-bandwidth
ipc_ reordered_semi_bandwidth
the semi-bandwidth following reordering (if any)
ipc_ out_of_range
number of indices out-of-range
ipc_ duplicates
number of duplicates
ipc_ upper
number of entries from the strict upper triangle
ipc_ missing_diagonals
number of missing diagonal entries for an allegedly-definite matrix
ipc_ semi_bandwidth_used
the semi-bandwidth used
ipc_ neg1
number of 1 by 1 pivots in the factorization
ipc_ neg2
number of 2 by 2 pivots in the factorization
bool perturbed
has the preconditioner been perturbed during the fctorization?
rpc_ fill_in_ratio
ratio of fill in to original nonzeros
rpc_ norm_residual
the norm of the solution residual
char bad_alloc[81]
name of array which provoked an allocate failure
ipc_ mc61_info[10]
the integer and real output arrays from mc61
rpc_ mc61_rinfo[15]
see mc61_info
struct psls_time_type time
times for various stages
struct sls_inform_type sls_inform
inform values from SLS
struct mi28_info mi28_info
the output info structure from HSL’s mi28
example calls#
This is an example of how to use the package to solve a definite linear system; the code is available in $GALAHAD/src/psls/C/pslst.c . A variety of supported matrix storage formats are shown.
Notice that C-style indexing is used, and that this is flagged by setting
control.f_indexing to false. The floating-point type rpc_
is set in galahad_precision.h to double by default, but to float
if the preprocessor variable SINGLE is defined. Similarly, the integer
type ipc_ from galahad_precision.h is set to int by default,
but to int64_t if the preprocessor variable INTEGER_64 is defined.
/* pslst.c */
/* Full test for the PSLS C interface using C sparse matrix indexing */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_psls.h"
int main(void) {
// Derived types
void *data;
struct psls_control_type control;
struct psls_inform_type inform;
// Set problem data
ipc_ n = 5; // dimension of A
ipc_ ne = 7; // number of elements of A
ipc_ dense_ne = n * ( n + 1 ) / 2; // number of elements of dense A
ipc_ row[] = {0, 1, 1, 2, 2, 3, 4}; // A indices & values, NB lower triangle
ipc_ col[] = {0, 0, 4, 1, 2, 2, 4};
ipc_ ptr[] = {0, 1, 3, 5, 6, 7};
rpc_ val[] = {2.0, 3.0, 6.0, 4.0, 1.0, 5.0, 1.0};
rpc_ dense[] = {2.0, 3.0, 0.0, 0.0, 4.0, 1.0, 0.0,
0.0, 5.0, 0.0, 0.0, 6.0, 0.0, 0.0, 1.0};
char st = ' ';
ipc_ status;
ipc_ status_apply;
printf(" C sparse matrix indexing\n\n");
printf(" basic tests of storage formats\n\n");
for( ipc_ d=1; d <= 3; d++){
// Initialize PSLS
psls_initialize( &data, &control, &status );
control.preconditioner = 2; // band preconditioner
control.semi_bandwidth = 1; // semibandwidth
strcpy( control.definite_linear_solver, "sils" );
// Set user-defined control options
control.f_indexing = false; // C sparse matrix indexing
switch(d){
case 1: // sparse co-ordinate storage
st = 'C';
psls_import( &control, &data, &status, n,
"coordinate", ne, row, col, NULL );
psls_form_preconditioner( &data, &status, ne, val );
break;
printf(" case %1" i_ipc_ " break\n",d);
case 2: // sparse by rows
st = 'R';
psls_import( &control, &data, &status, n,
"sparse_by_rows", ne, NULL, col, ptr );
psls_form_preconditioner( &data, &status, ne, val );
break;
case 3: // dense
st = 'D';
psls_import( &control, &data, &status, n,
"dense", ne, NULL, NULL, NULL );
psls_form_preconditioner( &data, &status, dense_ne, dense );
break;
}
// Set right-hand side b in x
rpc_ x[] = {8.0, 45.0, 31.0, 15.0, 17.0}; // values
if(status == 0){
psls_information( &data, &inform, &status );
psls_apply_preconditioner( &data, &status_apply, n, x );
}else{
status_apply = - 1;
}
printf("%c storage: status from form & factorize = %" i_ipc_
" apply = %" i_ipc_ "\n", st, status, status_apply );
//printf("x: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", x[i]);
// Delete internal workspace
psls_terminate( &data, &control, &inform );
}
}This is the same example, but now fortran-style indexing is used; the code is available in $GALAHAD/src/psls/C/pslstf.c .
/* pslstf.c */
/* Full test for the PSLS C interface using Fortran sparse matrix indexing */
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_psls.h"
int main(void) {
// Derived types
void *data;
struct psls_control_type control;
struct psls_inform_type inform;
// Set problem data
ipc_ n = 5; // dimension of A
ipc_ ne = 7; // number of elements of A
ipc_ dense_ne = n * ( n + 1 ) / 2; // number of elements of dense A
ipc_ row[] = {1, 2, 2, 3, 3, 4, 5}; // A indices & values, NB lower triangle
ipc_ col[] = {1, 1, 5, 2, 3, 3, 5};
ipc_ ptr[] = {1, 2, 4, 6, 7, 8};
rpc_ val[] = {2.0, 3.0, 6.0, 4.0, 1.0, 5.0, 1.0};
rpc_ dense[] = {2.0, 3.0, 0.0, 0.0, 4.0, 1.0, 0.0,
0.0, 5.0, 0.0, 0.0, 6.0, 0.0, 0.0, 1.0};
char st = ' ';
ipc_ status;
ipc_ status_apply;
printf(" Fortran sparse matrix indexing\n\n");
printf(" basic tests of storage formats\n\n");
for( ipc_ d=1; d <= 3; d++){
// Initialize PSLS
psls_initialize( &data, &control, &status );
control.preconditioner = 2; // band preconditioner
control.semi_bandwidth = 1; // semibandwidth
strcpy( control.definite_linear_solver, "sils" );
// Set user-defined control options
control.f_indexing = true; // fortran sparse matrix indexing
switch(d){
case 1: // sparse co-ordinate storage
st = 'C';
psls_import( &control, &data, &status, n,
"coordinate", ne, row, col, NULL );
psls_form_preconditioner( &data, &status, ne, val );
break;
printf(" case %1" i_ipc_ " break\n",d);
case 2: // sparse by rows
st = 'R';
psls_import( &control, &data, &status, n,
"sparse_by_rows", ne, NULL, col, ptr );
psls_form_preconditioner( &data, &status, ne, val );
break;
case 3: // dense
st = 'D';
psls_import( &control, &data, &status, n,
"dense", ne, NULL, NULL, NULL );
psls_form_preconditioner( &data, &status, dense_ne, dense );
break;
}
// Set right-hand side b in x
rpc_ x[] = {8.0, 45.0, 31.0, 15.0, 17.0}; // values
if(status == 0){
psls_information( &data, &inform, &status );
psls_apply_preconditioner( &data, &status_apply, n, x );
}else{
status_apply = - 1;
}
printf("%c storage: status from form & factorize = %" i_ipc_
" apply = %" i_ipc_ "\n", st, status, status_apply );
//printf("x: ");
//for( ipc_ i = 0; i < n; i++) printf("%f ", x[i]);
// Delete internal workspace
psls_terminate( &data, &control, &inform );
}
}