h02ce obtains integer solutions to sparse linear programming and quadratic programming problems.

Syntax

C#
public static void h02ce( int n, int m, int nnz, int iobj, int ncolh, H..::..H02CE_QPHX qphx, double[] a, int[] ha, int[] ka, double[] bl, double[] bu, string start, string[] names, int nname, string[] crname, ref int ns, double[] xs, int[] intvar, int lintvr, int mdepth, int[] istate, out int miniz, out int minz, out double obj, double[] clamda, int strtgy, int[] iz, double[] z, H..::..H02CE_MONIT monit, H..::..h02ceOptions options, out int ifail )

public static void h02ce(
	int n,
	int m,
	int nnz,
	int iobj,
	int ncolh,
	H..::..H02CE_QPHX qphx,
	double[] a,
	int[] ha,
	int[] ka,
	double[] bl,
	double[] bu,
	string start,
	string[] names,
	int nname,
	string[] crname,
	ref int ns,
	double[] xs,
	int[] intvar,
	int lintvr,
	int mdepth,
	int[] istate,
	out int miniz,
	out int minz,
	out double obj,
	double[] clamda,
	int strtgy,
	int[] iz,
	double[] z,
	H..::..H02CE_MONIT monit,
	H..::..h02ceOptions options,
	out int ifail
)

Visual Basic
Public Shared Sub h02ce ( _ n As Integer, _ m As Integer, _ nnz As Integer, _ iobj As Integer, _ ncolh As Integer, _ qphx As H..::..H02CE_QPHX, _ a As Double(), _ ha As Integer(), _ ka As Integer(), _ bl As Double(), _ bu As Double(), _ start As String, _ names As String(), _ nname As Integer, _ crname As String(), _ ByRef ns As Integer, _ xs As Double(), _ intvar As Integer(), _ lintvr As Integer, _ mdepth As Integer, _ istate As Integer(), _ <OutAttribute> ByRef miniz As Integer, _ <OutAttribute> ByRef minz As Integer, _ <OutAttribute> ByRef obj As Double, _ clamda As Double(), _ strtgy As Integer, _ iz As Integer(), _ z As Double(), _ monit As H..::..H02CE_MONIT, _ options As H..::..h02ceOptions, _ <OutAttribute> ByRef ifail As Integer _ )

Visual Basic

Public Shared Sub h02ce ( _
	n As Integer, _
	m As Integer, _
	nnz As Integer, _
	iobj As Integer, _
	ncolh As Integer, _
	qphx As H..::..H02CE_QPHX, _
	a As Double(), _
	ha As Integer(), _
	ka As Integer(), _
	bl As Double(), _
	bu As Double(), _
	start As String, _
	names As String(), _
	nname As Integer, _
	crname As String(), _
	ByRef ns As Integer, _
	xs As Double(), _
	intvar As Integer(), _
	lintvr As Integer, _
	mdepth As Integer, _
	istate As Integer(), _
	<OutAttribute> ByRef miniz As Integer, _
	<OutAttribute> ByRef minz As Integer, _
	<OutAttribute> ByRef obj As Double, _
	clamda As Double(), _
	strtgy As Integer, _
	iz As Integer(), _
	z As Double(), _
	monit As H..::..H02CE_MONIT, _
	options As H..::..h02ceOptions, _
	<OutAttribute> ByRef ifail As Integer _
)

Visual C++
public: static void h02ce( int n, int m, int nnz, int iobj, int ncolh, H..::..H02CE_QPHX^ qphx, array<double>^ a, array<int>^ ha, array<int>^ ka, array<double>^ bl, array<double>^ bu, String^ start, array<String^>^ names, int nname, array<String^>^ crname, int% ns, array<double>^ xs, array<int>^ intvar, int lintvr, int mdepth, array<int>^ istate, [OutAttribute] int% miniz, [OutAttribute] int% minz, [OutAttribute] double% obj, array<double>^ clamda, int strtgy, array<int>^ iz, array<double>^ z, H..::..H02CE_MONIT^ monit, H..::..h02ceOptions^ options, [OutAttribute] int% ifail )

Visual C++

public:
static void h02ce(
	int n, 
	int m, 
	int nnz, 
	int iobj, 
	int ncolh, 
	H..::..H02CE_QPHX^ qphx, 
	array<double>^ a, 
	array<int>^ ha, 
	array<int>^ ka, 
	array<double>^ bl, 
	array<double>^ bu, 
	String^ start, 
	array<String^>^ names, 
	int nname, 
	array<String^>^ crname, 
	int% ns, 
	array<double>^ xs, 
	array<int>^ intvar, 
	int lintvr, 
	int mdepth, 
	array<int>^ istate, 
	[OutAttribute] int% miniz, 
	[OutAttribute] int% minz, 
	[OutAttribute] double% obj, 
	array<double>^ clamda, 
	int strtgy, 
	array<int>^ iz, 
	array<double>^ z, 
	H..::..H02CE_MONIT^ monit, 
	H..::..h02ceOptions^ options, 
	[OutAttribute] int% ifail
)

F#
static member h02ce : n : int * m : int * nnz : int * iobj : int * ncolh : int * qphx : H..::..H02CE_QPHX * a : float[] * ha : int[] * ka : int[] * bl : float[] * bu : float[] * start : string * names : string[] * nname : int * crname : string[] * ns : int byref * xs : float[] * intvar : int[] * lintvr : int * mdepth : int * istate : int[] * miniz : int byref * minz : int byref * obj : float byref * clamda : float[] * strtgy : int * iz : int[] * z : float[] * monit : H..::..H02CE_MONIT * options : H..::..h02ceOptions * ifail : int byref -> unit

static member h02ce : 
        n : int * 
        m : int * 
        nnz : int * 
        iobj : int * 
        ncolh : int * 
        qphx : H..::..H02CE_QPHX * 
        a : float[] * 
        ha : int[] * 
        ka : int[] * 
        bl : float[] * 
        bu : float[] * 
        start : string * 
        names : string[] * 
        nname : int * 
        crname : string[] * 
        ns : int byref * 
        xs : float[] * 
        intvar : int[] * 
        lintvr : int * 
        mdepth : int * 
        istate : int[] * 
        miniz : int byref * 
        minz : int byref * 
        obj : float byref * 
        clamda : float[] * 
        strtgy : int * 
        iz : int[] * 
        z : float[] * 
        monit : H..::..H02CE_MONIT * 
        options : H..::..h02ceOptions * 
        ifail : int byref -> unit

Parameters

n: Type: System..::..Int32
On entry: $n$ , the number of variables (excluding slacks). This is the number of columns in the linear constraint matrix $A$ .

Constraint: $n \geq 1$ .

m: Type: System..::..Int32
On entry: $m$ , the number of general linear constraints (or slacks). This is the number of rows in $A$ , including the free row (if any; see iobj).

Constraint: $m \geq 1$ .

nnz: Type: System..::..Int32
On entry: the number of nonzero elements in $A$ .

Constraint: $1 \leq nnz \leq n \times m$ .

iobj: Type: System..::..Int32
On entry: if $iobj > 0$ , row iobj of $A$ is a free row containing the nonzero elements of the vector $c$ appearing in the linear objective term $c^{T} x$ .
If $iobj = 0$ , there is no free row, i.e., the problem is either an FP problem (in which case iobj must be set to zero), or a QP problem with $c = 0$ .

Constraint: $0 \leq iobj \leq m$ .

ncolh: Type: System..::..Int32
On entry: $n_{H}$ , the number of leading nonzero columns of the Hessian matrix $H$ . For FP and LP problems, ncolh must be set to zero.

Constraint: $0 \leq ncolh \leq n$ .

qphx

Type: NagLibrary..::..H..::..H02CE_QPHX

For QP problems, you must supply a version of qphx to compute the matrix product

H x

. If

H

has rows and columns consisting entirely of zeros, it is most efficient to order the variables

x = {(y z)}^{T}

so that

H x = (\begin{matrix} H_{1} & 0 \\ 0 & 0 \end{matrix}) (\begin{matrix} y \\ z \end{matrix}) = (\begin{matrix} H_{1} y \\ 0 \end{matrix}),

where the nonlinear variables

y

appear first as shown. For LP problems, qphx will never be called by h02ce.

A delegate of type H02CE_QPHX.

a: Type: array<System..::..Double>[]()[][]
An array of size [nnz]
On entry: the nonzero elements of $A$ , ordered by increasing column index. Note that multiple elements with the same row and column indices are not allowed.
On exit: used as internal workspace prior to being restored and hence is unchanged.

ha: Type: array<System..::..Int32>[]()[][]
An array of size [nnz]
On entry: $ha [i]$ must contain the row index of the nonzero element stored in $a [i]$ , for $i = 1, 2, \dots, nnz$ . Note that the row indices for a column may be supplied in any order.

Constraint: $1 \leq ha [i] \leq m$ , for $i = 0, 1, \dots, nnz - 1$ .

ka

Type: array<System..::..Int32>[]()[][]

An array of size [

n + 1

]

On entry:

ka [j]

must contain the index in a of the start of the

j

th column, for

j = 1, 2, \dots, n

. To specify the

j

th column as empty, set

ka [j] = ka [j + 1]

. Note that the first and last elements of ka must be such that

ka [0] = 1

and

ka [n] = nnz + 1

Constraints:

$ka [0] = 1$ ;
$ka [j] \geq 1$ , for $j = 1, 2, \dots, n - 1$ ;
$ka [n] = nnz + 1$ ;
$0 \leq ka [j + 1] - ka [j] \leq m$ , for $j = 0, 1, \dots, n - 1$ .

bl: Type: array<System..::..Double>[]()[][]
An array of size [ $n + m$ ]
On entry: $l$ , the lower bounds for all the variables and general constraints, in the following order. The first n elements of bl must contain the bounds on the variables $x$ , and the next m elements the bounds for the general linear constraints $A x$ (or slacks $s$ ) and the free row (if any). To specify a nonexistent lower bound (i.e., $l_{j} = - \infty$ ), set $bl [j] \leq - bigbnd$ , where $bigbnd$ is the value of the optional parameter Infinite Bound Size ( $default value = 10^{20}$ ). To specify the $j$ th constraint as an equality, set $bl [j] = bu [j] = β$ , say, where $|β| < bigbnd$ . Note that the lower bound corresponding to the free row must be set to $- \infty$ and stored in $bl [n + iobj - 1]$ .

Constraint: if $iobj > 0$ , $bl [n + iobj - 1] \leq - bigbnd$
(See also the description for bu.)

bu

Type: array<System..::..Double>[]()[][]

An array of size [

n + m

]

On entry:

u

, the upper bounds for all the variables and general constraints, in the following order. The first n elements of bl must contain the bounds on the variables

x

, and the next m elements the bounds for the general linear constraints

A x

(or slacks

s

) and the free row (if any). To specify a nonexistent upper bound (i.e.,

u_{j} = + \infty

), set

bu [j] \geq bigbnd

. Note that the upper bound corresponding to the free row must be set to

+ \infty

and stored in

bu [n + iobj - 1]

On exit: used as internal workspace prior to being restored and hence is unchanged.

Constraints:

if $iobj > 0$ , $bu [n + iobj - 1] \geq bigbnd$ ;
$bl [j] \leq bu [j]$ , for $j = 0, 1, \dots, n + m - 1$ ;
if $bl [j] = bu [j] = β$ , $|β| < bigbnd$ .

start

Type: System..::..String

On entry: indicates how a starting basis is to be obtained.

$start = "C"$: An internal crash procedure will be used to choose an initial basis matrix $B$ .
$start = "W"$: A basis is already defined in istate (probably from a previous call).

Constraint:

start = "C"

"W"

names

Type: array<System..::..String>[]()[][]

An array of size [

5

]

On entry: a set of names associated with the so-called MPSX form of the problem.

$names [0]$: Must contain the name for the problem (or be blank).
$names [1]$: Must contain the name for the free row (or be blank).
$names [2]$: Must contain the name for the constraint right-hand side (or be blank).
$names [3]$: Must contain the name for the ranges (or be blank).
$names [4]$: Must contain the name for the bounds (or be blank).

(These names are used in the monitoring file output; see [Description of Monitoring Information].)

nname

Type: System..::..Int32

On entry: the number of column (i.e., variable) and row names supplied in the array names.

$nname = 1$: There are no names. Default names will be used in the printed output.
$nname = n + m$: All names must be supplied.

Constraint:

nname = 1

n + m

crname: Type: array<System..::..String>[]()[][]
An array of size [nname]
On entry: the optional column and row names.
If $nname = 1$ , crname is not referenced and the printed output will use default names for the columns and rows.

If $nname = n + m$ , the first n elements must contain the names for the columns and the next m elements must contain the names for the rows. Note that the name for the free row (if any) must be stored in $crname [n + iobj]$ .

ns: Type: System..::..Int32%
On entry: $n_{S}$ , the number of superbasics. For QP problems, ns need not be specified if $start = "C"$ , but must retain its value from a previous call when $start = "W"$ . For FP and LP problems, ns need not be initialized.
On exit: the final number of superbasics. This will be zero for FP and LP problems.

xs: Type: array<System..::..Double>[]()[][]
An array of size [ $n + m$ ]
On entry: the initial values of the variables and slacks $(x, s)$ . (See the description for istate.)
On exit: $xs [i - 1]$ contains the final value of $x_{i}$ , for $i = 1, 2, \dots, n$ .

intvar: Type: array<System..::..Int32>[]()[][]
An array of size [lintvr]
On entry: specifies which components of the solution vector $x$ are constrained to be integer. Specifically, if $k$ elements of $x$ are required to take integer values then $intvar [i] = l_{i}$ , for $i = 1, 2, \dots, k$ , where $l_{i}$ is the integer index such that $x_{l_{i}}$ is integer. If $k < lintvr$ then $intvar [k + 1]$ must be set to $- 1$ to signal the end of the integer variable indices.
The order in which the indices of those components of $x$ required to be integer is presented determines the order in which the sub-problems are treated and solved. As such it can be a powerful tool to assist the method in achieving a solution efficiently. The general advice is to enter the important integer variables in the model early in intvar; secondary or less important variables should be entered near the end of the list. However some experimentation might be required to find the optimal order.

lintvr: Type: System..::..Int32
On entry: $k$ , the number of components of $x$ required to be integer. If $k = 0$ , then lintvr must be set to $1$ and $intvar [0]$ set to $- 1$ .

mdepth: Type: System..::..Int32
On entry: specifies the maximum depth the tree of sub-problems may be developed.
Suggested value: $mdepth = 2 \times n + 20$ .

Constraint: $mdepth > 0$ .

istate

Type: array<System..::..Int32>[]()[][]

An array of size [

n + m

]

On entry: if

start = "C"

, the first n elements of istate and xs must specify the initial states and values, respectively, of the variables

x

. (The slacks

s

need not be initialized.) An internal crash procedure is then used to select an initial basis matrix

B

. The initial basis matrix will be triangular (neglecting certain small elements in each column). It is chosen from various rows and columns of columns of

(A - I)

. Possible values for

istate [j]

are as follows:

$istate [j]$	State of $xs [j]$ during crash procedure
0 or $1$	Eligible for the basis
2	Ignored
3	Eligible for the basis (given preference over $0$ or $1$ )
4 or $5$	Ignored

If nothing special is known about the problem, or there is no wish to provide special information, you may set

istate [j] = 0

and

xs [j] = 0.0

, for

j = 1, 2, \dots, n

. All variables will then be eligible for the initial basis. Less trivially, to say that the

j

th variable will probably be equal to one of its bounds, set

istate [j] = 4

and

xs [j] = bl [j]

istate [j] = 5

and

xs [j] = bu [j]

as appropriate.

Following the crash procedure, variables for which

istate [j] = 2

are made superbasic. Other variables not selected for the basis are then made nonbasic at the value

xs [j]

bl [j] \leq xs [j] \leq bu [j]

, or at the value

bl [j]

bu [j]

closest to

xs [j]

start = "W"

, istate and xs must specify the initial states and values, respectively, of the variables and slacks

(x, s)

. If h02ce has been called previously with the same values of n and m, istate already contains satisfactory information.

Constraints:

if $start = "C"$ , $0 \leq istate [j] \leq 5$ , for $j = 0, 1, \dots, n - 1$ ;
if $start = "W"$ , $0 \leq istate [j] \leq 3$ , for $j = 0, 1, \dots, n + m - 1$ .

On exit: the final states of the variables and slacks

(x, s)

from the solution of the last sub-problem tackled. The significance of each possible value of

istate [j]

is as follows:

$istate [j]$	State of variable $j$	Normal value of $xs [j]$
$0$	Nonbasic	$bl [j]$
$1$	Nonbasic	$bu [j]$
$2$	Superbasic	Between $bl [j]$ and $bu [j]$
$3$	Basic	Between $bl [j]$ and $bu [j]$

Ninf = 0

(see [Description of the Printed Output]), basic and superbasic variables may be outside their bounds by as much as the value of the optional parameter Feasibility Tolerance (

default value = \max (10^{- 6}, \sqrt{ε})

, where

ε

is the machine precision). Note that unless the optional parameter

Scale Option = 0

(

default value = 2

) is specified, the Feasibility Tolerance applies to the variables of the scaled problem. In this case, the variables of the original problem may be as much as

0.1

outside their bounds, but this is unlikely unless the problem is very badly scaled.

Very occasionally some nonbasic variables may be outside their bounds by as much as the Feasibility Tolerance, and there may be some nonbasic variables for which

xs [j]

lies strictly between its bounds.

Ninf > 0

, some basic and superbasic variables may be outside their bounds by an arbitrary amount (bounded by Sinf (see [Description of the Printed Output]) if

Scale Option = 0

miniz: Type: System..::..Int32%
On exit: the minimum value of leniz required to start solving the problem. If $ifail = 14$ , h02ce may be called again with leniz suitably larger than miniz. (The bigger the better, since it is not certain how much workspace the basis factors need.)

minz: Type: System..::..Int32%
On exit: the minimum value of lenz required to start solving the problem. If $ifail = 15$ , h02ce may be called again with lenz suitably larger than minz. (The bigger the better, since it is not certain how much workspace the basis factors need.)

obj: Type: System..::..Double%
On exit: the value of the objective function.
If $Ninf = 0$ , obj includes the quadratic objective term $\frac{1}{2} x^{T} H x$ (if any).

If $Ninf > 0$ , obj is just the linear objective term $c^{T} x$ (if any). For FP problems, obj is set to zero.

clamda: Type: array<System..::..Double>[]()[][]
An array of size [ $n + m$ ]
On exit: a set of Lagrange-multipliers for the bounds on the variables and the general constraints. More precisely, the first n elements contain the multipliers (reduced costs) for the bounds on the variables, and the next m elements contain the multipliers (shadow prices) for the general linear constraints.

strtgy

Type: System..::..Int32

On entry: defines the branching strategy adopted by the method.

$strtgy = 0$: Each sub-problem first explored imposes a tighter upper bound on the component of $x$ .
$strtgy = 1$: Each sub-problem first explored imposes a tighter lower bound on the component of $x$ .
$strtgy = 2$: Each branch explored imposes a tighter upper bound on the component of $x$ if its fractional part is less than $0.5$ , otherwise it imposes a tighter lower bound.
$strtgy = 3$: Random choice is made between first exploring a tighter lower bound or a tighter upper bound sub-problem.

Constraint:

strtgy = 0

1

2

3

iz: Type: array<System..::..Int32>[]()[][]
An array of size [leniz]
the dimension of the array iz.
On entry: the dimension of the array iz as declared in the (sub)program from which h02ce is called.

Constraint: $leniz \geq 1$ .

z: Type: array<System..::..Double>[]()[][]
An array of size [lenz]
the dimension of the array z.
On entry: the dimension of the array z as declared in the (sub)program from which h02ce is called.

Constraint: $lenz \geq 1$ .
The amounts of workspace provided (i.e., leniz and lenz) and required (i.e., miniz and minz) are (by default) output on the current advisory message unit (as defined by (X04ABF not in this release)). Since the minimum values of leniz and lenz required to start solving the problem are returned in miniz and minz, respectively, you may prefer to obtain appropriate values from the output of a preliminary run with leniz and lenz set to $1$ . (h02ce will then terminate with $ifail = 14$ .)

monit: Type: NagLibrary..::..H..::..H02CE_MONIT
To provide feed-back on the progress of the branch and bound process. Additionally monit provides, via its parameter halt, the ability to terminate the process. (You might choose to do this when an integer solution is found, rather than search for a better solution.) If you do not require any intermediate output then monit may be the dummy method (H02CEY not in this release).
A delegate of type H02CE_MONIT.

options: Type: NagLibrary..::..H..::..h02ceOptions
An Object of type H.h02ceOptions. Used to configure optional parameters to this method.

ifail: Type: System..::..Int32%
On exit: $ifail = 0$ unless the method detects an error or a warning has been flagged (see [Error Indicators and Warnings]).

Description

h02ce is designed to obtain integer solutions to a class of quadratic programming problems addressed by e04nk. Specifically it solves the following problem:

\begin{array}{l} \underset{x \in R^{n}}{minimize} f (x) subject to l \leq \{\begin{matrix} x \\ A x \end{matrix}\} \leq u, \end{array}

(1)

where

x = {(x_{1}, x_{2}, \dots, x_{n})}^{T}

is a set of variables (some of which may be required to be integer),

A

is an

m

n

matrix and the objective function

f (x)

may be specified in a variety of ways depending upon the particular problem to be solved. The optional parameter Maximize may be used to specify an alternative problem in which

f (x)

is maximized. The possible forms for

f (x)

are listed in Table 1, in which the prefixes LP and QP stand for ‘linear programming’ and ‘quadratic programming’ respectively,

c

is an

n

-element vector and

H

is the

n

n

second-derivative matrix

\nabla^{2} f (x)

(the Hessian matrix).

Problem type	Objective function $f (x)$	Hessian matrix $H$
LP	$c^{T} x$	Not applicable
QP	$c^{T} x + \frac{1}{2} x^{T} H x$	Symmetric positive semidefinite

Table 1

For LP and QP problems, the unique global minimum value of

f (x)

is found. For QP problems, you must also provide a method that computes

H x

for any given vector

x

. (

H

need not be stored explicitly.)

(It is not expected that the feasibility problem of e04nk would be relevant here.)

The method employs a ‘Branch and Bound’ technique to enforce the integer constraints. In this technique the problem is first solved without the integer constraints. If a variable is found to be non-integral when it is required to have an integer value then two additional problems are constructed. One bounds the variable above by the nearest integer value below the optimal value previously obtained. The second problem is formed by bounding the variable below by the nearest integer value above the optimal value. This process is continued until an integer solution is found. At this point you may elect to stop, or may continue to search for better integer solutions by examining any other sub-problems that remain to be explained.

In practice the method tries to compute an integer solution as quickly as possible using a depth-first approach, since this helps determine a realistic cut-off value. If we have a cut-off value, say the value of the function at this first integer solution, and any sub-problem,

W

say, has a solution value greater than this cut-off value, then subsequent sub-problems of

W

must have solutions greater than the value of the solution at

W

and therefore need not be computed. Thus a knowledge of a good cut-off value can result in fewer sub-problems being solved and thus speed up the operation of the method. (See the description of monit in [Parameters] for details of how you can supply your own cut-off value.)

Each sub-problem is solved using (E04NKA not in this release). You are referred to the method document for e04nk for details of the algorithm used.

References

Gill P E, Hammarling S, Murray W, Saunders M A and Wright M H (1986) Users' guide for LSSOL (Version 1.0) Report SOL 86-1 Department of Operations Research, Stanford University

Gill P E and Murray W (1978) Numerically stable methods for quadratic programming Math. Programming 14 349–372

Gill P E, Murray W, Saunders M A and Wright M H (1986) Some theoretical properties of an augmented Lagrangian merit function Report SOL 86–6R Department of Operations Research, Stanford University

Gill P E, Murray W, Saunders M A and Wright M H (1989) A practical anti-cycling procedure for linearly constrained optimization Math. Programming 45 437–474

Gill P E, Murray W, Saunders M A and Wright M H (1991) Inertia-controlling methods for general quadratic programming SIAM Rev. 33 1–36

Hock W and Schittkowski K (1981) Test Examples for Nonlinear Programming Codes. Lecture Notes in Economics and Mathematical Systems 187 Springer–Verlag

Lawson C L, Hanson R J, Kincaid D R and Krogh F T (1979) Basic linear algebra supbrograms for Fortran usage ACM Trans. Math. Software 5 308–325

Murtagh B A and Saunders M A (1983) MINOS 5.0 user's guide Report SOL 83-20 Department of Operations Research, Stanford University

Error Indicators and Warnings

Errors or warnings detected by the method:

$ifail = - 1$: Halted at your request.

$ifail = 0$: Successful exit.

$ifail = 1$: Input parameter error immediately detected.

$ifail = 2$: No integer solution found.

$ifail = 3$: mdepth is too small.

$ifail = 4$: The problem is unbounded (or badly scaled). The objective function is not bounded below in the feasible region.

$ifail = 5$: The problem is infeasible. The general constraints cannot all be satisfied simultaneously to within the value of the optional parameter Feasibility Tolerance ( $default value = \max (10^{- 6}, \sqrt{ε})$ , where $ε$ is the machine precision).

$ifail = 6$: Too many iterations. The value of the optional parameter Iteration Limit ( $default value = \max (50, 5 (n + m))$ ) is too small.

$ifail = 7$: The reduced Hessian matrix $z^{T} H Z$ (see [Definition of the Working Set and Search Direction]) exceeds its assigned dimension. The value of the optional parameter Superbasics Limit ( $default value = \min (n_{H} + 1, n)$ ) is too small.

$ifail = 8$: The Hessian matrix $H$ appears to be indefinite. Check that qphx has been coded correctly and that all relevant elements of $H x$ have been assigned their correct values.

$ifail = 9$: An input parameter is invalid for an internal call to e04nk.

$ifail = 10$: Numerical error in trying to satisfy the general constraints. The basis is very ill-conditioned.

$ifail = 11$: Not enough integer workspace for the basis factors. Increase leniz and rerun h02ce.

$ifail = 12$: Not enough real workspace for the basis factors. Increase lenz and rerun h02ce.

$ifail = 13$: The basis is singular after $15$ attempts to factorize it (adding slacks where necessary). Either the problem is badly scaled or the value of the optional parameter LU Factor Tolerance ( $default value = 100.0$ ) is too large.

$ifail = 14$: Not enough integer workspace to start solving the problem. Increase leniz to at least miniz and rerun h02ce.

$ifail = 15$: Not enough real workspace to start solving the problem. Increase lenz to at least minz and rerun h02ce.

$ifail = 16$: An internal error has occurred. Contact NAG with details of your program.

$ifail = -9000$: An error occured, see message report.
$ifail = -8000$: Negative dimension for array $〈value〉$
$ifail = -6000$: Invalid Parameters $〈value〉$

Accuracy

h02ce implements a numerically stable active-set strategy and returns solutions that are as accurate as the condition of the problem warrants on the machine.

Parallelism and Performance

None.

Further Comments

This section contains a description of the printed output.

Description of the Printed Output

This section describes the (default) intermediate printout and final printout produced by h02ce. The intermediate printout is a subset of the monitoring information produced by the method at every iteration (see [Description of Monitoring Information]). You can control the level of printed output (see the description of the optional parameter Print Level in [Description of the Optional Parameters]). Note that the intermediate printout and final printout are produced only if

Print Level \geq 10

(the default).

The following line of summary output (

< 80

characters) is produced at every iteration. In all cases, the values of the quantities printed are those in effect on completion of the given iteration.

Itn	is the iteration count.
Step	is the step taken along the computed search direction. If a constraint is added during the current iteration, Step will be the step to the nearest constraint. When the problem is of type LP, the step can be greater than one during the optimality phase.
Ninf	is the number of violated constraints (infeasibilities). This will be zero during the optimality phase.
Sinf/Objective	is the value of the current objective function. If $x$ is not feasible, Sinf gives a weighted sum of the magnitudes of constraint violations. If $x$ is feasible, Objective is the value of the objective function. The output line for the final iteration of the feasibility phase (i.e., the first iteration for which Ninf is zero) will give the value of the true objective at the first feasible point. During the optimality phase, the value of the objective function will be nonincreasing. During the feasibility phase, the number of constraint infeasibilities will not increase until either a feasible point is found, or the optimality of the multipliers implies that no feasible point exists. Once optimal multipliers are obtained, the number of infeasibilities can increase, but the sum of infeasibilities will either remain constant or be reduced until the minimum sum of infeasibilities is found.
Norm rg	is $‖d_{S}‖$ , the Euclidean norm of the reduced gradient (see [The Main Iteration]). During the optimality phase, this norm will be approximately zero after a unit step. For FP and LP problems, Norm rg is not printed.

The final printout includes a listing of the status of every variable and constraint.

The following describes the printout for each variable. A full stop (.) is printed for any numerical value that is zero.

Variable

gives the name of the variable. If

nname = 1

, a default name is assigned to the

j

th variable, for

j = 1, 2, \dots, n

. If

nname = n + m

, the name supplied in

crname [j]

is assigned to the

j

th variable.

State

gives the state of the variable (LL if nonbasic on its lower bound, UL if nonbasic on its upper bound, EQ if nonbasic and fixed, FR if nonbasic and strictly between its bounds, BS if basic and SBS if superbasic).

A key is sometimes printed before State to give some additional information about the state of a variable. Note that unless the optional parameter

Scale Option = 0

(

default value = 2

) is specified, the tests for assigning a key are applied to the variables of the scaled problem.

A	Alternative optimum possible. The variable is nonbasic, but its reduced gradient is essentially zero. This means that if the variable were allowed to start moving away from its bound, there would be no change to the objective function. The values of the other free variables might change, giving a genuine alternative solution. However, if there are any degenerate variables (labelled D), the actual change might prove to be zero, since one of them could encounter a bound immediately. In either case the values of the Lagrange-multipliers might also change.
D	Degenerate. The variable is basic or superbasic, but it is equal to (or very close to) one of its bounds.
I	Infeasible. The variable is basic or superbasic and is currently violating one of its bounds by more than the value of the optional parameter Feasibility Tolerance ( $default value = \max (10^{- 6}, \sqrt{ε})$ , where $ε$ is the machine precision).
N	Not precisely optimal. The variable is nonbasic or superbasic. If the value of the reduced gradient for the variable exceeds the value of the optional parameter Optimality Tolerance ( $default value = \max (10^{- 6}, \sqrt{ε})$ ), the solution would not be declared optimal because the reduced gradient for the variable would not be considered negligible.

Value

is the value of the variable at the final iterate.

Lower Bound

is the lower bound specified for the variable. None indicates that

bl [j] \leq - bigbnd

Upper Bound

is the upper bound specified for the variable. None indicates that

bu [j] \geq bigbnd

Lagr Mult

is the Lagrange-multiplier for the associated bound. This will be zero if State is FR. If

x

is optimal, the multiplier should be non-negative if State is LL, non-positive if State is UL, and zero if State is BS or SBS.

Residual

is the difference between the variable Value and the nearer of its (finite) bounds

bl [j]

and

bu [j]

. A blank entry indicates that the associated variable is not bounded (i.e.,

bl [j] \leq - bigbnd

and

bu [j] \geq bigbnd

The meaning of the printout for linear constraints is the same as that given above for variables, with ‘variable’ replaced by ‘constraint’,

n

replaced by

m

crname [j]

replaced by

crname [n + j]

bl [j]

and

bu [j]

are replaced by

bl [n + j]

and

bu [n + j]

respectively, and with the following change in the heading.

Constrnt

gives the name of the linear constraint.

Note that movement off a constraint (as opposed to a variable moving away from its bound) can be interpreted as allowing the entry in the Residual column to become positive.

Numerical values are output with a fixed number of digits; they are not guaranteed to be accurate to this precision.

Example

This example minimizes the quadratic function

f (x) = c^{T} x + \frac{1}{2} x^{T} H x

, where

c = {(- 200.0, - 2000.0, - 2000.0, - 2000.0, - 2000.0, 400.0, 400.0)}^{T}

H = (\begin{array}{r} 2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 2 & 2 & 0 & 0 & 0 \\ 0 & 0 & 2 & 2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 2 & 2 \\ 0 & 0 & 0 & 0 & 0 & 2 & 2 \end{array})

subject to the bounds

\begin{array}{l} 0 \leq x_{1} \leq 200 \\ 0 \leq x_{2} \leq 2500 \\ 400 \leq x_{3} \leq 800 \\ 100 \leq x_{4} \leq 700 \\ 0 \leq x_{5} \leq 1500 \\ 0 \leq x_{6} \\ 0 \leq x_{7} \end{array}

to the linear constraints

\begin{array}{r} x_{1} & + & x_{2} & + & x_{3} & + & x_{4} & + & x_{5} & + & x_{6} & + & x_{7} & = & 2000 \\ 0.15 x_{1} & + & 0.04 x_{2} & + & 0.02 x_{3} & + & 0.04 x_{4} & + & 0.02 x_{5} & + & 0.01 x_{6} & + & 0.03 x_{7} & \leq & 60 \\ 0.03 x_{1} & + & 0.05 x_{2} & + & 0.08 x_{3} & + & 0.02 x_{4} & + & 0.06 x_{5} & + & 0.01 x_{6} & \leq & 100 \\ 0.02 x_{1} & + & 0.04 x_{2} & + & 0.01 x_{3} & + & 0.02 x_{4} & + & 0.02 x_{5} & \leq & 40 \\ 0.02 x_{1} & + & 0.03 x_{2} & + & 0.01 x_{5} & \leq & 30 \\ 1500 & \leq & 0.70 x_{1} & + & 0.75 x_{2} & + & 0.80 x_{3} & + & 0.75 x_{4} & + & 0.80 x_{5} & + & 0.97 x_{6} \\ 250 & \leq & 0.02 x_{1} & + & 0.06 x_{2} & + & 0.08 x_{3} & + & 0.12 x_{4} & + & 0.02 x_{5} & + & 0.01 x_{6} & + & 0.97 x_{7} & \leq & 300 \end{array}

and the variables

x_{2}

x_{3}

x_{4}

x_{5}

x_{6}

x_{7}

, are constrained to be integer.

The initial point, which is infeasible, is

x_{0} = {(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)}^{T} .

The optimal solution (to five figures) is

x^{*} = {(0.0, 355.0, 645.0, 164.0, 410.0, 275.0, 151.0)}^{T} .

One bound constraint and one linear constraint are active at the solution. Note that the Hessian matrix

H

is positive semidefinite.

Example program (C#): h02cee.cs

Example program data: h02cee.d

Example program results: h02cee.r

Algorithmic Details

This section contains a detailed description of the method used by h02ce.

Overview

h02ce employs a Branch and Bound technique (see [Description]) based on an inertia-controlling method to solve the sub-problems that maintains a Cholesky factorization of the reduced Hessian (see below). The method is similar to that of Gill and Murray (1978), and is described in detail by Gill et al. (1991). Here we briefly summarise the main features of the method. Where possible, explicit reference is made to the names of variables that are parameters of the method or appear in the printed output.

The method used has two distinct phases: finding an initial feasible point by minimizing the sum of infeasibilities (the feasibility phase), and minimizing the quadratic objective function within the feasible region (the optimality phase). The computations in both phases are performed by the same methods. The two-phase nature of the algorithm is reflected by changing the function being minimized from the sum of infeasibilities (the printed quantity Sinf; see [Description of Monitoring Information]) to the quadratic objective function (the printed quantity Objective; see [Description of Monitoring Information]).

In general, an iterative process is required to solve a quadratic program. Given an iterate

(x, s)

in both the original variables

x

and the slack variables

s

, a new iterate

(\bar{x}, \bar{s})

is defined by

(\begin{matrix} \bar{x} \\ \bar{s} \end{matrix}) = (\begin{matrix} x \\ s \end{matrix}) + α p,

(2)

where the step length

α

is a non-negative scalar (the printed quantity Step; see [Description of Monitoring Information]), and

p

is called the search direction. (For simplicity, we shall consider a typical iteration and avoid reference to the index of the iteration.) Once an iterate is feasible (i.e., satisfies the constraints), all subsequent iterates remain feasible.

Definition of the Working Set and Search Direction

At each iterate

(x, s)

, a working set of constraints is defined to be a linearly independent subset of the constraints that are satisfied ‘exactly’ (to within the value of the optional parameter Feasibility Tolerance; see [Description of the Optional Parameters]). The working set is the current prediction of the constraints that hold with equality at a solution of the LP or QP problem. Let

m_{W}

denote the number of constraints in the working set (including bounds), and let

W

denote the associated

m_{W}

(n + m)

working set matrix consisting of the

m_{W}

gradients of the working set constraints.

The search direction is defined so that constraints in the working set remain unaltered for any value of the step length. It follows that

p

must satisfy the identity

W p = 0 .

(3)

This characterisation allows

p

to be computed using any

n

n_{z}

full-rank matrix

Z

that spans the null space of

W

. (Thus,

n_{z} = n - m_{W}

and

W Z = 0

.) The null space matrix

Z

is defined from a sparse

L U

factorization of part of

W

(see (6) and (7) below). The direction

p

will satisfy (3) if

p = Z p_{Z},

(4)

where

p_{z}

is any

n_{z}

-vector.

The working set contains the constraints

A x - s = 0

and a subset of the upper and lower bounds on the variables

(x, s)

. Since the gradient of a bound constraint

x_{j} \geq l_{j}

x_{j} \leq u_{j}

is a vector of all zeros except for

\pm 1

in position

j

, it follows that the working set matrix contains the rows of

(A - I)

and the unit rows associated with the upper and lower bounds in the working set.

The working set matrix

W

can be represented in terms of a certain column partition of the matrix

(A - I)

. As in [Description] we partition the constraints

A x - s = 0

so that

B x_{B} + S x_{S} + N x_{N} = 0,

(5)

where

B

is a square nonsingular basis and

x_{B}

x_{S}

and

x_{n}

are the basic, superbasic and nonbasic variables respectively. The nonbasic variables are equal to their upper or lower bounds at

(x, s)

, and the superbasic variables are independent variables that are chosen to improve the value of the current objective function. The number of superbasic variables is

n_{S}

(the printed quantity Ns; see [Description of Monitoring Information]). Given values of

x_{N}

and

x_{S}

, the basic variables

x_{B}

are adjusted so that

(x, s)

satisfies (5).

P

is a permutation matrix such that

(A - I) P = (B S N)

, then the working set matrix

W

satisfies

W P = (\begin{matrix} B & S & N \\ 0 & 0 & I_{N} \end{matrix}),

(6)

where

I_{N}

is the identity matrix with the same number of columns as

N

The null space matrix

Z

is defined from a sparse

L U

factorization of part of

W

. In particular,

z

is maintained in ‘reduced gradient’ form, using the LUSOL package (see Gill et al. (1986)) to maintain sparse

L U

factors of the basis matrix

B

that alters as the working set

W

changes. Given the permutation

P

, the null space basis is given by

Z = P (\begin{matrix} - B^{- 1} S \\ I \\ 0 \end{matrix}) .

(7)

This matrix is used only as an operator, i.e., it is never computed explicitly. Products of the form

Z v

and

Z^{T} g

are obtained by solving with

B

B^{T}

. This choice of

Z

implies that

n_{Z}

, the number of ‘degrees of freedom’ at

(x, s)

, is the same as

n_{S}

, the number of superbasic variables.

Let

g_{Z}

and

H_{Z}

denote the reduced gradient and reduced Hessian of the objective function:

g_{Z} = Z^{T} g and H_{Z} = Z^{T} H Z,

(8)

where

g

is the objective gradient at

(x, s)

. Roughly speaking,

g_{z}

and

H_{z}

describe the first and second derivatives of an

n_{S}

-dimensional unconstrained problem for the calculation of

p_{Z}

. (The condition estimator of

H_{Z}

is the quantity Cond Hz in the monitoring file output; see [Description of Monitoring Information].)

At each iteration, an upper triangular factor

R

is available such that

H_{Z} = R^{T} R

. Normally,

R

is computed from

R^{T} R = Z^{T} H Z

at the start of the optimality phase and then updated as the QP working set changes. For efficiency, the dimension of

R

should not be excessive (say,

n_{S} \leq 1000

). This is guaranteed if the number of nonlinear variables is ‘moderate’.

If the QP problem contains linear variables,

H

is positive semidefinite and

R

may be singular with at least one zero diagonal element. However, an inertia-controlling strategy is used to ensure that only the last diagonal element of

R

can be zero. (See Gill et al. (1991) for a discussion of a similar strategy for indefinite quadratic programming.)

If the initial

R

is singular, enough variables are fixed at their current value to give a nonsingular

R

. This is equivalent to including temporary bound constraints in the working set. Thereafter,

R

can become singular only when a constraint is deleted from the working set (in which case no further constraints are deleted until

R

becomes nonsingular).

The Main Iteration

If the reduced gradient is zero,

(x, s)

is a constrained stationary point on the working set. During the feasibility phase, the reduced gradient will usually be zero only at a vertex (although it may be zero elsewhere in the presence of constraint dependencies). During the optimality phase, a zero reduced gradient implies that

x

minimizes the quadratic objective function when the constraints in the working set are treated as equalities. At a constrained stationary point, Lagrange-multipliers

λ

are defined from the equations

W^{T} λ = g (x) .

(9)

A Lagrange-multiplier

λ_{j}

corresponding to an inequality constraint in the working set is said to be optimal if

λ_{j} \leq σ

when the associated constraint is at its upper bound, or if

λ_{j} \geq - σ

when the associated constraint is at its lower bound, where

σ

depends on the value of the optional parameter Optimality Tolerance (see [Description of the Optional Parameters]). If a multiplier is nonoptimal, the objective function (either the true objective or the sum of infeasibilities) can be reduced by continuing the minimization with the corresponding constraint excluded from the working set. (This step is sometimes referred to as ‘deleting’ a constraint from the working set.) If optimal multipliers occur during the feasibility phase but the sum of infeasibilities is nonzero, there is no feasible point and the method terminates immediately with

ifail = 3

(see [Error Indicators and Warnings]).

The special form (6) of the working set allows the multiplier vector

λ

, the solution of (9), to be written in terms of the vector

d = (\begin{matrix} g \\ 0 \end{matrix}) - {(\begin{matrix} A & - I \end{matrix})}^{T} π = (\begin{matrix} g - A^{T} π \\ π \end{matrix}),

(10)

where

π

satisfies the equations

B^{T} π = g_{B}

, and

g_{B}

denotes the basic elements of

g

. The elements of

π

are the Lagrange-multipliers

λ_{j}

associated with the equality constraints

A x - s = 0

. The vector

d_{N}

of nonbasic elements of

d

consists of the Lagrange-multipliers

λ_{j}

associated with the upper and lower bound constraints in the working set. The vector

d_{S}

of superbasic elements of

d

is the reduced gradient

g_{z}

in (8). The vector

d_{B}

of basic elements of

d

is zero, by construction. (The Euclidean norm of

d_{S}

and the final values of

d_{S}

g

and

π

are the quantities Norm rg, Reduced Gradnt, Obj Gradient and Dual Activity in the monitoring file output; see [Description of Monitoring Information].)

If the reduced gradient is not zero, Lagrange-multipliers need not be computed and the search direction is given by

p = Z p_{z}

(see (7) and (11)). The step length is chosen to maintain feasibility with respect to the satisfied constraints.

There are two possible choices for

p_{z}

, depending on whether or not

H_{z}

is singular. If

H_{z}

is nonsingular,

R

is nonsingular and

p_{z}

in (4) is computed from the equations

R^{T} R p_{Z} = - g_{Z},

(11)

where

g_{z}

is the reduced gradient at

x

. In this case,

(x, s) + p

is the minimizer of the objective function subject to the working set constraints being treated as equalities. If

(x, s) + p

is feasible,

α

is defined to be unity. In this case, the reduced gradient at

(\bar{x}, \bar{s})

will be zero, and Lagrange-multipliers are computed at the next iteration. Otherwise,

α

is set to

α_{m}

, the step to the ‘nearest’ constraint along

p

. This constraint is added to the working set at the next iteration.

H_{z}

is singular, then

R

must also be singular, and an inertia-controlling strategy is used to ensure that only the last diagonal element of

R

is zero. (See Gill et al. (1991) for a discussion of a similar strategy for indefinite quadratic programming.) In this case,

p_{z}

satisfies

p_{Z}^{T} H_{Z} p_{Z} = 0 and g_{Z}^{T} p_{Z} \leq 0,

(12)

which allows the objective function to be reduced by any step of the form

(x, s) + α p

, where

α > 0

. The vector

p = Z p_{Z}

is a direction of unbounded descent for the QP problem in the sense that the QP objective is linear and decreases without bound along

p

. If no finite step of the form

(x, s) + α p

(where

α > 0

) reaches a constraint not in the working set, the QP problem is unbounded and the method terminates immediately with

ifail = 2

(see [Error Indicators and Warnings]). Otherwise,

α

is defined as the maximum feasible step along

p

and a constraint active at

(x, s) + α p

is added to the working set for the next iteration.

Miscellaneous

If the basis matrix is not chosen carefully, the condition of the null space matrix

Z

in (7) could be arbitrarily high. To guard against this, the method implements a ‘basis repair’ feature in which the LUSOL package (see Gill et al. (1986)) is used to compute the rectangular factorization

{(\begin{matrix} B & S \end{matrix})}^{T} = L U,

(13)

returning just the permutation

P

that makes

P L P^{T}

unit lower triangular. The pivot tolerance is set to require

{|P L P^{T}|}_{i j} \leq 2

, and the permutation is used to define

P

in (6). It can be shown that

‖Z‖

is likely to be little more than unity. Hence,

z

should be well-conditioned regardless of the condition of

W

. This feature is applied at the beginning of the optimality phase if a potential

B - S

ordering is known.

The EXPAND procedure (see Gill et al. (1989)) is used to reduce the possibility of cycling at a point where the active constraints are nearly linearly dependent. Although there is no absolute guarantee that cycling will not occur, the probability of cycling is extremely small (see Gill et al. (1986)). The main feature of EXPAND is that the feasibility tolerance is increased at the start of every iteration. This allows a positive step to be taken at every iteration, perhaps at the expense of violating the bounds on

(x, s)

by a small amount.

Suppose that the value of the optional parameter Feasibility Tolerance is

δ

. Over a period of

K

iterations (where

K

is the value of the optional parameter Expand Frequency; see [Description of the Optional Parameters]), the feasibility tolerance actually used by h02ce (i.e., the working feasibility tolerance) increases from

0.5 δ

δ

(in steps of

0.5 δ / K

At certain stages the following ‘resetting procedure’ is used to remove small constraint infeasibilities. First, all nonbasic variables are moved exactly onto their bounds. A count is kept of the number of nontrivial adjustments made. If the count is nonzero, the basic variables are recomputed. Finally, the working feasibility tolerance is reinitialized to

0.5 δ

If a problem requires more than

K

iterations, the resetting procedure is invoked and a new cycle of iterations is started. (The decision to resume the feasibility phase or optimality phase is based on comparing any constraint infeasibilities with

δ

The resetting procedure is also invoked when h02ce reaches an apparently optimal, infeasible or unbounded solution, unless this situation has already occurred twice. If any nontrivial adjustments are made, iterations are continued.

The EXPAND procedure not only allows a positive step to be taken at every iteration, but also provides a potential choice of constraints to be added to the working set. All constraints at a distance

α

(where

α \leq α_{m}

) along

p

from the current point are then viewed as acceptable candidates for inclusion in the working set. The constraint whose normal makes the largest angle with the search direction is added to the working set. This strategy helps keep the basis matrix

B

well-conditioned.

Description of Monitoring Information

This section describes the intermediate printout and final printout which constitutes the monitoring information produced by h02ce. (See also the description of the optional parameters Monitoring File and Print Level in [Description of the Optional Parameters].) You can control the level of printed output.

When

Print Level = 5

\geq 10

and

Monitoring File \geq 0

, the following line of intermediate printout (

< 120

characters) is produced at every iteration on the unit number specified by Monitoring File. Unless stated otherwise, the values of the quantities printed are those in effect on completion of the given iteration.

Itn	is the iteration count.
pp	is the partial price indicator. The variable selected by the last pricing operation came from the ppth partition of $A$ and $- I$ . Note that pp is reset to zero whenever the basis is refactorized.
dj	is the value of the reduced gradient (or reduced cost) for the variable selected by the pricing operation at the start of the current iteration.
+S	is the variable selected by the pricing operation to be added to the superbasic set.
-S	is the variable chosen to leave the superbasic set.
-B	is the variable removed from the basis (if any) to become nonbasic.
-B	is the variable chosen to leave the set of basics (if any) in a special basic $\leftrightarrow$ superbasic swap. The entry under -S has become basic if this entry is nonzero, and nonbasic otherwise. The swap is done to ensure that there are no superbasic slacks.
Step	is the value of the step length $α$ taken along the computed search direction $p$ . The variables $x$ have been changed to $x + α p$ . If a variable is made superbasic during the current iteration (i.e., +S is positive), Step will be the step to the nearest bound. During the optimality phase, the step can be greater than unity only if the reduced Hessian is not positive definite.
Pivot	is the $r$ th element of a vector $y$ satisfying $B y = a_{q}$ whenever $a_{q}$ (the $q$ th column of the constraint matrix $(A - I)$ ) replaces the $r$ th column of the basis matrix $B$ . Wherever possible, Step is chosen so as to avoid extremely small values of Pivot (since they may cause the basis to be nearly singular). In extreme cases, it may be necessary to increase the value of the optional parameter Pivot Tolerance ( $default value = ε^{0.67}$ , where $ε$ is the machine precision) to exclude very small elements of $y$ from consideration during the computation of Step.
Ninf	is the number of violated constraints (infeasibilities). This will be zero during the optimality phase.
Sinf/Objective	is the value of the current objective function. If $x$ is not feasible, Sinf gives a weighted sum of the magnitudes of constraint violations. If $x$ is feasible, Objective is the value of the objective function. The output line for the final iteration of the feasibility phase (i.e., the first iteration for which Ninf is zero) will give the value of the true objective at the first feasible point. During the optimality phase, the value of the objective function will be nonincreasing. During the feasibility phase, the number of constraint infeasibilities will not increase until either a feasible point is found, or the optimality of the multipliers implies that no feasible point exists. Once optimal multipliers are obtained, the number of infeasibilities can increase, but the sum of infeasibilities will either remain constant or be reduced until the minimum sum of infeasibilities is found.
L	is the number of nonzeros in the basis factor $L$ . Immediately after a basis factorization $B = L U$ , this is lenL, the number of subdiagonal elements in the columns of a lower triangular matrix. Further nonzeros are added to L when various columns of $B$ are later replaced. (Thus, L increases monotonically.)
U	is the number of nonzeros in the basis factor $U$ . Immediately after a basis factorization, this is lenU, the number of diagonal and superdiagonal elements in the rows of an upper triangular matrix. As columns of $B$ are replaced, the matrix $U$ is maintained explicitly (in sparse form). The value of U may fluctuate up or down; in general, it will tend to increase.
Ncp	is the number of compressions required to recover workspace in the data structure for $U$ . This includes the number of compressions needed during the previous basis factorization. Normally, Ncp should increase very slowly. If it does not, increase lenz by at least $L + U$ and rerun h02ce (possibly using $start = "W"$ ; see [Parameters]).
Norm rg	is $‖d_{S}‖$ , the Euclidean norm of the reduced gradient (see [The Main Iteration]). During the optimality phase, this norm will be approximately zero after a unit step. For FP and LP problems, Norm rg is not printed.
Ns	is the current number of superbasic variables. For FP and LP problems, Ns is not printed.
Cond Hz	is a lower bound on the condition number of the reduced Hessian (see [Definition of the Working Set and Search Direction]). The larger this number, the more difficult the problem. For FP and LP problems, Cond Hz is not printed.

When

Print Level \geq 20

and

Monitoring File \geq 0

, the following lines of intermediate printout (

< 120

characters) are produced on the unit number specified by Monitoring File whenever the matrix

B

B_{S} = {(\begin{matrix} B & S \end{matrix})}^{T}

is factorized. Gaussian elimination is used to compute an

L U

factorization of

B

B_{S}

, where

P L P^{T}

is a lower triangular matrix and

P U Q

is an upper triangular matrix for some permutation matrices

P

and

Q

. The factorization is stabilized in the manner described under the LU Factor Tolerance (

default value = 100.0

; see [Description of the Optional Parameters]).

Factorize

is the factorization count.

Demand

is a code giving the reason for the present factorization as follows:

Code	Meaning
$0$	First $L U$ factorization.
$1$	Number of updates reached the value of the optional parameter Factorization Frequency ( $default value = 100$ ).
$2$	Excessive nonzeros in updated factors.
$7$	Not enough storage to update factors.
$10$	Row residuals too large (see the description for the optional parameter Check Frequency).
$11$	Ill-conditioning has caused inconsistent results.

Iteration

is the iteration count.

Nonlinear

is the number of nonlinear variables in

B

(not printed if

B_{S}

is factorized).

Linear

is the number of linear variables in

B

(not printed if

B_{S}

is factorized).

Slacks

is the number of slack variables in

B

(not printed if

B_{S}

is factorized).

Elems

is the number of nonzeros in

B

(not printed if

B_{S}

is factorized).

Density

is the percentage nonzero density of

B

(not printed if

B_{S}

is factorized). More precisely,

Density = 100 \times Elems / {(Nonlinear + Linear + Slacks)}^{2}

Compressns

is the number of times the data structure holding the partially factorized matrix needed to be compressed, in order to recover unused workspace. Ideally, it should be zero. If it is more than

3

4

, increase leniz and lenz and rerun h02ce (possibly using

start = "W"

; see [Parameters]).

Merit

is the average Markowitz merit count for the elements chosen to be the diagonals of

P U Q

. Each merit count is defined to be

(c - 1) (r - 1)

, where

c

and

r

are the number of nonzeros in the column and row containing the element at the time it is selected to be the next diagonal. Merit is the average of m such quantities. It gives an indication of how much work was required to preserve sparsity during the factorization.

lenL

is the number of nonzeros in

L

lenU

is the number of nonzeros in

U

Increase

is the percentage increase in the number of nonzeros in

L

and

U

relative to the number of nonzeros in

B

. More precisely,

Increase = 100 \times (lenL + lenU - Elems) / Elems

is the number of rows in the problem. Note that

m = Ut + Lt + bp

is the number of triangular rows of

B

at the top of

U

is the number of columns remaining when the density of the basis matrix being factorized reached

0.3

Lmax

is the maximum subdiagonal element in the columns of

L

(not printed if

B_{S}

is factorized). This will not exceed the value of the LU Factor Tolerance.

Bmax

is the maximum nonzero element in

B

(not printed if

B_{S}

is factorized).

BSmax

is the maximum nonzero element in

B_{S}

(not printed if

B

is factorized).

Umax

is the maximum nonzero element in

U

, excluding elements of

B

that remain in

U

unchanged. (For example, if a slack variable is in the basis, the corresponding row of

B

will become a row of

U

without modification. Elements in such rows will not contribute to Umax. If the basis is strictly triangular, none of the elements of

B

will contribute, and Umax will be zero.)

Ideally, Umax should not be significantly larger than Bmax. If it is several orders of magnitude larger, it may be advisable to reset the LU Factor Tolerance to a value near

1.0

Umax is not printed if

B_{S}

is factorized.

Umin

is the magnitude of the smallest diagonal element of

P U Q

(not printed if

B_{S}

is factorized).

Growth

is the value of the ratio

Umax / Bmax

, which should not be too large.

Providing Lmax is not large (say

< 10.0

), the ratio

\max (Bmax, Umax) / Umin

is an estimate of the condition number of

B

. If this number is extremely large, the basis is nearly singular and some numerical difficulties could occur in subsequent computations. (However, an effort is made to avoid near singularity by using slacks to replace columns of

B

that would have made Umin extremely small, and the modified basis is refactorized.)

Growth is not printed if

B_{S}

is factorized.

is the number of triangular columns of

B

at the beginning of

L

is the size of the ‘bump’ or block to be factorized nontrivially after the triangular rows and columns have been removed.

is the number of columns remaining when the density of the basis matrix being factorized reached

0.6

When

Print Level \geq 20

and

Monitoring File \geq 0

, the following lines of intermediate printout (

< 80

characters) are produced on the unit number specified by Monitoring File whenever

start = "C"

(see [Parameters]). They refer to the number of columns selected by the crash procedure during each of several passes through

A

, whilst searching for a triangular basis matrix.

Slacks	is the number of slacks selected initially.
Free cols	is the number of free columns in the basis.
Preferred	is the number of ‘preferred’ columns in the basis (i.e., $istate [j] = 3$ for some $j \leq n$ ).
Unit	is the number of unit columns in the basis.
Double	is the number of double columns in the basis.
Triangle	is the number of triangular columns in the basis.
Pad	is the number of slacks used to pad the basis.

When

Print Level \geq 20

and

Monitoring File \geq 0

, the following lines of intermediate printout (

< 80

characters) are produced on the unit number specified by Monitoring File. They refer to the elements of the names array (see [Parameters]).

Name	gives the name for the problem (blank if none).
Objective	gives the name of the free row for the problem (blank if none).
RHS	gives the name of the constraint right-hand side for the problem (blank if none).
Ranges	gives the name of the ranges for the problem (blank if none).
Bounds	gives the name of the bounds for the problem (blank if none).

When

Print Level = 1

\geq 10

and

Monitoring File \geq 0

, the following lines of final printout (

< 120

characters) are produced on the unit number specified by Monitoring File.

Let

a_{j}

denote the

j

th column of

A

, for

j = 1, 2, \dots, n

. The following describes the printout for each column (or variable). A full stop (.) is printed for any numerical value that is zero.

Number

is the column number

j

. (This is used internally to refer to

x_{j}

in the intermediate output.)

Column

gives the name of

x_{j}

State

A key is sometimes printed before State to give some additional information about the state of

x_{j}

. Note that unless the optional parameter

Scale Option = 0

(

default value = 2

) is specified, the tests for assigning a key are applied to the variables of the scaled problem.

A	Alternative optimum possible. The variable is nonbasic, but its reduced gradient is essentially zero. This means that if the variable were allowed to start moving away from its bound, there would be no change to the objective function. The values of the other free variables might change, giving a genuine alternative solution. However, if there are any degenerate variables (labelled D), the actual change might prove to be zero, since one of them could encounter a bound immediately. In either case the values of the Lagrange-multipliers might also change.
D	Degenerate. The variable is basic or superbasic, but it is equal to (or very close to) one of its bounds.
I	Infeasible. The variable is basic or superbasic and is currently violating one of its bounds by more than the value of the optional parameter Feasibility Tolerance ( $default value = \max (10^{- 6}, \sqrt{ε})$ , where $ε$ is the machine precision).
N	Not precisely optimal. The variable is nonbasic or superbasic. If the value of the reduced gradient for the variable exceeds the value of the optional parameter Optimality Tolerance ( $default value = \max (10^{- 6}, \sqrt{ε})$ ), the solution would not be declared optimal because the reduced gradient for the variable would not be considered negligible.

Activity

is the value of

x_{j}

at the final iterate.

Obj Gradient

is the value of

g_{j}

at the final iterate. For FP problems,

g_{j}

is set to zero.

Lower Bound

is the lower bound specified for the variable. None indicates that

bl [j] \leq - bigbnd

Upper Bound

is the upper bound specified for the variable. None indicates that

bu [j] \geq bigbnd

Reduced Gradnt

is the value of

d_{j}

at the final iterate (see [The Main Iteration]). For FP problems,

d_{j}

is set to zero.

m + j

is the value of

m + j

Let

v_{i}

denote the

i

th row of

A

, for

i = 1, 2, \dots, m

. The following describes the printout for each row (or constraint). A full stop (.) is printed for any numerical value that is zero.

Number

is the value of

n + i

. (This is used internally to refer to

s_{i}

in the intermediate output.)

Row

gives the name of

ν_{i}

State

gives the state of the variable (LL if active on its lower bound, UL if active on its upper bound, EQ if active and fixed, BS if inactive when

s_{i}

is basic and SBS if inactive when

s_{i}

is superbasic).

A key is sometimes printed before State to give some additional information about the state of

s_{i}

. Note that unless the optional parameter

Scale Option = 0

(

default value = 2

) is specified, the tests for assigning a key are applied to the variables of the scaled problem.

A	Alternative optimum possible. The variable is nonbasic, but its reduced gradient is essentially zero. This means that if the variable were allowed to start moving away from its bound, there would be no change to the objective function. The values of the other free variables might change, giving a genuine alternative solution. However, if there are any degenerate variables (labelled D), the actual change might prove to be zero, since one of them could encounter a bound immediately. In either case the values of the Lagrange-multipliers might also change.
D	Degenerate. The variable is basic or superbasic, but it is equal to (or very close to) one of its bounds.
I	Infeasible. The variable is basic or superbasic and is currently violating one of its bounds by more than the value of the optional parameter Feasibility Tolerance ( $default value = \max (10^{- 6}, \sqrt{ε})$ , where $ε$ is the machine precision).
N	Not precisely optimal. The variable is nonbasic or superbasic. If the value of the reduced gradient for the variable exceeds the value of the optional parameter Optimality Tolerance ( $default value = \max (10^{- 6}, \sqrt{ε})$ ), the solution would not be declared optimal because the reduced gradient for the variable would not be considered negligible.

Activity

is the value of

v_{i}

at the final iterate.

Slack Activity

is the value by which

v_{i}

differs from its nearest bound. (For the free row (if any), it is set to Activity.)

Lower Bound

is the lower bound specified for the variable. None indicates that

bl [j] \leq - bigbnd

Upper Bound

is the upper bound specified for the variable. None indicates that

bu [j] \geq bigbnd

gives the index

i

v_{i}

Numerical values are output with a fixed number of digits; they are not guaranteed to be accurate to this precision.

Syntax

Parameters

Description

References

Error Indicators and Warnings

Accuracy

Parallelism and Performance

Further Comments

Description of the Printed Output

Example

Algorithmic Details

Overview

Definition of the Working Set and Search Direction

The Main Iteration

Miscellaneous

Description of Monitoring Information

See Also