NAG Library Routine Document

e04nff  (qp_dense_solve_old)
e04nfa (qp_dense_solve)

1

Purpose

2

Specification

2.1

Specification for e04nff

2.2

Specification for e04nfa

3

Description

4

References

5

Arguments

1:     $\mathbf{n}$ – IntegerInput

On entry: $n$, the number of variables.

Constraint: $\mathbf{n}>0$.

2:     $\mathbf{nclin}$ – IntegerInput

On entry: $m_L$, the number of general linear constraints.

Constraint: $\mathbf{nclin}\ge 0$.

3:     $\mathbf{a}\left(\mathbf{lda},*\right)$ – Real (Kind=nag_wp) arrayInput

Note: the second dimension of the array a must be at least $\mathbf{n}$ if $\mathbf{nclin}>0$ and at least $1$ if $\mathbf{nclin}=0$.

On entry: the $\mathit{i}$th row of a must contain the coefficients of the $\mathit{i}$th general linear constraint, for $\mathit{i}=1,2,\dots ,m_L$.

If $\mathbf{nclin}=0$, a is not referenced.

4:     $\mathbf{lda}$ – IntegerInput

On entry: the first dimension of the array a as declared in the (sub)program from which e04nff/e04nfa is called.

Constraint: $\mathbf{lda}\ge \max \left(1,\mathbf{nclin}\right)$.

5:     $\mathbf{bl}\left(\mathbf{n}+\mathbf{nclin}\right)$ – Real (Kind=nag_wp) arrayInput

6:     $\mathbf{bu}\left(\mathbf{n}+\mathbf{nclin}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: bl must contain the lower bounds and bu the upper bounds, for all the constraints in the following order. The first $n$ elements of each array must contain the bounds on the variables, and the next $m_L$ elements the bounds for the general linear constraints (if any). To specify a nonexistent lower bound (i.e., $l_j=-\infty$), set $\mathbf{bl}\left(j\right)\le -\mathit{bigbnd}$, and to specify a nonexistent upper bound (i.e., $u_j=+\infty$), set $\mathbf{bu}\left(j\right)\ge \mathit{bigbnd}$; the default value of $\mathit{bigbnd}$ is ${10}^{20}$, but this may be changed by the optional parameter Infinite Bound Size. To specify the $j$th constraint as an equality, set $\mathbf{bl}\left(j\right)=\mathbf{bu}\left(j\right)=\beta$, say, where $\left|\beta \right|<\mathit{bigbnd}$.

Constraints:

• $\mathbf{bl}\left(\mathit{j}\right)\le \mathbf{bu}\left(\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,\mathbf{n}+\mathbf{nclin}$;
• if $\mathbf{bl}\left(j\right)=\mathbf{bu}\left(j\right)=\beta$, $\left|\beta \right|<\mathit{bigbnd}$.

7:     $\mathbf{cvec}\left(*\right)$ – Real (Kind=nag_wp) arrayInput

Note: the dimension of the array cvec must be at least $\mathbf{n}$ if the problem is of type LP, QP2 (the default) or QP4, and at least $1$ otherwise.

On entry: the coefficients of the explicit linear term of the objective function when the problem is of type LP, QP2 (the default) and QP4.

If the problem is of type FP, QP1, or QP3, cvec is not referenced.

8:     $\mathbf{h}\left(\mathbf{ldh},*\right)$ – Real (Kind=nag_wp) arrayInput

Note: the second dimension of the array h must be at least $\mathbf{n}$ if it is to be used to store $H$ explicitly, and at least $1$ otherwise.

On entry: may be used to store the quadratic term $H$ of the QP objective function if desired. In some cases, you need not use h to store $H$ explicitly (see the specification of subroutine qphess). The elements of h are referenced only by subroutine qphess. The number of rows of $H$ is denoted by $m$, whose default value is $n$. (The optional parameter Hessian Rows may be used to specify a value of $m<n$.)

If the default version of qphess is used and the problem is of type QP1 or QP2 (the default), the first $m$ rows and columns of h must contain the leading $m$ by $m$ rows and columns of the symmetric Hessian matrix $H$. Only the diagonal and upper triangular elements of the leading $m$ rows and columns of h are referenced. The remaining elements need not be assigned.

If the default version of qphess is used and the problem is of type QP3 or QP4, the first $m$ rows of h must contain an $m$ by $n$ upper trapezoidal factor of the symmetric Hessian matrix $H^{\mathrm{T}}H$. The factor need not be of full rank, i.e., some of the diagonal elements may be zero. However, as a general rule, the larger the dimension of the leading nonsingular sub-matrix of h, the fewer iterations will be required. Elements outside the upper trapezoidal part of the first $m$ rows of h need not be assigned.

If a non-default version of qphess is supplied, then in some cases it may be desirable to use a one-dimensional array to transmit data to qphess. (This is illustrated in the example program in Section 10 in e04ngf/e04nga.) h is then declared as an ldh by $1$ array, where $\mathbf{ldh}\ge \mathbf{n}\times \left(\mathbf{n}+1\right)/2$.

In other situations, it may be desirable to compute $Hx$ or $H^{\mathrm{T}}Hx$ without accessing h – for example, if $H$ or $H^{\mathrm{T}}H$ is sparse or has special structure. The arguments h and ldh may then refer to any convenient array.

If the problem is of type FP or LP, h is not referenced.

9:     $\mathbf{ldh}$ – IntegerInput

On entry: the first dimension of the array h as declared in the (sub)program from which e04nff/e04nfa is called.

Constraints:

• if the problem is of type QP1, QP2 (the default), QP3 or QP4, $\mathbf{ldh}\ge \mathbf{n}$ or at least the value of the optional parameter Hessian Rows;
• if the problem is of type FP or LP, $\mathbf{ldh}\ge 1$.

10:   $\mathbf{qphess}$ – Subroutine, supplied by the NAG Library or the user.External Procedure

In general, you need not provide a version of qphess, because a ‘default’ subroutine with name e04nfu/e54nfu is included in the Library. However, the algorithm of e04nff/e04nfa requires only the product of $H$ or $H^{\mathrm{T}}H$ and a vector $x$; and in some cases you may obtain increased efficiency by providing a version of qphess that avoids the need to define the elements of the matrices $H$ or $H^{\mathrm{T}}H$ explicitly.

qphess is not referenced if the problem is of type FP or LP, in which case qphess may be the routine e04nfu/e54nfu.

The specification of qphess for e04nff is:

Fortran Interface

Subroutine qphess ( n, jthcol, h, ldh, x, hx) Integer, Intent (In) :: n, jthcol, ldh Real (Kind=nag_wp), Intent (In) :: h(ldh,*), x(n) Real (Kind=nag_wp), Intent (Out) :: hx(n)

C Header Interface

#include nagmk26.h void  qphess ( const Integer *n, const Integer *jthcol, const double h[], const Integer *ldh, const double x[], double hx[])

The specification of qphess for e04nfa is:

Fortran Interface

Subroutine qphess ( n, jthcol, h, ldh, x, hx, iuser, ruser, iwsav) Integer, Intent (In) :: n, jthcol, ldh Integer, Intent (Inout) :: iuser(*), iwsav(610) Real (Kind=nag_wp), Intent (In) :: h(ldh,*), x(n) Real (Kind=nag_wp), Intent (Inout) :: ruser(*) Real (Kind=nag_wp), Intent (Out) :: hx(n)

C Header Interface

#include nagmk26.h void  qphess ( const Integer *n, const Integer *jthcol, const double h[], const Integer *ldh, const double x[], double hx[], Integer iuser[], double ruser[], Integer iwsav[])

1:     $\mathbf{n}$ – IntegerInput

On entry: this is the same argument as supplied to this routine. See the description for the top level argument n.

2:     $\mathbf{jthcol}$ – IntegerInput

On entry: specifies whether or not the vector $x$ is a column of the identity matrix.

$\mathbf{jthcol}=j>0$

The vector $x$ is the $j$th column of the identity matrix, and hence $Hx$ or $H^{\mathrm{T}}Hx$ is the $j$th column of $H$ or $H^{\mathrm{T}}H$, respectively. This may in some cases require very little computation and qphess may be coded to take advantage of this. However special code is not necessary because $x$ is always stored explicitly in the array x.

$\mathbf{jthcol}=0$

$x$ has no special form.

3:     $\mathbf{h}\left(\mathbf{ldh},*\right)$ – Real (Kind=nag_wp) arrayInput

On entry: this is the same argument as supplied to this routine. See the description for the top level argument h.

4:     $\mathbf{ldh}$ – IntegerInput

On entry: this is the same argument as supplied to this routine. See the description for the top level argument ldh.

5:     $\mathbf{x}\left(\mathbf{n}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: the vector $x$.

6:     $\mathbf{hx}\left(\mathbf{n}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the product $Hx$ if the problem is of type QP1 or QP2 (the default), or the product $H^{\mathrm{T}}Hx$ if the problem is of type QP3 or QP4.

Note: the following are additional arguments for specific use with e04nfa. Users of e04nff therefore need not read the remainder of this description.

7:     $\mathbf{iuser}\left(*\right)$ – Integer arrayUser Workspace

8:     $\mathbf{ruser}\left(*\right)$ – Real (Kind=nag_wp) arrayUser Workspace

qphess is called with the arguments iuser and ruser as supplied to e04nff/e04nfa. You should use the arrays iuser and ruser to supply information to qphess.

9:     $\mathbf{iwsav}\left(610\right)$ – Integer arrayCommunication Array

iwsav contains information that is required by the default routine e54nfu.

qphess must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which e04nff/e04nfa is called. Arguments denoted as Input must not be changed by this procedure.

Note: qphess should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by e04nff/e04nfa. If your code inadvertently does return any NaNs or infinities, e04nff/e04nfa is likely to produce unexpected results.

11:   $\mathbf{istate}\left(\mathbf{n}+\mathbf{nclin}\right)$ – Integer arrayInput/Output

On entry: need not be set if the (default) optional parameter Cold Start is used.

If the optional parameter Warm Start has been chosen, istate specifies the desired status of the constraints at the start of the feasibility phase. More precisely, the first $n$ elements of istate refer to the upper and lower bounds on the variables, and the next $m_L$ elements refer to the general linear constraints (if any). Possible values for $\mathbf{istate}\left(j\right)$ are as follows:

$\mathbf{istate}\left(j\right)$	Meaning
0	The corresponding constraint should not be in the initial working set.
1	The constraint should be in the initial working set at its lower bound.
2	The constraint should be in the initial working set at its upper bound.
3	The constraint should be in the initial working set as an equality. This value must not be specified unless $\mathbf{bl}\left(j\right)=\mathbf{bu}\left(j\right)$.

The values $-2$, $-1$ and $4$ are also acceptable but will be reset to zero by the routine. If e04nff/e04nfa has been called previously with the same values of n and nclin, istate already contains satisfactory information. (See also the description of the optional parameter Warm Start.) The routine also adjusts (if necessary) the values supplied in x to be consistent with istate.

Constraint: $-2\le \mathbf{istate}\left(\mathit{j}\right)\le 4$, for $\mathit{j}=1,2,\dots ,\mathbf{n}+\mathbf{nclin}$.

On exit: the status of the constraints in the working set at the point returned in x. The significance of each possible value of $\mathbf{istate}\left(j\right)$ is as follows:

$\mathbf{istate}\left(j\right)$	Meaning
$-2$	The constraint violates its lower bound by more than the feasibility tolerance.
$-1$	The constraint violates its upper bound by more than the feasibility tolerance.
$\phantom{-}0$	The constraint is satisfied to within the feasibility tolerance, but is not in the working set.
$\phantom{-}1$	This inequality constraint is included in the working set at its lower bound.
$\phantom{-}2$	This inequality constraint is included in the working set at its upper bound.
$\phantom{-}3$	This constraint is included in the working set as an equality. This value of istate can occur only when $\mathbf{bl}\left(j\right)=\mathbf{bu}\left(j\right)$.
$\phantom{-}4$	This corresponds to optimality being declared with $\mathbf{x}\left(j\right)$ being temporarily fixed at its current value. This value of istate can occur only when $\mathbf{ifail}=\mathbf{1}$ on exit.

12:   $\mathbf{x}\left(\mathbf{n}\right)$ – Real (Kind=nag_wp) arrayInput/Output

On entry: an initial estimate of the solution.

On exit: the point at which e04nff/e04nfa terminated. If $\mathbf{ifail}=\mathbf{0}$, $\mathbf{1}$ or $\mathbf{4}$, x contains an estimate of the solution.

13:   $\mathbf{iter}$ – IntegerOutput

On exit: the total number of iterations performed.

14:   $\mathbf{obj}$ – Real (Kind=nag_wp)Output

On exit: the value of the objective function at $x$ if $x$ is feasible, or the sum of infeasibilities at $x$ otherwise. If the problem is of type FP and $x$ is feasible, obj is set to zero.

15:   $\mathbf{ax}\left(\max \left(1,\mathbf{nclin}\right)\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the final values of the linear constraints $Ax$.

If $\mathbf{nclin}=0$, ax is not referenced.

16:   $\mathbf{clamda}\left(\mathbf{n}+\mathbf{nclin}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the values of the Lagrange multipliers for each constraint with respect to the current working set. The first $n$ elements contain the multipliers for the bound constraints on the variables, and the next $m_L$ elements contain the multipliers for the general linear constraints (if any). If $\mathbf{istate}\left(j\right)=0$ (i.e., constraint $j$ is not in the working set), $\mathbf{clamda}\left(j\right)$ is zero. If $x$ is optimal, $\mathbf{clamda}\left(j\right)$ should be non-negative if $\mathbf{istate}\left(j\right)=1$, non-positive if $\mathbf{istate}\left(j\right)=2$ and zero if $\mathbf{istate}\left(j\right)=4$.

17:   $\mathbf{iwork}\left(\mathbf{liwork}\right)$ – Integer arrayWorkspace

18:   $\mathbf{liwork}$ – IntegerInput

On entry: the dimension of the array iwork as declared in the (sub)program from which e04nff/e04nfa is called.

Constraint: $\mathbf{liwork}\ge 2\times \mathbf{n}+3$.

19:   $\mathbf{work}\left(\mathbf{lwork}\right)$ – Real (Kind=nag_wp) arrayWorkspace

20:   $\mathbf{lwork}$ – IntegerInput

On entry: the dimension of the array work as declared in the (sub)program from which e04nff/e04nfa is called.

Constraints:

• if the problem is of type QP2 (the default) or QP4,
  • if $\mathbf{nclin}>0$, $\mathbf{lwork}\ge 2\times {\mathbf{n}}^2+8\times \mathbf{n}+5\times \mathbf{nclin}$;
  • if $\mathbf{nclin}=0$, $\mathbf{lwork}\ge {\mathbf{n}}^2+8\times \mathbf{n}$;
• if the problem is of type QP1 or QP3,
  • if $\mathbf{nclin}>0$, $\mathbf{lwork}\ge 2\times {\mathbf{n}}^2+7\times \mathbf{n}+5\times \mathbf{nclin}$;
  • if $\mathbf{nclin}=0$, $\mathbf{lwork}\ge {\mathbf{n}}^2+7\times \mathbf{n}$;
• if the problem is of type LP,
  • if $\mathbf{nclin}=0$, $\mathbf{lwork}\ge 8\times \mathbf{n}+1$;
  • if $\mathbf{nclin}\ge \mathbf{n}$, $\mathbf{lwork}\ge 2\times {\mathbf{n}}^2+8\times \mathbf{n}+5\times \mathbf{nclin}$;
  • otherwise $\mathbf{lwork}\ge 2\times {\left(\mathbf{nclin}+1\right)}^2+8\times \mathbf{n}+5\times \mathbf{nclin}$;
• if the problem is of type FP,
  • if $\mathbf{nclin}=0$, $\mathbf{lwork}\ge 7\times \mathbf{n}+1$;
  • if $\mathbf{nclin}\ge \mathbf{n}$, $\mathbf{lwork}\ge 2\times {\mathbf{n}}^2+7\times \mathbf{n}+5\times \mathbf{nclin}$;
  • otherwise $\mathbf{lwork}\ge 2\times {\left(\mathbf{nclin}+1\right)}^2+7\times \mathbf{n}+5\times \mathbf{nclin}$.

The amounts of workspace provided and required are (by default) output on the current advisory message unit (as defined by x04abf). As an alternative to computing liwork and lwork from the formulas given above, you may prefer to obtain appropriate values from the output of a preliminary run with liwork and lwork set to $1$. (e04nff/e04nfa will then terminate with $\mathbf{ifail}=\mathbf{6}$.)

21:   $\mathbf{ifail}$ – IntegerInput/Output

Note: for e04nfa, ifail does not occur in this position in the argument list. See the additional arguments described below.

On entry: ifail must be set to $0$, $-1\text{ or }1$. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.

For environments where it might be inappropriate to halt program execution when an error is detected, the value $-1\text{ or }1$ is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, because for this routine the values of the output arguments may be useful even if $\mathbf{ifail}\ne \mathbf{0}$ on exit, the recommended value is $-1$. When the value $-\mathbf{1}\text{ or }1$ is used it is essential to test the value of ifail on exit.

On exit: $\mathbf{ifail}=\mathbf{0}$ unless the routine detects an error or a warning has been flagged (see Section 6).

e04nff/e04nfa returns with $\mathbf{ifail}=\mathbf{0}$ if $x$ is a strong local minimizer, i.e., the reduced gradient (Norm Gz; see Section 9.2) is negligible, the Lagrange multipliers (Lagr Mult; see Section 9.2) are optimal and $H_R$ (the reduced Hessian of $f\left(x\right)$; see Section 11.2) is positive semidefinite.

Note: the following are additional arguments for specific use with e04nfa. Users of e04nff therefore need not read the remainder of this description.

21:   $\mathbf{iuser}\left(*\right)$ – Integer arrayUser Workspace

22:   $\mathbf{ruser}\left(*\right)$ – Real (Kind=nag_wp) arrayUser Workspace

iuser and ruser are not used by e04nff/e04nfa, but are passed directly to qphess and may be used to pass information to this routine.

23:   $\mathbf{lwsav}\left(120\right)$ – Logical arrayCommunication Array

24:   $\mathbf{iwsav}\left(610\right)$ – Integer arrayCommunication Array

25:   $\mathbf{rwsav}\left(475\right)$ – Real (Kind=nag_wp) arrayCommunication Array

The arrays lwsav, iwsav and rwsav must not be altered between calls to any of the routines e04nfa, e04nga or e04nha.

26:   $\mathbf{ifail}$ – IntegerInput/Output

Note: see the argument description for ifail above.

6

Error Indicators and Warnings

$\mathbf{ifail}=1$

The iterations were terminated at a dead-point. The necessary conditions for optimality are satisfied but the sufficient conditions are not. (The reduced gradient is negligible, the Lagrange multipliers are optimal, but $H_R$ is singular or there are some very small multipliers.) If ${\nabla }^{2}f\left(x\right)$ is not positive definite, $x$ is not necessarily a local solution of the problem and verification of optimality requires further information. If ${\nabla }^{2}f\left(x\right)$ is positive semidefinite or the problem is of type LP, $x$ gives the global minimum value of the objective function, but the final $x$ is not unique.

$\mathbf{ifail}=2$

The solution appears to be unbounded, i.e., the objective function is not bounded below in the feasible region. This value of ifail occurs if a step larger than Infinite Step Size ($\text{default value}={10}^{20}$) would have to be taken in order to continue the algorithm, or the next step would result in an element of $x$ having magnitude larger than Infinite Bound Size ($\text{default value}={10}^{20}$).

$\mathbf{ifail}=3$

No feasible point was found, i.e., it was not possible to satisfy all the constraints to within the feasibility tolerance. In this case, the constraint violations at the final $x$ will reveal a value of the tolerance for which a feasible point will exist – for example, when the feasibility tolerance for each violated constraint exceeds its Slack (see Section 9.2) at the final point. The modified problem (with an altered feasibility tolerance) may then be solved using a Warm Start. You should check that there are no constraint redundancies. If the data for the constraints are accurate only to the absolute precision $\sigma$, you should ensure that the value of the optional parameter Feasibility Tolerance ($\text{default value}=\sqrt{\epsilon }$, where $\epsilon$ is the machine precision) is greater than $\sigma$. For example, if all elements of $A$ are of order unity and are accurate only to three decimal places, the Feasibility Tolerance should be at least ${10}^{-3}$.

$\mathbf{ifail}=4$

The limiting number of iterations was reached before normal termination occurred.

The values of the optional parameters Feasibility Phase Iteration Limit ($\text{default value}=\max \left(50,5\left(n+m_L\right)\right)$) and Optimality Phase Iteration Limit ($\text{default value}=\max \left(50,5\left(n+m_L\right)\right)$) may be too small. If the method appears to be making progress (e.g., the objective function is being satisfactorily reduced), either increase the iterations limit and rerun e04nff/e04nfa or, alternatively, rerun e04nff/e04nfa using the Warm Start facility to specify the initial working set.

$\mathbf{ifail}=5$

The reduced Hessian exceeds its assigned dimension. The algorithm needed to expand the reduced Hessian when it was already at its maximum dimension, as specified by the optional parameter Maximum Degrees of Freedom ($\text{default value}=n$).

The value of the optional parameter Maximum Degrees of Freedom is too small. Rerun e04nff/e04nfa with a larger value (possibly using the Warm Start facility to specify the initial working set).

$\mathbf{ifail}=6$

An input argument is invalid.

$\mathbf{ifail}=7$

The designated problem type was not FP, LP, QP1, QP2, QP3 or QP4. Rerun e04nff/e04nfa with the optional parameter Problem Type set to one of these values.

$\mathbf{\text{Overflow}}$

If the printed output before the overflow error contains a warning about serious ill-conditioning in the working set when adding the $j$th constraint, it may be possible to avoid the difficulty by increasing the magnitude of the Feasibility Tolerance ($\text{default value}=\sqrt{\epsilon }$, where $\epsilon$ is the machine precision) and rerunning the program. If the message recurs even after this change, the offending linearly dependent constraint (with index ‘$j$’) must be removed from the problem.

$\mathbf{ifail}=-99$

An unexpected error has been triggered by this routine. Please contact NAG.

See Section 3.9 in How to Use the NAG Library and its Documentation for further information.

$\mathbf{ifail}=-399$

Your licence key may have expired or may not have been installed correctly.

See Section 3.8 in How to Use the NAG Library and its Documentation for further information.

$\mathbf{ifail}=-999$

Dynamic memory allocation failed.

See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

7

Accuracy

8

Parallelism and Performance

9

Further Comments

9.1

Scaling

9.2

Description of the Printed Output

10

Example

10.1

Program Text

Program Text (e04nffe.f90)

Program Text (e04nfae.f90)

10.2

Program Data

Program Data (e04nffe.d)

Program Data (e04nfae.d)

10.3

Program Results

Program Results (e04nffe.r)

Program Results (e04nfae.r)

11

Algorithmic Details

11.1

Overview

11.2

Definition of Search Direction

11.3

Main Iteration

11.4

Choosing the Initial Working Set

12

Optional Parameters

• Check Frequency
• Cold Start
• Crash Tolerance
• Defaults
• Expand Frequency
• Feasibility Phase Iteration Limit
• Feasibility Tolerance
• Hessian Rows
• Infinite Bound Size
• Infinite Step Size
• Iteration Limit
• Iters
• Itns
• List
• Maximum Degrees of Freedom
• Minimum Sum of Infeasibilities
• Monitoring File
• Nolist
• Optimality Phase Iteration Limit
• Optimality Tolerance
• Print Level
• Problem Type
• Rank Tolerance
• Warm Start

Begin 
  Print Level = 5 
End

Call e04ngf (ioptns, inform)

Call e04nhf ('Print Level = 5')

12.1

Description of the Optional Parameters

• the keywords, where the minimum abbreviation of each keyword is underlined (if no characters of an optional qualifier are underlined, the qualifier may be omitted);
• a parameter value, where the letters $a$, $i$ and $r$ denote options that take character, integer and real values respectively;
• the default value, where the symbol $\epsilon$ is a generic notation for machine precision (see x02ajf).

13

Description of Monitoring Information