GALAHAD DPS package#

purpose#

The dps package constructs a symmetric, positive definite matrix $M$ from a given $H$ so that $H$ is is diagonal in the norm $\|v\|_M = \sqrt{v^T M v}$ induced by $M$, and consequently minimizers of trust-region and regularized quadratic subproblems may be computed efficiently. The aim is either to minimize the quadratic objective function

\[q(x) = f + g^T x + \frac{1}{2} x^T H x,\]

where the vector $x$ is required to satisfy the ellipsoidal trust-region constraint $\|x\|_{M} \leq \Delta$, or to minimize the regularized quadratic objective

\[r(x) = q(x) + \frac{\sigma}{p} \|x\|_M^p,\]

where the radius $\Delta > 0$, the weight $\sigma > 0$, and the power $p \geq 2$. A factorization of the matrix $H$ will be required, so this package is most suited for the case where such a factorization, either dense or sparse, may be found efficiently.

See Section 4 of $GALAHAD/doc/dps.pdf for additional details.

method#

The required solution $x_*$ necessarily satisfies the optimality condition $H x_* + \lambda_* M x_* + g = 0$, where $\lambda_* \geq 0$ is a Lagrange multiplier that corresponds to the constraint $\|x\|_M \leq \Delta$ in the trust-region case, and is given by $\lambda_* = \sigma \|x_*\|^{p-2}$ for the regularization problem involve $r(x)$. In addition $H + \lambda_* M$ will be positive semi-definite; in most instances it will actually be positive definite, but in special “hard” cases singularity is a possibility.

The matrix $H$ is decomposed as

\[H = P L D L^T P^T\]

by calling the matrix factorization package SLS. Here $P$ is a permutation matrix, $L$ is unit lower triangular and $D$ is block diagonal, with blocks of dimension at most two. The spectral decomposition of each diagonal block of $D$ is computed, and each eigenvalue $\theta$ is replaced by $\max ( | \theta | , \theta_{\min} ) $, where $\theta_{\min}$ is a positive user-supplied value. The resulting block diagonal matrix is $B$, from which we define the modified-absolute-value

\[M = P L B L^T P^T;\]

an alternative due to Goldfarb uses instead the simpler

\[M = P L L^T P^T.\]

Given the factors of $H$ (and $M$), the required solution is found by making the change of variables $y = B^{1/2} L^T P^T x$ (or $y = L^T P^T x$ in the Goldfarb case) which results in ``diagonal’’ trust-region and regularization subproblems, whose solution may be easily obtained suing a Newton or higher-order iteration of a resulting “secular” equation. If subsequent problems, for which $H$ and $g$ are unchanged, are to be attempted, the existing factorization and solution may easily be exploited.

The dominant cost is that for the factorization of the symmetric, but potentially indefinite, matrix $H$ using the package SLS.

references#

The method is described in detail for the trust-region case in

N. I. M. Gould and J. Nocedal, ``The modified absolute-value factorization for trust-region minimization’’. In ``High Performance Algorithms and Software in Nonlinear Optimization’’ (R. De Leone, A. Murli, P. M. Pardalos and G. Toraldo, eds.), Kluwer Academic Publishers, (1998) 225–241,

while the adaptation for the regularization case is obvious. The method used to solve the diagonal trust-region and regularization subproblems are as given by

H. S. Dollar, N. I. M. Gould and D. P. Robinson, ``On solving trust-region and other regularised subproblems in optimization’’. Mathematical Programming Computation 2(1) (2010) 21–57

with simplifications due to the diagonal Hessian.

matrix storage#

The symmetric $n$ by $n$ matrix $H$ 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 $H$ 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 $H$ is symmetric, only the lower triangular part (that is the part $H_{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 H_val will hold the value $H_{ij}$ (and, by symmetry, $H_{ji}$) for $0 \leq j \leq i \leq n-1$. The string H_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 $H$, its row index i, column index j and value $H_{ij}$, $0 \leq j \leq i \leq n-1$, are stored as the $l$-th components of the integer arrays H_row and H_col and real array H_val, respectively, while the number of nonzeros is recorded as H_ne = $ne$. Note that only the entries in the lower triangle should be stored. The string H_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 $H$ the i-th component of the integer array H_ptr holds the position of the first entry in this row, while H_ptr(n) holds the total number of entries. The column indices j, $0 \leq j \leq i$, and values $H_{ij}$ of the entries in the i-th row are stored in components l = H_ptr(i), …, H_ptr(i+1)-1 of the integer array H_col, and real array H_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 H_type = ‘sparse_by_rows’ should be specified.

Diagonal storage format: If $H$ is diagonal (i.e., $H_{ij} = 0$ for all $0 \leq i \neq j \leq n-1$) only the diagonals entries $H_{ii}$, $0 \leq i \leq n-1$ need be stored, and the first n components of the array H_val may be used for the purpose. The string H_type = ‘diagonal’ should be specified.

Multiples of the identity storage format: If $H$ 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 H_val. The string H_type = ‘scaled_identity’ should be specified.

The identity matrix format: If $H$ is the identity matrix, no values need be stored. The string H_type = ‘identity’ should be specified.

The zero matrix format: The same is true if $H$ is the zero matrix, but now the string H_type = ‘zero’ or ‘none’ should be specified.

introduction to function calls#

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

dps_initialize - provide default control parameters and set up initial data structures
dps_read_specfile (optional) - override control values by reading replacement values from a file
dps_import - import control and matrix data structures
dps_reset_control (optional) - possibly change control parameters if a sequence of problems are being solved
one of
- dps_solve_tr_problem - solve the trust-region problem (1)
- dps_solve_rq_problem - solve the regularized-quadratic problem (2)
optionally one of
- dps_resolve_tr_problem - resolve the trust-region problem (1) when the non-matrix data has changed
- dps_resolve_rq_problem - resolve the regularized-quadratic problem (2) when the non-matrix data has changed
dps_information (optional) - recover information about the solution and solution process
dps_terminate - deallocate data structures

See the examples section for illustrations of use.

callable functions#

overview of functions provided#

// typedefs

typedef float spc_;
typedef double rpc_;
typedef int ipc_;

// structs

struct dps_control_type;
struct dps_inform_type;
struct dps_time_type;

// global functions

void dps_initialize(void **data, struct dps_control_type* control, ipc_ *status);
void dps_read_specfile(struct dps_control_type* control, const char specfile[]);

void dps_import(
    struct dps_control_type* control,
    void **data,
    ipc_ *status,
    ipc_ n,
    const char H_type[],
    ipc_ ne,
    const ipc_ H_row[],
    const ipc_ H_col[],
    const ipc_ H_ptr[]
);

void dps_reset_control(
    struct dps_control_type* control,
    void **data,
    ipc_ *status
);

void dps_solve_tr_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    ipc_ ne,
    rpc_ H_val[],
    rpc_ c[],
    rpc_ f,
    rpc_ radius,
    rpc_ x[]
);

void dps_solve_rq_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    ipc_ ne,
    rpc_ H_val[],
    rpc_ c[],
    rpc_ f,
    rpc_ power,
    rpc_ weight,
    rpc_ x[]
);

void dps_resolve_tr_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    rpc_ c[],
    rpc_ f,
    rpc_ radius,
    rpc_ x[]
);

void dps_resolve_rq_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    rpc_ c[],
    rpc_ f,
    rpc_ power,
    rpc_ weight,
    rpc_ x[]
);

void dps_information(void **data, struct dps_inform_type* inform, ipc_ *status);

void dps_terminate(
    void **data,
    struct dps_control_type* control,
    struct dps_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 dps_initialize(void **data, struct dps_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 dps_control_type)

status

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

0
The initialization was successful.

void dps_read_specfile(struct dps_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/dps/DPS.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/dps.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 dps_control_type)
specfile	is a character string containing the name of the specification file

void dps_import(
    struct dps_control_type* control,
    void **data,
    ipc_ *status,
    ipc_ n,
    const char H_type[],
    ipc_ ne,
    const ipc_ H_row[],
    const ipc_ H_col[],
    const ipc_ H_ptr[]
)

Import problem data into internal storage prior to solution.

Parameters:

control	is a struct whose members provide control paramters for the remaining prcedures (see dps_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: 1 The import was successful, and the package is ready for the solve phase -1 An allocation error occurred. A message indicating the offending array is written on unit control.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 control.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 requirement that type contains its relevant string ‘dense’, ‘coordinate’ or ‘sparse_by_rows’ has been violated.
n	is a scalar variable of type ipc_, that holds the number of variables
H_type	is a one-dimensional array of type char that specifies the symmetric storage scheme used for the Hessian. 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 H in the sparse co-ordinate storage scheme. It need not be set for any of the other schemes.
H_row	is a one-dimensional array of size ne and type ipc_, that holds the row indices of the lower triangular part of H 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
H_col	is a one-dimensional array of size ne and type ipc_, that holds the column indices of the lower triangular part of H in either the sparse co-ordinate, or the sparse row-wise storage scheme. It need not be set when the dense or diagonal storage schemes are used, and in this case can be NULL
H_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 H, 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 dps_reset_control(
    struct dps_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 dps_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:

1. The import was successful, and the package is ready for the solve phase

void dps_solve_tr_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    ipc_ ne,
    rpc_ H_val[],
    rpc_ c[],
    rpc_ f,
    rpc_ radius,
    rpc_ x[]
)

Find the global minimizer of the trust-region problem (1).

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: 0 The run was successful -1 An allocation error occurred. A message indicating the offending array is written on unit control.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 control.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 requirement that type contains its relevant string ‘dense’, ‘coordinate’ or ‘sparse_by_rows’ has been violated. -9 The analysis phase of the factorization failed; the return status from the factorization package is given in the component inform.factor_status -10 The factorization failed; the return status from the factorization package is given in the component inform.factor_status. -16 The problem is so ill-conditioned that further progress is impossible. -40 An error has occured when building the preconditioner.
n	is a scalar variable of type ipc_, that holds the number of variables
ne	is a scalar variable of type ipc_, that holds the number of entries in the lower triangular part of the Hessian matrix $H$.
H_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 Hessian matrix $H$ in any of the available storage schemes.
c	is a one-dimensional array of size n and type rpc_, that holds the linear term $c$ in the objective function. The j-th component of c, j = 0, … , n-1, contains $c_j$.
f	is a scalar variable pointer of type rpc_, that holds the value of the holds the constant term $f$ in the objective function.
radius	is a scalar variable pointer of type rpc_, that holds the value of the trust-region radius, $\Delta > 0$.
x	is a one-dimensional array of size n and type rpc_, that holds the values $x$ of the optimization variables. The j-th component of x, j = 0, … , n-1, contains $x_j$.

void dps_solve_rq_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    ipc_ ne,
    rpc_ H_val[],
    rpc_ c[],
    rpc_ f,
    rpc_ power,
    rpc_ weight,
    rpc_ x[]
)

Find the global minimizer of the regularized-quadartic problem (2).

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: 0 The run was successful -1 An allocation error occurred. A message indicating the offending array is written on unit control.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 control.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 requirement that type contains its relevant string ‘dense’, ‘coordinate’ or ‘sparse_by_rows’ has been violated. -9 The analysis phase of the factorization failed; the return status from the factorization package is given in the component inform.factor_status -10 The factorization failed; the return status from the factorization package is given in the component inform.factor_status. -16 The problem is so ill-conditioned that further progress is impossible. -40 An error has occured when building the preconditioner.
n	is a scalar variable of type ipc_, that holds the number of variables
ne	is a scalar variable of type ipc_, that holds the number of entries in the lower triangular part of the Hessian matrix $H$.
H_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 Hessian matrix $H$ in any of the available storage schemes.
c	is a one-dimensional array of size n and type rpc_, that holds the linear term $c$ in the objective function. The j-th component of c, j = 0, … , n-1, contains $c_j$.
f	is a scalar variable pointer of type rpc_, that holds the value of the holds the constant term $f$ in the objective function.
weight	is a scalar variable pointer of type rpc_, that holds the value of the regularization weight, $\sigma > 0$.
power	is a scalar variable pointer of type rpc_, that holds the value of the regularization power, $p \geq 2$.
x	is a one-dimensional array of size n and type rpc_, that holds the values $x$ of the optimization variables. The j-th component of x, j = 0, … , n-1, contains $x_j$.

void dps_resolve_tr_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    rpc_ c[],
    rpc_ f,
    rpc_ radius,
    rpc_ x[]
)

Find the global minimizer of the trust-region problem (1) if some non-matrix components have changed since a call to dps_solve_tr_problem.

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: 0 The run was successful -1 An allocation error occurred. A message indicating the offending array is written on unit control.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 control.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 requirement that type contains its relevant string ‘dense’, ‘coordinate’ or ‘sparse_by_rows’ has been violated. -16 The problem is so ill-conditioned that further progress is impossible.
n	is a scalar variable of type ipc_, that holds the number of variables
c	is a one-dimensional array of size n and type rpc_, that holds the linear term $c$ in the objective function. The j-th component of c, j = 0, … , n-1, contains $c_j$.
f	is a scalar variable pointer of type rpc_, that holds the value of the constant term $f$ in the objective function.
radius	is a scalar variable pointer of type rpc_, that holds the value of the trust-region radius, $\Delta > 0$.
x	is a one-dimensional array of size n and type rpc_, that holds the values $x$ of the optimization variables. The j-th component of x, j = 0, … , n-1, contains $x_j$.

void dps_resolve_rq_problem(
    void **data,
    ipc_ *status,
    ipc_ n,
    rpc_ c[],
    rpc_ f,
    rpc_ power,
    rpc_ weight,
    rpc_ x[]
)

Find the global minimizer of the regularized-quadartic problem (2) if some non-matrix components have changed since a call to dps_solve_rq_problem.

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: 0 The run was successful -1 An allocation error occurred. A message indicating the offending array is written on unit control.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 control.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. -16 The problem is so ill-conditioned that further progress is impossible.
n	is a scalar variable of type ipc_, that holds the number of variables
c	is a one-dimensional array of size n and type rpc_, that holds the linear term $c$ in the objective function. The j-th component of c, j = 0, … , n-1, contains $c_j$.
f	is a scalar variable pointer of type rpc_, that holds the value of the holds the constant term $f$ in the objective function.
weight	is a scalar variable pointer of type rpc_, that holds the value of the regularization weight, $\sigma > 0$.
power	is a scalar variable pointer of type rpc_, that holds the value of the regularization power, $p \geq 2$.
x	is a one-dimensional array of size n and type rpc_, that holds the values $x$ of the optimization variables. The j-th component of x, j = 0, … , n-1, contains $x_j$.

void dps_information(void **data, struct dps_inform_type* inform, ipc_ *status)

Provides output information

Parameters:

data

holds private internal data

inform

is a struct containing output information (see dps_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 dps_terminate(
    void **data,
    struct dps_control_type* control,
    struct dps_inform_type* inform
)

Deallocate all internal private storage

Parameters:

data	holds private internal data
control	is a struct containing control information (see dps_control_type)
inform	is a struct containing output information (see dps_inform_type)

available structures#

dps_control_type structure#

#include <galahad_dps.h>

struct dps_control_type {
    // fields

    bool f_indexing;
    ipc_ error;
    ipc_ out;
    ipc_ problem;
    ipc_ print_level;
    ipc_ new_h;
    ipc_ taylor_max_degree;
    rpc_ eigen_min;
    rpc_ lower;
    rpc_ upper;
    rpc_ stop_normal;
    rpc_ stop_absolute_normal;
    bool goldfarb;
    bool space_critical;
    bool deallocate_error_fatal;
    char problem_file[31];
    char symmetric_linear_solver[31];
    char prefix[31];
    struct sls_control_type sls_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_ problem

unit to write problem data into file problem_file

ipc_ print_level

controls level of diagnostic output

ipc_ new_h

how much of $H$ has changed since the previous call. Possible values are

0 unchanged
1 values but not indices have changed
2 values and indices have changed

ipc_ taylor_max_degree

maximum degree of Taylor approximant allowed

rpc_ eigen_min

smallest allowable value of an eigenvalue of the block diagonal factor of $H$

rpc_ lower

lower and upper bounds on the multiplier, if known

rpc_ upper

see lower

rpc_ stop_normal

stop trust-region solution when $| ||x||_M - \delta | \leq$ max( .stop_normal \* delta, .stop_absolute_normal )

rpc_ stop_absolute_normal

see stop_normal

bool goldfarb

use the Goldfarb variant of the trust-region/regularization norm rather than the modified absolute-value version

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 problem_file[31]

name of file into which to write problem data

char symmetric_linear_solver[31]

the name of the symmetric-indefinite linear equation solver used. Possible choices are currently: ‘sils’, ‘ma27’, ‘ma57’, ‘ma77’, ‘ma86’, ‘ma97’, ‘ssids’, ‘mumps’, ‘pardiso’, ‘mkl_pardiso’, ‘pastix’, ‘wsmp’, and ‘sytr’, although only ‘sytr’ 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 the Cholesky factorization and solution

dps_time_type structure#

#include <galahad_dps.h>

struct dps_time_type {
    // fields

    rpc_ total;
    rpc_ analyse;
    rpc_ factorize;
    rpc_ solve;
    rpc_ clock_total;
    rpc_ clock_analyse;
    rpc_ clock_factorize;
    rpc_ clock_solve;
};

detailed documentation#

time derived type as a C struct

components#

rpc_ total

total CPU time spent in the package

rpc_ analyse

CPU time spent reordering H prior to factorization.

rpc_ factorize

CPU time spent factorizing H.

rpc_ solve

CPU time spent solving the diagonal model system.

rpc_ clock_total

total clock time spent in the package

rpc_ clock_analyse

clock time spent reordering H prior to factorization

rpc_ clock_factorize

clock time spent factorizing H

rpc_ clock_solve

clock time spent solving the diagonal model system

dps_inform_type structure#

#include <galahad_dps.h>

struct dps_inform_type {
    // fields

    ipc_ status;
    ipc_ alloc_status;
    ipc_ mod_1by1;
    ipc_ mod_2by2;
    rpc_ obj;
    rpc_ obj_regularized;
    rpc_ x_norm;
    rpc_ multiplier;
    rpc_ pole;
    bool hard_case;
    char bad_alloc[81];
    struct dps_time_type time;
    struct sls_inform_type sls_inform;
};

detailed documentation#

inform derived type as a C struct

components#

ipc_ status

return status. See DPS_solve for details

ipc_ alloc_status

STAT value after allocate failure.

ipc_ mod_1by1

the number of 1 by 1 blocks from the factorization of H that were modified when constructing $M$

ipc_ mod_2by2

the number of 2 by 2 blocks from the factorization of H that were modified when constructing $M$

rpc_ obj

the value of the quadratic function

rpc_ obj_regularized

the value of the regularized quadratic function

rpc_ x_norm

the M-norm of the solution

rpc_ multiplier

the Lagrange multiplier associated with the constraint/regularization

rpc_ pole

a lower bound max(0,-lambda_1), where lambda_1 is the left-most eigenvalue of $(H,M)$

bool hard_case

has the hard case occurred?

char bad_alloc[81]

name of array that provoked an allocate failure

struct dps_time_type time

time information

struct sls_inform_type sls_inform

information from SLS

Index

example calls#

This is an example of how to use the package to solve a trust-region subproblem; the code is available in $GALAHAD/src/dps/C/dpst.c . A variety of supported Hessian and constraint 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 preproccesor 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 preproccesor variable INTEGER_64 is defined.

/* dpst.c */
/* Full test for the DPS C interface using C sparse matrix indexing */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_dps.h"
#ifdef REAL_128
#include <quadmath.h>
#endif

int main(void) {

    // Derived types
    void *data;
    struct dps_control_type control;
    struct dps_inform_type inform;

    // Set problem data
    ipc_ n = 3; // dimension of H
    ipc_ H_ne = 4; // number of elements of H
    ipc_ H_dense_ne = 6; // number of elements of H
    ipc_ H_row[] = {0, 1, 2, 2}; // row indices, NB lower triangle
    ipc_ H_col[] = {0, 1, 2, 0};
    ipc_ H_ptr[] = {0, 1, 2, 4};
    rpc_ H_val[] = {1.0, 2.0, 3.0, 4.0};
    rpc_ H_dense[] = {1.0, 0.0, 2.0, 4.0, 0.0, 3.0};
    rpc_ f = 0.96;
    rpc_ radius = 1.0;
    rpc_ half_radius = 0.5;
    rpc_ c[] = {0.0, 2.0, 0.0};

    char st = ' ';
    ipc_ status;
    rpc_ x[n];

    printf(" C sparse matrix indexing\n\n");

    printf(" basic tests of storage formats\n\n");

    for( ipc_ storage_type=1; storage_type <= 3; storage_type++){

    // Initialize DPS
      dps_initialize( &data, &control, &status );

    // Set user-defined control options
      control.f_indexing = false; // C sparse matrix indexing
      strcpy(control.symmetric_linear_solver,"sytr ") ;

      switch(storage_type){
        case 1: // sparse co-ordinate storage
            st = 'C';
            // import the control parameters and structural data
            dps_import( &control, &data, &status, n,
                       "coordinate", H_ne, H_row, H_col, NULL );
            // solve the problem
            dps_solve_tr_problem( &data, &status, n, H_ne, H_val,
                                  c, f, radius, x );
            break;
        case 2: // sparse by rows
            st = 'R';
            // import the control parameters and structural data
            dps_import( &control, &data, &status, n,
                        "sparse_by_rows", H_ne, NULL, H_col, H_ptr );
            dps_solve_tr_problem( &data, &status, n, H_ne, H_val,
                                  c, f, radius, x );
            break;
        case 3: // dense
            st = 'D';
            // import the control parameters and structural data
            dps_import( &control, &data, &status, n,
                        "dense", H_ne, NULL, NULL, NULL );
            dps_solve_tr_problem( &data, &status, n, H_dense_ne, H_dense,
                                  c, f, radius, x );
            break;
        }

      dps_information( &data, &inform, &status );
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_dps.h
#include "galahad_pquad_dps.h"
#else
      printf("format %c: DPS_solve_problem exit status   = %1" i_ipc_
             ", f = %.2f\n", st, inform.status, inform.obj );
#endif

      switch(storage_type){
        case 1: // sparse co-ordinate storage
            st = 'C';
            // solve the problem
            dps_resolve_tr_problem( &data, &status, n,
                                    c, f, half_radius, x );
            break;
        case 2: // sparse by rows
            st = 'R';
            dps_resolve_tr_problem( &data, &status, n,
                                    c, f, half_radius, x );
            break;
        case 3: // dense
            st = 'D';
            dps_resolve_tr_problem( &data, &status, n,
                                    c, f, half_radius, x );
            break;
        }

      dps_information( &data, &inform, &status );
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_dpsr.h
#include "galahad_pquad_dpsr.h"
#else
      printf("format %c: DPS_resolve_problem exit status = %1" i_ipc_
             ", f = %.2f\n", st, inform.status, inform.obj );
#endif
      //printf("x: ");
      //for( ipc_ i = 0; i < n+m; i++) printf("%f ", x[i]);

      // Delete internal workspace
      dps_terminate( &data, &control, &inform );
   }
}

This is the same example, but now fortran-style indexing is used; the code is available in $GALAHAD/src/dps/C/dpstf.c .

/* dpstf.c */
/* Full test for the DPS C interface using Fortran sparse matrix indexing */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "galahad_precision.h"
#include "galahad_cfunctions.h"
#include "galahad_dps.h"
#ifdef REAL_128
#include <quadmath.h>
#endif

int main(void) {

    // Derived types
    void *data;
    struct dps_control_type control;
    struct dps_inform_type inform;

    // Set problem data
    ipc_ n = 3; // dimension of H
    ipc_ H_ne = 4; // number of elements of H
    ipc_ H_dense_ne = 6; // number of elements of H
    ipc_ H_row[] = {1, 2, 3, 3}; // row indices, NB lower triangle
    ipc_ H_col[] = {1, 2, 3, 1};
    ipc_ H_ptr[] = {1, 2, 3, 5};
    rpc_ H_val[] = {1.0, 2.0, 3.0, 4.0};
    rpc_ H_dense[] = {1.0, 0.0, 2.0, 4.0, 0.0, 3.0};
    rpc_ f = 0.96;
    rpc_ radius = 1.0;
    rpc_ half_radius = 0.5;
    rpc_ c[] = {0.0, 2.0, 0.0};

    char st = ' ';
    ipc_ status;
    rpc_ x[n];

    printf(" Fortran sparse matrix indexing\n\n");

    printf(" basic tests of storage formats\n\n");

    for( ipc_ storage_type=1; storage_type <= 3; storage_type++){

    // Initialize DPS
      dps_initialize( &data, &control, &status );

    // Set user-defined control options
      control.f_indexing = true; // fortran sparse matrix indexing
      strcpy(control.symmetric_linear_solver,"sytr ") ;

      switch(storage_type){
        case 1: // sparse co-ordinate storage
            st = 'C';
            // import the control parameters and structural data
            dps_import( &control, &data, &status, n,
                       "coordinate", H_ne, H_row, H_col, NULL );
            // solve the problem
            dps_solve_tr_problem( &data, &status, n, H_ne, H_val,
                                  c, f, radius, x );
            break;
        case 2: // sparse by rows
            st = 'R';
            // import the control parameters and structural data
            dps_import( &control, &data, &status, n,
                        "sparse_by_rows", H_ne, NULL, H_col, H_ptr );
            dps_solve_tr_problem( &data, &status, n, H_ne, H_val,
                                  c, f, radius, x );
            break;
        case 3: // dense
            st = 'D';
            // import the control parameters and structural data
            dps_import( &control, &data, &status, n,
                        "dense", H_ne, NULL, NULL, NULL );
            dps_solve_tr_problem( &data, &status, n, H_dense_ne, H_dense,
                                  c, f, radius, x );
            break;
        }

      dps_information( &data, &inform, &status );
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_dps.h
#include "galahad_pquad_dps.h"
#else
      printf("format %c: DPS_solve_problem exit status   = %1" i_ipc_
             ", f = %.2f\n", st, inform.status, inform.obj );
#endif

      switch(storage_type){
        case 1: // sparse co-ordinate storage
            st = 'C';
            // solve the problem
            dps_resolve_tr_problem( &data, &status, n,
                                    c, f, half_radius, x );
            break;
        case 2: // sparse by rows
            st = 'R';
            dps_resolve_tr_problem( &data, &status, n,
                                    c, f, half_radius, x );
            break;
        case 3: // dense
            st = 'D';
            dps_resolve_tr_problem( &data, &status, n,
                                    c, f, half_radius, x );
            break;
        }

      dps_information( &data, &inform, &status );
#ifdef REAL_128
// interim replacement for quad output: $GALAHAD/include/galahad_pquad_dpsr.h
#include "galahad_pquad_dpsr.h"
#else
      printf("format %c: DPS_resolve_problem exit status = %1" i_ipc_
             ", f = %.2f\n", st, inform.status, inform.obj );
#endif
      //printf("x: ");
      //for( ipc_ i = 0; i < n+m; i++) printf("%f ", x[i]);

      // Delete internal workspace
      dps_terminate( &data, &control, &inform );
   }
}