***********************************************************************
* *
* PIND - A DISCRETE NEWTON ALGORITHM FOR LARGE-SCALE *
* EQUALITY CONSTRAINED OPTIMIZATION BASED ON THE *
* USE OF INDEFINITE PRECONDITIONERS. *
* *
***********************************************************************
1. Introduction:
----------------
The double-precision FORTRAN 77 basic subroutine PIND is designed
to find a close approximation to a local minimum of a twice continuously
differentiable function F(X) on the set defined by the nonlinear
equality constraints
FC_1(X) = 0, FC_2(X) = 0, ..., FC_NC(X)=0.
Here X is a vector of NF variables and FC_I(X), 1 <= I <= NC, are twice
continuously differentiable functions. We assume that the sparsity
pattern of the objective Hessian matrix is known and that NF and NC are
large, but constraint functions FC_I(X), 1 <= I <= NC depend on a small
number of variables. The latter assumption implies that the mapping
CF(X) = [FC_1(X), FC_2(X), ..., FC_NC(X)] has a sparse Jacobian matrix,
which will be denoted by CG(X) (it has NC rows and NF columns). The
sparsity pattern of the objective Hessian matrix (only the upper part) is
stored in the standard compressed row format using arrays IH and JH. For
example, if the Hessian matrix has the following pattern
H = | * * * 0 * |
| * * 0 * 0 |
| * 0 * 0 * |
| 0 * 0 * 0 |
| * 0 * 0 * |
(asterisks denote nonzero elements) then arrays IH and JH contain
elements
IH(1)=1, IH(2)=5, IH(3)=7, IH(4)=9, IH(5)=10, IH(6)=11,
JH(1)=1, JH(2)=2, JH(3)=3, JH(4)=5, JH(5)=2, JH(6)=4, JH(7)=3,
JH(8)=5, JH(9)=4, JH(10)=5,
i.e., IH contains pointers of the diagonal elements in the upper part
of the Hessian matrix and JH contains column indices of the nonzero
elements. Note that IH has NF+1 elements, the last element is equal
MH+1, where MH is the number of nonzero elements stored. The sparsity
pattern of the constraint Jacobian matrix is stored in the standard
compressed row format using arrays ICG and JCG. For example, if the
Jacobian matrix has the following pattern
CG = | * * 0 * |
| * * * 0 |
| * 0 0 * |
| 0 * * 0 |
| * 0 * 0 |
(asterisks denote nonzero elements) then arrays ICG and JCG contain
elements
ICG(1)=1, ICG(2)=4, ICG(3)=7, ICG(4)=9, ICG(5)=11, ICG(6)=13,
JCG(1)=1, JCG(2)=2, JCG(3)=4, JCG(4)=1, JCG(5)=2, JCG(6)=3,
JCG(7)=1, JCG(8)=4, JCG(9)=2, JCG(10)=3, JCG(11)=1, JCG(12)=3,
i.e., ICG contains pointers of the first elements in rows of the
Jacobian matrix and JCG contains column indices of the nonzero elements.
Note that ICG has NC+1 elements, the last element is equal MC+1, where
MC is the number of nonzero elements.
To simplify user's work, an additional easy to use subroutine
PINDU is added. It calls the basic general subroutine PIND. All
subroutines contain a description of formal parameters and extensive
comments. Furthermore, test program TINDU is included, which contains
several test problems (see e.g. [2]). This test program serve as an
example for using the subroutine PINDU, verify its correctness and
demonstrate its efficiency.
In this short guide, we describe all subroutines which can be
called from the user's program. A detailed description of the method is
given in [1]. In the description of formal parameters, we introduce a
type of the argument that specifies whether the argument must have a
value defined on entry to the subroutine (I), whether it is a value
which will be returned (O), or both (U), or whether it is an auxiliary
value (A). Note that the arguments of the type I can be changed on
output under some circumstances, especially if improper input values
were given. Besides formal parameters, we can use a COMMON /STAT/ block
containing statistical information. This block, used in each subroutine
has the following form:
COMMON /STAT/ NRES,NDEC,NIIT,NIT,NFV,NFG,NFH
The arguments have the following meaning:
Argument Type Significance
----------------------------------------------------------------------
NRES O Positive INTEGER variable that indicates the number of
restarts.
NDEC O Positive INTEGER variable that indicates the number of
matrix decompositions.
NIIT O Positive INTEGER variable that indicates the number of
inner iterations (for solving linear systems).
NIT O Positive INTEGER variable that indicates the number of
iterations.
NFV O Positive INTEGER variable that indicates the number of
function evaluations.
NFG O Positive INTEGER variable that specifies the number of
gradient evaluations.
NFH O Positive INTEGER variable that specifies the number of
Hessian evaluations.
2. Subroutine PINDU:
--------------------
The calling sequence is
CALL PINDU(NF,NC,X,IH,JH,CF,ICG,JCG,IPAR,RPAR,F,GMAX,CMAX,
& IPRNT,ITERM)
The arguments have the following meaning.
Argument Type Significance
----------------------------------------------------------------------
NF I Positive INTEGER variable that specifies the number of
variables of the partially separable function.
NC I Positive INTEGER variable that specifies the number of
constraints.
X(NF) U On input, DOUBLE PRECISION vector with the initial
estimate to the solution. On output, the approximation
to the minimum.
IH(NF+1) I INTEGER array which contains pointers of the diagonal
elements in the upper part of the Lagrangian function
Hessian matrix.
JH(NH) I INTEGER array which contains column indices of the
nonzero elements and additional working space for the
Choleski factor.
CF(NC+1) O DOUBLE PRECISION vector which contains values of
constraint functions.
ICG(NC+1) I INTEGER array which contains pointers of the first
elements in rows of the constraint Jacobian matrix.
JCG(MC) I INTEGER array which contains column indices of the
nonzero elements.
IPAR(7) I INTEGER parameters:
IPAR(1)=MIT, IPAR(2)=MFV, IPAR(3)=MFG,
IPAR(4)=MOS, IPAR(5)=MD, IPAR(6)=MDE,
IPAR(7)=IFIL.
Parameters MIT, MFV, MFG, MOS, MD, MDE are described
in Section 3 together with other parameters of the
subroutine PIND. Parameter IFIL specifies a relative
size of the space reserved for fill-in. The choice
IFIL=0 causes that the default value IFIL=1 will be
taken.
RPAR(5) I DOUBLE PRECISION parameters:
RPAR(1)=XMAX, RPAR(2)=TOLX, RPAR(3)=TOLC,
RPAR(4)=TOLG, RPAR(5)=RPF1.
Parameters XMAX, TOLX, TOLC, TOLG, RPF1 are described
in Section 3 together with other parameters of the
subroutine PIND.
F O DOUBLE PRECISION value of the objective function at the
solution X.
GMAX O DOUBLE PRECISION maximum absolute value of a partial
derivative of the Lagrangian function.
CMAX O maximum constraint violation.
IPRNT I INTEGER variable that specifies PRINT:
IPRNT= 0 - print is suppressed,
IPRNT= 1 - basic print of final results,
IPRNT=-1 - extended print of final results,
IPRNT= 2 - basic print of intermediate and final
results,
IPRNT=-2 - extended print of intermediate and final
results.
ITERM O INTEGER variable that indicates the cause of termination:
ITERM= 1 - if |X - XO| was less than or equal to TOLX
in two subsequent iterations,
ITERM= 2 - if |F - FO| was less than or equal to TOLF
in two subsequent iterations,
ITERM= 3 - if F is less than or equal to TOLB,
ITERM= 4 - if GMAX is less than or equal to TOLG,
ITERM=11 - if NIT exceeded MIT,
ITERM=12 - if NFV exceeded MFV,
ITERM=13 - if NFG exceeded MFG,
ITERM< 0 - if the method failed. If ITERM=-6, then the
termination criterion has not been satisfied,
but the point obtained is usually acceptable.
Values ITERM<=-40 detect a lack of space. In
this case, parameter IPAR(7)=IFIL has to be
increased (IFIL=2, IFIL=3, etc.).
The subroutine PINDU requires the user supplied subroutines OBJ,
DOBJ that define the objective function and its gradient and CON, DCON
that define constraint functions and their gradients. These subroutines
have the form
SUBROUTINE OBJ(NF,X,F)
SUBROUTINE DOBJ(NF,X,G)
SUBROUTINE CON(NF,KC,X,FC)
SUBROUTINE DCON(NF,KC,X,GC)
The arguments of the user supplied subroutines have the following
meaning.
Argument Type Significance
----------------------------------------------------------------------
NF I Positive INTEGER variable that specifies the number of
variables of the objective function.
X(NF) I DOUBLE PRECISION an estimate to the solution.
F O DOUBLE PRECISION value of the objective function at the
point X.
G(NF) O DOUBLE PRECISION gradient of the objective function
at the point X.
KC I INTEGER index of the partial function.
FC O DOUBLE PRECISION value of the KC-th partial function at
the point X.
GC(NF) O DOUBLE PRECISION gradient of the KC-th partial function
at the point X.
3. Subroutine PIND:
-------------------
This general subroutine is called from all subroutines described
in Section 2. The calling sequence is
CALL PIND(NF,NC,MAXH,MAXB,X,GF,H,IH,JH,CG,ICG,JCG,KCG,CF,CFO,
& GC,G,S,XO,GO,XS,CZ,CZO,CP,B,IB,JB,CGD,DD,BD,COL,PSR,PERC,INVP,
& WN11,WN12,WN13,WN14,XMAX,TOLX,TOLC,TOLG,RFF1,F,GMAX,CMAX,MIT,
& MFV,MFG,MOS,MD,MDE,IPRNT,ITERM)
The arguments NF, NC, X, IH, JH, CF, ICG, JCG, F, GMAX, CMAX, IPRNT,
ITERM have the same meaning as in Section 2. Other arguments have the
following meaning:
Argument Type Significance
----------------------------------------------------------------------
MAXH I INTEGER size of array H.
MAXB I INTEGER size of array B.
GF(NF) A DOUBLE PRECISION gradient of the objective function.
H(MAXH) A DOUBLE PRECISION nonzero elements of the approximation
of the Hessian matrix and nonzero elements of the
Choleski factor.
CG(MC) A DOUBLE PRECISION nonzero elements of the Jacobian
matrix.
KCG A INTEGER auxiliary array.
CFO(NC+1) A DOUBLE PRECISION vector which contains old values of
constraint functions.
GC(NF) A DOUBLE PRECISION gradient of the constraint function.
G(NF) A DOUBLE PRECISION gradient of the merit function.
S(NF+NC) A DOUBLE PRECISION direction vector.
XO(NF+NC) A DOUBLE PRECISION array which contains increments of
variables.
GO(NF+NC) A DOUBLE PRECISION array which contains increments of
gradients.
XS(NF+NC) A DOUBLE PRECISION auxiliary array.
CZ(NC) A DOUBLE PRECISION vector of Lagrange multipliers.
CZO(NC) A DOUBLE PRECISION old vector of Lagrange multipliers.
CP(NF+NC) A DOUBLE PRECISION auxiliary array.
B(MAXB) A nonzero elements of the preconditioner together with
the nonzero elements of its Choleski factor.
IB(NC+1) A INTEGER array which contains pointers of the diagonal
elements in the upper part of the preconditioner.
JB(MAXB) A INTEGER array which contains column indices of the
nonzero elements.
GCD(NF*ND) A DOUBLE PRECISION array which contains dense columns.
DD(ND) A DOUBLE PRECISION auxiliary array.
BD(ND*ND) A DOUBLE PRECISION auxiliary array.
COL(NF) A INTEGER auxiliary array.
PSR(NC+1) A INTEGER pointer vector in the compact form of the
Choleski factor.
PERC(NC) A INTEGER permutation vector
INVP(NC) A INTEGER inverse permutation vector
WN11(NC+1) A INTEGER auxiliary array.
WN12(NC+1) A INTEGER auxiliary array.
WN13(NC+1) A INTEGER auxiliary array.
WN14(NC+1) A INTEGER auxiliary array.
XMAX I DOUBLE PRECISION maximum stepsize; the choice XMAX=0
causes that the default value 1.0D+16 will be taken.
TOLX I DOUBLE PRECISION tolerance for the change of the
coordinate vector X; the choice TOLX=0 causes that the
default value TOLX=1.0D-16 will be taken.
TOLC I DOUBLE PRECISION tolerance for the constraint violation;
the choice TOLC=0 causes that the default
value TOLC=1.0D-6 will be taken.
TOLG I DOUBLE PRECISION tolerance for the Lagrangian function
gradient; the choice TOLG=0 causes that the default
value TOLG=1.0D-6 will be taken.
RPF1 I DOUBLE PRECISION value of the penalty parameter in the
merit function; the choice RPF1=0 causes that the
default value 1.0D-4 will be taken.
MIT I INTEGER variable that specifies the maximum number of
iterations; the choice MIT=0 causes that the default
value 1000 will be taken.
MFV I INTEGER variable that specifies the maximum number of
function evaluations; the choice |MFV|=0 causes that
the default value 1000 will be taken.
MFG I INTEGER variable that specifies the maximum number of
gradient evaluations; the choice |MFV|=0 causes that
the default value 10000 will be taken.
MOS I INTEGER choice of preconditioning strategy:
MOS= 1 - preconditioning by the constraint
preconditioner with complete
Gill-Murray decomposition.
MOS=-1 - preconditioning by the constraint
preconditioner with incomplete
Gill-Murray decomposition.
The choice MOS=0 causes that the default value MOS=1
will be taken.
MD I INTEGER maximum number of the dense rows; the choice
MD=0 causes that the default value 10 will be taken.
MDE I INTEGER maximum number of nonzero elements in a sparse
rows; the choice MDE=0 causes that the default value
50 will be taken.
The choice of parameter XMAX can be sensitive in many cases. First, the
objective function can be evaluated only in a relatively small region
(if it contains exponentials) so that the maximum stepsize is necessary.
Secondly, the problem can be very ill-conditioned far from the solution
point so that large steps can be unsuitable. Finally, if the problem has
more local solutions, a suitably chosen maximum stepsize can lead to
obtaining a better local solution. A suitable choice of parameter RPF1
can improve the efficiency of the method. The default value RPF1=0.0001
is relatively small. Larger value (RPF1=1 say) can sometimes give better
result.
The subroutine PIND requires the user supplied subroutines OBJ
DOBJ, CON, DCON, which are described in Section 2.
4. Verification of the subroutines:
-----------------------------------
Subroutine PINDU can be verified and tested using the program
TINDU. This program calls the subroutines TIUS20, TINS20 (initiation),
TFFU20 (objective function evaluation), TFGU20 (objective gradient
evaluation), TCFU20 (constraint functions evaluation) and TFGU20
(constraint gradients evaluation) containing 18 equality constrained
test problems with at most 1000 variables [2]. The results obtained by
the program TINDU on a PC computer with Microsoft Power Station Fortran
compiler have the following form.
NIT= 9 NFV= 10 NFG= 60 C= 0.177636E-14 G= 0.202238E-11 ITERM= 4
NIT= 12 NFV= 13 NFG= 182 C= 0.425254E-10 G= 0.621097E-07 ITERM= 4
NIT= 10 NFV= 11 NFG= 66 C= 0.172595E-11 G= 0.842788E-08 ITERM= 4
NIT= 16 NFV= 19 NFG= 102 C= 0.141707E-10 G= 0.405805E-09 ITERM= 4
NIT= 18 NFV= 19 NFG= 190 C= 0.642577E-08 G= 0.311534E-06 ITERM= 4
NIT= 18 NFV= 19 NFG= 266 C= 0.125855E-11 G= 0.711327E-08 ITERM= 4
NIT= 9 NFV= 11 NFG= 70 C= 0.118305E-11 G= 0.221512E-11 ITERM= 4
NIT= 38 NFV= 51 NFG= 273 C= 0.444089E-15 G= 0.111759E-07 ITERM= 4
NIT= 28 NFV= 70 NFG= 203 C= 0.134587E-11 G= 0.125862E-08 ITERM= 4
NIT= 12 NFV= 15 NFG= 78 C= 0.140633E-07 G= 0.486837E-07 ITERM= 4
NIT= 9 NFV= 10 NFG= 60 C= 0.261324E-11 G= 0.299467E-11 ITERM= 4
NIT= 7 NFV= 8 NFG= 56 C= 0.143574E-11 G= 0.288525E-10 ITERM= 4
NIT= 15 NFV= 17 NFG= 128 C= 0.177636E-14 G= 0.641856E-06 ITERM= 4
NIT= 12 NFV= 13 NFG= 91 C= 0.229161E-07 G= 0.274480E-07 ITERM= 4
NIT= 19 NFV= 35 NFG= 120 C= 0.550671E-13 G= 0.127818E-11 ITERM= 4
NIT= 8 NFV= 9 NFG= 45 C= 0.679852E-09 G= 0.608118E-07 ITERM= 4
NIT= 9 NFV= 10 NFG= 50 C= 0.222045E-15 G= 0.177636E-14 ITERM= 4
NIT= 10 NFV= 13 NFG= 55 C= 0.368464E-08 G= 0.809397E-06 ITERM= 4
NITER = 259 NFVAL = 353 NSUCC = 18
TIME= 0:00:00.84
The rows corresponding to individual test problems contain the number of
iterations NIT, the number of function evaluations NFV, the number of
gradient evaluations NFG, the constraint violation C, the norm of the
Lagrangian function gradient G and the cause of termination ITERM.
References:
-----------
[1] Luksan L., Vlcek J.: LSA: Algorithms for large-scale unconstrained
and box constrained optimization Technical Report V-896. Prague,
ICS AS CR, 2004.
[2] Luksan L., Vlcek J.: Sparse and partially separable test problems
for unconstrained and equality constrained optimization. Research
Report V-767, Institute of Computer Science, Academy of Sciences
of the Czech Republic, Prague, Czech Republic, 1998.