GALAHAD WCP package#
purpose#
The wcp
package uses a primal-dual interior-point method to find a
well-centered point within a polyhedral set.
The aim is to find a point that lies interior to the boundary of the
polyhedron definied by the general linear constraints and simple bounds
See Section 4 of $GALAHAD/doc/wcp.pdf for a brief description of the method employed and other details.
terminolgy#
More specifically, if possible, the package finds a solution to the system of primal optimality equations
method#
The algorithm is iterative, and at each major iteration attempts to find a solution to the perturbed system (1), (2),
Given any solution to (1)–(2) and (10) satisfying (9),
the perturbations are reduced (sometimes to zero) so as to ensure that the
current solution is feasible for the next perturbed problem. Specifically,
the perturbation
Ultimately the intention is to drive all the perturbations to zero. It can be shown that if the original problem (1)–(4) and (6) has a solution, the perturbations will be zero after a finite number of major iterations. Equally, if there is no interior solution (6), the sets of (primal and dual) variables that are necessarily at (one of) their bounds for all feasible points—we refer to these as implicit equalities—will be identified, as will the possibility that there is no point (interior or otherwise) in the primal and/or dual feasible regions.
Each major iteration requires the solution
SBLS
.
In order to make the solution as efficient as possible, the
variables and constraints are reordered internally
by the package QPP
prior to solution.
In particular, fixed variables, and
free (unbounded on both sides) constraints are temporarily removed.
In addition, an attempt to identify and remove linearly dependent
equality constraints may be made by factorizing
SBLS
, and examining small pivot blocks.
reference#
The basic algorithm, its convergence analysis and results of numerical experiments are given in
C. Cartis and N. I. M. Gould, Finding a point n the relative interior of a polyhedron. Technical Report TR-2006-016, Rutherford Appleton Laboratory (2006).
matrix storage#
The unsymmetric
Dense storage format:
The matrix
Dense by columns storage format:
The matrix
Sparse co-ordinate storage format:
Only the nonzero entries of the matrices are stored.
For the
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
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
introduction to function calls#
To solve a given problem, functions from the wcp package must be called in the following order:
wcp_initialize - provide default control parameters and set up initial data structures
wcp_read_specfile (optional) - override control values by reading replacement values from a file
wcp_import - set up problem data structures and fixed values
wcp_reset_control (optional) - possibly change control parameters if a sequence of problems are being solved
wcp_find_wcp - find a well-centered point
wcp_information (optional) - recover information about the solution and solution process
wcp_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), Float64 (double precision) or, if supported, Float128 (quadruple precision).
callable functions#
function wcp_initialize(T, data, control, status)
Set default control values and initialize private data
Parameters:
data |
holds private internal data |
control |
is a structure containing control information (see wcp_control_type) |
status |
is a scalar variable of type Int32 that gives the exit status from the package. Possible values are (currently):
|
function wcp_read_specfile(T, control, 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/wcp/WCP.template. See also Table 2.1 in the Fortran documentation provided in $GALAHAD/doc/wcp.pdf for a list of how these keywords relate to the components of the control structure.
Parameters:
control |
is a structure containing control information (see wcp_control_type) |
specfile |
is a one-dimensional array of type Vararg{Cchar} that must give the name of the specification file |
function wcp_import(T, control, data, status, n, m, A_type, A_ne, A_row, A_col, A_ptr)
Import problem data into internal storage prior to solution.
Parameters:
control |
is a structure whose members provide control parameters for the remaining procedures (see wcp_control_type) |
data |
holds private internal data |
status |
is a scalar variable of type Int32 that gives the exit status from the package. Possible values are:
|
n |
is a scalar variable of type Int32 that holds the number of variables. |
m |
is a scalar variable of type Int32 that holds the number of general linear constraints. |
A_type |
is a one-dimensional array of type Vararg{Cchar} that specifies the unsymmetric storage scheme used for the constraint Jacobian, |
A_ne |
is a scalar variable of type Int32 that holds the number of entries in |
A_row |
is a one-dimensional array of size A_ne and type Int32 that holds the row indices of |
A_col |
is a one-dimensional array of size A_ne and type Int32 that holds the column indices of |
A_ptr |
is a one-dimensional array of size n+1 and type Int32 that holds the starting position of each row of |
function wcp_reset_control(T, control, data, status)
Reset control parameters after import if required.
Parameters:
control |
is a structure whose members provide control parameters for the remaining procedures (see wcp_control_type) |
data |
holds private internal data |
status |
is a scalar variable of type Int32 that gives the exit status from the package. Possible values are:
|
function wcp_find_wcp(T, data, status, n, m, g, a_ne, A_val, c_l, c_u, x_l, x_u, x, c, y_l, y_u, z_l, z_u, x_stat, c_stat)
Find a well-centered point in the feasible region
Parameters:
data |
holds private internal data |
status |
is a scalar variable of type Int32 that gives the entry and exit status from the package. Possible exit values are:
|
n |
is a scalar variable of type Int32 that holds the number of variables |
m |
is a scalar variable of type Int32 that holds the number of general linear constraints. |
g |
is a one-dimensional array of size n and type T that holds the target vector |
a_ne |
is a scalar variable of type Int32 that holds the number of entries in the constraint Jacobian matrix |
A_val |
is a one-dimensional array of size a_ne and type T that holds the values of the entries of the constraint Jacobian matrix |
c_l |
is a one-dimensional array of size m and type T that holds the lower bounds |
c_u |
is a one-dimensional array of size m and type T that holds the upper bounds |
x_l |
is a one-dimensional array of size n and type T that holds the lower bounds |
x_u |
is a one-dimensional array of size n and type T that holds the upper bounds |
x |
is a one-dimensional array of size n and type T that holds the values |
c |
is a one-dimensional array of size m and type T that holds the residual |
y_l |
is a one-dimensional array of size n and type T that holds the values |
y_u |
is a one-dimensional array of size n and type T that holds the values |
z_l |
is a one-dimensional array of size n and type T that holds the values |
z_u |
is a one-dimensional array of size n and type T that holds the values |
x_stat |
is a one-dimensional array of size n and type Int32 that gives the optimal status of the problem variables. If x_stat(j) is negative, the variable |
c_stat |
is a one-dimensional array of size m and type Int32 that gives the optimal status of the general linear constraints. If c_stat(i) is negative, the constraint value |
function wcp_information(T, data, inform, status)
Provides output information.
Parameters:
data |
holds private internal data |
inform |
is a structure containing output information (see wcp_inform_type) |
status |
is a scalar variable of type Int32 that gives the exit status from the package. Possible values are (currently):
|
function wcp_terminate_s(data, control, inform)
Deallocate all internal private storage.
Parameters:
data |
holds private internal data |
control |
is a structure containing control information (see wcp_control_type) |
inform |
is a structure containing output information (see wcp_inform_type) |
available structures#
wcp_control_type structure#
struct wcp_control_type{T} f_indexing::Bool error::Int32 out::Int32 print_level::Int32 start_print::Int32 stop_print::Int32 maxit::Int32 initial_point::Int32 factor::Int32 max_col::Int32 indmin::Int32 valmin::Int32 itref_max::Int32 infeas_max::Int32 perturbation_strategy::Int32 restore_problem::Int32 infinity::T stop_p::T stop_d::T stop_c::T prfeas::T dufeas::T mu_target::T mu_accept_fraction::T mu_increase_factor::T required_infeas_reduction::T implicit_tol::T pivot_tol::T pivot_tol_for_dependencies::T zero_pivot::T perturb_start::T alpha_scale::T identical_bounds_tol::T reduce_perturb_factor::T reduce_perturb_multiplier::T insufficiently_feasible::T perturbation_small::T cpu_time_limit::T clock_time_limit::T remove_dependencies::Bool treat_zero_bounds_as_general::Bool just_feasible::Bool balance_initial_complementarity::Bool use_corrector::Bool space_critical::Bool deallocate_error_fatal::Bool record_x_status::Bool record_c_status::Bool prefix::NTuple{31,Cchar} fdc_control::fdc_control_type{T} sbls_control::sbls_control_type{T}
detailed documentation#
control derived type as a Julia structure
components#
Bool f_indexing
use C or Fortran sparse matrix indexing
Int32 error
error and warning diagnostics occur on stream error
Int32 out
general output occurs on stream out
Int32 print_level
the level of output required is specified by print_level
Int32 start_print
any printing will start on this iteration
Int32 stop_print
any printing will stop on this iteration
Int32 maxit
at most maxit inner iterations are allowed
Int32 initial_point
how to choose the initial point. Possible values are
0 the values input in X, shifted to be at least prfeas from their nearest bound, will be used
1 the nearest point to the “bound average” 0.5(X_l+X_u) that satisfies the linear constraints will be used
Int32 factor
the factorization to be used. Possible values are
0 automatic
1 Schur-complement factorization
2 augmented-system factorization
Int32 max_col
the maximum number of nonzeros in a column of A which is permitted with the Schur-complement factorization
Int32 indmin
an initial guess as to the integer workspace required by SBLS
Int32 valmin
an initial guess as to the real workspace required by SBLS
Int32 itref_max
the maximum number of iterative refinements allowed
Int32 infeas_max
the number of iterations for which the overall infeasibility of the problem is not reduced by at least a factor .required_infeas_reduction before the problem is flagged as infeasible (see required_infeas_reducti
Int32 perturbation_strategy
the strategy used to reduce relaxed constraint bounds. Possible values are
0 do not perturb the constraints
1 reduce all perturbations by the same amount with linear reduction
2 reduce each perturbation as much as possible with linear reduction
3 reduce all perturbations by the same amount with superlinear reduction
4 reduce each perturbation as much as possible with superlinear reduction
Int32 restore_problem
indicate whether and how much of the input problem should be restored on output. Possible values are
0 nothing restored
1 scalar and vector parameters
2 all parameters
T infinity
any bound larger than infinity in modulus will be regarded as infinite
T stop_p
the required accuracy for the primal infeasibility
T stop_d
the required accuracy for the dual infeasibility
T stop_c
the required accuracy for the complementarity
T prfeas
initial primal variables will not be closer than prfeas from their bound
T dufeas
initial dual variables will not be closer than dufeas from their bounds
T mu_target
the target value of the barrier parameter. If mu_target is not positive, it will be reset to an appropriate value
T mu_accept_fraction
the complemtary slackness x_i.z_i will be judged to lie within an acceptable margin around its target value mu as soon as mu_accept_fraction * mu <= x_i.z_i <= ( 1 / mu_accept_fraction ) * mu; the perturbations will be reduced as soon as all of the complemtary slacknesses x_i.z_i lie within acceptable bounds. mu_accept_fraction will be reset to ensure that it lies in the interval (0,1]
T mu_increase_factor
the target value of the barrier parameter will be increased by mu_increase_factor for infeasible constraints every time the perturbations are adjusted
T required_infeas_reduction
if the overall infeasibility of the problem is not reduced by at least a factor required_infeas_reduction over .infeas_max iterations, the problem is flagged as infeasible (see infeas_max)
T implicit_tol
any primal or dual variable that is less feasible than implicit_tol will be regarded as defining an implicit constraint
T pivot_tol
the threshold pivot used by the matrix factorization. See the documentation for SBLS for details (obsolete)
T pivot_tol_for_dependencies
the threshold pivot used by the matrix factorization when attempting to detect linearly dependent constraints. See the documentation for SBLS for details (obsolete)
T zero_pivot
any pivots smaller than zero_pivot in absolute value will be regarded to zero when attempting to detect linearly dependent constraints (obsolete)
T perturb_start
the constraint bounds will initially be relaxed by .perturb_start; this perturbation will subsequently be reduced to zero. If perturb_start < 0, the amount by which the bounds are relaxed will be computed automatically
T alpha_scale
the test for rank defficiency will be to factorize ( alpha_scale I A^T ) ( A 0 )
T identical_bounds_tol
any pair of constraint bounds (c_l,c_u) or (x_l,x_u) that are closer tha identical_bounds_tol will be reset to the average of their values
T reduce_perturb_factor
the constraint perturbation will be reduced as follows:
if the variable lies outside a bound, the corresponding perturbation will be reduced to reduce_perturb_factor * current pertubation
( 1 - reduce_perturb_factor ) * violation
otherwise, if the variable lies within insufficiently_feasible of its bound the pertubation will be reduced to reduce_perturb_multiplier * current pertubation
otherwise if will be set to zero
T reduce_perturb_multiplier
see reduce_perturb_factor
T insufficiently_feasible
see reduce_perturb_factor
T perturbation_small
if the maximum constraint pertubation is smaller than perturbation_small and the violation is smaller than implicit_tol, the method will deduce that there is a feasible point but no interior
T cpu_time_limit
the maximum CPU time allowed (-ve means infinite)
T clock_time_limit
the maximum elapsed clock time allowed (-ve means infinite)
Bool remove_dependencies
the equality constraints will be preprocessed to remove any linear dependencies if true
Bool treat_zero_bounds_as_general
any problem bound with the value zero will be treated as if it were a general value if true
Bool just_feasible
if .just_feasible is true, the algorithm will stop as soon as a feasible point is found. Otherwise, the optimal solution to the problem will be found
Bool balance_initial_complementarity
if .balance_initial_complementarity is .true. the initial complemetarity will be balanced
Bool use_corrector
if .use_corrector, a corrector step will be used
Bool space_critical
if .space_critical true, every effort will be made to use as little space as possible. This may result in longer computation time
Bool deallocate_error_fatal
if .deallocate_error_fatal is true, any array/pointer deallocation error will terminate execution. Otherwise, computation will continue
Bool record_x_status
if .record_x_status is true, the array inform.X_status will be allocated and the status of the bound constraints will be reported on exit.
Bool record_c_status
if .record_c_status is true, the array inform.C_status will be allocated and the status of the general constraints will be reported on exit.
NTuple{31,Cchar} prefix
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’
fdc_control_type{T} fdc_control_type fdc_control
control parameters for FDC
sbls_control_type{T} sbls_control_type sbls_control
control parameters for SBLS
wcp_time_type structure#
struct wcp_time_type{T} total::T preprocess::T find_dependent::T analyse::T factorize::T solve::T clock_total::T clock_preprocess::T clock_find_dependent::T clock_analyse::T clock_factorize::T clock_solve::T
detailed documentation#
time derived type as a Julia structure
components#
T total
the total CPU time spent in the package
T preprocess
the CPU time spent preprocessing the problem
T find_dependent
the CPU time spent detecting linear dependencies
T analyse
the CPU time spent analysing the required matrices prior to factorization
T factorize
the CPU time spent factorizing the required matrices
T solve
the CPU time spent computing the search direction
T clock_total
the total clock time spent in the package
T clock_preprocess
the clock time spent preprocessing the problem
T clock_find_dependent
the clock time spent detecting linear dependencies
T clock_analyse
the clock time spent analysing the required matrices prior to factorization
T clock_factorize
the clock time spent factorizing the required matrices
T clock_solve
the clock time spent computing the search direction
wcp_inform_type structure#
struct wcp_inform_type{T} status::Int32 alloc_status::Int32 bad_alloc::NTuple{81,Cchar} iter::Int32 factorization_status::Int32 factorization_integer::Int64 factorization_real::Int64 nfacts::Int32 c_implicit::Int32 x_implicit::Int32 y_implicit::Int32 z_implicit::Int32 obj::T mu_final_target_max::T non_negligible_pivot::T feasible::Bool time::wcp_time_type{T} fdc_inform::fdc_inform_type{T} sbls_inform::sbls_inform_type{T}
detailed documentation#
inform derived type as a Julia structure
components#
Int32 status
return status. See WCP_solve for details
Int32 alloc_status
the status of the last attempted allocation/deallocation
NTuple{81,Cchar} bad_alloc
the name of the array for which an allocation/deallocation error occurred
Int32 iter
the total number of iterations required
Int32 factorization_status
the return status from the factorization
Int64 factorization_integer
the total integer workspace required for the factorization
Int64 factorization_real
the total real workspace required for the factorization
Int32 nfacts
the total number of factorizations performed
Int32 c_implicit
the number of general constraints that lie on (one) of their bounds for feasible solutions
Int32 x_implicit
the number of variables that lie on (one) of their bounds for all feasible solutions
Int32 y_implicit
the number of Lagrange multipliers for general constraints that lie on (one) of their bounds for all feasible solutions
Int32 z_implicit
the number of dual variables that lie on (one) of their bounds for all feasible solutions
T obj
the value of the objective function at the best estimate of the solution determined by WCP_solve
T mu_final_target_max
the largest target value on termination
T non_negligible_pivot
the smallest pivot which was not judged to be zero when detecting linear dependent constraints
Bool feasible
is the returned primal-dual “solution” strictly feasible?
wcp_time_type{T} wcp_time_type time
timings (see above)
fdc_inform_type{T} fdc_inform_type fdc_inform
inform parameters for FDC
sbls_inform_type{T} sbls_inform_type sbls_inform
inform parameters for SBLS
example calls#
This is an example of how to use the package to find a well-centred point; the code is available in $GALAHAD/src/wcp/Julia/test_wcp.jl . A variety of supported Hessian and constraint matrix storage formats are shown.
# test_wcp.jl
# Simple code to test the Julia interface to WCP
using GALAHAD
using Test
using Printf
using Accessors
using Quadmath
function test_wcp(::Type{T}) where T
# Derived types
data = Ref{Ptr{Cvoid}}()
control = Ref{wcp_control_type{T}}()
inform = Ref{wcp_inform_type{T}}()
# Set problem data
n = 3 # dimension
m = 2 # number of general constraints
g = T[0.0, 2.0, 0.0] # linear term in the objective
A_ne = 4 # Jacobian elements
A_row = Cint[1, 1, 2, 2] # row indices
A_col = Cint[1, 2, 2, 3] # column indices
A_ptr = Cint[1, 3, 5] # row pointers
A_val = T[2.0, 1.0, 1.0, 1.0] # values
c_l = T[1.0, 2.0] # constraint lower bound
c_u = T[2.0, 2.0] # constraint upper bound
x_l = T[-1.0, -Inf, -Inf] # variable lower bound
x_u = T[1.0, Inf, 2.0] # variable upper bound
# Set output storage
c = zeros(T, m) # constraint values
x_stat = zeros(Cint, n) # variable status
c_stat = zeros(Cint, m) # constraint status
st = ' '
status = Ref{Cint}()
@printf(" Fortran sparse matrix indexing\n\n")
@printf(" basic tests of wcp storage formats\n\n")
for d in 1:3
# Initialize WCP
wcp_initialize(T, data, control, status)
# Set user-defined control options
@reset control[].f_indexing = true # Fortran sparse matrix indexing
# Start from 0
x = T[0.0, 0.0, 0.0]
y_l = T[0.0, 0.0]
y_u = T[0.0, 0.0]
z_l = T[0.0, 0.0, 0.0]
z_u = T[0.0, 0.0, 0.0]
# sparse co-ordinate storage
if d == 1
st = 'C'
wcp_import(T, control, data, status, n, m, "coordinate", A_ne, A_row, A_col, C_NULL)
wcp_find_wcp(T, data, status, n, m, g, A_ne, A_val,
c_l, c_u, x_l, x_u, x, c, y_l, y_u, z_l, z_u,
x_stat, c_stat)
end
# sparse by rows
if d == 2
st = 'R'
wcp_import(T, control, data, status, n, m,
"sparse_by_rows", A_ne, C_NULL, A_col, A_ptr)
wcp_find_wcp(T, data, status, n, m, g, A_ne, A_val,
c_l, c_u, x_l, x_u, x, c, y_l, y_u, z_l, z_u,
x_stat, c_stat)
end
# dense
if d == 3
st = 'D'
A_dense_ne = 6 # number of elements of A
A_dense = T[2.0, 1.0, 0.0, 0.0, 1.0, 1.0]
wcp_import(T, control, data, status, n, m,
"dense", A_dense_ne, C_NULL, C_NULL, C_NULL)
wcp_find_wcp(T, data, status, n, m, g, A_dense_ne, A_dense,
c_l, c_u, x_l, x_u, x, c, y_l, y_u, z_l, z_u,
x_stat, c_stat)
end
wcp_information(T, data, inform, status)
if inform[].status == 0
@printf("%c:%6i iterations. Optimal objective value = %5.2f status = %1i\n", st,
inform[].iter, inform[].obj, inform[].status)
else
@printf("%c: WCP_solve exit status = %1i\n", st, inform[].status)
end
# @printf("x: ")
# for i = 1:n
# @printf("%f ", x[i])
# end
# @printf("\n")
# @printf("gradient: ")
# for i = 1:n
# @printf("%f ", g[i])
# end
# @printf("\n")
# Delete internal workspace
wcp_terminate(T, data, control, inform)
end
return 0
end
@testset "WCP" begin
@test test_wcp(Float32) == 0
@test test_wcp(Float64) == 0
@test test_wcp(Float128) == 0
end